Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
fit_vs_directcount.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file fit_vs_directcount.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 double SandB_CBFunction(double *x, double *p)
58 {
59  double norm1 = p[0]; // amplitude of Y(1S) peak
60  double alpha = p[1]; // tail start
61  double n = p[2]; // tail shape
62  double mu1 = p[3]; // upsilon (1S) mass
63  double sigma = p[4]; // upsilon width
64  double norm2 = p[5];
65  double norm3 = p[6];
66  double mu2 = mu1*1.0595;
67  double mu3 = mu1*1.0946;
68 
69  double A = pow(n/fabs(alpha),n)*TMath::Exp(-pow(fabs(alpha),2)/2.);
70  double B = n/fabs(alpha) - fabs(alpha);
71  double k1 = (x[0]-mu1)/sigma;
72  double k2 = (x[0]-mu2)/sigma;
73  double k3 = (x[0]-mu3)/sigma;
74 
75  double val,val1,val2,val3;
76 
77  if( k1 > -alpha ) { val1 = norm1*TMath::Exp(-0.5*pow(k1,2)); }
78  else { val1 = norm1*A*pow(B-k1,-n); }
79  if( k2 > -alpha ) { val2 = norm2*TMath::Exp(-0.5*pow(k2,2)); }
80  else { val2 = norm2*A*pow(B-k2,-n); }
81  if( k3 > -alpha ) { val3 = norm3*TMath::Exp(-0.5*pow(k3,2)); }
82  else { val3 = norm3*A*pow(B-k3,-n); }
83 
84  double bgnorm1 = p[7];
85  double bgslope1 = p[8];
86 
87  double bg = exp(bgnorm1+x[0]*bgslope1);
88 
89  val = val1 + val2 + val3 + bg;
90  if( TMath::IsNaN(val) ) val = 0.0;
91 
92  return val;
93 }
94 
95 //=======================================================================
96 
97 void fit_vs_directcount() {
98 
99 const int nbins = 15;
100 
101 // TRUE numbers
102 //Number of upsilons in AuAu plot = 11752 971 111
103 //Number of upsilons in pp plot = 9755 2537 1414
104 double Nups1[nbins+1],Nups2[nbins+1],Nups3[nbins+1];
105 Nups1[0]=1020.82;
106 Nups2[0]=84.3446;
107 Nups3[0]=9.64186;
108 Nups1[1]=2519.56;
109 Nups2[1]=208.177;
110 Nups3[1]=23.7977;
111 Nups1[2]=2870.95;
112 Nups2[2]=237.21;
113 Nups3[2]=27.1167;
114 Nups1[3]=2326.04;
115 Nups2[3]=192.187;
116 Nups3[3]=21.9699;
117 Nups1[4]=1500.87;
118 Nups2[4]=124.008;
119 Nups3[4]=14.176;
120 Nups1[5]=820.126;
121 Nups2[5]=67.7622;
122 Nups3[5]=7.74625;
123 Nups1[6]=396.538;
124 Nups2[6]=32.7637;
125 Nups3[6]=3.74538;
126 Nups1[7]=175.598;
127 Nups2[7]=14.5086;
128 Nups3[7]=1.65856;
129 Nups1[8]=73.2069;
130 Nups2[8]=6.04867;
131 Nups3[8]=0.691454;
132 Nups1[9]=29.3682;
133 Nups2[9]=2.42653;
134 Nups3[9]=0.277389;
135 Nups1[10]=11.5318;
136 Nups2[10]=0.952802;
137 Nups3[10]=0.10892;
138 Nups1[11]=4.49014;
139 Nups2[11]=0.370994;
140 Nups3[11]=0.0424103;
141 Nups1[12]=1.75068;
142 Nups2[12]=0.144648;
143 Nups3[12]=0.0165355;
144 Nups1[13]=0.688396;
145 Nups2[13]=0.0568782;
146 Nups3[13]=0.00650204;
147 Nups1[14]=0.274397;
148 Nups2[14]=0.0226719;
149 Nups3[14]=0.00259174;
150 
151 double Nups1pp[nbins+1],Nups2pp[nbins+1],Nups3pp[nbins+1];
152 Nups1pp[0]=847.355;
153 Nups2pp[0]=220.373;
154 Nups3pp[0]=122.825;
155 Nups1pp[1]=2091.41;
156 Nups2pp[1]=543.917;
157 Nups3pp[1]=303.153;
158 Nups1pp[2]=2383.1;
159 Nups2pp[2]=619.776;
160 Nups3pp[2]=345.433;
161 Nups1pp[3]=1930.78;
162 Nups2pp[3]=502.14;
163 Nups3pp[3]=279.869;
164 Nups1pp[4]=1245.83;
165 Nups2pp[4]=324.005;
166 Nups3pp[4]=180.585;
167 Nups1pp[5]=680.763;
168 Nups2pp[5]=177.047;
169 Nups3pp[5]=98.6775;
170 Nups1pp[6]=329.155;
171 Nups2pp[6]=85.6039;
172 Nups3pp[6]=47.7115;
173 Nups1pp[7]=145.759;
174 Nups2pp[7]=37.9077;
175 Nups3pp[7]=21.1279;
176 Nups1pp[8]=60.767;
177 Nups2pp[8]=15.8038;
178 Nups3pp[8]=8.80825;
179 Nups1pp[9]=24.3777;
180 Nups2pp[9]=6.33996;
181 Nups3pp[9]=3.53358;
182 Nups1pp[10]=9.57218;
183 Nups2pp[10]=2.48945;
184 Nups3pp[10]=1.3875;
185 Nups1pp[11]=3.72714;
186 Nups2pp[11]=0.969323;
187 Nups3pp[11]=0.540253;
188 Nups1pp[12]=1.45319;
189 Nups2pp[12]=0.377933;
190 Nups3pp[12]=0.210641;
191 Nups1pp[13]=0.571418;
192 Nups2pp[13]=0.14861;
193 Nups3pp[13]=0.0828277;
194 Nups1pp[14]=0.227769;
195 Nups2pp[14]=0.0592364;
196 Nups3pp[14]=0.0330155;
197 
198 // new suppression
199 double raapt[9],raa1s[9],raa2s[9],oldraa1s[9],oldraa2s[9];
200 raapt[0] = 1.5;
201 raapt[1] = 4.5;
202 raapt[2] = 7.5;
203 raapt[3] = 10.5;
204 raapt[4] = 13.5;
205 oldraa1s[0] = 0.535;
206 oldraa1s[1] = 0.535;
207 oldraa1s[2] = 0.535;
208 oldraa1s[3] = 0.535;
209 oldraa1s[4] = 0.535;
210 oldraa2s[0] = 0.170;
211 oldraa2s[1] = 0.170;
212 oldraa2s[2] = 0.170;
213 oldraa2s[3] = 0.170;
214 oldraa2s[4] = 0.170;
215 raa1s[0] = 0.4960;
216 raa1s[1] = 0.4960;
217 raa1s[2] = 0.4955;
218 raa1s[3] = 0.4968;
219 raa1s[4] = 0.4743;
220 raa2s[0] = 0.1710;
221 raa2s[1] = 0.1629;
222 raa2s[2] = 0.1326;
223 raa2s[3] = 0.1232;
224 raa2s[4] = 0.0928;
225 TGraph* grRAA1S = new TGraph(5,raapt,raa1s);
226 TGraph* grRAA2S = new TGraph(5,raapt,raa2s);
227 TGraph* groldRAA1S = new TGraph(5,raapt,oldraa1s);
228 TGraph* groldRAA2S = new TGraph(5,raapt,oldraa2s);
229 for(int i=0; i<nbins; i++) {
230  Nups1[i] = Nups1[i] * grRAA1S->Eval(0.5+i*1.0)/groldRAA1S->Eval(0.5+i*1.0);
231  Nups2[i] = Nups2[i] * grRAA2S->Eval(0.5+i*1.0)/groldRAA2S->Eval(0.5+i*1.0);
232 }
233 grRAA1S->SetLineColor(kBlue);
234 grRAA1S->SetLineStyle(3);
235 grRAA1S->SetLineWidth(3);
236 grRAA2S->SetLineColor(kMagenta);
237 grRAA2S->SetLineStyle(3);
238 grRAA2S->SetLineWidth(3);
239 
240 
241 
242 //gROOT->LoadMacro("sPHENIXStyle/sPhenixStyle.C");
243 //SetsPhenixStyle();
244 gStyle->SetOptStat(0);
245 gStyle->SetOptFit(0);
246 
247 char tmpname[999];
248 TLatex* tl[nbins+1];
249 char tlchar[999];
250 
251 TH1D* hhfit[nbins+1];
252 TH1D* hhcorrbg_scaled[nbins+1];
253 TH1D* hhups1pp[nbins+1];
254 TH1D* hhups2pp[nbins+1];
255 TH1D* hhups3pp[nbins+1];
256 TH1D* hhupspp[nbins+1];
257 TH1D* hhups1[nbins+1];
258 TH1D* hhups2[nbins+1];
259 TH1D* hhups3[nbins+1];
260 TH1D* hhups[nbins+1];
261 TH1D* hhall_pp[nbins+1];
262 TH1D* hhcorrbg_pp[nbins+1];
263 
264 TF1* fCBups1s = new TF1("fCBups1s",CBFunction,5.,14.,5);
265 TF1* fCBups2s = new TF1("fCBups2s",CBFunction,5.,14.,5);
266 TF1* fCBups1spp = new TF1("fCBups1spp",CBFunction,5.,14.,5);
267 TF1* fCBups2spp = new TF1("fCBups2spp",CBFunction,5.,14.,5);
268 TF1* fTCB = new TF1("fTCB",TripleCBFunction,5.,14.,7);
269 TF1* fSandB = new TF1("fSandB",SandB_CBFunction,5.,14.,9);
270 TF1* fSandBpp = new TF1("fSandBpp",SandB_CBFunction,5.,14.,9);
271 TF1* fSandBauau = new TF1("fSandBauau",SandB_CBFunction,5.,14.,9);
272 
273 double tonypar1 = 0.98; // Tony's 100 pion simulation April 2019
274 double tonypar2 = 0.93; // Tony's 100 pion simulation April 2019
275 //double tonypar3 = 9.437; // Tony's 100 pion simulation April 2019
276 double tonypar3 = 9.448; // benchmark
277 double tonypar4 = 0.100; // benchmark
278 
279 double u1start = 9.25;
280 double u1stop = 9.65;
281 double u2start = 9.80;
282 double u2stop = 10.20;
283 double u3start = 10.20;
284 double u3stop = 10.55;
285 double sum1[nbins+1] = {0.};
286 double truesum1[nbins+1] = {0.};
287 double truesum1pp[nbins+1] = {0.};
288 double ersum1[nbins+1] = {0.};
289 double sum2[nbins+1] = {0.};
290 double truesum2[nbins+1] = {0.};
291 double truesum2pp[nbins+1] = {0.};
292 double ersum2[nbins+1] = {0.};
293 double sum1pp[nbins+1] = {0.};
294 double ersum1pp[nbins+1] = {0.};
295 double sum2pp[nbins+1] = {0.};
296 double ersum2pp[nbins+1] = {0.};
297 double sumfit1[nbins+1] = {0.};
298 double ersumfit1[nbins+1] = {0.};
299 double sumfit2[nbins+1] = {0.};
300 double ersumfit2[nbins+1] = {0.};
301 double sumfit1pp[nbins+1] = {0.};
302 double ersumfit1pp[nbins+1] = {0.};
303 double sumfit2pp[nbins+1] = {0.};
304 double ersumfit2pp[nbins+1] = {0.};
305 double sum1ppbg[nbins+1] = {0.};
306 double ersum1ppbg[nbins+1] = {0.};
307 double sum2ppbg[nbins+1] = {0.};
308 double ersum2ppbg[nbins+1] = {0.};
309 
310 double xx1[nbins+1]; for(int i=0; i<nbins+1; i++) {xx1[i] = 0.50 + double(i);}
311 double xx2[nbins+1]; for(int i=0; i<nbins+1; i++) {xx2[i] = 0.40 + double(i);}
312 double xx3[nbins+1]; for(int i=0; i<nbins+1; i++) {xx3[i] = 0.60 + double(i);}
313 
314 //=====================================================================================================
315 
316 // read histograms created by newplotbg.C
317 
318 TFile* f=new TFile("ups_corrbg_24b_auau.root");
319 for(int i=0; i<nbins+1; i++) {
320 
321  sprintf(tmpname,"hhfit_%d",i);
322  hhfit[i] = (TH1D*)f->Get(tmpname);
323  hhfit[i]->SetDirectory(gROOT);
324  sprintf(tmpname,"hhcorrbg_scaled_%d",i);
325  hhcorrbg_scaled[i] = (TH1D*)f->Get(tmpname);
326  hhcorrbg_scaled[i]->SetDirectory(gROOT);
327 
328  sprintf(tmpname,"hhups1pp_%d",i);
329  hhups1pp[i] = (TH1D*)f->Get(tmpname);
330  hhups1pp[i]->SetDirectory(gROOT);
331  sprintf(tmpname,"hhups2pp_%d",i);
332  hhups2pp[i] = (TH1D*)f->Get(tmpname);
333  hhups2pp[i]->SetDirectory(gROOT);
334  sprintf(tmpname,"hhups3pp_%d",i);
335  hhups3pp[i] = (TH1D*)f->Get(tmpname);
336  hhups3pp[i]->SetDirectory(gROOT);
337  sprintf(tmpname,"hhupspp_%d",i);
338  hhupspp[i] = (TH1D*)f->Get(tmpname);
339  hhupspp[i]->SetDirectory(gROOT);
340 
341  sprintf(tmpname,"hhups1_%d",i);
342  hhups1[i] = (TH1D*)f->Get(tmpname);
343  hhups1[i]->SetDirectory(gROOT);
344  sprintf(tmpname,"hhups2_%d",i);
345  hhups2[i] = (TH1D*)f->Get(tmpname);
346  hhups2[i]->SetDirectory(gROOT);
347  sprintf(tmpname,"hhups3_%d",i);
348  hhups3[i] = (TH1D*)f->Get(tmpname);
349  hhups3[i]->SetDirectory(gROOT);
350  sprintf(tmpname,"hhups_%d",i);
351  hhups[i] = (TH1D*)f->Get(tmpname);
352  hhups[i]->SetDirectory(gROOT);
353 
354  sprintf(tmpname,"hhall_pp_%d",i);
355  hhall_pp[i] = (TH1D*)f->Get(tmpname);
356  hhall_pp[i]->SetDirectory(gROOT);
357  sprintf(tmpname,"hhcorrbg_pp_%d",i);
358  hhcorrbg_pp[i] = (TH1D*)f->Get(tmpname);
359  hhcorrbg_pp[i]->SetDirectory(gROOT);
360 
361 }
362 f->Close();
363 
364  // fg and bg parameters from newplotbg.C
365  double bgpar0[nbins+1],bgpar1[nbins+1];
366  for(int i=0; i<nbins+1; i++) {
367  bgpar0[i] = hhcorrbg_scaled[i]->GetFunction("expo")->GetParameter(0);
368  bgpar1[i] = hhcorrbg_scaled[i]->GetFunction("expo")->GetParameter(1);
369  }
370  double fgpar0 = hhups[nbins]->GetFunction("fTCBauau")->GetParameter(0);
371  double fgpar5 = hhups[nbins]->GetFunction("fTCBauau")->GetParameter(5);
372  double fgpar6 = hhups[nbins]->GetFunction("fTCBauau")->GetParameter(6);
373 
374 //-------------------------------------------------------------------------------------
375 
376 TCanvas* cups = new TCanvas("cups","Upsilons and correlated BG",150,150,700,700);
377 
378  // "true" distribution
379  fSandBauau->SetParameter(0,fgpar0);
380  fSandBauau->SetParameter(1,tonypar1);
381  fSandBauau->SetParameter(2,tonypar2);
382  fSandBauau->SetParameter(3,tonypar3);
383  fSandBauau->SetParameter(4,tonypar4);
384  fSandBauau->SetParameter(5,fgpar5);
385  fSandBauau->SetParameter(6,fgpar6);
386  fSandBauau->SetParameter(7,bgpar0[nbins]);
387  fSandBauau->SetParameter(8,bgpar1[nbins]);
388  fSandBauau->SetLineColor(kRed);
389 // fSandBauau->Draw("same");
390 
391  // "true" Upsilon1S
392  fCBups1s->SetParameter(0,fSandBauau->GetParameter(0));
393  fCBups1s->SetParameter(1,fSandBauau->GetParameter(1));
394  fCBups1s->SetParameter(2,fSandBauau->GetParameter(2));
395  fCBups1s->SetParameter(3,fSandBauau->GetParameter(3));
396  fCBups1s->SetParameter(4,fSandBauau->GetParameter(4));
397  double true_ups1s_integral = fCBups1s->Integral(5.,14.);
398  double true_ups1s_ampl = fSandBauau->GetParameter(0);
399  double true_ups1s_amplerr = fSandBauau->GetParError(0);
400  double binsize = hhfit[nbins]->GetBinWidth(1);
401  cout << "TRUE Integral = " << true_ups1s_integral/binsize << " +- " << true_ups1s_integral*(true_ups1s_amplerr/true_ups1s_ampl)/binsize << " ( " << true_ups1s_amplerr/true_ups1s_ampl*100. << "% )" << endl;
402 
403  // "true" corr bg
404  TF1* fcorrbg[nbins+1];
405  fcorrbg[nbins] = new TF1("fcorrbg_15","exp([0]+[1]*x)",7.,14.);
406  fcorrbg[nbins]->SetParameters(bgpar0[nbins],bgpar1[nbins]);
407  //fcorrbg[nbins]->SetLineStyle(2);
408  fcorrbg[nbins]->SetLineWidth(1);
409  fcorrbg[nbins]->SetLineColor(kRed);
410 
411 //--- FIT all pT -------------------------------------------------------------------------
412 
413  fSandB->SetParameter(0,fgpar0);
414  fSandB->FixParameter(1,tonypar1);
415  fSandB->FixParameter(2,tonypar2);
416  fSandB->FixParameter(3,tonypar3);
417  fSandB->FixParameter(4,tonypar4);
418  fSandB->SetParameter(5,fgpar5);
419  fSandB->SetParameter(6,fgpar6);
420  fSandB->SetParameter(7,bgpar0[nbins]);
421  fSandB->SetParameter(8,bgpar1[nbins]);
422  hhfit[nbins]->Fit(fSandB,"qrl","",8.,11.);
423 
424  hhfit[nbins]->GetXaxis()->SetTitleOffset(1.0);
425  hhfit[nbins]->GetYaxis()->SetTitleOffset(1.6);
426  hhfit[nbins]->SetAxisRange(8.0,11.0);
427  hhfit[nbins]->Draw();
428 
429  TF1* fcorrbg_fromfit[nbins+1];
430  sprintf(tmpname,"fcorrbg_fromfit_%d",15);
431  fcorrbg_fromfit[nbins] = new TF1(tmpname,"exp([0]+[1]*x)",7.,14.);
432  fcorrbg_fromfit[nbins]->SetParameters(fSandB->GetParameter(7),fSandB->GetParameter(8));
433  fcorrbg_fromfit[nbins]->SetLineStyle(2);
434  fcorrbg_fromfit[nbins]->SetLineWidth(3);
435  fcorrbg_fromfit[nbins]->Draw("same");
436  fcorrbg[nbins]->Draw("same"); // true correlated bg
437 
438  // counting range
439  TLine* line1 = new TLine(u1start,0.,u1start,hhfit[nbins]->GetMaximum()*1.1); line1->SetLineStyle(2); line1->Draw();
440  TLine* line2 = new TLine(u1stop,0.,u1stop,hhfit[nbins]->GetMaximum()*1.1); line2->SetLineStyle(2); line2->Draw();
441 
442  // Upsilon1S from fit
443  fCBups1s->SetParameter(0,fSandB->GetParameter(0));
444  fCBups1s->SetParameter(1,fSandB->GetParameter(1));
445  fCBups1s->SetParameter(2,fSandB->GetParameter(2));
446  fCBups1s->SetParameter(3,fSandB->GetParameter(3));
447  fCBups1s->SetParameter(4,fSandB->GetParameter(4));
448  double ups1s_integral = fCBups1s->Integral(5.,14.);
449  double ups1s_ampl = fSandB->GetParameter(0);
450  double ups1s_amplerr = fSandB->GetParError(0);
451  cout << "Integral = " << ups1s_integral/binsize << " +- " << ups1s_integral*(ups1s_amplerr/ups1s_ampl)/binsize << " ( " << ups1s_amplerr/ups1s_ampl*100. << "% )" << endl;
452 
453 //--- Direct Counting all pT ------------------------------------------------------
454 
455  int fbin1 = hhfit[nbins]->FindBin(u1start + 0.001);
456  int lbin1 = hhfit[nbins]->FindBin(u1stop - 0.001);
457  int fbin2 = hhfit[nbins]->FindBin(u2start + 0.001);
458  int lbin2 = hhfit[nbins]->FindBin(u2stop - 0.001);
459  int fbin3 = hhfit[nbins]->FindBin(u3start + 0.001);
460  int lbin3 = hhfit[nbins]->FindBin(u3stop - 0.001);
461 
462  for(int j=fbin1; j<=lbin1; j++) {
463  //sum1[nbins] += (hhfit[nbins]->GetBinContent(j) - fcorrbg_fromfit[nbins]->Eval(hhfit[nbins]->GetBinCenter(j)));
464  sum1[nbins] += (hhfit[nbins]->GetBinContent(j) - fcorrbg[nbins]->Eval(hhfit[nbins]->GetBinCenter(j)));
465  ersum1[nbins] += hhfit[nbins]->GetBinError(j)*hhfit[nbins]->GetBinError(j);
466  truesum1[nbins] += hhups1[nbins]->GetBinContent(j);
467  }
468  ersum1[nbins] = sqrt(ersum1[nbins]);
469  cout << "Direct count = " << truesum1[nbins] << " +- " << ersum1[nbins] << " ( " << ersum1[nbins]/truesum1[nbins]*100. << "% )" << endl;
470 
471 //============================================================================================
472 // VS pT
473 //============================================================================================
474 
475 int npts=10;
476 
477 TCanvas* dummy = new TCanvas("dummy","dummy",0,0,500,500);
478 
479  for(int i=0; i<npts; i++) {
480 
481  fSandB->SetParameter(0,hhfit[i]->GetMaximum()/9.);
482  fSandB->SetParLimits(0,0.,99999.);
483  fSandB->FixParameter(1,tonypar1);
484  fSandB->FixParameter(2,tonypar2);
485  fSandB->FixParameter(3,tonypar3);
486  fSandB->FixParameter(4,tonypar4);
487  fSandB->SetParameter(5,fgpar5/1.);
488  fSandB->SetParLimits(5,0.5,9999.);
489  fSandB->SetParameter(6,fgpar6/1.);
490  fSandB->SetParLimits(6,0.,999.);
491  fSandB->SetParameter(7,bgpar0[i]/9.);
492  fSandB->SetParLimits(7,0.,99.);
493  fSandB->FixParameter(8,bgpar1[i]);
494  hhfit[i]->Fit(fSandB,"qrl","",8.,11.);
495 
496  fCBups1s->SetParameter(0,fSandB->GetParameter(0));
497  fCBups1s->SetParameter(1,fSandB->GetParameter(1));
498  fCBups1s->SetParameter(2,fSandB->GetParameter(2));
499  fCBups1s->SetParameter(3,fSandB->GetParameter(3));
500  fCBups1s->SetParameter(4,fSandB->GetParameter(4));
501  double ups1s_integral = fCBups1s->Integral(5.,14.);
502  double ups1s_ampl = fSandB->GetParameter(0);
503  double ups1s_amplerr = fSandB->GetParError(0);
504 // cout << i << " Integral1 = " << ups1s_integral/binsize << " +- " << ups1s_integral*(ups1s_amplerr/ups1s_ampl)/binsize << " ( " << ups1s_amplerr/ups1s_ampl*100. << "% )" << endl;
505  sumfit1[i] = ups1s_integral/binsize;
506  ersumfit1[i] = ups1s_integral*(ups1s_amplerr/ups1s_ampl)/binsize;
507 
508  fCBups2s->SetParameter(0,fSandB->GetParameter(5));
509  fCBups2s->SetParameter(1,fSandB->GetParameter(1));
510  fCBups2s->SetParameter(2,fSandB->GetParameter(2));
511  fCBups2s->SetParameter(3,fSandB->GetParameter(3));
512  fCBups2s->SetParameter(4,fSandB->GetParameter(4));
513  double ups2s_integral = fCBups2s->Integral(5.,14.);
514  double ups2s_ampl = fSandB->GetParameter(5);
515  double ups2s_amplerr = fSandB->GetParError(5);
516  cout << i << " Integral2 = " << ups2s_integral/binsize << " +- " << ups2s_integral*(ups2s_amplerr/ups2s_ampl)/binsize << " ( " << ups2s_amplerr/ups2s_ampl*100. << "% )" << endl;
517  sumfit2[i] = ups2s_integral/binsize;
518  ersumfit2[i] = ups2s_integral*(ups2s_amplerr/ups2s_ampl)/binsize;
519 
520  sprintf(tmpname,"fcorrbg_fromfit_%d",i); // correlated bg from fit
521  fcorrbg_fromfit[i] = new TF1(tmpname,"exp([0]+[1]*x)",7.,14.);
522  fcorrbg_fromfit[i]->SetLineStyle(2);
523  fcorrbg_fromfit[i]->SetParameters(fSandB->GetParameter(7),fSandB->GetParameter(8));
524 
525  sprintf(tmpname,"fcorrbg_%d",i); // true correlated bg
526  fcorrbg[i] = new TF1(tmpname,"exp([0]+[1]*x)",7.,14.);
527  fcorrbg[i]->SetLineWidth(1);
528  fcorrbg[i]->SetLineColor(kRed);
529  fcorrbg[i]->SetParameters(bgpar0[i],bgpar1[i]);
530 
531 //--- vs pT direct counting ----------------------------------------------------
532 
533  sum1[i]=0.;
534  truesum1[i]=0.;
535  ersum1[i]=0.;
536  sum1pp[i]=0.;
537  ersum1pp[i]=0.;
538  for(int j=fbin1; j<=lbin1; j++) {
539  //sum1[i] += (hhfit[i]->GetBinContent(j) - fcorrbg_fromfit[i]->Eval(hhfit[i]->GetBinCenter(j)));
540  sum1[i] += (hhfit[i]->GetBinContent(j) - fcorrbg[i]->Eval(hhfit[i]->GetBinCenter(j)));
541  ersum1[i] += hhfit[i]->GetBinError(j)*hhfit[i]->GetBinError(j);
542  sum1pp[i] += hhupspp[i]->GetBinContent(j);
543  ersum1pp[i] += hhupspp[i]->GetBinError(j)*hhupspp[i]->GetBinError(j);
544  truesum1[i] += hhups1[i]->GetBinContent(j);
545  }
546  ersum1[i] = sqrt(ersum1[i]);
547  ersum1pp[i] = sqrt(ersum1pp[i]);
548  //cout << " Direct count1 = " << truesum1[i] << " +- " << ersum1[i] << " ( " << ersum1[i]/sum1[i]*100. << "% )" << endl;
549 
550  sum2[i]=0.;
551  truesum2[i]=0.;
552  ersum2[i]=0.;
553  sum2pp[i]=0.;
554  ersum2pp[i]=0.;
555  for(int j=fbin2; j<=lbin2; j++) {
556  //sum2[i] += (hhfit[i]->GetBinContent(j) - fcorrbg_fromfit[i]->Eval(hhfit[i]->GetBinCenter(j)));
557  sum2[i] += (hhfit[i]->GetBinContent(j) - fcorrbg[i]->Eval(hhfit[i]->GetBinCenter(j)));
558  ersum2[i] += hhfit[i]->GetBinError(j)*hhfit[i]->GetBinError(j);
559  sum2pp[i] += hhupspp[i]->GetBinContent(j);
560  ersum2pp[i] += hhupspp[i]->GetBinError(j)*hhupspp[i]->GetBinError(j);
561  truesum2[i] += hhups2[i]->GetBinContent(j);
562  }
563  ersum2[i] = sqrt(ersum2[i]);
564  ersum2pp[i] = sqrt(ersum2pp[i]);
565  cout << " Direct count2 = " << truesum2[i] << " +- " << ersum2[i] << " ( " << ersum2[i]/sum2[i]*100. << "% )" << endl;
566 
567 
568  } // end loop over pT bins
569 
570 delete dummy;
571 
572 TCanvas* cupsvspt_fit = new TCanvas("cupsvspt_fit","Upsilons vs. p_{T} after fit",100,100,1100,900);
573 cupsvspt_fit->Divide(3,3);
574 for(int i=0; i<9; i++) {
575 if(i==0) {cupsvspt_fit->cd(0);}
576 if(i==1) {cupsvspt_fit->cd(1);}
577 if(i==2) {cupsvspt_fit->cd(2);}
578 if(i==3) {cupsvspt_fit->cd(3);}
579 if(i==4) {cupsvspt_fit->cd(4);}
580 if(i==5) {cupsvspt_fit->cd(5);}
581 if(i==6) {cupsvspt_fit->cd(6);}
582 if(i==7) {cupsvspt_fit->cd(7);}
583 if(i==8) {cupsvspt_fit->cd(8);}
584  hhfit[i]->SetAxisRange(8.0,11.0);
585  hhfit[i]->Draw();
586  fcorrbg_fromfit[i]->Draw("same");
587  fcorrbg[i]->Draw("same");
588  sprintf(tlchar,"%d-%d GeV",i,i+1); tl[i] = new TLatex(8.2,hhfit[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
589 }
590 
591 
592 
593 
594 //==========================================================================
595 // p+p
596 //==========================================================================
597 cout << endl << "---------- p+p ---------------------------------------" << endl << endl;
598 
599 double bgpar0pp[nbins+1],bgpar1pp[nbins+1];
600 for(int i=0; i<nbins+1; i++) {
601  bgpar0pp[i] = hhcorrbg_pp[i]->GetFunction("expo")->GetParameter(0);
602  bgpar1pp[i] = hhcorrbg_pp[i]->GetFunction("expo")->GetParameter(1);
603 }
604 TF1* fcorrbg_pp[nbins+1];
605 TF1* fcorrbg_pp_fromfit[nbins+1];
606 
607 TCanvas* cpp = new TCanvas("cpp","p+p",100,100,700,700);
608 
609  double tmpmax = hhall_pp[nbins]->GetMaximum();
610 
611  fSandBpp->SetParameter(0,tmpmax/9.);
612  fSandBpp->SetParLimits(0,0.,99999.);
613  fSandBpp->FixParameter(1,tonypar1);
614  fSandBpp->FixParameter(2,tonypar2);
615  fSandBpp->FixParameter(3,tonypar3);
616  fSandBpp->FixParameter(4,0.089); // p+p mass resolution
617  fSandBpp->SetParameter(5,tmpmax/99.);
618  fSandBpp->SetParLimits(5,0.5,9999.);
619  fSandBpp->SetParameter(6,tmpmax/99.);
620  fSandBpp->SetParLimits(6,0.,999.);
621  fSandBpp->SetParameter(7,bgpar0pp[nbins]/9.);
622  fSandBpp->SetParLimits(7,0.,99.);
623  fSandBpp->SetParameter(8,bgpar1pp[nbins]);
624  hhall_pp[nbins]->Fit(fSandBpp,"qrl","",7.,12.);
625  hhall_pp[nbins]->SetAxisRange(8.,11.);
626  hhall_pp[nbins]->Draw();
627 // hhcorrbg_pp[nbins]->GetFunction("expo")->SetLineWidth(1);
628 // hhcorrbg_pp[nbins]->GetFunction("expo")->SetLineColor(kRed);
629 // hhcorrbg_pp[nbins]->GetFunction("expo")->Draw("same"); // "true" correlated bg
630  sprintf(tmpname,"fcorrbg_pp_%d",nbins); // "true" correlated bg copy
631  fcorrbg_pp[nbins] = new TF1(tmpname,"exp([0]+[1]*x)",7.,14.);
632  fcorrbg_pp[nbins]->SetLineColor(kRed);
633  fcorrbg_pp[nbins]->SetLineWidth(1);
634  fcorrbg_pp[nbins]->SetParameters(bgpar0[nbins], bgpar1pp[nbins]);
635  fcorrbg_pp[nbins]->Draw("same");
636 
637  sprintf(tmpname,"fcorrbg_pp_fromfit_%d",nbins); // correlated bg from fit
638  fcorrbg_pp_fromfit[nbins] = new TF1(tmpname,"exp([0]+[1]*x)",7.,14.);
639  fcorrbg_pp_fromfit[nbins]->SetLineStyle(2);
640  fcorrbg_pp_fromfit[nbins]->SetParameters(hhall_pp[nbins]->GetFunction("fSandBpp")->GetParameter(7), hhall_pp[nbins]->GetFunction("fSandBpp")->GetParameter(8));
641  fcorrbg_pp_fromfit[nbins]->Draw("same");
642 
643 
644  fCBups1spp->SetParameter(0,fSandBpp->GetParameter(0));
645  fCBups1spp->SetParameter(1,fSandBpp->GetParameter(1));
646  fCBups1spp->SetParameter(2,fSandBpp->GetParameter(2));
647  fCBups1spp->SetParameter(3,fSandBpp->GetParameter(3));
648  fCBups1spp->SetParameter(4,fSandBpp->GetParameter(4));
649  double ups1spp_integral = fCBups1spp->Integral(5.,14.);
650  double ups1spp_ampl = fSandBpp->GetParameter(0);
651  double ups1spp_amplerr = fSandBpp->GetParError(0);
652  cout << "p+p Integral1 = " << ups1spp_integral/binsize << " +- " << ups1spp_integral*(ups1spp_amplerr/ups1spp_ampl)/binsize << " ( " << ups1spp_amplerr/ups1spp_ampl*100. << "% )" << endl;
653  fCBups2spp->SetParameter(0,fSandBpp->GetParameter(5));
654  fCBups2spp->SetParameter(1,fSandBpp->GetParameter(1));
655  fCBups2spp->SetParameter(2,fSandBpp->GetParameter(2));
656  fCBups2spp->SetParameter(3,fSandBpp->GetParameter(3));
657  fCBups2spp->SetParameter(4,fSandBpp->GetParameter(4));
658  double ups2spp_integral = fCBups2spp->Integral(5.,14.);
659  double ups2spp_ampl = fSandBpp->GetParameter(5);
660  double ups2spp_amplerr = fSandBpp->GetParError(5);
661  cout << "p+p Integral2 = " << ups2spp_integral/binsize << " +- " << ups2spp_integral*(ups2spp_amplerr/ups2spp_ampl)/binsize << " ( " << ups2spp_amplerr/ups2spp_ampl*100. << "% )" << endl;
662 
663 //------ p+p vs. pT -------------------------------------------------------------------------------
664 
665 TCanvas* dummy2 = new TCanvas("dummy2","dummy2",0,0,500,500);
666 
667  for(int i=0; i<npts; i++) {
668 
669  double tmpmax = hhall_pp[i]->GetMaximum()/9.;
670 
671  fSandBpp->SetParameter(0,tmpmax);
672  fSandBpp->SetParLimits(0,0.,99999.);
673  fSandBpp->FixParameter(1,tonypar1);
674  fSandBpp->FixParameter(2,tonypar2);
675  fSandBpp->FixParameter(3,tonypar3);
676  fSandBpp->FixParameter(4,0.089); // p+p mass resolution
677  fSandBpp->SetParameter(5,tmpmax);
678  fSandBpp->SetParLimits(5,0.0,9999.);
679  fSandBpp->SetParameter(6,tmpmax);
680  fSandBpp->SetParLimits(6,0.0,999.);
681  double tmppar7 = bgpar0pp[i]; if(tmppar7<0.) {tmppar7=0.;}
682  fSandBpp->SetParameter(7,tmppar7);
683  fSandBpp->SetParLimits(7,0.,99.);
684  fSandBpp->FixParameter(8,bgpar1pp[i]);
685  hhall_pp[i]->Fit(fSandBpp,"qrl","",8.,11.);
686 
687  sprintf(tmpname,"fcorrbg_pp_%d",i); // "true" correlated bg copy
688  fcorrbg_pp[i] = new TF1(tmpname,"exp([0]+[1]*x)",7.,14.);
689  fcorrbg_pp[i]->SetLineColor(kRed);
690  fcorrbg_pp[i]->SetLineWidth(1);
691  fcorrbg_pp[i]->SetParameters(bgpar0pp[i], bgpar1pp[i]);
692 
693  sprintf(tmpname,"fcorrbg_pp_fromfit_%d",i); // correlated bg from fit
694  fcorrbg_pp_fromfit[i] = new TF1(tmpname,"exp([0]+[1]*x)",7.,14.);
695  fcorrbg_pp_fromfit[i]->SetLineStyle(2);
696  fcorrbg_pp_fromfit[i]->SetParameters(hhall_pp[i]->GetFunction("fSandBpp")->GetParameter(7), hhall_pp[i]->GetFunction("fSandBpp")->GetParameter(8));
697 
698  fCBups1spp->SetParameter(0,fSandBpp->GetParameter(0));
699  fCBups1spp->SetParameter(1,fSandBpp->GetParameter(1));
700  fCBups1spp->SetParameter(2,fSandBpp->GetParameter(2));
701  fCBups1spp->SetParameter(3,fSandBpp->GetParameter(3));
702  fCBups1spp->SetParameter(4,fSandBpp->GetParameter(4));
703  double ups1spp_integral = fCBups1spp->Integral(5.,14.);
704  double ups1spp_ampl = fSandBpp->GetParameter(0);
705  double ups1spp_amplerr = fSandBpp->GetParError(0);
706  fCBups2spp->SetParameter(0,fSandBpp->GetParameter(5));
707  fCBups2spp->SetParameter(1,fSandBpp->GetParameter(1));
708  fCBups2spp->SetParameter(2,fSandBpp->GetParameter(2));
709  fCBups2spp->SetParameter(3,fSandBpp->GetParameter(3));
710  fCBups2spp->SetParameter(4,fSandBpp->GetParameter(4));
711  double ups2spp_integral = fCBups2spp->Integral(5.,14.);
712  double ups2spp_ampl = fSandBpp->GetParameter(5);
713  double ups2spp_amplerr = fSandBpp->GetParError(5);
714  cout << "p+p Integral1 vs pT = " << i << " " << ups1spp_integral/binsize << " +- " << ups1spp_integral*(ups1spp_amplerr/ups1spp_ampl)/binsize << " ( " << ups1spp_amplerr/ups1spp_ampl*100. << "% )" << endl;
715 
716 
717  sumfit1pp[i] = ups1spp_integral/binsize;
718  ersumfit1pp[i] = ups1spp_integral*(ups1spp_amplerr/ups1spp_ampl)/binsize;
719  sumfit2pp[i] = ups2spp_integral/binsize;
720  ersumfit2pp[i] = ups2spp_integral*(ups2spp_amplerr/ups2spp_ampl)/binsize;
721 
722  truesum1pp[i]=0.;
723  truesum2pp[i]=0.;
724  sum1ppbg[i]=0.;
725  ersum1ppbg[i]=0.;
726  sum2ppbg[i]=0.;
727  ersum2ppbg[i]=0.;
728  for(int j=fbin1; j<=lbin1; j++) {
729  //sum1ppbg[i] += (hhall_pp[i]->GetBinContent(j) - fcorrbg_pp_fromfit[i]->Eval(hhall_pp[i]->GetBinCenter(j)));
730  sum1ppbg[i] += hhall_pp[i]->GetBinContent(j) - fcorrbg_pp[i]->Eval(hhall_pp[i]->GetBinCenter(j));
731  ersum1ppbg[i] += hhall_pp[i]->GetBinError(j)*hhall_pp[i]->GetBinError(j);
732  truesum1pp[i] += hhups1pp[i]->GetBinContent(j);
733  }
734  ersum1ppbg[i] = sqrt(ersum1ppbg[i]);
735  for(int j=fbin2; j<=lbin2; j++) {
736  sum2ppbg[i] += hhall_pp[i]->GetBinContent(j) - fcorrbg_pp[i]->Eval(hhall_pp[i]->GetBinCenter(j));
737  ersum2ppbg[i] += hhall_pp[i]->GetBinError(j)*hhall_pp[i]->GetBinError(j);
738  truesum2pp[i] += hhups2pp[i]->GetBinContent(j);
739  }
740  ersum2ppbg[i] = sqrt(ersum2ppbg[i]);
741 
742 
743  } // end loop over p+p pT bins
744 
745 delete dummy;
746 
747 TCanvas* cupsvspt_pp = new TCanvas("cupsvspt_pp","p+p Upsilons vs. p_{T}",100,100,1100,900);
748 cupsvspt_pp->Divide(3,3);
749 for(int i=0; i<9; i++) {
750 if(i==0) {cupsvspt_pp->cd(0);}
751 if(i==1) {cupsvspt_pp->cd(1);}
752 if(i==2) {cupsvspt_pp->cd(2);}
753 if(i==3) {cupsvspt_pp->cd(3);}
754 if(i==4) {cupsvspt_pp->cd(4);}
755 if(i==5) {cupsvspt_pp->cd(5);}
756 if(i==6) {cupsvspt_pp->cd(6);}
757 if(i==7) {cupsvspt_pp->cd(7);}
758 if(i==8) {cupsvspt_pp->cd(8);}
759  hhall_pp[i]->SetAxisRange(8.0,11.0);
760  hhall_pp[i]->Draw();
761  fcorrbg_pp[i]->Draw("same");
762  fcorrbg_pp_fromfit[i]->Draw("same");
763  sprintf(tlchar,"%d-%d GeV",i,i+1); tl[i] = new TLatex(8.2,hhall_pp[i]->GetMaximum()*0.95,tlchar); tl[i]->Draw();
764 }
765 
766 //
767 //=== compare errors from direct counting and fit
768 //
769 
770 TCanvas* cercomp = new TCanvas("cercomp","error comparizon",50,50,700,700);
771 TH2F* hhcomp = new TH2F("hhcomp"," ",10,0.,10.,10,0.,1.0);
772 hhcomp->GetXaxis()->SetTitle("Transverse momentum [GeV/c]");
773 hhcomp->GetXaxis()->SetTitleOffset(1.0);
774 hhcomp->Draw();
775 
776 double reler1[nbins+1],reler2[nbins+1];
777 double relerfit1[nbins+1],relerfit2[nbins+1];
778 double reler1pp[nbins+1],reler2pp[nbins+1];
779 double reler1ppbg[nbins+1],reler2ppbg[nbins+1];
780 double relerfit1pp[nbins+1],relerfit2pp[nbins+1];
781 for(int i=0; i<npts; i++) {
782  reler1[i] = ersum1[i]/truesum1[i]; reler2[i] = ersum2[i]/truesum2[i];
783  relerfit1[i] = ersumfit1[i]/Nups1[i]; relerfit2[i] = ersumfit2[i]/Nups2[i];
784  reler1pp[i] = ersum1pp[i]/truesum1pp[i]; reler2pp[i] = ersum2pp[i]/truesum2pp[i]; // direct count no bg
785  reler1ppbg[i] = ersum1ppbg[i]/truesum1pp[i]; reler2ppbg[i] = ersum2ppbg[i]/truesum2pp[i]; // direct count with bg
786  relerfit1pp[i] = ersumfit1pp[i]/Nups1pp[i]; relerfit2pp[i] = ersumfit2pp[i]/Nups2pp[i]; // fit
787 }
788 
789 TGraphErrors* grcomp1 = new TGraphErrors(npts,xx1,reler1,0,0);
790 grcomp1->SetMarkerStyle(20);
791 grcomp1->SetMarkerColor(kBlack);
792 grcomp1->SetLineColor(kBlack);
793 grcomp1->SetLineWidth(2);
794 grcomp1->SetMarkerSize(1.5);
795 grcomp1->Draw("p");
796 
797 TGraphErrors* grcompfit1 = new TGraphErrors(npts,xx1,relerfit1,0,0);
798 grcompfit1->SetMarkerStyle(24);
799 grcompfit1->SetMarkerColor(kBlack);
800 grcompfit1->SetLineColor(kBlack);
801 grcompfit1->SetLineWidth(2);
802 grcompfit1->SetMarkerSize(1.5);
803 grcompfit1->Draw("p");
804 
805 TCanvas* cercomppp = new TCanvas("cercomppp","p+p error comparizon",50,50,700,700);
806 TH2F* hhcomppp = new TH2F("hhcomppp"," ",10,0.,10.,10,0.,1.0);
807 hhcomppp->GetXaxis()->SetTitle("Transverse momentum [GeV/c]");
808 hhcomppp->GetXaxis()->SetTitleOffset(1.0);
809 hhcomppp->Draw();
810 
811 TGraphErrors* grcomp1pp = new TGraphErrors(npts,xx1,reler1pp,0,0);
812 grcomp1pp->SetMarkerStyle(20);
813 grcomp1pp->SetMarkerColor(kBlue);
814 grcomp1pp->SetLineColor(kBlue);
815 grcomp1pp->SetLineWidth(2);
816 grcomp1pp->SetMarkerSize(1.5);
817 grcomp1pp->Draw("p");
818 
819 TGraphErrors* grcomp1ppbg = new TGraphErrors(npts,xx1,reler1ppbg,0,0);
820 grcomp1ppbg->SetMarkerStyle(21);
821 grcomp1ppbg->SetMarkerColor(kBlue);
822 grcomp1ppbg->SetLineColor(kBlue);
823 grcomp1ppbg->SetLineWidth(2);
824 grcomp1ppbg->SetMarkerSize(1.5);
825 grcomp1ppbg->Draw("p");
826 
827 TGraphErrors* grcompfit1pp = new TGraphErrors(npts,xx1,relerfit1pp,0,0);
828 grcompfit1pp->SetMarkerStyle(24);
829 grcompfit1pp->SetMarkerColor(kBlue);
830 grcompfit1pp->SetLineColor(kBlue);
831 grcompfit1pp->SetLineWidth(2);
832 grcompfit1pp->SetMarkerSize(1.5);
833 grcompfit1pp->Draw("p");
834 
835 
836 //==========================================================================
837 // RAA
838 //==========================================================================
839 
840 double raa1[nbins+1],raa2[nbins+1],raa3[nbins+1],erraa1[nbins+1],erraa2[nbins+1],erraa3[nbins+1];
841 double raa1fit[nbins+1],raa2fit[nbins+1],raa3fit[nbins+1],erraa1fit[nbins+1],erraa2fit[nbins+1],erraa3fit[nbins+1];
842 
843 for(int i=0; i<nbins; i++) {
844 
845  raa1[i] = 0.535 * grRAA1S->Eval(0.5+i*1.0)/groldRAA1S->Eval(0.5+i*1.0);
846  raa2[i] = 0.170 * grRAA2S->Eval(0.5+i*1.0)/groldRAA2S->Eval(0.5+i*1.0);
847  raa3[i] = 0.035;
848  raa1fit[i] = 0.535 * grRAA1S->Eval(0.5+i*1.0)/groldRAA1S->Eval(0.5+i*1.0);
849  raa2fit[i] = 0.170 * grRAA2S->Eval(0.5+i*1.0)/groldRAA2S->Eval(0.5+i*1.0);
850  raa3fit[i] = 0.035;
851 
852  erraa1[i] = raa1[i]*sqrt(pow(ersum1[i]/truesum1[i],2)+pow(ersum1pp[i]/truesum1pp[i],2));
853  erraa2[i] = raa2[i]*sqrt(pow(ersum2[i]/truesum2[i],2)+pow(ersum2pp[i]/truesum2pp[i],2));
854  erraa1fit[i] = raa1fit[i]*sqrt(pow(ersumfit1[i]/Nups1[i],2)+pow(ersumfit1pp[i]/Nups1pp[i],2));
855  erraa2fit[i] = raa2fit[i]*sqrt(pow(ersumfit2[i]/Nups2[i],2)+pow(ersumfit2pp[i]/Nups2pp[i],2));
856 
857 }
858 
859 // plot icomparison direct counting vs. fit
860 
861 TCanvas* craa = new TCanvas("craa","R_{AA}",120,120,600,600);
862 TH2F* hh2 = new TH2F("hh2"," ",10,0.,10.,10,0.,0.8);
863 hh2->GetXaxis()->SetTitle("Transverse momentum [GeV/c]");
864 hh2->GetXaxis()->SetTitleOffset(1.0);
865 hh2->GetXaxis()->SetTitleColor(1);
866 hh2->GetXaxis()->SetTitleSize(0.040);
867 hh2->GetXaxis()->SetLabelSize(0.040);
868 hh2->GetYaxis()->SetTitle("R_{AA}");
869 hh2->GetYaxis()->SetTitleOffset(1.3);
870 hh2->GetYaxis()->SetTitleSize(0.040);
871 hh2->GetYaxis()->SetLabelSize(0.040);
872 hh2->Draw();
873 
874 TGraphErrors* gr1 = new TGraphErrors(npts,xx2,raa1,0,erraa1);
875 gr1->SetMarkerStyle(20);
876 gr1->SetMarkerColor(kBlack);
877 gr1->SetLineColor(kBlack);
878 gr1->SetLineWidth(2);
879 gr1->SetMarkerSize(1.5);
880 gr1->Draw("p");
881 
882 TGraphErrors* gr1fit = new TGraphErrors(npts,xx3,raa1fit,0,erraa1fit);
883 gr1fit->SetMarkerStyle(24);
884 gr1fit->SetMarkerColor(kBlack);
885 gr1fit->SetLineColor(kBlack);
886 gr1fit->SetLineWidth(2);
887 gr1fit->SetMarkerSize(1.5);
888 gr1fit->Draw("p");
889 
890 TGraphErrors* gr2 = new TGraphErrors(npts-2,xx2,raa2,0,erraa2);
891 gr2->SetMarkerStyle(20);
892 gr2->SetMarkerColor(kRed);
893 gr2->SetLineColor(kRed);
894 gr2->SetLineWidth(2);
895 gr2->SetMarkerSize(1.5);
896 gr2->Draw("p");
897 
898 TGraphErrors* gr2fit = new TGraphErrors(npts-2,xx3,raa2fit,0,erraa2fit);
899 gr2fit->SetMarkerStyle(24);
900 gr2fit->SetMarkerColor(kRed);
901 gr2fit->SetLineColor(kRed);
902 gr2fit->SetLineWidth(2);
903 gr2fit->SetMarkerSize(1.5);
904 gr2fit->Draw("p");
905 
906 delete dummy;
907 
908 // plot fit
909 
910 TCanvas* craafit = new TCanvas("craafit","R_{AA}",20,20,600,600);
911 TH2F* hh2f = new TH2F("hh2f"," ",10,0.,10.,10,0.,0.8);
912 hh2f->GetXaxis()->SetTitle("Transverse momentum [GeV/c]");
913 hh2f->GetXaxis()->SetTitleOffset(1.0);
914 hh2f->GetXaxis()->SetTitleColor(1);
915 hh2f->GetXaxis()->SetTitleSize(0.040);
916 hh2f->GetXaxis()->SetLabelSize(0.040);
917 hh2f->GetYaxis()->SetTitle("R_{AA}");
918 hh2f->GetYaxis()->SetTitleOffset(1.3);
919 hh2f->GetYaxis()->SetTitleSize(0.040);
920 hh2f->GetYaxis()->SetLabelSize(0.040);
921 hh2f->Draw();
922 
923 TGraphErrors* gr1f = new TGraphErrors(npts,xx1,raa1fit,0,erraa1fit);
924 gr1f->SetMarkerStyle(20);
925 gr1f->SetMarkerColor(kBlack);
926 gr1f->SetLineColor(kBlack);
927 gr1f->SetLineWidth(2);
928 gr1f->SetMarkerSize(1.5);
929 gr1f->SetName("gr2c");
930 gr1f->Draw("p");
931 
932 TGraphErrors* gr2f = new TGraphErrors(npts-2,xx1,raa2fit,0,erraa2fit);
933 gr2f->SetMarkerStyle(20);
934 gr2f->SetMarkerColor(kRed);
935 gr2f->SetLineColor(kRed);
936 gr2f->SetLineWidth(2);
937 gr2f->SetMarkerSize(1.5);
938 gr2f->SetName("gr2f");
939 gr2f->Draw("p");
940 
941 TLegend *leg = new TLegend(0.2,0.77,0.5,0.92);
942 leg->SetFillColor(10);
943 leg->SetFillStyle(1001);
944 //TLegendEntry *entry=leg->AddEntry("gr1f","Y(1S)","p");
945 //entry=leg->AddEntry("gr2f","Y(2S)","p");
946 //entry=leg->AddEntry("","","");
947 leg->AddEntry("gr1f","Y(1S)","p");
948 leg->AddEntry("gr2f","Y(2S)","p");
949 leg->AddEntry("","","");
950 leg->Draw();
951 TLatex* l1 = new TLatex(3.5,0.73,"sPHENIX Simulation");
952 l1->Draw();
953 
954 grRAA1S->Draw("l");
955 grRAA2S->Draw("l");
956 
957 // plot direct count
958 
959 TCanvas* craacount = new TCanvas("craacount","R_{AA}",20,20,600,600);
960 TH2F* hh2c = new TH2F("hh2c"," ",10,0.,10.,10,0.,0.8);
961 hh2c->GetXaxis()->SetTitle("Transverse momentum [GeV/c]");
962 hh2c->GetXaxis()->SetTitleOffset(1.0);
963 hh2c->GetXaxis()->SetTitleColor(1);
964 hh2c->GetXaxis()->SetTitleSize(0.040);
965 hh2c->GetXaxis()->SetLabelSize(0.040);
966 hh2c->GetYaxis()->SetTitle("R_{AA}");
967 hh2c->GetYaxis()->SetTitleOffset(1.3);
968 hh2c->GetYaxis()->SetTitleSize(0.040);
969 hh2c->GetYaxis()->SetLabelSize(0.040);
970 hh2c->Draw();
971 
972 TGraphErrors* gr1c = new TGraphErrors(npts,xx1,raa1,0,erraa1);
973 gr1c->SetMarkerStyle(20);
974 gr1c->SetMarkerColor(kBlack);
975 gr1c->SetLineColor(kBlack);
976 gr1c->SetLineWidth(2);
977 gr1c->SetMarkerSize(1.5);
978 gr1c->SetName("gr1c");
979 gr1c->Draw("p");
980 
981 TGraphErrors* gr2c = new TGraphErrors(npts-3,xx1,raa2,0,erraa2);
982 gr2c->SetMarkerStyle(20);
983 gr2c->SetMarkerColor(kRed);
984 gr2c->SetLineColor(kRed);
985 gr2c->SetLineWidth(2);
986 gr2c->SetMarkerSize(1.5);
987 gr2c->SetName("gr2c");
988 gr2c->Draw("p");
989 
990 TLegend *legc = new TLegend(0.2,0.77,0.5,0.92);
991 legc->SetFillColor(10);
992 legc->SetFillStyle(1001);
993 legc->AddEntry("gr1c","Y(1S)","p");
994 legc->AddEntry("gr2c","Y(2S)","p");
995 legc->AddEntry("","","");
996 //TLegendEntry *entry2=legc->AddEntry("gr1c","Y(1S)","p");
997 //entry2=legc->AddEntry("gr2c","Y(2S)","p");
998 //entry2=legc->AddEntry("","","");
999 legc->Draw();
1000 TLatex* l1c = new TLatex(3.5,0.73,"sPHENIX Simulation");
1001 l1c->Draw();
1002 
1003 grRAA1S->Draw("l");
1004 grRAA2S->Draw("l");
1005 
1006 }
1007