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