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