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