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