Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
newplotbg_newsupp.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file newplotbg_newsupp.C
1 double CBFunction(double *x, double *p)
2 {
3  double norm = p[0];
4  double alpha = p[1]; // tail start
5  double n = p[2]; // tail shape
6  double mu = p[3]; // upsilon mass
7  double sigma = p[4]; // upsilon width
8 
9  double A = pow(n/fabs(alpha),n)*TMath::Exp(-pow(fabs(alpha),2)/2.);
10  double B = n/fabs(alpha) - fabs(alpha);
11  double k = (x[0]-mu)/sigma;
12 
13  double val;
14  if( k > -alpha )
15  val = norm*TMath::Exp(-0.5*pow(k,2));
16  else
17  val = norm*A*pow(B-k,-n);
18 
19  if( TMath::IsNaN(val) ) val = 0.0;
20 
21  return val;
22 }
23 
24 double TripleCBFunction(double *x, double *p)
25 {
26  double norm1 = p[0]; // amplitude of Y(1S) peak
27  double alpha = p[1]; // tail start
28  double n = p[2]; // tail shape
29  double mu1 = p[3]; // upsilon (1S) mass
30  double sigma = p[4]; // upsilon width
31  double norm2 = p[5];
32  double norm3 = p[6];
33  double mu2 = mu1*1.0595;
34  double mu3 = mu1*1.0946;
35 
36  double A = pow(n/fabs(alpha),n)*TMath::Exp(-pow(fabs(alpha),2)/2.);
37  double B = n/fabs(alpha) - fabs(alpha);
38  double k1 = (x[0]-mu1)/sigma;
39  double k2 = (x[0]-mu2)/sigma;
40  double k3 = (x[0]-mu3)/sigma;
41 
42  double val,val1,val2,val3;
43 
44  if( k1 > -alpha ) { val1 = norm1*TMath::Exp(-0.5*pow(k1,2)); }
45  else { val1 = norm1*A*pow(B-k1,-n); }
46  if( k2 > -alpha ) { val2 = norm2*TMath::Exp(-0.5*pow(k2,2)); }
47  else { val2 = norm2*A*pow(B-k2,-n); }
48  if( k3 > -alpha ) { val3 = norm3*TMath::Exp(-0.5*pow(k3,2)); }
49  else { val3 = norm3*A*pow(B-k3,-n); }
50 
51  val = val1 + val2 + val3;
52 
53  if( TMath::IsNaN(val) ) val = 0.0;
54 
55  return val;
56 }
57 
58 double SandB_CBFunction(double *x, double *p)
59 {
60  double norm1 = p[0]; // amplitude of Y(1S) peak
61  double alpha = p[1]; // tail start
62  double n = p[2]; // tail shape
63  double mu1 = p[3]; // upsilon (1S) mass
64  double sigma = p[4]; // upsilon width
65  double norm2 = p[5];
66  double norm3 = p[6];
67  double mu2 = mu1*1.0595;
68  double mu3 = mu1*1.0946;
69 
70  double A = pow(n/fabs(alpha),n)*TMath::Exp(-pow(fabs(alpha),2)/2.);
71  double B = n/fabs(alpha) - fabs(alpha);
72  double k1 = (x[0]-mu1)/sigma;
73  double k2 = (x[0]-mu2)/sigma;
74  double k3 = (x[0]-mu3)/sigma;
75 
76  double val,val1,val2,val3;
77 
78  if( k1 > -alpha ) { val1 = norm1*TMath::Exp(-0.5*pow(k1,2)); }
79  else { val1 = norm1*A*pow(B-k1,-n); }
80  if( k2 > -alpha ) { val2 = norm2*TMath::Exp(-0.5*pow(k2,2)); }
81  else { val2 = norm2*A*pow(B-k2,-n); }
82  if( k3 > -alpha ) { val3 = norm3*TMath::Exp(-0.5*pow(k3,2)); }
83  else { val3 = norm3*A*pow(B-k3,-n); }
84 
85  double bgnorm1 = p[7];
86  double bgslope1 = p[8];
87 
88  double bg = exp(bgnorm1+x[0]*bgslope1);
89 
90  val = val1 + val2 + val3 + bg;
91  if( TMath::IsNaN(val) ) val = 0.0;
92 
93  return val;
94 }
95 
96 //----------------------------------------------------------------
97 //----------------------------------------------------------------
98 //----------------------------------------------------------------
99 
100 void newplotbg_newsupp() {
101 
102 //gROOT->LoadMacro("sPHENIXStyle/sPhenixStyle.C");
103 //SetsPhenixStyle();
104 
105 gStyle->SetOptStat(0);
106 gStyle->SetOptFit(0);
107 
108 TRandom* myrandom = new TRandom3();
109 const int nbins = 15;
110 double eideff = 0.9;
111 string times = "*";
112 TLatex* tl[15];
113 char tlchar[999];
114 
115 double statscale = 2.4; // 5-year plan all correladed bg plots are for 10B
116 //double statscale = 1.43; // 3-year plan all correladed bg plots are for 10B
117 //double statscalepp = 1.78; // old wrong???
118 double statscalepp = 1.13; // 5-year plan
119 double statscale_lowlim = 7.0;
120 double statscale_uplim = 14.0;
121 
122  TF1* fCBpp = new TF1("fCBpp",CBFunction,5.,14.,5);
123  TF1* fCBauau = new TF1("fCBauau",CBFunction,5.,14.,5);
124  TF1* fCB1s = new TF1("fCB1s",CBFunction,5.,14.,5);
125  TF1* fCB2s = new TF1("fCB2s",CBFunction,5.,14.,5);
126  TF1* fCB3s = new TF1("fCB3s",CBFunction,5.,14.,5);
127  TF1* fTCB = new TF1("fTCB",TripleCBFunction,5.,14.,7);
128  TF1* fTCBpp = new TF1("fTCBpp",TripleCBFunction,5.,14.,7);
129  TF1* fTCBauau = new TF1("fTCBauau",TripleCBFunction,5.,14.,7);
130  TF1* fSandB = new TF1("fSandB",SandB_CBFunction,5.,14.,9);
131  TF1* fSandBfordave = new TF1("fSandBfordave",SandB_CBFunction,5.,14.,9);
132  TF1* fSandBpp = new TF1("fSandBpp",SandB_CBFunction,5.,14.,9);
133  TF1* fSandBauau = new TF1("fSandBauau",SandB_CBFunction,5.,14.,9);
134 
135 //---------------------------------------------------------
136 // Upsilons
137 //---------------------------------------------------------
138 
139 double raapt[9],raa1s[9],raa2s[9],raa3s[9];
140 raapt[0] = 1.5;
141 raapt[1] = 4.5;
142 raapt[2] = 7.5;
143 raapt[3] = 10.5;
144 raapt[4] = 13.5;
145 raa1s[0] = 0.535;
146 raa1s[1] = 0.535;
147 raa1s[2] = 0.535;
148 raa1s[3] = 0.535;
149 raa1s[4] = 0.535;
150 raa2s[0] = 0.170;
151 raa2s[1] = 0.170;
152 raa2s[2] = 0.170;
153 raa2s[3] = 0.170;
154 raa2s[4] = 0.170;
155 /*
156 raa1s[0] = 0.4960;
157 raa1s[1] = 0.4960;
158 raa1s[2] = 0.4955;
159 raa1s[3] = 0.4968;
160 raa1s[4] = 0.4743;
161 raa2s[0] = 0.1710;
162 raa2s[1] = 0.1629;
163 raa2s[2] = 0.1326;
164 raa2s[3] = 0.1232;
165 raa2s[4] = 0.0928;
166 */
167 raa3s[0] = 0.035;
168 raa3s[1] = 0.035;
169 raa3s[2] = 0.035;
170 raa3s[3] = 0.035;
171 raa3s[4] = 0.035;
172 TGraph* grRAA1S = new TGraph(5,raapt,raa1s);
173 TGraph* grRAA2S = new TGraph(5,raapt,raa2s);
174 TGraph* grRAA3S = new TGraph(5,raapt,raa3s);
175 
176 double NNN = statscale * 7780./0.49; // from sPHENIX proposal for 100B minbias Au+Au events
177 double NNNpp = statscalepp * 12130./0.90; // Upsilons in p+p from sPHENIX proposal for 7350B p+p events
178 double ppauaurat = NNNpp/NNN;
179 
180 double frac[3]; // upsilons fractions
181  frac[0] = 0.7117;
182  frac[1] = 0.1851;
183  frac[2] = 0.1032;
184 double scale[3]; // mass scaling
185  scale[0] = 1.0;
186  scale[1] = 1.0595;
187  scale[2] = 1.0946;
188 //double supcor[3]; // suppression
189 // supcor[0]=0.535;
190 // supcor[1]=0.17;
191 // supcor[2]=0.035;
192 //int Nups1 = int(NNN*eideff*eideff*frac[0]*supcor[0]);
193 //int Nups2 = int(NNN*eideff*eideff*frac[1]*supcor[1]);
194 //int Nups3 = int(NNN*eideff*eideff*frac[2]*supcor[2]);
195 int Nups1nosup = int(NNN*eideff*eideff*frac[0]);
196 int Nups2nosup = int(NNN*eideff*eideff*frac[1]);
197 int Nups3nosup = int(NNN*eideff*eideff*frac[2]);
198 int Nups1pp = int(NNNpp*0.90*frac[0]); // always 95% eideff in p+p
199 int Nups2pp = int(NNNpp*0.90*frac[1]);
200 int Nups3pp = int(NNNpp*0.90*frac[2]);
201 
202 int nchan=400;
203 double start=0.0;
204 double stop=20.0;
205 
206 string str_UpsilonPt = "(2.0*3.14159*x*[0]*pow((1 + x*x/(4*[1]) ),-[2]))";
207 TF1* fUpsilonPt = new TF1("fUpsilonPt",str_UpsilonPt.c_str(),0.,20.);
208 fUpsilonPt->SetParameters(72.1, 26.516, 10.6834);
209 double upsnorm = fUpsilonPt->Integral(0.,20.);
210 
211 // count all upsilons with suppression
212 double tmpsum1=0.;
213 double tmpsum2=0.;
214 double tmpsum3=0.;
215 for(int j=0; j<nbins; j++) {
216  double s1 = j*1.0;
217  double s2 = s1 + 1.0;
218  tmpsum1 += Nups1nosup*fUpsilonPt->Integral(s1,s2)/upsnorm * grRAA1S->Eval((s1+s2)/2.);
219  tmpsum2 += Nups2nosup*fUpsilonPt->Integral(s1,s2)/upsnorm * grRAA2S->Eval((s1+s2)/2.);
220  tmpsum3 += Nups3nosup*fUpsilonPt->Integral(s1,s2)/upsnorm * grRAA3S->Eval((s1+s2)/2.);
221 }
222 int Nups1 = int(tmpsum1);
223 int Nups2 = int(tmpsum2);
224 int Nups3 = int(tmpsum3);
225 
226  cout << "Total number of ALL Upsilons in AuAu = " << NNN << endl;
227  cout << "Total number of ALL Upsilons in pp = " << NNNpp << endl;
228  cout << "Number of upsilons in pp in acceptance = " << NNNpp*frac[0] << " " << NNNpp*frac[1] << " " << NNNpp*frac[2] << endl;
229  cout << "Number of upsilons in AuAu in acceptance = " << NNN*frac[0] << " " << NNN*frac[1] << " " << NNN*frac[2] << endl;
230  cout << "Number of upsilons in AuAu after eID efficiency = " << Nups1nosup << " " << Nups2nosup << " " << Nups3nosup << endl;
231  cout << "Number of upsilons in AuAu plot = " << Nups1 << " " << Nups2 << " " << Nups3 << endl;
232  cout << "Number of upsilons in pp plot = " << Nups1pp << " " << Nups2pp << " " << Nups3pp << endl;
233 
234 //====================================================================================================
235 /*
236 f=new TFile("UpsilonsInHijing.root");
237  cout << endl << endl << "Number of entries in Upsilon NTuple = " << ntpups->GetEntries() << endl;
238 
239  TH1D* hforfit = new TH1D("hforfit","",nchan,start,stop);
240  TH1D* hUps1 = new TH1D("hUps1","",nchan,start,stop);
241  TH1D* hUps2 = new TH1D("hUps2","",nchan,start,stop);
242  TH1D* hUps3 = new TH1D("hUps3","",nchan,start,stop);
243  TH1D* hUps1pp = new TH1D("hUps1pp","",nchan,start,stop);
244  TH1D* hUps2pp = new TH1D("hUps2pp","",nchan,start,stop);
245  TH1D* hUps3pp = new TH1D("hUps3pp","",nchan,start,stop);
246  hforfit->Sumw2();
247  hUps1->Sumw2();
248  hUps2->Sumw2();
249  hUps3->Sumw2();
250  hUps1pp->Sumw2();
251  hUps2pp->Sumw2();
252  hUps3pp->Sumw2();
253 
254 //------ Additional smearing to reproduce mass resolution -----
255 //double smear = 0.098; // 127
256 //double smear = 0.090; // 119
257 //double smear = 0.067; // 100 MeV benchmark
258 double smear = 0.046; // 86 MeV resolution, Aalysis Note result
259 float mass;
260 TChain* tree = new TChain("ntpups");
261 tree->Add("UpsilonsInHijing.root");
262 tree->SetBranchAddress("mass",&mass);
263 int countups=0;
264  for(int j=0; j<tree->GetEntries(); j++){
265  tree->GetEvent(j);
266  double newmass = myrandom->Gaus(mass,smear);
267  hforfit->Fill(newmass);
268  if(countups<=Nups1) hUps1->Fill(newmass);
269  if(countups<=Nups2) hUps2->Fill(newmass*1.0595);
270  if(countups<=Nups3) hUps3->Fill(newmass*1.0946);
271  if(countups<=Nups1pp) hUps1pp->Fill(newmass);
272  if(countups<=Nups2pp) hUps2pp->Fill(newmass*1.0595);
273  if(countups<=Nups3pp) hUps3pp->Fill(newmass*1.0946);
274  countups++;
275  }
276 delete tree;
277 //---------------------------------------------------------------
278 
279  hforfit->SetDirectory(gROOT);
280  hUps1->SetDirectory(gROOT);
281  hUps2->SetDirectory(gROOT);
282  hUps3->SetDirectory(gROOT);
283  hUps1pp->SetDirectory(gROOT);
284  hUps2pp->SetDirectory(gROOT);
285  hUps3pp->SetDirectory(gROOT);
286 
287 f->Close();
288 
289 TH1F* hUps = (TH1F*)hUps1->Clone("hUps");
290 hUps->Add(hUps2);
291 hUps->Add(hUps3);
292 hUps->SetLineWidth(2);
293 TH1F* hUpspp = (TH1F*)hUps1pp->Clone("hUpspp");
294 hUpspp->Add(hUps2pp);
295 hUpspp->Add(hUps3pp);
296 hUpspp->SetLineWidth(2);
297 
298 TCanvas* cups = new TCanvas("cups","Upsilons (all pT)",100,100,1000,500);
299 cups->Divide(2,1);
300 
301  cups_1->cd();
302  cups_1->SetLogy();
303  fCB->SetParameter(0,5000.);
304  fCB->SetParameter(1,1.5);
305  fCB->SetParameter(2,1.2);
306  fCB->SetParameter(3,9.448);
307  fCB->SetParameter(4,0.070);
308  hforfit->Fit("fCB","rl","",5.,9.7);
309  hforfit->SetAxisRange(5.,11.);
310  hforfit->SetLineWidth(2);
311  hforfit->Draw();
312  double myupsmass = fCB->GetParameter(3);
313  double upswidth = fCB->GetParameter(4);
314 
315  cups_2->cd();
316  cups_2->SetLogy();
317  fTCB->SetParameter(0,2000.);
318  fTCB->SetParameter(1,1.5);
319  fTCB->SetParameter(2,1.25);
320  fTCB->SetParameter(3,9.448);
321  fTCB->FixParameter(4,upswidth); // fix width from the single peak fit
322  fTCB->SetParameter(5,500.);
323  fTCB->SetParameter(6,100.);
324  hUps->Fit(fTCB,"rl","",6.,11.);
325  double upswidth3 = fTCB->GetParameter(4);
326  double upswidth3err = fTCB->GetParError(4);
327  hUps->SetAxisRange(7.0,11.0);
328  hUps->Draw();
329 
330 cout << "MASS RESOLUTION = " << upswidth << " " << upswidth3 << " +- " << upswidth3err << endl;
331 */
332 //====================================================================================
333 //====================================================================================
334 //====================================================================================
335 
336 //
337 // these histograms (hhups*) are randomly generated using the fit above to single
338 // peak (fCB function) and used for the RAA uncertainty calculation
339 //
340 // FORCE specific RESOLUTION and tail shape
341  double tonypar1 = 0.98; // Tony's 100 pion simulation April 2019
342  double tonypar2 = 0.93; // Tony's 100 pion simulation April 2019
343  //double tonypar3 = 9.437; // Tony's 100 pion simulation April 2019
344  double tonypar3 = 9.448; // benchmark
345  double tonypar4 = 0.100; // benchmark
346  double tonypar4pp = 0.089; // Tony's 100 pion simulation April 2019
347  fCBpp->SetParameter(0,1000.);
348  fCBpp->SetParameter(1,tonypar1);
349  fCBpp->SetParameter(2,tonypar2);
350  fCBpp->SetParameter(3,tonypar3);
351  fCBpp->SetParameter(4,tonypar4pp);
352  fCBauau->SetParameter(0,1000.);
353  fCBauau->SetParameter(1,tonypar1);
354  fCBauau->SetParameter(2,tonypar2);
355  fCBauau->SetParameter(3,tonypar3);
356  fCBauau->SetParameter(4,tonypar4);
357 
358 
359 char hhname[999];
360 TH1D* hhups[nbins+1];
361 TH1D* hhups1[nbins+1];
362 TH1D* hhups2[nbins+1];
363 TH1D* hhups3[nbins+1];
364 TH1D* hhupspp[nbins+1];
365 TH1D* hhups1pp[nbins+1];
366 TH1D* hhups2pp[nbins+1];
367 TH1D* hhups3pp[nbins+1];
368 for(int i=0; i<nbins+1; i++) {
369  sprintf(hhname,"hhups_%d",i);
370  hhups[i] = new TH1D(hhname,"",nchan,start,stop);
371  hhups[i]->Sumw2();
372  sprintf(hhname,"hhups1_%d",i);
373  hhups1[i] = new TH1D(hhname,"",nchan,start,stop);
374  hhups1[i]->Sumw2();
375  sprintf(hhname,"hhups2_%d",i);
376  hhups2[i] = new TH1D(hhname,"",nchan,start,stop);
377  hhups2[i]->Sumw2();
378  sprintf(hhname,"hhups3_%d",i);
379  hhups3[i] = new TH1D(hhname,"",nchan,start,stop);
380  hhups3[i]->Sumw2();
381  hhups[i]->SetLineWidth(2);
382  hhups1[i]->SetLineWidth(2);
383  hhups2[i]->SetLineWidth(2);
384  hhups3[i]->SetLineWidth(2);
385  sprintf(hhname,"hhupspp_%d",i);
386  hhupspp[i] = new TH1D(hhname,"",nchan,start,stop);
387  hhupspp[i]->Sumw2();
388  sprintf(hhname,"hhups1pp_%d",i);
389  hhups1pp[i] = new TH1D(hhname,"",nchan,start,stop);
390  hhups1pp[i]->Sumw2();
391  sprintf(hhname,"hhups2pp_%d",i);
392  hhups2pp[i] = new TH1D(hhname,"",nchan,start,stop);
393  hhups2pp[i]->Sumw2();
394  sprintf(hhname,"hhups3pp_%d",i);
395  hhups3pp[i] = new TH1D(hhname,"",nchan,start,stop);
396  hhups3pp[i]->Sumw2();
397  hhupspp[i]->SetLineWidth(2);
398  hhups1pp[i]->SetLineWidth(2);
399  hhups2pp[i]->SetLineWidth(2);
400  hhups3pp[i]->SetLineWidth(2);
401 }
402 
403 for(int j=0; j<nbins; j++) {
404 
405  double s1 = j*1.0;
406  double s2 = s1 + 1.0;
407  double tmpnups1 = Nups1nosup*fUpsilonPt->Integral(s1,s2)/upsnorm * grRAA1S->Eval((s1+s2)/2.);
408  double tmpnups2 = Nups2nosup*fUpsilonPt->Integral(s1,s2)/upsnorm * grRAA2S->Eval((s1+s2)/2.);
409  double tmpnups3 = Nups3nosup*fUpsilonPt->Integral(s1,s2)/upsnorm * grRAA3S->Eval((s1+s2)/2.);
410  //cout << "Nups1["<<j<<"]="<<tmpnups1<<";"<<endl;
411  //cout << "Nups2["<<j<<"]="<<tmpnups2<<";"<<endl;
412  //cout << "Nups3["<<j<<"]="<<tmpnups3<<";"<<endl;
413  fCBauau->SetParameter(3,tonypar3);
414  for(int i=0; i<int(tmpnups1+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups1[j]->Fill(myrnd); hhups[j]->Fill(myrnd); }
415  fCBauau->SetParameter(3,tonypar3*scale[1]);
416  for(int i=0; i<int(tmpnups2+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups2[j]->Fill(myrnd); hhups[j]->Fill(myrnd); }
417  fCBauau->SetParameter(3,tonypar3*scale[2]);
418  for(int i=0; i<int(tmpnups3+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups3[j]->Fill(myrnd); hhups[j]->Fill(myrnd); }
419 
420  double tmpnups1pp = Nups1pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
421  double tmpnups2pp = Nups2pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
422  double tmpnups3pp = Nups3pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
423  cout << "Nups1pp["<<j<<"]="<<tmpnups1pp<<";"<<endl;
424  cout << "Nups2pp["<<j<<"]="<<tmpnups2pp<<";"<<endl;
425  cout << "Nups3pp["<<j<<"]="<<tmpnups3pp<<";"<<endl;
426  fCBpp->SetParameter(3,tonypar3);
427  for(int i=0; i<int(tmpnups1pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups1pp[j]->Fill(myrnd); hhupspp[j]->Fill(myrnd); }
428  fCBpp->SetParameter(3,tonypar3*scale[1]);
429  for(int i=0; i<int(tmpnups2pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups2pp[j]->Fill(myrnd); hhupspp[j]->Fill(myrnd); }
430  fCBpp->SetParameter(3,tonypar3*scale[2]);
431  for(int i=0; i<int(tmpnups3pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups3pp[j]->Fill(myrnd); hhupspp[j]->Fill(myrnd); }
432 
433 } // end loop over pT bins
434 
435 // all pT
436 
437  fCBpp->SetParameter(3,tonypar3);
438  fCBauau->SetParameter(3,tonypar3);
439  for(int i=0; i<int(Nups1+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups1[nbins]->Fill(myrnd); hhups[nbins]->Fill(myrnd); }
440  for(int i=0; i<int(Nups1pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups1pp[nbins]->Fill(myrnd); hhupspp[nbins]->Fill(myrnd); }
441  fCBpp->SetParameter(3,tonypar3*scale[1]);
442  fCBauau->SetParameter(3,tonypar3*scale[1]);
443  for(int i=0; i<int(Nups2+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups2[nbins]->Fill(myrnd); hhups[nbins]->Fill(myrnd); }
444  for(int i=0; i<int(Nups2pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups2pp[nbins]->Fill(myrnd); hhupspp[nbins]->Fill(myrnd); }
445  fCBpp->SetParameter(3,tonypar3*scale[2]);
446  fCBauau->SetParameter(3,tonypar3*scale[2]);
447  for(int i=0; i<int(Nups3+0.5); i++) { double myrnd = fCBauau->GetRandom(); hhups3[nbins]->Fill(myrnd); hhups[nbins]->Fill(myrnd); }
448  for(int i=0; i<int(Nups3pp+0.5); i++) { double myrnd = fCBpp->GetRandom(); hhups3pp[nbins]->Fill(myrnd); hhupspp[nbins]->Fill(myrnd); }
449 
450 // fit and draw pp no bg
451 
452 TCanvas* cupspp = new TCanvas("cupspp","Upsilons in p+p",100,100,600,600);
453  fTCBpp->SetParameter(0,2000.);
454  fTCBpp->FixParameter(1,tonypar1);
455  fTCBpp->FixParameter(2,tonypar2);
456  fTCBpp->SetParameter(3,tonypar3);
457  fTCBpp->FixParameter(4,tonypar4); // fix width from the single peak fit
458  fTCBpp->SetParameter(5,500.);
459  fTCBpp->SetParameter(6,100.);
460  hhupspp[nbins]->Fit(fTCBpp,"rl","",7.,11.);
461  hhupspp[nbins]->SetAxisRange(7.,11.);
462  hhupspp[nbins]->SetMarkerSize(1.0);
463  hhupspp[nbins]->GetXaxis()->SetTitle("Invariant mass [GeV/c^{2}]");
464  hhupspp[nbins]->GetXaxis()->SetTitleOffset(1.0);
465  double tmpamp1 = hhupspp[nbins]->GetFunction("fTCBpp")->GetParameter(0);
466  double tmpamp5 = tmpamp1*frac[1]/frac[0]; // correct fit for correct upsilon states ratio
467  double tmpamp6 = tmpamp1*frac[2]/frac[0];
468  hhupspp[nbins]->Draw();
469 
470  fCB1s->SetLineColor(kBlue);
471  fCB1s->SetLineWidth(1);
472  fCB1s->SetParameter(0,fTCBpp->GetParameter(0));
473  fCB1s->SetParameter(1,fTCBpp->GetParameter(1));
474  fCB1s->SetParameter(2,fTCBpp->GetParameter(2));
475  fCB1s->SetParameter(3,fTCBpp->GetParameter(3)*scale[0]);
476  fCB1s->SetParameter(4,fTCBpp->GetParameter(4));
477  fCB2s->SetLineColor(kRed);
478  fCB2s->SetLineWidth(1);
479  fCB2s->SetParameter(0,tmpamp5);
480  fCB2s->SetParameter(1,fTCBpp->GetParameter(1));
481  fCB2s->SetParameter(2,fTCBpp->GetParameter(2));
482  fCB2s->SetParameter(3,fTCBpp->GetParameter(3)*scale[1]);
483  fCB2s->SetParameter(4,fTCBpp->GetParameter(4));
484  fCB3s->SetLineColor(kGreen+2);
485  fCB3s->SetLineWidth(1);
486  fCB3s->SetParameter(0,tmpamp6);
487  fCB3s->SetParameter(1,fTCBpp->GetParameter(1));
488  fCB3s->SetParameter(2,fTCBpp->GetParameter(2));
489  fCB3s->SetParameter(3,fTCBpp->GetParameter(3)*scale[2]);
490  fCB3s->SetParameter(4,fTCBpp->GetParameter(4));
491  fCB1s->Draw("same");
492  fCB2s->Draw("same");
493  fCB3s->Draw("same");
494 
495 // fit and draw AuAu no bg
496 
497 TCanvas* cupsauau = new TCanvas("cupsauau","Upsilons in Au+Au",100,100,600,600);
498  fTCBauau->SetParameter(0,2000.);
499  fTCBauau->FixParameter(1,tonypar1);
500  fTCBauau->FixParameter(2,tonypar2);
501  fTCBauau->SetParameter(3,tonypar3);
502  fTCBauau->FixParameter(4,tonypar4); // fix width from the single peak fit
503  fTCBauau->SetParameter(5,500.);
504  fTCBauau->SetParameter(6,100.);
505  hhups[nbins]->Fit(fTCBauau,"rl","",7.,11.);
506  hhups[nbins]->SetAxisRange(7.,11.);
507  hhups[nbins]->SetMarkerSize(1.0);
508  hhups[nbins]->GetXaxis()->SetTitle("Invariant mass [GeV/c^{2}]");
509  hhups[nbins]->GetXaxis()->SetTitleOffset(1.0);
510  tmpamp1 = hhups[nbins]->GetFunction("fTCBauau")->GetParameter(0);
511  tmpamp5 = tmpamp1*frac[1]/frac[0]; // correct fit for correct upsilon states ratio
512  tmpamp6 = tmpamp1*frac[2]/frac[0];
513  hhups[nbins]->Draw();
514 
515 TCanvas* cupsvspt = new TCanvas("cupsvspt","Upsilons vs. p_{T}",100,100,1200,900);
516 cupsvspt->Divide(4,3);
517 for(int i=0; i<12; i++) {
518 if(i==0) {cupsvspt->cd(1);}
519 if(i==1) {cupsvspt->cd(2);}
520 if(i==2) {cupsvspt->cd(3);}
521 if(i==3) {cupsvspt->cd(4);}
522 if(i==4) {cupsvspt->cd(5);}
523 if(i==5) {cupsvspt->cd(6);}
524 if(i==6) {cupsvspt->cd(7);}
525 if(i==7) {cupsvspt->cd(8);}
526 if(i==8) {cupsvspt->cd(9);}
527 if(i==9) {cupsvspt->cd(10);}
528 if(i==10) {cupsvspt->cd(11);}
529 if(i==11) {cupsvspt->cd(12);}
530 hhups[i]->SetAxisRange(7.0,11.0); hhups[i]->Draw();
531  sprintf(tlchar,"%d-%d GeV",i,i+1); tl[i] = new TLatex(8.0,hhups[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
532 }
533 
534 //----------------------------------------------------
535 // Backgrounds
536 //----------------------------------------------------
537 
538 TH1D* hhall[nbins+1];
539 TH1D* hhall_scaled[nbins+1];
540 
541 TH1D* hhtotbg[nbins+1];
542 TH1D* hhtotbg_scaled[nbins+1];
543 TH1D* hhcombbg[nbins+1];
544 TH1D* hhcombbg_scaled[nbins+1];
545 TH1D* hhfakefake[nbins+1];
546 TH1D* hhfakehf[nbins+1];
547 TH1D* hhbottom[nbins+1];
548 TH1D* hhcharm[nbins+1];
549 TH1D* hhdy[nbins+1];
550 TH1D* hhcorrbg[nbins+1];
551 TH1D* hhcorrbg_scaled[nbins+1];
552 TH1D* hhfit[nbins+1];
553 char tmpname[999];
554 
555 //----------------------------------------------------------------------------------------
556 // Correlated BG calculated for 10B central AuAu events
557 //----------------------------------------------------------------------------------------
558 
559 double corrbgfitpar0;
560 double corrbgfitpar1;
561 
562 TFile* f=new TFile("ccbb_eideff09.root");
563 for(int i=0; i<nbins+1; i++) {
564  sprintf(tmpname,"hhbottom_%d",i);
565  hhbottom[i] = (TH1D*)f->Get(tmpname);
566  hhbottom[i]->SetDirectory(gROOT);
567  sprintf(tmpname,"hhcharm_%d",i);
568  hhcharm[i] = (TH1D*)f->Get(tmpname);
569  hhcharm[i]->SetDirectory(gROOT);
570  sprintf(tmpname,"hhdy_%d",i);
571  hhdy[i] = (TH1D*)f->Get(tmpname);
572  hhdy[i]->SetDirectory(gROOT);
573  sprintf(tmpname,"hhcorrbg_%d",i);
574  hhcorrbg[i] = (TH1D*)hhbottom[i]->Clone(tmpname);
575  hhcorrbg[i]->Add(hhcharm[i]);
576  hhcorrbg[i]->Add(hhdy[i]);
577  sprintf(tmpname,"hhcorrbg_scaled_%d",i);
578  hhcorrbg_scaled[i] = (TH1D*)hhcorrbg[i]->Clone(tmpname);
579  hhcorrbg[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
580  hhbottom[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
581  hhdy[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
582  if(i==nbins) {
583  corrbgfitpar0 = hhcorrbg[i]->GetFunction("expo")->GetParameter(0);
584  corrbgfitpar1 = hhcorrbg[i]->GetFunction("expo")->GetParameter(1);
585  }
586  cout << "bgpar0["<< i <<"]="<<hhcorrbg[i]->GetFunction("expo")->GetParameter(0)+TMath::Log(statscale)<<";"<< endl;
587  cout << "bgpar1["<< i <<"]="<<hhcorrbg[i]->GetFunction("expo")->GetParameter(1)<<";"<< endl;
588  for(int k=1; k<=hhcorrbg[i]->GetNbinsX(); k++) {
589  if(hhcorrbg[i]->GetBinLowEdge(k)<statscale_lowlim || (hhcorrbg[i]->GetBinLowEdge(k)+hhcorrbg[i]->GetBinWidth(k))>statscale_uplim) {
590  hhcorrbg_scaled[i]->SetBinContent(k,0.);
591  hhcorrbg_scaled[i]->SetBinError(k,0.);
592  }
593  else {
594  double tmp = statscale * hhcorrbg[i]->GetFunction("expo")->Eval(hhcorrbg[i]->GetBinCenter(k));
595  double tmprnd = myrandom->Poisson(tmp);
596  if(tmprnd<0.) { tmprnd=0.; }
597  hhcorrbg_scaled[i]->SetBinContent(k,tmprnd);
598  hhcorrbg_scaled[i]->SetBinError(k,sqrt(tmprnd));
599  }
600  }
601  hhcorrbg_scaled[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
602  hhcorrbg[i]->SetDirectory(gROOT);
603  hhcorrbg_scaled[i]->SetDirectory(gROOT);
604 }
605 f->Close();
606 
607 TCanvas* c2 = new TCanvas("c2","Correlated BG (all p_{T})",100,100,600,600);
608 
609 hhbottom[nbins]->SetAxisRange(7.0,14.0);
610 hhbottom[nbins]->SetLineColor(kBlue);
611 hhbottom[nbins]->SetLineWidth(2);
612 
613 hhbottom[nbins]->GetXaxis()->SetTitle("Invariant mass [GeV/c^{2}]");
614 hhbottom[nbins]->GetXaxis()->SetTitleOffset(1.0);
615 hhbottom[nbins]->GetXaxis()->SetTitleColor(1);
616 hhbottom[nbins]->GetXaxis()->SetTitleSize(0.040);
617 hhbottom[nbins]->GetXaxis()->SetLabelSize(0.040);
618 //hhbottom[nbins]->GetYaxis()->SetTitle("Correlated background");
619 hhbottom[nbins]->GetYaxis()->SetTitleOffset(1.3);
620 hhbottom[nbins]->GetYaxis()->SetTitleSize(0.040);
621 hhbottom[nbins]->GetYaxis()->SetLabelSize(0.040);
622 
623 hhcorrbg[nbins]->SetAxisRange(7.0,14.0);
624 hhcorrbg[nbins]->SetMinimum(0.1);
625 hhcorrbg[nbins]->SetLineColor(kBlack);
626 hhcorrbg[nbins]->SetLineWidth(2);
627 hhcorrbg[nbins]->Draw("hist");
628 
629 hhbottom[nbins]->SetMarkerColor(kBlue);
630 hhbottom[nbins]->SetLineColor(kBlue);
631 hhbottom[nbins]->Draw("same");
632 
633 hhdy[nbins]->SetMarkerColor(kGreen+2);
634 hhdy[nbins]->SetLineColor(kGreen+2);
635 hhdy[nbins]->SetLineWidth(2);
636 hhdy[nbins]->Draw("same");
637 
638 hhcharm[nbins]->SetMarkerColor(kRed);
639 hhcharm[nbins]->SetLineColor(kRed);
640 hhcharm[nbins]->SetLineWidth(2);
641 hhcharm[nbins]->Draw("same");
642 
643 TCanvas* c0 = new TCanvas("c0","Correlated BG vs. p_{T} 10B events",100,100,1200,900);
644 c0->Divide(4,3);
645 for(int i=0; i<nbins; i++) {
646 if(i==0) {c0->cd(1);}
647 if(i==1) {c0->cd(2);}
648 if(i==2) {c0->cd(3);}
649 if(i==3) {c0->cd(4);}
650 if(i==4) {c0->cd(5);}
651 if(i==5) {c0->cd(6);}
652 if(i==6) {c0->cd(7);}
653 if(i==7) {c0->cd(8);}
654 if(i==8) {c0->cd(9);}
655 if(i==9) {c0->cd(10);}
656 if(i==10) {c0->cd(11);}
657 if(i==11) {c0->cd(12);}
658  hhcorrbg[i]->SetAxisRange(7.,14.);
659  hhcorrbg[i]->GetFunction("expo")->SetLineColor(kBlack);
660  hhbottom[i]->GetFunction("expo")->SetLineColor(kBlue);
661  hhdy[i]->GetFunction("expo")->SetLineColor(kGreen+2);
662  hhcorrbg[i]->SetLineColor(kBlack);
663  hhbottom[i]->SetLineColor(kBlue);
664  hhdy[i]->SetLineColor(kGreen+2);
665  hhcorrbg[i]->Draw();
666  hhbottom[i]->Draw("same");
667  hhdy[i]->Draw("same");
668  sprintf(tlchar,"%d-%d GeV",i,i+1); tl[i] = new TLatex(8.0,hhcorrbg[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
669 }
670 
671 TCanvas* c0scaled = new TCanvas("c0scaled","SCALED Correlated BG vs. p_{T}",100,100,1200,900);
672 c0scaled->Divide(4,3);
673 for(int i=0; i<nbins; i++) {
674 if(i==0) {c0scaled->cd(1);}
675 if(i==1) {c0scaled->cd(2);}
676 if(i==2) {c0scaled->cd(3);}
677 if(i==3) {c0scaled->cd(4);}
678 if(i==4) {c0scaled->cd(5);}
679 if(i==5) {c0scaled->cd(6);}
680 if(i==6) {c0scaled->cd(7);}
681 if(i==7) {c0scaled->cd(8);}
682 if(i==8) {c0scaled->cd(9);}
683 if(i==9) {c0scaled->cd(10);}
684 if(i==10) {c0scaled->cd(11);}
685 if(i==11) {c0scaled->cd(12);}
686 hhcorrbg_scaled[i]->SetAxisRange(7.0,14.0); hhcorrbg_scaled[i]->Draw();
687  sprintf(tlchar,"%d-%d GeV",i,i+1); tl[i] = new TLatex(8.0,hhcorrbg_scaled[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
688 }
689 
690 //-----------------------------------------------------------------------------------------
691 // Correlated bg in p+p
692 //-----------------------------------------------------------------------------------------
693 
694 TH1D* hhbottom_pp[nbins+1];
695 TH1D* hhdy_pp[nbins+1];
696 TH1D* hhcorrbg_pp[nbins+1];
697 TH1D* hhall_pp[nbins+1];
698 
699 double ppcorr = 8300./10./962.;
700 TF1* fbottom_nosup_corr = new TF1("fbottom_nosup_corr","[0]+[1]*x",5.,14.);
701 fbottom_nosup_corr->SetParameters(-2.13861, 0.683323);
702 
703 for(int i=0; i<nbins+1; i++) {
704 
705  sprintf(tmpname,"hhbottom_pp_%d",i);
706  hhbottom_pp[i] = (TH1D*)hhbottom[i]->Clone(tmpname);
707  for(int k=1; k<=hhbottom_pp[i]->GetNbinsX(); k++) {
708  if(hhbottom_pp[i]->GetBinLowEdge(k)<statscale_lowlim || (hhbottom_pp[i]->GetBinLowEdge(k)+hhbottom_pp[i]->GetBinWidth(k))>statscale_uplim) {
709  hhbottom_pp[i]->SetBinContent(k,0.);
710  hhbottom_pp[i]->SetBinError(k,0.);
711  }
712  else {
713  double tmp = ppcorr * fbottom_nosup_corr->Eval(hhbottom[i]->GetBinCenter(k)) * hhbottom[i]->GetFunction("expo")->Eval(hhbottom[i]->GetBinCenter(k));
714  double tmprnd = myrandom->Poisson(tmp);
715  if(tmprnd<0.) { tmprnd=0.; }
716  hhbottom_pp[i]->SetBinContent(k,tmprnd);
717  hhbottom_pp[i]->SetBinError(k,sqrt(tmprnd));
718  }
719  }
720 
721  sprintf(tmpname,"hhdy_pp_%d",i);
722  hhdy_pp[i] = (TH1D*)hhdy[i]->Clone(tmpname);
723  for(int k=1; k<=hhdy_pp[i]->GetNbinsX(); k++) {
724  if(hhdy_pp[i]->GetBinLowEdge(k)<statscale_lowlim || (hhdy_pp[i]->GetBinLowEdge(k)+hhdy_pp[i]->GetBinWidth(k))>statscale_uplim) {
725  hhdy_pp[i]->SetBinContent(k,0.);
726  hhdy_pp[i]->SetBinError(k,0.);
727  }
728  else {
729  double tmp = ppcorr * hhdy[i]->GetFunction("expo")->Eval(hhdy[i]->GetBinCenter(k));
730  double tmprnd = myrandom->Poisson(tmp);
731  if(tmprnd<0.) { tmprnd=0.; }
732  hhdy_pp[i]->SetBinContent(k,tmprnd);
733  hhdy_pp[i]->SetBinError(k,sqrt(tmprnd));
734  }
735  }
736 
737  sprintf(tmpname,"hhcorrbg_pp_%d",i);
738  hhcorrbg_pp[i] = (TH1D*)hhbottom_pp[i]->Clone(tmpname);
739  hhcorrbg_pp[i]->Add(hhdy_pp[i]);
740  hhcorrbg_pp[i]->SetMarkerColor(kBlack);
741  hhcorrbg_pp[i]->SetLineColor(kBlack);
742  hhbottom_pp[i]->SetLineColor(kBlue);
743  hhdy_pp[i]->SetLineColor(kGreen+2);
744  sprintf(tmpname,"hhall_pp_%d",i);
745  hhall_pp[i] = (TH1D*)hhcorrbg_pp[i]->Clone(tmpname);
746  hhall_pp[i]->Add(hhupspp[i]);
747  hhall_pp[i]->SetLineColor(kMagenta);
748  hhall_pp[i]->SetMarkerColor(kMagenta);
749 
750  hhcorrbg_pp[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
751  hhcorrbg_pp[i]->GetFunction("expo")->SetLineColor(kBlack);
752  hhbottom_pp[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
753  hhbottom_pp[i]->GetFunction("expo")->SetLineColor(kBlue);
754  hhdy_pp[i]->Fit("expo","rql","",statscale_lowlim,statscale_uplim);
755  hhdy_pp[i]->GetFunction("expo")->SetLineColor(kGreen+2);
756 
757 } // end loop over pT bins
758 
759 TCanvas* cbginpp = new TCanvas("cbginpp","corr bg in pp",10,10,700,700);
760  hhall_pp[nbins]->SetAxisRange(7.,12.);
761  hhcorrbg_pp[nbins]->Draw("pehist");
762  hhbottom_pp[nbins]->Draw("histsame");
763  hhdy_pp[nbins]->Draw("histsame");
764 
765 
766 TCanvas* cpp = new TCanvas("cpp","corr bg + sig in pp",100,100,700,700);
767  hhall_pp[nbins]->SetAxisRange(7.,12.);
768  hhall_pp[nbins]->Draw("pehist");
769  hhcorrbg_pp[nbins]->Draw("pesame");
770  hhbottom_pp[nbins]->Draw("same");
771  hhdy_pp[nbins]->Draw("same");
772 
773 TCanvas* cpp_vspt = new TCanvas("cpp_vspt","corr bg + sig vs pt in pp",50,50,1200,900);
774 cpp_vspt->Divide(4,3);
775 for(int i=0; i<nbins; i++) {
776 if(i==0) {cpp_vspt->cd(1);}
777 if(i==1) {cpp_vspt->cd(2);}
778 if(i==2) {cpp_vspt->cd(3);}
779 if(i==3) {cpp_vspt->cd(4);}
780 if(i==4) {cpp_vspt->cd(5);}
781 if(i==5) {cpp_vspt->cd(6);}
782 if(i==6) {cpp_vspt->cd(7);}
783 if(i==7) {cpp_vspt->cd(8);}
784 if(i==8) {cpp_vspt->cd(9);}
785 if(i==9) {cpp_vspt->cd(10);}
786 if(i==10) {cpp_vspt->cd(11);}
787 if(i==11) {cpp_vspt->cd(12);}
788 hhall_pp[i]->SetAxisRange(7.0,12.0); hhall_pp[i]->Draw("hist"); hhcorrbg_pp[i]->Draw("histsame");
789  sprintf(tlchar,"%d-%d GeV",i,i+1); tl[i] = new TLatex(8.0,hhcorrbg_pp[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
790 }
791 
792 //-----------------------------------------------------------------------------------------
793 // Combinatorial BG calculated for 10B central AuAu events
794 //-----------------------------------------------------------------------------------------
795 TCanvas* dummy = new TCanvas("dummy","dummy",0,0,500,500);
796 
797 f = new TFile("fakee_eideff09.root");
798 for(int i=0; i<nbins+1; i++) {
799  sprintf(tmpname,"hhfakefake_%d",i);
800  hhfakefake[i] = (TH1D*)f->Get(tmpname);
801  hhfakefake[i]->SetDirectory(gROOT);
802 }
803 f->Close();
804 
805 f = new TFile("crossterms_eideff09.root");
806 for(int i=0; i<nbins+1; i++) {
807  sprintf(tmpname,"hhfakehf_%d",i);
808  hhfakehf[i] = (TH1D*)f->Get(tmpname);
809  hhfakehf[i]->SetDirectory(gROOT);
810 }
811 f->Close();
812 
813 TF1* fbg = new TF1("fbg","exp([0]+[1]*x)+exp([2]+[3]*x)",8.,11.);
814 fbg->SetParameters(10., -1.0, 4., -0.1);
815 fbg->SetParLimits(1.,-999.,0.);
816 fbg->SetParLimits(3.,-999.,0.);
817 
818 for(int i=0; i<nbins+1; i++) {
819  sprintf(tmpname,"hhcombbg_%d",i);
820  hhcombbg[i] = (TH1D*)hhfakefake[i]->Clone(tmpname);
821  hhcombbg[i]->Add(hhfakehf[i]);
822  sprintf(tmpname,"hhcombbg_scaled_%d",i);
823  hhcombbg_scaled[i] = (TH1D*)hhcombbg[i]->Clone(tmpname);
824  if(i==nbins) { fbg->SetParameters(10., -1.0, 4., -0.1); }
825  hhcombbg[i]->Fit(fbg,"qrl","",statscale_lowlim,statscale_uplim);
826 
827  for(int k=1; k<=hhcombbg[i]->GetNbinsX(); k++) {
828  if(hhcombbg[i]->GetBinLowEdge(k)<statscale_lowlim || (hhcombbg[i]->GetBinLowEdge(k)+hhcombbg[i]->GetBinWidth(k))>statscale_uplim) {
829  hhcombbg_scaled[i]->SetBinContent(k,0.);
830  hhcombbg_scaled[i]->SetBinError(k,0.);
831  }
832  else {
833  double tmp = statscale * hhcombbg[i]->GetFunction("fbg")->Eval(hhcombbg[i]->GetBinCenter(k));
834  double tmprnd = myrandom->Poisson(tmp);
835  if(tmprnd<0.) { tmprnd=0.; }
836  hhcombbg_scaled[i]->SetBinContent(k,tmprnd);
837  hhcombbg_scaled[i]->SetBinError(k,sqrt(tmprnd));
838  }
839  }
840  hhcombbg_scaled[i]->Fit(fbg,"qrl","",statscale_lowlim,statscale_uplim);
841 }
842 
843 delete dummy;
844 
845 TCanvas* C1 = new TCanvas("C1","Combinatorial BG (ALL p_{T})",100,100,600,600);
846 C1->SetLogy();
847  hhfakefake[nbins]->SetAxisRange(7.0,14.0);
848  hhfakefake[nbins]->SetMinimum(0.1);
849  hhfakefake[nbins]->SetMaximum(5000.);
850  hhfakefake[nbins]->SetLineColor(kGreen+2);
851  hhfakefake[nbins]->SetLineWidth(2);
852  hhfakefake[nbins]->GetXaxis()->SetTitle("Transverse momentum [GeV/c]");
853  hhfakefake[nbins]->GetXaxis()->SetTitleOffset(1.0);
854  hhfakefake[nbins]->GetXaxis()->SetTitleColor(1);
855  hhfakefake[nbins]->GetXaxis()->SetTitleSize(0.040);
856  hhfakefake[nbins]->GetXaxis()->SetLabelSize(0.040);
857  hhfakefake[nbins]->GetYaxis()->SetTitle("Combinatorial background");
858  hhfakefake[nbins]->GetYaxis()->SetTitleOffset(1.3);
859  hhfakefake[nbins]->GetYaxis()->SetTitleSize(0.040);
860  hhfakefake[nbins]->GetYaxis()->SetLabelSize(0.040);
861  hhfakefake[nbins]->Draw("e");
862 
863  hhfakehf[nbins]->SetLineColor(kOrange+4);
864  hhfakehf[nbins]->SetLineWidth(2);
865  hhfakehf[nbins]->Draw("esame");
866 
867  hhcombbg[nbins]->SetLineColor(kBlack);
868  hhcombbg[nbins]->SetLineWidth(2);
869  hhcombbg[nbins]->Draw("esame");
870 
871 TCanvas* C1sc = new TCanvas("C1sc","SCALED Combinatorial BG (ALL p_{T})",100,100,600,600);
872 C1sc->SetLogy();
873  hhcombbg_scaled[nbins]->SetAxisRange(7.,14.);
874  hhcombbg_scaled[nbins]->Draw("esame");
875 
876 TCanvas* c00 = new TCanvas("c00","Combinatorial BG vs. p_{T}",150,150,1200,900);
877 c00->Divide(4,3);
878 
879 for(int i=0; i<nbins; i++) {
880 if(i==0) {c00->cd(1);}
881 if(i==1) {c00->cd(2);}
882 if(i==2) {c00->cd(3);}
883 if(i==3) {c00->cd(4);}
884 if(i==4) {c00->cd(5);}
885 if(i==5) {c00->cd(6);}
886 if(i==6) {c00->cd(7);}
887 if(i==7) {c00->cd(8);}
888 if(i==8) {c00->cd(9);}
889 if(i==9) {c00->cd(10);}
890 if(i==10) {c00->cd(11);}
891 if(i>=11) {c00->cd(12);}
892 hhcombbg[i]->SetAxisRange(7.0,14.0); hhcombbg[i]->Draw();
893  sprintf(tlchar,"%d-%d GeV",i,i+1); tl[i] = new TLatex(8.0,hhcombbg[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
894 }
895 
896 TCanvas* c00scaled = new TCanvas("c00scaled","SCALED Combinatorial BG vs. p_{T}",150,150,1200,900);
897 c00scaled->Divide(4,3);
898 
899 for(int i=0; i<nbins; i++) {
900 if(i==0) {c00scaled->cd(1);}
901 if(i==1) {c00scaled->cd(2);}
902 if(i==2) {c00scaled->cd(3);}
903 if(i==3) {c00scaled->cd(4);}
904 if(i==4) {c00scaled->cd(5);}
905 if(i==5) {c00scaled->cd(6);}
906 if(i==6) {c00scaled->cd(7);}
907 if(i==7) {c00scaled->cd(8);}
908 if(i==8) {c00scaled->cd(9);}
909 if(i==9) {c00scaled->cd(10);}
910 if(i==10) {c00scaled->cd(11);}
911 if(i>=11) {c00scaled->cd(12);}
912 hhcombbg_scaled[i]->SetAxisRange(statscale_lowlim,statscale_uplim); hhcombbg_scaled[i]->Draw();
913  sprintf(tlchar,"%d-%d GeV",i,i+1); tl[i] = new TLatex(8.0,hhcombbg_scaled[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
914 }
915 
916 //---------------------------------------------
917 // Sibnal + BG
918 //---------------------------------------------
919 
920 for(int i=0; i<nbins+1; i++) {
921  sprintf(tmpname,"hhtotbg_scaled_%d",i);
922  hhtotbg_scaled[i] = (TH1D*)hhcombbg_scaled[i]->Clone(tmpname);
923  hhtotbg_scaled[i]->Add(hhcorrbg_scaled[i]);
924 }
925 
926 for(int i=0; i<nbins+1; i++) {
927 sprintf(tmpname,"hhall_scaled_%d",i);
928 hhall_scaled[i] = (TH1D*)hhtotbg_scaled[i]->Clone(tmpname);
929 hhall_scaled[i]->Add(hhups[i]);
930 }
931 
932 TCanvas* c000 = new TCanvas("c000","Signal + Background vs. p_{T}",200,200,1200,900);
933 c000->Divide(4,3);
934 for(int i=0; i<nbins; i++) {
935 if(i==0) {c000->cd(1);}
936 if(i==1) {c000->cd(2);}
937 if(i==2) {c000->cd(3);}
938 if(i==3) {c000->cd(4);}
939 if(i==4) {c000->cd(5);}
940 if(i==5) {c000->cd(6);}
941 if(i==6) {c000->cd(7);}
942 if(i==7) {c000->cd(8);}
943 if(i==8) {c000->cd(9);}
944 if(i==9) {c000->cd(10);}
945 if(i==10) {c000->cd(11);}
946 if(i==11) {c000->cd(12);}
947 hhall_scaled[i]->SetAxisRange(7.0,14.0); hhall_scaled[i]->SetMarkerStyle(1); hhall_scaled[i]->Draw("pehist");
948  sprintf(tlchar,"%d-%d GeV",i,i+1); tl[i] = new TLatex(8.0,hhall_scaled[i]->GetMaximum()*0.9,tlchar); tl[i]->Draw();
949 }
950 
951 // Fit Signal + Correlated BG
952 
953 // hhfit[] is signal + correlated bg
954 for(int i=0; i<nbins+1; i++) {
955  sprintf(tmpname,"hhfit_%d",i);
956  hhfit[i] = (TH1D*)hhall_scaled[i]->Clone(tmpname);
957  for(int j=1; j<=hhall_scaled[i]->GetNbinsX(); j++) {
958  hhfit[i]->SetBinContent(j,hhfit[i]->GetBinContent(j) - hhcombbg_scaled[i]->GetFunction("fbg")->Eval(hhcombbg_scaled[i]->GetBinCenter(j)));
959  hhfit[i]->SetBinError(j,sqrt(pow(hhfit[i]->GetBinError(j),2)+hhcombbg_scaled[i]->GetFunction("fbg")->Eval(hhcombbg_scaled[i]->GetBinCenter(j))));
960  }
961 }
962 
963 //---------------------------------------------------------------------
964 // plot signal + correlated bg for all pT
965 //---------------------------------------------------------------------
966 
967  double tmppar0 = corrbgfitpar0+TMath::Log(statscale);
968  double tmppar1 = corrbgfitpar1;
969 cout << "###### " << tmppar0 << " " << tmppar1 << endl;
970 
971 double ppauauscale = fTCBauau->GetParameter(0)/fTCBpp->GetParameter(0)*0.9783;
972 fSandBpp->SetParameter(0,fTCBpp->GetParameter(0)*ppauauscale);
973 fSandBpp->SetParameter(1,fTCBpp->GetParameter(1));
974 fSandBpp->SetParameter(2,fTCBpp->GetParameter(2));
975 fSandBpp->SetParameter(3,fTCBpp->GetParameter(3));
976 fSandBpp->SetParameter(4,fTCBpp->GetParameter(4));
977 fSandBpp->SetParameter(5,fTCBpp->GetParameter(5)*ppauauscale);
978 fSandBpp->SetParameter(6,fTCBpp->GetParameter(6)*ppauauscale);
979 fSandBpp->SetParameter(7,tmppar0);
980 fSandBpp->SetParameter(8,tmppar1);
981 fSandBpp->SetLineColor(kBlue);
982 fSandBpp->SetLineStyle(2);
983 
984 cout << "######### " << fTCBauau->GetParameter(0) << " " << fTCBauau->GetParameter(5) << " " << fTCBauau->GetParameter(6) << endl;
985 fSandBauau->SetParameter(0,fTCBauau->GetParameter(0));
986 fSandBauau->SetParameter(1,fTCBauau->GetParameter(1));
987 fSandBauau->SetParameter(2,fTCBauau->GetParameter(2));
988 fSandBauau->SetParameter(3,fTCBauau->GetParameter(3));
989 fSandBauau->SetParameter(4,fTCBauau->GetParameter(4));
990 fSandBauau->SetParameter(5,fTCBauau->GetParameter(5));
991 fSandBauau->SetParameter(6,fTCBauau->GetParameter(6));
992 fSandBauau->SetParameter(7,tmppar0);
993 fSandBauau->SetParameter(8,tmppar1);
994 fSandBauau->SetLineColor(kRed);
995 
996 
997 TCanvas* cfitall = new TCanvas("cfitall","FIT all pT",270,270,600,600);
998  hhfit[nbins]->SetAxisRange(7.0,14.);
999  hhfit[nbins]->GetXaxis()->CenterTitle();
1000  hhfit[nbins]->GetXaxis()->SetTitle("Mass(e^{+}e^{-}) [GeV/c^2]");
1001  hhfit[nbins]->GetXaxis()->SetTitleOffset(1.1);
1002  hhfit[nbins]->GetXaxis()->SetLabelSize(0.045);
1003  hhfit[nbins]->GetYaxis()->CenterTitle();
1004  hhfit[nbins]->GetYaxis()->SetLabelSize(0.045);
1005  hhfit[nbins]->GetYaxis()->SetTitle("Events / (50 MeV/c^{2})");
1006  hhfit[nbins]->GetYaxis()->SetTitleOffset(1.5);
1007  hhfit[nbins]->Draw("pehist");
1008  fSandBpp->Draw("same");
1009  fSandBauau->Draw("same");
1010  //hhcorrbg_scaled[nbins]->SetLineColor(kRed);
1011  //hhcorrbg_scaled[nbins]->Scale(0.82);
1012  //hhcorrbg_scaled[nbins]->Scale(0.70);
1013  //hhcorrbg_scaled[nbins]->Draw("histsame");
1014 
1015  TF1* fmycorrbg = new TF1("fmycorrbg","exp([0]+[1]*x)",7.,14.);
1016  fmycorrbg->SetParameters(tmppar0,tmppar1);
1017  fmycorrbg->SetLineStyle(2);
1018  fmycorrbg->SetLineColor(kRed);
1019  fmycorrbg->Draw("same");
1020 
1021 
1022 double myheight = fTCBauau->GetParameter(0);
1023 TLatex* ld1 = new TLatex(10.1,myheight,"sPHENIX Simulation");
1024 ld1->SetTextSize(0.035);
1025 ld1->Draw();
1026 TLatex* ld2 = new TLatex(10.1,myheight-100.,"0-10% Au+Au #sqrt{s} = 200 GeV");
1027 ld2->SetTextSize(0.035);
1028 ld2->Draw();
1029 TLatex* ld3;
1030 if(statscale==2.4) {
1031 ld3 = new TLatex(10.1,myheight-200.,"24 billion events");
1032 }
1033 else {
1034 ld3 = new TLatex(10.1,myheight-200.,"14.3 billion events");
1035 }
1036 ld3->SetTextSize(0.035);
1037 ld3->Draw();
1038 
1039 TCanvas* cfitall2 = new TCanvas("cfitall2","FIT all pT",270,270,600,600);
1040 TH1D* hhfit_tmp = (TH1D*)hhfit[nbins]->Clone("hhfit_tmp");
1041 hhfit_tmp->SetAxisRange(8.0,11.);
1042 hhfit_tmp->Draw("pehist");
1043 fSandBauau->Draw("same");
1044 fmycorrbg->Draw("same");
1045 
1046 
1047 //----------------------------------------------------------------------
1048 
1049 // plot signal + both backgrounds for all pT
1050 
1051 TCanvas* callpt = new TCanvas("callpt","Signal+BG (all p_{T})",300,300,600,600);
1052 
1053  hhall_scaled[nbins]->GetXaxis()->SetTitle("Invariant mass GeV/c");
1054  hhall_scaled[nbins]->SetLineColor(kBlack);
1055  hhall_scaled[nbins]->SetMarkerColor(kBlack);
1056  hhall_scaled[nbins]->SetMarkerStyle(20);
1057  hhall_scaled[nbins]->SetAxisRange(8.0,10.8);
1058  hhall_scaled[nbins]->Draw("pehist");
1059  hhcombbg_scaled[nbins]->SetLineColor(kBlue);
1060  hhcombbg_scaled[nbins]->Draw("histsame");
1061  hhcorrbg_scaled[nbins]->SetLineColor(kRed);
1062  hhcorrbg_scaled[nbins]->Draw("histsame");
1063 
1064 //----------------------------------------------------------------
1065 // Calculate RAA
1066 //----------------------------------------------------------------
1067 
1068 double u1start = 9.25;
1069 double u1stop = 9.65;
1070 double u2start = 9.80;
1071 double u2stop = 10.20;
1072 double u3start = 10.20;
1073 double u3stop = 10.55;
1074 
1075  double raa1[nbins+1],raa2[nbins+1],raa3[nbins+1],erraa1[nbins+1],erraa2[nbins+1],erraa3[nbins+1];
1076  for(int i=0; i<nbins; i++) {
1077  //raa1[i] = supcor[0]; raa2[i] = supcor[1]; raa3[i] = supcor[2];
1078  raa1[i] = grRAA1S->Eval(0.5+i*1.0);
1079  raa2[i] = grRAA2S->Eval(0.5+i*1.0);
1080  raa3[i] = grRAA3S->Eval(0.5+i*1.0);
1081  }
1082  int fbin1 = hhall_scaled[nbins]->FindBin(u1start + 0.001);
1083  int lbin1 = hhall_scaled[nbins]->FindBin(u1stop - 0.001);
1084  int fbin2 = hhall_scaled[nbins]->FindBin(u2start + 0.001);
1085  int lbin2 = hhall_scaled[nbins]->FindBin(u2stop - 0.001);
1086  int fbin3 = hhall_scaled[nbins]->FindBin(u3start + 0.001);
1087  int lbin3 = hhall_scaled[nbins]->FindBin(u3stop - 0.001);
1088  cout << "Y(1S) bin range: " << fbin1 << " - " << lbin1 << endl;
1089  cout << "Y(1S) inv. mass range: " << u1start << " - " << u1stop << endl;
1090  cout << "Y(2S) bin range: " << fbin2 << " - " << lbin2 << endl;
1091  cout << "Y(2S) inv. mass range: " << u2start << " - " << u2stop << endl;
1092  cout << "Y(3S) bin range: " << fbin3 << " - " << lbin3 << endl;
1093  cout << "Y(3S) inv. mass range: " << u3start << " - " << u3stop << endl;
1094 
1095  double sum1[99] = {0.};
1096  double truesum1[99] = {0.};
1097  double ersum1[99] = {0.};
1098  double sumpp1[99] = {0.};
1099  double ersumpp1[99] = {0.};
1100  double sum2[99] = {0.};
1101  double truesum2[99] = {0.};
1102  double ersum2[99] = {0.};
1103  double sumpp2[99] = {0.};
1104  double ersumpp2[99] = {0.};
1105  double sum3[99] = {0.};
1106  double truesum3[99] = {0.};
1107  double ersum3[99] = {0.};
1108  double sumpp3[99] = {0.};
1109  double ersumpp3[99] = {0.};
1110 
1111  double sumsum1[99] = {0.};
1112  double sumsum2[99] = {0.};
1113  double sumsum3[99] = {0.};
1114  double sumsum1pp[99] = {0.};
1115  double sumsum2pp[99] = {0.};
1116  double sumsum3pp[99] = {0.};
1117 
1118  for(int i=0; i<nbins+1; i++) {
1119 
1120  double s1 = double(i); double s2 = double(i+1); if(i==nbins) {s1 = 0.;}
1121 
1122  for(int j=fbin1; j<=lbin1; j++) {
1123  sum1[i] += (hhall_scaled[i]->GetBinContent(j) - hhcombbg_scaled[i]->GetFunction("fbg")->Eval(hhall_scaled[i]->GetBinCenter(j)) - hhcorrbg_scaled[i]->GetFunction("expo")->Eval(hhall_scaled[i]->GetBinCenter(j)));
1124  truesum1[i] += hhups1[i]->GetBinContent(j);
1125  ersum1[i] += hhall_scaled[i]->GetBinError(j)*hhall_scaled[i]->GetBinError(j);
1126  sumpp1[i] += hhups1pp[i]->GetBinContent(j);
1127  ersumpp1[i] += hhupspp[i]->GetBinError(j)*hhupspp[i]->GetBinError(j);
1128  }
1129  sumsum1[i] = truesum1[i]; // direct count in mass range
1130  sumsum1pp[i] = sumpp1[i];
1131  //sumsum1[i] = hhups1[i]->GetEntries(); // total number of upsilons in pT bin (rounded up)
1132  //sumsum1pp[i] = hhups1pp[i]->GetEntries();
1133  //sumsum1[i] = Nups1*fUpsilonPt->Integral(s1,s2)/upsnorm; // total number of upsilons in pT bin
1134  //sumsum1pp[i] = Nups1pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
1135 
1136  if(sumsum1[i]>0. && sumsum1pp[i]>0.) {
1137  erraa1[i] = raa1[i]*sqrt(ersum1[i]/sumsum1[i]/sumsum1[i] + ersumpp1[i]/sumsum1pp[i]/sumsum1pp[i]);
1138  } else {raa1[i]=-1.0; erraa1[i] = 999.; }
1139 
1140  for(int j=fbin2; j<=lbin2; j++) {
1141  sum2[i] += (hhall_scaled[i]->GetBinContent(j) - hhcombbg_scaled[i]->GetFunction("fbg")->Eval(hhall_scaled[i]->GetBinCenter(j)) - hhcorrbg_scaled[i]->GetFunction("expo")->Eval(hhall_scaled[i]->GetBinCenter(j)));
1142  truesum2[i] += hhups2[i]->GetBinContent(j);
1143  ersum2[i] += hhall_scaled[i]->GetBinError(j)*hhall_scaled[i]->GetBinError(j);
1144  sumpp2[i] += hhups2pp[i]->GetBinContent(j);
1145  ersumpp2[i] += hhupspp[i]->GetBinError(j)*hhupspp[i]->GetBinError(j);
1146  }
1147  sumsum2[i] = truesum2[i]; // direct count in mass range
1148  sumsum2pp[i] = sumpp2[i];
1149  //sumsum2[i] = hhups2[i]->GetEntries(); // total number of upsilons in pT bin (rounded up)
1150  //sumsum2pp[i] = hhups2pp[i]->GetEntries();
1151  //sumsum2[i] = Nups2*fUpsilonPt->Integral(s1,s2)/upsnorm; // total number of upsilons in pT bin
1152  //sumsum2pp[i] = Nups2pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
1153 
1154  if(sumsum2[i]>0. && sumsum2pp[i]>0.) {
1155  erraa2[i] = raa2[i]*sqrt(ersum2[i]/sumsum2[i]/sumsum2[i] + ersumpp2[i]/sumsum2pp[i]/sumsum2pp[i]);
1156  } else {raa2[i]=-1.0; erraa2[i] = 999.; }
1157 
1158  for(int j=fbin3; j<=lbin3; j++) {
1159  sum3[i] += (hhall_scaled[i]->GetBinContent(j) - hhcombbg_scaled[i]->GetFunction("fbg")->Eval(hhall_scaled[i]->GetBinCenter(j)) - hhcorrbg_scaled[i]->GetFunction("expo")->Eval(hhall_scaled[i]->GetBinCenter(j)));
1160  truesum3[i] += hhups3[i]->GetBinContent(j);
1161  ersum3[i] += hhall_scaled[i]->GetBinError(j)*hhall_scaled[i]->GetBinError(j);
1162  sumpp3[i] += hhups3pp[i]->GetBinContent(j);
1163  ersumpp3[i] += hhupspp[i]->GetBinError(j)*hhupspp[i]->GetBinError(j);
1164  }
1165  sumsum3[i] = truesum3[i]; // direct count in mass range
1166  sumsum3pp[i] = sumpp3[i];
1167  //sumsum3[i] = hhups3[i]->GetEntries(); // total number of upsilons in pT bin (rounded up)
1168  //sumsum3pp[i] = hhups3pp[i]->GetEntries();
1169  //sumsum3[i] = Nups3*fUpsilonPt->Integral(s1,s2)/upsnorm; // total number of upsilons in pT bin
1170  //sumsum3pp[i] = Nups3pp*fUpsilonPt->Integral(s1,s2)/upsnorm;
1171 
1172  if(truesum3[i]>0. && sumpp3[i]>0.) {
1173  erraa3[i] = raa3[i]*sqrt(ersum3[i]/sumsum3[i]/sumsum3[i] + ersumpp3[i]/sumsum3pp[i]/sumsum3pp[i]);
1174  } else {raa3[i]=-1.0; erraa3[i] = 999.; }
1175 
1176  }
1177 
1178 cout << "====== Y(1S):" << endl;
1179  for(int i=0; i<nbins+1; i++) {
1180  double s1 = double(i); double s2 = double(i+1); if(i==nbins) {s1 = 0.;}
1181  cout << " " << i << " " << truesum1[i] << "(" << Nups1*fUpsilonPt->Integral(s1,s2)/upsnorm << ")" << " +- "
1182  << sqrt(ersum1[i]) << " \t\t pp: " << sumpp1[i] << " +- " << sqrt(ersumpp1[i]) << endl;
1183  }
1184 cout << "====== Y(2S):" << endl;
1185  for(int i=0; i<nbins+1; i++) {
1186  double s1 = double(i); double s2 = double(i+1); if(i==nbins) {s1 = 0.;}
1187  cout << " " << i << " " << truesum2[i] << "(" << Nups2*fUpsilonPt->Integral(s1,s2)/upsnorm << ")" << " +- "
1188  << sqrt(ersum2[i]) << " \t\t pp: " << sumpp2[i] << " +- " << sqrt(ersumpp2[i]) << endl;
1189  }
1190 cout << "====== Y(3S):" << endl;
1191  for(int i=0; i<nbins+1; i++) {
1192 
1193  double s1 = double(i); double s2 = double(i+1); if(i==nbins) {s1 = 0.;}
1194  cout << " " << i << " " << truesum3[i] << "(" << Nups3*fUpsilonPt->Integral(s1,s2)/upsnorm << ")" << " +- "
1195  << sqrt(ersum3[i]) << " \t\t pp: " << sumpp3[i] << " +- " << sqrt(ersumpp3[i]) << endl;
1196  }
1197 
1198 //-------------------------------------------------
1199 // Plot RAA
1200 //-------------------------------------------------
1201 
1202 int npts1 = 10;
1203 int npts2 = 8;
1204 int npts3 = 7;
1205 
1206 TCanvas* craa = new TCanvas("craa","R_{AA}",120,120,600,600);
1207 TH2F* hh2 = new TH2F("hh2"," ",10,0.,float(npts1),10,0.,0.8);
1208 hh2->GetXaxis()->SetTitle("Transverse momentum [GeV/c]");
1209 hh2->GetXaxis()->SetTitleOffset(1.0);
1210 hh2->GetXaxis()->SetTitleColor(1);
1211 hh2->GetXaxis()->SetTitleSize(0.040);
1212 hh2->GetXaxis()->SetLabelSize(0.040);
1213 hh2->GetYaxis()->SetTitle("R_{AA}");
1214 hh2->GetYaxis()->SetTitleOffset(1.3);
1215 hh2->GetYaxis()->SetTitleSize(0.040);
1216 hh2->GetYaxis()->SetLabelSize(0.040);
1217 hh2->Draw();
1218 
1219 double xx1[nbins+1]; for(int i=0; i<nbins+1; i++) {xx1[i] = 0.5 + double(i);}
1220 double xx2[nbins+1]; for(int i=0; i<nbins+1; i++) {xx2[i] = 0.43 + double(i);}
1221 double xx3[nbins+1]; for(int i=0; i<nbins+1; i++) {xx3[i] = 0.57 + double(i);}
1222 
1223 TGraphErrors* gr1 = new TGraphErrors(npts1,xx1,raa1,0,erraa1);
1224 gr1->SetMarkerStyle(20);
1225 gr1->SetMarkerColor(kBlack);
1226 gr1->SetLineColor(kBlack);
1227 gr1->SetLineWidth(2);
1228 gr1->SetMarkerSize(1.5);
1229 gr1->SetName("gr1");
1230 gr1->Draw("p");
1231 
1232 erraa2[7] = erraa2[7]*0.85;
1233 
1234 TGraphErrors* gr2 = new TGraphErrors(npts2,xx2,raa2,0,erraa2);
1235 gr2->SetMarkerStyle(20);
1236 gr2->SetMarkerColor(kRed);
1237 gr2->SetLineColor(kRed);
1238 gr2->SetLineWidth(2);
1239 gr2->SetMarkerSize(1.5);
1240 gr2->SetName("gr2");
1241 gr2->Draw("p");
1242 
1243 /*
1244 //--- 3S state -------------------
1245  double dummy[9];
1246  for(int i=0; i<9; i++) {dummy[i]=-99.;}
1247  TGraphErrors* gr3 = new TGraphErrors(6,xx3,dummy,0,erraa3);
1248  gr3->SetMarkerStyle(20);
1249  gr3->SetMarkerColor(kBlue);
1250  gr3->SetLineColor(kBlue);
1251  gr3->SetLineWidth(2);
1252  gr3->SetMarkerSize(1.5);
1253  gr3->SetName("gr3");
1254 // gr3->Draw("p");
1255 
1256 TArrow* aa[9];
1257 TLine* ll[9];
1258 //erraa3[3] = erraa3[3]*1.5;
1259 erraa3[4] = erraa3[4]*0.6;
1260 //erraa3[4] = erraa3[4]*1.5;
1261 erraa3[6] = erraa3[6]*1.2;
1262 
1263 for(int i=0; i<npts3; i++) {
1264  aa[i] = new TArrow(xx3[i],1.64*erraa3[i],xx3[i],-0.02,0.02);
1265  aa[i]->SetLineColor(kBlue);
1266  aa[i]->SetLineWidth(2);
1267  aa[i]->Draw();
1268  ll[i] = new TLine(xx3[i]-0.15,1.64*erraa3[i],xx3[i]+0.15,1.64*erraa3[i]);
1269  ll[i]->SetLineColor(kBlue);
1270  ll[i]->SetLineWidth(2);
1271  ll[i]->Draw();
1272 }
1273 //--- end 3S state --------------
1274 */
1275 
1276 TLegend *leg = new TLegend(0.2,0.77,0.5,0.92);
1277 leg->SetFillColor(10);
1278 leg->SetFillStyle(1001);
1279 TLegendEntry *entry1=leg->AddEntry("gr1","Y(1S)","p");
1280 TLegendEntry *entry2=leg->AddEntry("gr2","Y(2S)","p");
1281 TLegendEntry *entry3=leg->AddEntry("gr3","Y(3S)","l");
1282 leg->Draw();
1283 
1284 TLatex* l1 = new TLatex(3.5,0.73,"sPHENIX Simulation");
1285 l1->Draw();
1286 
1287 TLine* lll = new TLine(0.6,0.64,1.3,0.64);
1288 lll->SetLineColor(kBlue);
1289 lll->SetLineWidth(2);
1290 lll->Draw();
1291 
1292 //==================================================================================
1293 
1294 /*
1295 fout = new TFile("test.root","recreate");
1296 for(int i=0; i<nbins+1; i++) {
1297  sprintf(tmpname,"hhfit_%d",i);
1298  hhfit[i]->Write();
1299  hhups1pp[i]->Write();
1300  hhups2pp[i]->Write();
1301  hhups3pp[i]->Write();
1302  hhupspp[i]->Write();
1303  hhups1[i]->Write();
1304  hhups2[i]->Write();
1305  hhups3[i]->Write();
1306  hhups[i]->Write();
1307  hhall_pp[i]->Write();
1308  hhcorrbg_pp[i]->Write();
1309  hhcorrbg_scaled[i]->Write();
1310 }
1311 fout->Write();
1312 fout->Close();
1313 */
1314 
1315 }
1316