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