Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ccbb.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file ccbb.cc
1 
2 #include "TFile.h"
3 #include "TF1.h"
4 #include "TH1.h"
5 #include "TH2.h"
6 #include "TGraph.h"
7 #include "TChain.h"
8 #include "TRandom3.h"
9 #include "TGraphAsymmErrors.h"
10 #include "TGraphErrors.h"
11 
12 #include <fstream>
13 #include <iostream>
14 //#include <iomanip>
15 //#include <vector>
16 
17 using namespace std;
18 
19 int ccbb(double eideff, string ofname);
20 
21 TH1D* hcharm_pt;
22 TH1D* hbottom_pt;
23 TH1D* hhf_pt;
24 TH1D* hcharm_pt_nosupp;
25 TH1D* hbottom_pt_nosupp;
26 TH1D* hhf_pt_nosupp;
27 
28 double norm_charm(double* x, double* par) {
29  double pt = x[0];
30  double norm = par[0];
31  int bin = hcharm_pt->GetXaxis()->FindBin(pt);
32  double value = hcharm_pt->GetBinContent(bin);
33  return norm*value;
34 }
35 double norm_charm_nosupp(double* x, double* par) {
36  double pt = x[0];
37  double norm = par[0];
38  int bin = hcharm_pt_nosupp->GetXaxis()->FindBin(pt);
39  double value = hcharm_pt_nosupp->GetBinContent(bin);
40  return norm*value;
41 }
42 
43 double norm_bottom(double* x, double* par) {
44  double pt = x[0];
45  double norm = par[0];
46  int bin = hbottom_pt->GetXaxis()->FindBin(pt);
47  double value = hbottom_pt->GetBinContent(bin);
48  return norm*value;
49 }
50 double norm_bottom_nosupp(double* x, double* par) {
51  double pt = x[0];
52  double norm = par[0];
53  int bin = hbottom_pt_nosupp->GetXaxis()->FindBin(pt);
54  double value = hbottom_pt_nosupp->GetBinContent(bin);
55  return norm*value;
56 }
57 
58 double norm_hf(double* x, double* par) {
59  double pt = x[0];
60  double norm = par[0];
61  int bin = hhf_pt->GetXaxis()->FindBin(pt);
62  double value = hhf_pt->GetBinContent(bin);
63  return norm*value;
64 }
65 double norm_hf_nosupp(double* x, double* par) {
66  double pt = x[0];
67  double norm = par[0];
68  int bin = hhf_pt_nosupp->GetXaxis()->FindBin(pt);
69  double value = hhf_pt_nosupp->GetBinContent(bin);
70  return norm*value;
71 }
72 
73 
74 
75 //====================================================================================
76 
77 int main(int argc, char* argv[])
78 {
79  double eideff = 0.9;
80  string ofname="ccbb.root";
81  if(argc==1) cout << argv[0] << " is running with standard parameters..." << endl;
82  if(argc==2) { eideff = atoi(argv[1]); }
83  if(argc==3) { eideff = atoi(argv[1]); ofname=argv[2]; }
84  cout << "eID efficiency = " << eideff << endl;
85  cout << "output file name: " << ofname << endl;
86  return ccbb(eideff, ofname);
87 }
88 
89 //------------------------------------------------------------------------------------
90 
91 int ccbb(double eideff, string ofname) {
92 
93 // 0 = ppg182, min. bias AuAu
94 // 1 = FONLL, p+p
95 // 2 = my own function
96 // 3 = Kazuya, 0-10% central AuAu
97 int which_botfrac = 3;
98 
99 // 0 = Min.Bias from ppg182
100 // 1 = 0-10% central from Kazuya
101 // 2 = no suppression
102 // 3 = my own suppression
103 int which_supp = 2;
104 
105 double ptcut = 2.0;
106 
107 double NeventsAuAu = 10.0e+09;
108 
109 TRandom* rnd = new TRandom3();
110 rnd->SetSeed(0);
111 
112 const int nbins = 15;
113 double binlim[nbins+1];
114 for(int i=0; i<=nbins; i++) {binlim[i]=double(i);}
115 
116 double startmass=0.;
117 double stopmass=20.;
118 int nchanmass = 400;
119 double startpt=0.;
120 double stoppt=10.;
121 int nchanpt = 100;
122 double bsizept = (stoppt-startpt)/float(nchanpt);
123 
124 char hname[999];
125 
126 double erpt[99], erpt2[99]; for(int i=0; i<99; i++) {erpt[i]=0.075; erpt2[i]=0.1;}
127 std::string tmp0,tmp1,tmp2,tmp3,tmp4,tmp5;
128 
129 // last histogram is integrated over pT
130 TH1D* hcharm[nbins+1];
131 TH1D* hcharm_ckin3_4[nbins+1];
132 TH1D* hbottom[nbins+1];
133 TH1D* hdy[nbins+1];
134 for(int i=0; i<nbins+1; i++) { sprintf(hname,"hcharm_%d",i); hcharm[i] = new TH1D(hname,"",nchanmass,startmass,stopmass); hcharm[i]->Sumw2(); }
135 for(int i=0; i<nbins+1; i++) { sprintf(hname,"hcharm_ckin3_3_%d",i); hcharm_ckin3_4[i] = new TH1D(hname,"",nchanmass,startmass,stopmass); hcharm_ckin3_4[i]->Sumw2(); }
136 for(int i=0; i<nbins+1; i++) { sprintf(hname,"hbottom_%d",i); hbottom[i] = new TH1D(hname,"",nchanmass,startmass,stopmass); hbottom[i]->Sumw2(); }
137 for(int i=0; i<nbins+1; i++) { sprintf(hname,"hdy_%d",i); hdy[i] = new TH1D(hname,"",nchanmass,startmass,stopmass); hdy[i]->Sumw2();}
138 hcharm_pt = new TH1D("hcharm_pt", "",nchanpt,startpt,stoppt);
139 hbottom_pt = new TH1D("hbottom_pt","",nchanpt,startpt,stoppt);
140 hcharm_pt->Sumw2();
141 hbottom_pt->Sumw2();
142 hcharm_pt_nosupp = new TH1D("hcharm_pt_nosupp", "",nchanpt,startpt,stoppt);
143 hbottom_pt_nosupp = new TH1D("hbottom_pt_nosupp","",nchanpt,startpt,stoppt);
144 hcharm_pt_nosupp->Sumw2();
145 hbottom_pt_nosupp->Sumw2();
146 
147 float mass,pt,eta,pt1,pt2,eta1,eta2;
148 float pt_sngl,eta_sngl;
149 
150 float fitstart=2.0;
151 float fitstop=8.0;
152 TF1* fnorm_charm = new TF1("fnorm_charm", norm_charm, fitstart,fitstop,1);
153 TF1* fnorm_bottom = new TF1("fnorm_bottom",norm_bottom,fitstart,fitstop,1);
154 TF1* fnorm_hf = new TF1("fnorm_hf" , norm_hf, fitstart,fitstop,1);
155 TF1* fnorm_charm_nosupp = new TF1("fnorm_charm_nosupp", norm_charm_nosupp, fitstart,fitstop,1);
156 TF1* fnorm_bottom_nosupp = new TF1("fnorm_bottom_nosupp",norm_bottom_nosupp,fitstart,fitstop,1);
157 TF1* fnorm_hf_nosupp = new TF1("fnorm_hf_nosupp" , norm_hf_nosupp, fitstart,fitstop,1);
158 
159 /*
160 //--------------------------------------------------------------------------------------
161 //--------------------------------------------------------------------------------------
162 TFile *fFONLL=new TFile("FONLL_ratio.root");
163  TGraph *gYieldFONLL_c = (TGraph*)fFONLL->Get("gRatio_c")->Clone();
164  TGraph *gYieldFONLL_b = (TGraph*)fFONLL->Get("gRatio_b")->Clone();
165  int Ndata=gYieldFONLL_c->GetN();
166  double *Xc=gYieldFONLL_c->GetX();
167  double *Yc=gYieldFONLL_c->GetY();
168  double *Xb=gYieldFONLL_b->GetX();
169  double *Yb=gYieldFONLL_b->GetY();
170  double XR[999], YR[999], ZR[999];
171  for(int i=0; i<Ndata; i++) {XR[i]=Xc[i]; YR[i]=Yb[i]/(Yb[i]+Yc[i]); ZR[i]=Yc[i]/(Yb[i]+Yc[i]);}
172  TGraph *gFracFONLL_b = new TGraph(Ndata,XR,YR);
173  TGraph *gFracFONLL_c = new TGraph(Ndata,XR,ZR);
174 fFONLL->Close();
175 */
176 //
177 // bottom/(charm+bottom) and charm/(charm+bottom) ratios in p+p from FONLL
178 // from Darren (ppg182)
179 //
180 ifstream fin_fonll;
181 double pt_fonll[99],fonll[99],fonllup[99],fonll_erup[99],fonll_erdn[99],dtmp0,dtmp1;
182 fin_fonll.open("FONLL_ratio.txt");
183 if(!fin_fonll) {cerr << " ERROR: Cannot open input text file. " << endl;}
184 for(int i=0; i<19; i++) {
185  fin_fonll >> pt_fonll[i] >> tmp0 >> fonll[i] >> tmp1 >> dtmp0 >> tmp2 >> dtmp1 >> tmp3;
186  fonllup[i] = dtmp0;
187  fonll_erup[i] = dtmp0 - fonll[i];
188  fonll_erdn[i] = fonll[i] - dtmp1;
189  cout << pt_fonll[i] << " " << fonll[i] << " " << fonll_erup[i] << " " << fonll_erdn[i] << endl;
190 }
191 fin_fonll.close();
192  TGraphAsymmErrors *gfonll = new TGraphAsymmErrors(19,pt_fonll,fonll,0,0,fonll_erdn,fonll_erup);
193  TGraph *gfonllup = new TGraph(19,pt_fonll,fonllup);
194 
195 
196 //----------------------------------------------------------------------
197 // my bottom/(charm+bottom) fraction just for testing
198 //----------------------------------------------------------------------
199 TF1* mybotfrac = new TF1("mybotfrac","(0.896512*pow(x+(1.050669),5.792579)/(2630.566406+pow(x+(1.050669),5.792579)))",0.,10.); // my FONLL-like function
200 //TF1* mybotfrac = new TF1("mybotfrac","0.98*pow(x-1.0,3)/(2+pow(x-1.0,3))",1.,10.); // myfunction close to Kazuya's upper limit
201 double myx[999];
202 double myy[999];
203 double pyy[999];
204 int myn = 150;
205 for(int i=0; i<myn; i++) { myx[i] = i*0.1; myy[i] = mybotfrac->Eval(myx[i]); }
206 //for(int i=0; i<myn; i++) { myx[i] = i*0.1; myy[i] = gfonll->Eval(myx[i])*mybotfrac->Eval(myx[i]); }
207 TGraph* gr_mybotfrac = new TGraph(myn,myx,myy);
208 
209 TF1* pythia_botfrac = new TF1("pythia_botfrac","(0.898794*pow(x+(1.035767),5.749620)/(2453.650635+pow(x+(1.035767),5.749620)))*1.0",0.,10.); // PYTHIA
210 for(int i=0; i<myn; i++) { pyy[i] = pythia_botfrac->Eval(myx[i]); }
211 TGraph* gr_pythia_botfrac = new TGraph(myn,myx,pyy);
212 
213 //--------------------------------------------------------------------------------------
214 // bottom/(charm+bottom) ratio in Min. Bias Au+Au from ppg182:
215 // arXiv:1509.04662
216 // A. Adare et al. (PHENIX Collaboration) Phys. Rev. C 93, 034904 (2016)
217 //--------------------------------------------------------------------------------------
218 double ratcb[99],ratcberup[99],ratcberdn[99],ratcbpt[99];
219 for(int i=0; i<15; i++) { ratcbpt[i] = 1.1+0.2*i; }
220 ratcbpt[15] = 4.25;
221 ratcbpt[16] = 4.75;
222 ratcbpt[17] = 5.5;
223 ratcbpt[18] = 6.5;
224 ratcbpt[19] = 7.5;
225 ratcbpt[20] = 8.5;
226 ratcb[0] = 3.99e-02;
227 ratcb[1] = 6.15e-02;
228 ratcb[2] = 9.41e-02;
229 ratcb[3] = 1.40e-01;
230 ratcb[4] = 2.01e-01;
231 ratcb[5] = 2.73e-01;
232 ratcb[6] = 3.49e-01;
233 ratcb[7] = 4.22e-01;
234 ratcb[8] = 4.84e-01;
235 ratcb[9] = 5.32e-01;
236 ratcb[10] = 5.64e-01;
237 ratcb[11] = 5.81e-01;
238 ratcb[12] = 5.89e-01;
239 ratcb[13] = 5.87e-01;
240 ratcb[14] = 5.79e-01;
241 ratcb[15] = 5.57e-01;
242 ratcb[16] = 5.13e-01;
243 ratcb[17] = 4.79e-01;
244 ratcb[18] = 4.71e-01;
245 ratcb[19] = 4.61e-01;
246 ratcb[20] = 4.36e-01;
247  TGraph *gFrac_b = new TGraph(21,ratcbpt,ratcb);
248 
249 double fracc[99];
250 for(int i=0; i<21; i++) { fracc[i] = 1.0-ratcb[i]; }
251  TGraph *gFrac_c = new TGraph(21,ratcbpt,fracc);
252 
253 
254 //-------------------------------------------------------------------------------
255 // bottom fraction from Kazuya Nagashima's talk at QM2017
256 // PHENIX preliminary for 0-10% central Au+Au
257 // http://www.phenix.bnl.gov/phenix/WWW/p/draft/nagasy/QM2017/slide/PHENIX_HF_KazuyaNagashima.pdf
258 //-------------------------------------------------------------------------------
259  TGraphAsymmErrors* bfrac_cent1;
260  TFile *fk = new TFile("kazuya/bfrac_cent1.root");
261  bfrac_cent1 = (TGraphAsymmErrors*)fk->Get("bfrac_cent1");
262  fk->Close();
263 
264 //-----------------------------------------------------------------------------------------
265 // R_AA of HF electrons according to Kazuya Nagashima; 0-10% central; 2004+2014 data
266 // Different from ppg077
267 // Probably because it's combined 2004 and 2014 results (ppg077 = 2004 only)
268 //-----------------------------------------------------------------------------------------
269  ifstream kfile;
270  kfile.open("kazuya/csv/ppg077_HFRAA_cent1.csv");
271  TGraphErrors *tRaadata = new TGraphErrors();
272  TGraphErrors *tRaadata_sys = new TGraphErrors();
273  double var1, var2, var3, var4, var5, var6;
274  for(int i=0;i<28;i++) {
275  kfile >> var1 >> var2 >> var3 >> var4 >> var5 >> var6;
276  tRaadata->SetPoint(i,var1,var2);
277  tRaadata->SetPointError(i,0.,var3);
278  tRaadata_sys->SetPoint(i,var1,var2);
279  tRaadata_sys->SetPointError(i,0.1,var5);
280  }
281  kfile.close();
282 
283 //--------------------------------------------------------------------------------------
284 // Charm and Bottom R_AA in 0-10% central Au+Au from Kazuya Nagashima's QM2017 presentation
285 //--------------------------------------------------------------------------------------
286  TGraphAsymmErrors* tcharm;
287  TGraphAsymmErrors* tbottom;
288  TFile *fk2 = new TFile("kazuya/Raa_cent1.root");
289  tcharm = (TGraphAsymmErrors*)fk2->Get("tcharm");
290  tbottom = (TGraphAsymmErrors*)fk2->Get("tbottom");
291  fk2->Close();
292 
293  TGraph* mytcharm = new TGraph();
294  TGraph* mytbottom = new TGraph();
295  mytcharm->SetPoint(0,0.0,1.0);
296  mytcharm->SetPoint(1,1.5,1.0);
297  mytbottom->SetPoint(0,0.0,1.0);
298  mytbottom->SetPoint(1,1.5,1.0);
299  double *tmpx = tcharm->GetX();
300  double *tmpy = tcharm->GetY();
301  for(int i=2; i<8; i++) { mytcharm->SetPoint(i,tmpx[i-2],tmpy[i-2]); }
302  mytcharm->SetPoint(8,15.,tmpy[5]);
303  tmpx = tbottom->GetX();
304  tmpy = tbottom->GetY();
305  for(int i=2; i<8; i++) { mytbottom->SetPoint(i,tmpx[i-2],tmpy[i-2]); }
306  mytbottom->SetPoint(2,tmpx[0],1.0);
307  mytbottom->SetPoint(8,15.,tmpy[5]);
308 
309 
310 //--------------------------------------------------------------------------------------
311 // Charm and Bottom R_AA in Min. Bias Au+Au from ppg182:
312 // arXiv:1509.04662
313 // A. Adare et al. (PHENIX Collaboration) Phys. Rev. C 93, 034904 (2016)
314 //--------------------------------------------------------------------------------------
315 
316 double pt182[9],raa182charm[9],raa182bottom[9],raa182charm_erdn[9],raa182charm_erup[9],raa182bottom_erdn[9],raa182bottom_erup[9];
317 double pt182orig[9],raa182charmorig[9],raa182bottomorig[9];
318 pt182[0]=0.0;
319 pt182[1]=1.5;
320 raa182charm[0] = 1.0;
321 raa182charm[1] = 1.0;
322 raa182bottom[0] = 1.0;
323 raa182bottom[1] = 1.0;
324 
325 ifstream fin182;
326 fin182.open("ppg182.txt");
327 if(!fin182) {cerr << " ERROR: Cannot open input text file. " << endl;}
328 fin182 >> tmp0 >> tmp1 >> tmp2 >> tmp3;
329 fin182 >> tmp0 >> tmp1 >> tmp2 >> tmp3;
330 for(int i=2; i<8; i++) {
331  fin182 >> pt182[i] >> raa182charm[i] >> raa182charm_erup[i-2] >> raa182charm_erdn[i-2];
332  pt182orig[i-2]=pt182[i];
333  raa182charmorig[i-2] = raa182charm[i];
334 }
335 fin182 >> tmp0 >> tmp1 >> tmp2 >> tmp3;
336 fin182 >> tmp0 >> tmp1 >> tmp2 >> tmp3;
337 for(int i=2; i<8; i++) {
338  fin182 >> pt182[i] >> raa182bottom[i] >> raa182bottom_erup[i-2] >> raa182bottom_erdn[i-2];
339  raa182bottomorig[i-2] = raa182bottom[i];
340 }
341 
342 pt182[8]=15.;
343 raa182charm[8] = raa182charm[7];
344 raa182bottom[2] = 1.0; // measured value is 1.15
345 raa182bottom[8] = raa182bottom[7];
346  TGraph *gRAA_charm = new TGraph(9,pt182,raa182charm);
347  TGraph *gRAA_bottom = new TGraph(9,pt182,raa182bottom);
348  TGraphAsymmErrors *gRAASysErr_charm = new TGraphAsymmErrors(6,pt182orig,raa182charmorig, 0,0,raa182charm_erdn, raa182charm_erup);
349  TGraphAsymmErrors *gRAASysErr_bottom = new TGraphAsymmErrors(6,pt182orig,raa182bottomorig,0,0,raa182bottom_erdn,raa182bottom_erup);
350 
351 
352 // My own RAA from comparison with data
353  //TF1* fmyraac = new TF1("fmyraac","0.303395+18.995424/(66.573738+pow(x,7.237049))",1.0,9.0); // fit above 3 GeV
354  //TF1* fmyraab = new TF1("fmyraab","0.347859+64550.024555/(140970.430941+pow(x,9.780939))",1.0,9.0); // fit above 3 Gev
355  TF1* fmyraac = new TF1("fmyraac","0.327075+20.482705/(66.589763+pow(x,7.237551))",1.0,9.0); // fit above 2 GeV
356  TF1* fmyraab = new TF1("fmyraab","0.370667+68699.897905/(140801.861045+pow(x,9.779860))",1.0,9.0); // fit above 2 Gev
357  double mysuppc[999],mysuppb[999],myxx[999];
358  int mycount=0;
359  for(int i=0; i<myn; i++) { if(myx[i]>1.0) { myxx[mycount] = myx[i]; mysuppc[mycount] = fmyraac->Eval(myx[i]); mycount++; } }
360  TGraph* gr_mysupp_charm = new TGraph(mycount,myxx,mysuppc);
361  mycount=0;
362  for(int i=0; i<myn; i++) { if(myx[i]>1.0) { mysuppb[mycount] = fmyraab->Eval(myx[i]); mycount++; } }
363  TGraph* gr_mysupp_bottom = new TGraph(mycount,myxx,mysuppb);
364 
365 
366 // Choose which suppression to use
367 TGraph* gCharmSuppression;
368 TGraph* gBottomSuppression;
369 if(which_supp==0) {
370  gCharmSuppression = new TGraph(*gRAA_charm);
371  gBottomSuppression = new TGraph(*gRAA_bottom);
372 } else if(which_supp==1) {
373  gCharmSuppression = new TGraph(*mytcharm);
374  gBottomSuppression = new TGraph(*mytbottom);
375 } else if(which_supp==2) {
376  double x1[4] = {0.,5.,10.,15.};
377  double y1[4] = {1.,1.,1., 1.};
378  gCharmSuppression = new TGraph(4,x1,y1);
379  gBottomSuppression = new TGraph(4,x1,y1);
380 } else if(which_supp==3) {
381  gCharmSuppression = new TGraph(*gr_mysupp_charm);
382  gBottomSuppression = new TGraph(*gr_mysupp_bottom);
383 } else {
384  cerr << "ERROR: Wrong suppression factors !!!" << endl;
385 }
386 
387 // Choose which Bottom Fraction to use
388 TGraph* gBottomFraction;
389 if(which_botfrac==0) { // ppg182, min. bias AuAu
390  gBottomFraction = new TGraph(*gFrac_b);
391 } else if(which_botfrac==1) { // FONLL p+p
392  //gBottomFraction = new TGraph(*gFracFONLL_b);
393  gBottomFraction = new TGraph(*gfonll);
394 } else if(which_botfrac==2) { // my own function
395  gBottomFraction = new TGraph(*gr_mybotfrac);
396 } else if(which_botfrac==3) { // Kazuya, 0-10% central AuAu
397  gBottomFraction = new TGraph(*bfrac_cent1);
398 }
399 
400 
401 
402 //-----------------------------------------------------------------------------------------
403 // R_AA of HF electrons from ppg077; 0-10% central Au+Au; 2004 data
404 //-----------------------------------------------------------------------------------------
405 
406 double ppg077_pt[99],ppg077_raa[99],ppg077_staterup[99],ppg077_staterdn[99],ppg077_syserup[99],ppg077_syserdn[99];
407 
408 ifstream finppg077;
409 finppg077.open("ppg077_0-10.txt");
410 if(!finppg077) {cerr << " ERROR: Cannot open input text file. " << endl;}
411 finppg077 >> tmp0 >> tmp1 >> tmp2 >> tmp3 >> tmp4 >> tmp5;
412 finppg077 >> tmp0 >> tmp1;
413 finppg077 >> tmp0 >> tmp1 >> tmp2 >> tmp3 >> tmp4 >> tmp5;
414 for(unsigned int i=0; i<28; i++) { finppg077 >> ppg077_pt[i] >> ppg077_raa[i] >> ppg077_staterup[i] >> ppg077_staterdn[i] >> ppg077_syserup[i] >> ppg077_syserdn[i]; }
415 
416 TGraphAsymmErrors* grRAA_ppg077 = new TGraphAsymmErrors(28,ppg077_pt,ppg077_raa,0,0,ppg077_staterdn,ppg077_staterup);
417 TGraphAsymmErrors* grRAA_ppg077SysErr = new TGraphAsymmErrors(28,ppg077_pt,ppg077_raa,erpt2,erpt2,ppg077_syserdn,ppg077_syserup);
418 
419 
420 //------------------------------------------------------------------------------------------
421 // Invariant yield of HF electrons in Au+Au from ppg066, units are (GeV/c)^-2
422 // Multiply it by pT, multiply by 2pi.
423 // This should be dN/dpT per event per unit of rapidity.
424 // Then multiply by 2 since we are looking at +-1 unit of rapidity.
425 // First two errors are statistical (up/down) third and fourth are systematic (up/down)
426 // 0-10% central
427 // ppg066 data are in form: 1/2pipT * dN/dpT/dy units are (GeV/c)^{-2}
428 // these are heavy flavor electrons (c+b)
429 //------------------------------------------------------------------------------------------
430 double xau1[99],exau1[99],yau1[99],eyau1up[99],eyau1dn[99],eyau1upsys[99],eyau1dnsys[99];
431 double ycau1[99],eycau1up[99],eycau1dn[99],eycau1upsys[99],eycau1dnsys[99];
432 double ybau1[99],eybau1up[99],eybau1dn[99],eybau1upsys[99],eybau1dnsys[99];
433 double xau2[99],exau2[99],yau2[99],eyau2up[99],eyau2dn[99],eyau2upsys[99],eyau2dnsys[99];
434 double ycau2[99],eycau2up[99],eycau2dn[99],eycau2upsys[99],eycau2dnsys[99];
435 double ybau2[99],eybau2up[99],eybau2dn[99],eybau2upsys[99],eybau2dnsys[99];
436 
437 ifstream finau1;
438 finau1.open("ppg066_0-10.txt");
439 if(!finau1) {cerr << " ERROR: Cannot open input text file. " << endl;}
440 for(int i=0; i<28; i++) {
441 finau1 >> xau1[i] >> yau1[i] >> eyau1up[i] >> eyau1dn[i] >> eyau1upsys[i] >> eyau1dnsys[i];
442 yau1[i] = yau1[i]* xau1[i]*2.*3.141592654*2. * NeventsAuAu; // multiply by pT, by 2pi and by 2 units of rapidity to get dN/dpT per event
443 eyau1up[i] = eyau1up[i]* xau1[i]*2.*3.141592654*2. * NeventsAuAu;
444 eyau1dn[i] = eyau1dn[i]* xau1[i]*2.*3.141592654*2. * NeventsAuAu;
445 eyau1upsys[i] = eyau1upsys[i]* xau1[i]*2.*3.141592654*2. * NeventsAuAu;
446 eyau1dnsys[i] = eyau1dnsys[i]* xau1[i]*2.*3.141592654*2. * NeventsAuAu;
447 
448 double mybottom = gBottomFraction->Eval(xau1[i]);
449 double mycharm = 1.-mybottom;
450 ycau1[i] = yau1[i]* mycharm; // multiply by c/(c+b) ratio to get charm yield
451 eycau1up[i] = eyau1up[i]* mycharm;
452 eycau1dn[i] = eyau1dn[i]* mycharm;
453 eycau1upsys[i] = eyau1upsys[i]* mycharm;
454 eycau1dnsys[i] = eyau1dnsys[i]* mycharm;
455 ybau1[i] = yau1[i]* mybottom; // multiply by b/(c+b) ratio to get bottom yield
456 eybau1up[i] = eyau1up[i]* mybottom;
457 eybau1dn[i] = eyau1dn[i]* mybottom;
458 eybau1upsys[i] = eyau1upsys[i]* mybottom;
459 eybau1dnsys[i] = eyau1dnsys[i]* mybottom;
460 
461 }
462 finau1.close();
463 
464  TGraphAsymmErrors *gHFAuAu = new TGraphAsymmErrors(28,xau1,yau1,0,0,eyau1dn,eyau1up);
465  TGraphAsymmErrors *gCharmAuAu = new TGraphAsymmErrors(28,xau1,ycau1,0,0,eycau1dn,eycau1up);
466  TGraphAsymmErrors *gBottomAuAu = new TGraphAsymmErrors(28,xau1,ybau1,0,0,eybau1dn,eybau1up);
467  TGraphAsymmErrors *gHFAuAuSysErr = new TGraphAsymmErrors(28,xau1,yau1, erpt,erpt,eyau1dnsys,eyau1upsys);
468  TGraphAsymmErrors *gCharmAuAuSysErr = new TGraphAsymmErrors(28,xau1,ycau1,erpt,erpt,eycau1dnsys,eycau1upsys);
469  TGraphAsymmErrors *gBottomAuAuSysErr = new TGraphAsymmErrors(28,xau1,ybau1,erpt,erpt,eybau1dnsys,eybau1upsys);
470 
471 
472 //-------------------------------------------------------------------------------------------
473 // Invariant cross-section of HF electrons from ppg066 (run5 data) in p+p
474 // Convert to invariant yield by dividing by 42mb (p+p total cross-section)
475 //-------------------------------------------------------------------------------------------
476  TGraphErrors *gsys_all4;
477  TGraphErrors *gxs_all4;
478  TFile *fpp = new TFile("/phenix/WWW/p/info/ppg/066/material/RUN5pp_best_fit.root");
479  gxs_all4 = (TGraphErrors*)fpp->Get("gxs_all4");
480  gsys_all4 = (TGraphErrors*)fpp->Get("gsys_all4");
481  fpp->Close();
482  int Ndatapp=gsys_all4->GetN();
483  double *X1 = gxs_all4->GetX();
484  double *Y1 = gxs_all4->GetY();
485  double XX[99], YY[99], Estat[99], Esys[99];
486  double YYC[99], YYB[99], ECstat[99], EBstat[99], ECsys[99], EBsys[99];
487  for(int i=0; i<Ndatapp; i++) {
488  XX[i] = X1[i];
489  YY[i] = Y1[i]/42.e-03 * XX[i]*2.*3.141592654*2. * NeventsAuAu; // multiply by pT, by 2pi and by 2 units of rapidity to get dN/dpT per event
490  Estat[i] = gxs_all4->GetErrorY(i)/42.e-03 * XX[i]*2.*3.141592654*2. * NeventsAuAu;
491  Esys[i] = gsys_all4->GetErrorY(i)/42.e-03 * XX[i]*2.*3.141592654*2. * NeventsAuAu;
492  double mybottom = pythia_botfrac->Eval(XX[i]); // PYTHIA
493  double mycharm = 1.-mybottom;
494  //double mybottom = gfonllup->Eval(XX[i]); // use FONLL upper limit prediction to separate charm and bottom
495  //double mycharm = 1.-mybottom;
496  //double mybottom = gfonll->Eval(XX[i]); // use FONLL prediction to separate charm and bottom
497  //double mycharm = 1.-mybottom;
498  YYC[i] = YY[i] * mycharm;
499  ECstat[i] = Estat[i] * mycharm;
500  ECsys[i] = Esys[i] * mycharm;
501  YYB[i] = YY[i] * mybottom;
502  EBstat[i] = Estat[i] * mybottom;
503  EBsys[i] = Esys[i] * mybottom;
504  }
505  TGraphErrors *ppgraph_stat = new TGraphErrors(Ndatapp,XX,YY,0,Estat);
506  TGraphErrors *ppgraph_stat_charm = new TGraphErrors(Ndatapp,XX,YYC,0,ECstat);
507  TGraphErrors *ppgraph_stat_bottom = new TGraphErrors(Ndatapp,XX,YYB,0,EBstat);
508  TGraphErrors *ppgraph_sys = new TGraphErrors(Ndatapp,XX,YY,0,Esys);
509 
510 
511 
512 
513 
514 
515 
516 //------------------------------------------------------------------------------------------
517 // CCbar from PYTHIA
518 //------------------------------------------------------------------------------------------
519 
520 TChain* chain_ccbar = new TChain("ntp2");
521 
522 // each of these files has 1.63386e+09 generated event
523 // cross-section = 0.187161 mb
524 // 10B 0-10% central Au+Au events before suppression corerespond to 9550B p+p events.
525 // each file corresponds to 1.63386e+09 * 40/0.187161 = 349.19e+09 Min.Bias p+p events
526 // need 9550/349.19 = 27.35 files.
527 // we can use those 27 files if we believe PYTHIA cross-section
528 // PYTHIA ccbar cross-section scale up by 3.31994 +- 0.0403619 (fit above 2GeV); (was 2.09),
529 // so we need 9080 runs (was 5716 runs)
530 // all files below are for 5716 runs
531 // 27 files (2700 runs)
532 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_10.root"); // ckin3_3 is really not ckin3
533 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_11.root");
534 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_12.root");
535 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_13.root");
536 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_14.root");
537 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_15.root");
538 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_16.root");
539 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_17.root");
540 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_18.root");
541 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_19.root");
542 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_20.root");
543 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_21.root");
544 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_22.root");
545 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_23.root");
546 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_24.root");
547 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_25.root");
548 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_26.root");
549 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_27.root");
550 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_28.root");
551 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_29.root");
552 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_30.root");
553 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_31.root");
554 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_32.root");
555 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_34.root");
556 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_35.root");
557 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_36.root");
558 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_37.root");
559 
560 // 17 more files
561 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_4.root");
562 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_5.root");
563 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_6.root");
564 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_7.root");
565 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_8.root");
566 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_38.root");
567 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_41.root");
568 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_42.root");
569 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_43.root");
570 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_44.root");
571 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_48.root");
572 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_49.root");
573 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_51.root");
574 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_53.root");
575 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_54.root");
576 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_55.root");
577 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_56.root");
578 
579 // six files below have 560 runs together (5.6 files)
580 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_9.root");
581 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_39.root");
582 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_46.root");
583 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_47.root");
584 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_50.root");
585 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_52.root");
586 
587 // two files below have together 388 runs (3.88 files)
588 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_1.root");
589 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_2.root");
590 
591 // ten files below are equivalent to 6 files above (six files = 3.6 files above)
592 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar/ccbar_0.root");
593 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar/ccbar_1.root");
594 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar/ccbar_2.root");
595 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar/ccbar_3.root");
596 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar/ccbar_4.root");
597 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar/ccbar_5.root");
598 /*
599 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar/ccbar_6.root");
600 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar/ccbar_7.root");
601 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar/ccbar_8.root");
602 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar/ccbar_9.root");
603 */
604 
605 // Single electron normalization.
606 // PYTHIA cross-section = 0.18713mb
607 // measured cs = 544 +/- 39(stat) +/- 142(syst) +/- 200(model) mkb; from ppg085, arXiv:0802.0050, PLB 670 iss. 4-5 p.313 (2009)
608 // Files 0 through 9 have 1.67471e+07 generated p+p events each
609 // Files 10 through 14 have 6.6809e+07 generated events each
610 // All files correspond to 40/0.18713 * (10*1.67471e+07 + 5*6.6809e+07) = 107.20B Min.Bias p+p events
611 // 10B 0-10% central Au+Au events (Ncoll = 955) correspond to 9550B p+p events (before suppression)
612 // Normalization factor = 9550/107.20 = 89.1
613 double charm_scale = 89.1;
614 TChain* chain_ccbar_sngl = new TChain("ntp1");
615 chain_ccbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_norm/ccbar_norm_0.root");
616 chain_ccbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_norm/ccbar_norm_1.root");
617 chain_ccbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_norm/ccbar_norm_2.root");
618 chain_ccbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_norm/ccbar_norm_3.root");
619 chain_ccbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_norm/ccbar_norm_4.root");
620 chain_ccbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_norm/ccbar_norm_5.root");
621 chain_ccbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_norm/ccbar_norm_6.root");
622 chain_ccbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_norm/ccbar_norm_7.root");
623 chain_ccbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_norm/ccbar_norm_8.root");
624 chain_ccbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_norm/ccbar_norm_9.root");
625 chain_ccbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_norm/ccbar_norm_10.root");
626 chain_ccbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_norm/ccbar_norm_11.root");
627 chain_ccbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_norm/ccbar_norm_12.root");
628 chain_ccbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_norm/ccbar_norm_13.root");
629 chain_ccbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_norm/ccbar_norm_14.root");
630 
631 cout << "Number of entries in CCbar chain: " << chain_ccbar->GetEntries() << endl;
632 cout << "Number of entries in CCbar SINGLES chain: " << chain_ccbar_sngl->GetEntries() << endl;
633 
634 chain_ccbar->SetBranchAddress("mass",&mass);
635 chain_ccbar->SetBranchAddress("pt", &pt);
636 chain_ccbar->SetBranchAddress("eta", &eta);
637 chain_ccbar->SetBranchAddress("pt1", &pt1);
638 chain_ccbar->SetBranchAddress("pt2", &pt2);
639 chain_ccbar->SetBranchAddress("eta1",&eta1);
640 chain_ccbar->SetBranchAddress("eta2",&eta2);
641 
642 chain_ccbar_sngl->SetBranchAddress("pt", &pt_sngl);
643 chain_ccbar_sngl->SetBranchAddress("eta", &eta_sngl);
644 
645 // Singles for normalization
646 for(int j=0; j<chain_ccbar_sngl->GetEntries(); j++){
647  chain_ccbar_sngl->GetEvent(j);
648  if(j%100000==0) cout << "entry # " << j << " " << endl;
649  if(fabs(eta_sngl)<1.0) { hcharm_pt_nosupp->Fill(pt_sngl); }
650  double random = rnd->Uniform(0.,1.);
651  double myraa = gCharmSuppression->Eval(pt_sngl);
652  if(random<=myraa) { if(fabs(eta_sngl)<1.0) { hcharm_pt->Fill(pt_sngl); } }
653 }
654 cout << "Number of entries in charm pT histogram = " << hcharm_pt->GetEntries() << endl;
655 
656 hcharm_pt->Scale(1.0/bsizept); // convert to dN/dpT
657 hcharm_pt->Scale(charm_scale); // scale up to 10B most central AuAu according to cross-section
658  // Normalize PYTHIA to measured charm
659  fnorm_charm->SetParameter(0,1.);
660  gCharmAuAu->Fit(fnorm_charm,"QRN0","",fitstart,fitstop);
661  double norm_charm = fnorm_charm->GetParameter(0);
662  cout << "charm normalization = " << norm_charm << endl;
663 
664 hcharm_pt_nosupp->Scale(1.0/bsizept); // convert to dN/dpT
665 hcharm_pt_nosupp->Scale(charm_scale); // scale up to 10B most central AuAu according to cross-section
666  fnorm_charm_nosupp->SetParameter(0,1.);
667  ppgraph_stat_charm->Fit(fnorm_charm_nosupp,"QRN0","",fitstart,fitstop);
668  double norm_charm_nosupp = fnorm_charm_nosupp->GetParameter(0);
669  cout << "charm normalization (no suppression) = " << norm_charm_nosupp << endl;
670 
671 
672 // invariant mass distributions for CCbar
673 for(int j=0; j<chain_ccbar->GetEntries()*eideff*eideff; j++){
674  chain_ccbar->GetEvent(j);
675  if(j%20000==0) cout << "entry # " << j << endl;
676  if(pt1>ptcut && pt2>ptcut && fabs(eta1)<1.0 && fabs(eta2)<1.0) {
677  double random1 = rnd->Uniform(0.,1.);
678  double myraa1 = gCharmSuppression->Eval(pt1);
679  double random2 = rnd->Uniform(0.,1.);
680  double myraa2 = gCharmSuppression->Eval(pt2);
681  if(random1<=myraa1 && random2<myraa2) {
682  int mybin = -1;
683  for(int i=0; i<nbins; i++) { if(pt<binlim[i+1]) { mybin = i; break; } }
684  if(mybin>=0 && mybin<nbins) { (hcharm[mybin])->Fill(mass); }
685  (hcharm[nbins])->Fill(mass); // all pT
686  }
687  }
688 }
689 
690 //------------------------------------------------------------------------------------------
691 // CCbar with ckin(3) = 4 from PYTHIA
692 //------------------------------------------------------------------------------------------
693 
694 TChain* chain_ccbar_ckin3_4 = new TChain("ntp2");
695 // Each "run" has 4.98345e+06 generated events
696 // Cross-section 0.825489e-03 mb
697 // ccbar_ckin3_4/ccbar cross-section normalization factor = 241478.0e+09 / 9428.13e+09 = 25.61
698 // Fit gives scaling factor = 21.6616 +/- 0.539134
699 // // corerection for acceptance = 1.1824 (we need 1.1824 times more ckin3=4 data because of lower acceptance)
700 // According to PYTHIA cross-section we need 9550/24147.8 * 1.1824 * 100 = 46.76 runs
701 // to get an equivalint of 10B 0-10% central Au+Au events.
702 // Additional cross-section scaling (from measured HF cross-section) = 3.31994 +- 0.0403619 (was 2.09)
703 // This, we need 155.2 runs (was 97.7 runs).
704 //
705 // The file below contains 49 runs, it corresponds to 10B 0-10% AuAu events (before suppression) if we use PYTHIA cross-section.
706 //chain_ccbar_ckin3_4->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_4/ccbar_ckin3_4_9550B_ppevents.root");
707 // The file below has 100 runs
708 chain_ccbar_ckin3_4->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_4/ccbar_ckin3_4_0.root");
709 // the file below has 55 runs
710 chain_ccbar_ckin3_4->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_4/ccbar_ckin3_4_88.root");
711 
712 
713 cout << "Number of entries in CCbar with ckin(3) = 4 chain: " << chain_ccbar_ckin3_4->GetEntries() << endl;
714 
715 chain_ccbar_ckin3_4->SetBranchAddress("mass",&mass);
716 chain_ccbar_ckin3_4->SetBranchAddress("pt", &pt);
717 chain_ccbar_ckin3_4->SetBranchAddress("eta", &eta);
718 chain_ccbar_ckin3_4->SetBranchAddress("pt1", &pt1);
719 chain_ccbar_ckin3_4->SetBranchAddress("pt2", &pt2);
720 chain_ccbar_ckin3_4->SetBranchAddress("eta1",&eta1);
721 chain_ccbar_ckin3_4->SetBranchAddress("eta2",&eta2);
722 
723 // invariant mass distributions for CCbar with ckin(3) = 4
724 for(int j=0; j<chain_ccbar_ckin3_4->GetEntries()*eideff*eideff; j++){
725  chain_ccbar_ckin3_4->GetEvent(j);
726  if(j%20000==0) cout << "entry # " << j << endl;
727  if(pt1>ptcut && pt2>ptcut && fabs(eta1)<1.0 && fabs(eta2)<1.0) {
728  double random1 = rnd->Uniform(0.,1.);
729  double myraa1 = gCharmSuppression->Eval(pt1);
730  double random2 = rnd->Uniform(0.,1.);
731  double myraa2 = gCharmSuppression->Eval(pt2);
732  if(random1<=myraa1 && random2<myraa2) {
733  int mybin = -1;
734  for(int i=0; i<nbins; i++) { if(pt<binlim[i+1]) { mybin = i; break; } }
735  if(mybin>=0 && mybin<nbins) { (hcharm_ckin3_4[mybin])->Fill(mass); }
736  (hcharm_ckin3_4[nbins])->Fill(mass); // all pT
737  }
738  }
739 }
740 
741 //------------------------------------------------------------------------------------------
742 // BBbar from PYTHIA
743 //------------------------------------------------------------------------------------------
744 
745 TChain* chain_bbbar = new TChain("ntp2");
746 // PYTHIA cross-section = 0.000734778 mb
747 // One "run" has 1.189e+06 generated events = 64.734e+09 MinBias p+p events.
748 // So we need 147.5 runs if we believe pythia cross-section
749 // But bbbar cross-section in PYTHIA seems to be 3.22271 +- 0.0388839 (was 2.66) times smaller than one measured by PHENIX
750 // So we need 475.5 runs
751 // The file below has 148 runs and can be used if we believe PYTHIA cross-section
752 chain_bbbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_bbbar/bbbar_9550B_ppevents.root");
753 // The file below has 114 runs
754 //chain_bbbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_bbbar/bbbar_7350B_ppevents.root");
755 //
756 // Files 5, 6, 7 contain 100 runs, file 8 contains 93 runs, file 88 contains 28 runs.
757 // One run has 1.19143e+06 generated events
758 // all files have 476 runs
759 chain_bbbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_bbbar/bbbar_5.root");
760 chain_bbbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_bbbar/bbbar_6.root");
761 chain_bbbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_bbbar/bbbar_7.root");
762 chain_bbbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_bbbar/bbbar_88.root");
763 
764 // Single bottom normalization.
765 // PYTHIA cross-section = 0.00073482mb
766 // FONLL cs = 1.87 +0.99 -0.67 mkb
767 // measured cs = 3.2 +1.2-1.1(stat) +1.4-1.3(sys) [A.Adare et al., Phys.Rev.Lett., 103.082002; arXiv:0903.4851]
768 // measured dSigma/dY at midrapidity = 0.92 +0.34-0.31(stat) +0.39-0.36(sys) mkb [same]
769 // measured from dielectroins cs = 3.9 +2.5(stat) +3-2(sys) mkb [A. Adare et al., Phys. Lett. B670, 313 (2009)]
770 // File 0 has 6.17663e+06 generated p+p events
771 // Corresponds to 40/0.00073482*6.17663e+06 = 336.23B minbias p+p events
772 // 10B 0-10% central Au+Au events (Ncoll = 955) correspond to 9550B p+p events (before suppression)
773 // Bottom normalization factor = 9550/336.23 = 28.4
774 double bottom_scale = 28.4;
775 TChain* chain_bbbar_sngl = new TChain("ntp1");
776 chain_bbbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_bbbar_norm/bbbar_norm_0.root");
777 
778 cout << "Number of entries in BBbar chain: " << chain_bbbar->GetEntries() << endl;
779 cout << "Number of entries in BBbar SINGLES chain: " << chain_bbbar_sngl->GetEntries() << endl;
780 
781 chain_bbbar->SetBranchAddress("mass",&mass);
782 chain_bbbar->SetBranchAddress("pt", &pt);
783 chain_bbbar->SetBranchAddress("eta", &eta);
784 chain_bbbar->SetBranchAddress("pt1", &pt1);
785 chain_bbbar->SetBranchAddress("pt2", &pt2);
786 chain_bbbar->SetBranchAddress("eta1",&eta1);
787 chain_bbbar->SetBranchAddress("eta2",&eta2);
788 
789 chain_bbbar_sngl->SetBranchAddress("pt", &pt_sngl);
790 chain_bbbar_sngl->SetBranchAddress("eta", &eta_sngl);
791 
792 // Singles for normalization
793 for(int j=0; j<chain_bbbar_sngl->GetEntries(); j++){
794  chain_bbbar_sngl->GetEvent(j);
795  if(j%100000==0) cout << "entry # " << j << " " << endl;
796  if(fabs(eta_sngl)<1.0) { hbottom_pt_nosupp->Fill(pt_sngl); }
797  double random = rnd->Uniform(0.,1.);
798  double myraa = gBottomSuppression->Eval(pt_sngl);
799  if(random<=myraa) { if(fabs(eta_sngl)<1.0) { hbottom_pt->Fill(pt_sngl); } }
800 }
801 cout << "Number of entries in bottom pT histogram = " << hbottom_pt->GetEntries() << endl;
802 
803 hbottom_pt->Scale(1.0/bsizept); // convert to dN/dpT
804 hbottom_pt->Scale(bottom_scale); // scale up to 10B most central AuAu according to cross-section
805  // normalize PYTHIA to measured bottom
806  fnorm_bottom->SetParameter(0,1.);
807  gBottomAuAu->Fit(fnorm_bottom,"QRN0","",fitstart,fitstop);
808  double norm_bottom = fnorm_bottom->GetParameter(0);
809  cout << "bottom normalization = " << norm_bottom << endl;
810 
811 hbottom_pt_nosupp->Scale(1.0/bsizept); // convert to dN/dpT
812 hbottom_pt_nosupp->Scale(bottom_scale); // scale up to 10B most central AuAu according to cross-section
813  fnorm_bottom_nosupp->SetParameter(0,1.);
814  ppgraph_stat_bottom->Fit(fnorm_bottom_nosupp,"QRN0","",fitstart,fitstop);
815  double norm_bottom_nosupp = fnorm_bottom_nosupp->GetParameter(0);
816  cout << "bottom normalization (no suppression) = " << norm_bottom_nosupp << endl;
817 
818  // normalize PYTHIA to measured Heavy Flavor
819  hhf_pt = (TH1D*)hbottom_pt->Clone("hhf_pt");
820  hhf_pt->SetName("hhf_pt");
821  hhf_pt->Add(hcharm_pt);
822  fnorm_hf->SetParameter(0,1.);
823  gHFAuAu->Fit(fnorm_hf,"QRN0","",fitstart,fitstop);
824  double norm_hf = fnorm_hf->GetParameter(0);
825  cout << "HF normalization = " << norm_hf << endl;
826 
827  // HF no suppression
828  hhf_pt_nosupp = (TH1D*)hbottom_pt_nosupp->Clone("hhf_pt_nosupp");
829  hhf_pt_nosupp->SetName("hhf_pt_nosupp");
830  hhf_pt_nosupp->Add(hcharm_pt_nosupp);
831  fnorm_hf_nosupp->SetParameter(0,1.);
832  ppgraph_stat->Fit(fnorm_hf_nosupp,"QRN0","",fitstart,fitstop);
833  double norm_hf_nosupp = fnorm_hf_nosupp->GetParameter(0);
834  cout << "HF normalization (no suppression) = " << norm_hf_nosupp << endl;
835 
836 // invariant mass distributions for BBbar
837 for(int j=0; j<chain_bbbar->GetEntries()*eideff*eideff; j++){
838  chain_bbbar->GetEvent(j);
839  if(j%20000==0) cout << "entry # " << j << endl;
840  if(pt1>ptcut && pt2>ptcut && fabs(eta1)<1.0 && fabs(eta2)<1.0) {
841  double random1 = rnd->Uniform(0.,1.);
842  double myraa1 = gBottomSuppression->Eval(pt1);
843  double random2 = rnd->Uniform(0.,1.);
844  double myraa2 = gBottomSuppression->Eval(pt2);
845  if(random1<=myraa1 && random2<=myraa2) {
846  int mybin = -1;
847  for(int i=0; i<nbins; i++) { if(pt<binlim[i+1]) { mybin = i; break; } }
848  if(mybin>=0 && mybin<nbins) { (hbottom[mybin])->Fill(mass); }
849  (hbottom[nbins])->Fill(mass);
850  }
851  }
852 }
853 
854 
855 //-----------------------------------------------------------------------
856 // DY from PYTHIA
857 //----------------------------------------------------------------------
858 TChain* chain_dy = new TChain("ntp2");
859 // PYTHIA cross-section = 0.000118202 mb
860 // so, one generated event corresponds to 40/0.000118202 = 338409.5 MinBias p+p events
861 // one run has 0.465274e+06 generated events
862 // so, 10B 0-10% central AuAu events (9550B p+p min.bias events) correspond to 60.65 runs;
863 // the file below has 61 runs.
864 // do we need to scale up like charm and bottom?
865 chain_dy->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_dy/dy_9550B_ppevents.root");
866 
867 cout << "Number of entries in DY chain: " << chain_dy->GetEntries() << endl;
868 
869 chain_dy->SetBranchAddress("mass",&mass);
870 chain_dy->SetBranchAddress("pt", &pt);
871 chain_dy->SetBranchAddress("eta", &eta);
872 chain_dy->SetBranchAddress("pt1", &pt1);
873 chain_dy->SetBranchAddress("pt2", &pt2);
874 chain_dy->SetBranchAddress("eta1",&eta1);
875 chain_dy->SetBranchAddress("eta2",&eta2);
876 
877 // Invariant mass distributions for DY
878 for(int j=0; j<chain_dy->GetEntries()*eideff*eideff; j++){
879  chain_dy->GetEvent(j);
880  if(j%20000==0) cout << "entry # " << j << endl;
881  if(pt1>ptcut && pt2>ptcut && fabs(eta1)<1.0 && fabs(eta2)<1.0) {
882  int mybin = -1;
883  for(int i=0; i<nbins; i++) { if(pt<binlim[i+1]) { mybin = i; break; } }
884  if(mybin>=0 && mybin<nbins) { (hdy[mybin])->Fill(mass); }
885  (hdy[nbins])->Fill(mass);
886  }
887 }
888 
889 
890 
891 //----------------------------------------------------------------------
892 // Write out
893 //----------------------------------------------------------------------
894 
895 TFile* fout = new TFile(ofname.c_str(),"RECREATE");
896 
897  TH1D* hhcharm[nbins];
898  TH1D* hhcharm_ckin3_4[nbins];
899  TH1D* hhbottom[nbins];
900  TH1D* hhdy[nbins];
901  for(int i=0; i<nbins+1; i++) {
902  sprintf(hname,"hhbottom_%d",i); hhbottom[i] = (TH1D*)(hbottom[i])->Clone(hname);
903  sprintf(hname,"hhcharm_%d",i); hhcharm[i] = (TH1D*)(hcharm[i])->Clone(hname);
904  sprintf(hname,"hhcharm_ckin3_4_%d",i); hhcharm_ckin3_4[i] = (TH1D*)(hcharm_ckin3_4[i])->Clone(hname);
905  sprintf(hname,"hhdy_%d",i); hhdy[i] = (TH1D*)(hdy[i])->Clone(hname);
906  }
907 
908  TH1D* hhbottom_pt = (TH1D*)hbottom_pt->Clone("hhbottom_pt");
909  TH1D* hhcharm_pt = (TH1D*)hcharm_pt->Clone("hhcharm_pt");
910  TH1D* hhhf_pt = (TH1D*)hhbottom_pt->Clone("hhhf_pt");
911  hhhf_pt->Add(hhcharm_pt);
912  hhcharm_pt->SetLineColor(kRed+2);
913  hhbottom_pt->SetLineColor(kBlue+1);
914  hhhf_pt->SetLineColor(kBlack);
915 
916  TH1D* hhbottom_pt_nosupp = (TH1D*)hbottom_pt_nosupp->Clone("hhbottom_pt_nosupp");
917  TH1D* hhcharm_pt_nosupp = (TH1D*)hcharm_pt_nosupp->Clone("hhcharm_pt_nosupp");
918  TH1D* hhhf_pt_nosupp = (TH1D*)hhbottom_pt_nosupp->Clone("hhhf_pt_nosupp");
919  hhhf_pt_nosupp->Add(hhcharm_pt_nosupp);
920  hhcharm_pt_nosupp->SetLineColor(kRed+2);
921  hhbottom_pt_nosupp->SetLineColor(kBlue+1);
922  hhhf_pt_nosupp->SetLineColor(kBlack);
923 
924 
925  //TGraph *ggFracFONLL_b = new TGraph(Ndata,XR,YR); ggFracFONLL_b->SetName("ggFracFONLL_b");
926  TGraph *ggFrac_c = new TGraph(21,ratcbpt,fracc); ggFrac_c->SetName("ggFrac_c");
927  TGraph *ggFrac_b = new TGraph(21,ratcbpt,ratcb); ggFrac_b->SetName("ggFrac_b");
928  TGraphAsymmErrors *ggCharmAuAu = new TGraphAsymmErrors(28,xau1,ycau1,0,0,eycau1dn,eycau1up); ggCharmAuAu->SetName("ggCharmAuAu");
929  TGraphAsymmErrors *ggBottomAuAu = new TGraphAsymmErrors(28,xau1,ybau1,0,0,eybau1dn,eybau1up); ggBottomAuAu->SetName("ggBottomAuAu");
930  TGraphAsymmErrors *ggHFAuAu = new TGraphAsymmErrors(28,xau1,yau1,0,0,eyau1dn,eyau1up); ggHFAuAu->SetName("ggHFAuAu");
931  TGraphAsymmErrors *ggCharmAuAuSysErr = new TGraphAsymmErrors(28,xau1,ycau1,erpt,erpt,eycau1dnsys,eycau1upsys); ggCharmAuAuSysErr->SetName("ggCharmAuAuSysErr");
932  TGraphAsymmErrors *ggBottomAuAuSysErr = new TGraphAsymmErrors(28,xau1,ybau1,erpt,erpt,eybau1dnsys,eybau1upsys); ggBottomAuAuSysErr->SetName("ggBottomAuAuSysErr");
933  TGraphAsymmErrors *ggHFAuAuSysErr = new TGraphAsymmErrors(28,xau1,yau1, erpt,erpt,eyau1dnsys,eyau1upsys); ggHFAuAuSysErr->SetName("ggHFAuAuSysErr");
934  TGraphAsymmErrors* ggrRAA_ppg077 = new TGraphAsymmErrors(28,ppg077_pt,ppg077_raa,0,0,ppg077_staterdn,ppg077_staterup); ggrRAA_ppg077->SetName("ggrRAA_ppg077");
935  TGraphAsymmErrors* ggrRAA_ppg077SysErr = new TGraphAsymmErrors(28,ppg077_pt,ppg077_raa,erpt2,erpt2,ppg077_syserdn,ppg077_syserup); ggrRAA_ppg077SysErr->SetName("ggrRAA_ppg077SysErr");
936 
937  TGraph *ggRAA_charm = new TGraph(9,pt182,raa182charm); ggRAA_charm->SetName("ggRAA_charm");
938  TGraph *ggRAA_bottom = new TGraph(9,pt182,raa182bottom); ggRAA_bottom->SetName("ggRAA_bottom");
939  TGraphAsymmErrors *ggRAASysErr_charm = new TGraphAsymmErrors(6,pt182orig,raa182charmorig, 0,0,raa182charm_erdn, raa182charm_erup); ggRAASysErr_charm->SetName("ggRAASysErr_charm");
940  TGraphAsymmErrors *ggRAASysErr_bottom = new TGraphAsymmErrors(6,pt182orig,raa182bottomorig,0,0,raa182bottom_erdn,raa182bottom_erup); ggRAASysErr_bottom->SetName("ggRAASysErr_bottom");
941 
942  ggFrac_b->SetLineColor(kMagenta);
943  ggFrac_b->SetLineWidth(2);
944 // ggFracFONLL_b->SetLineColor(kBlack);
945 // ggFracFONLL_b->SetLineWidth(2);
946 
947  ggHFAuAu->SetMarkerStyle(20);
948  ggHFAuAu->SetMarkerColor(kBlack);
949  ggBottomAuAu->SetMarkerStyle(20);
950  ggBottomAuAu->SetMarkerColor(kBlue);
951  ggBottomAuAu->SetLineColor(kBlue);
952  ggCharmAuAu->SetMarkerStyle(20);
953  ggCharmAuAu->SetMarkerColor(kRed);
954  ggCharmAuAu->SetLineColor(kRed);
955 
956  ggHFAuAuSysErr->SetFillColor(0);
957  ggHFAuAuSysErr->SetMarkerStyle(6);
958  ggHFAuAuSysErr->SetMarkerSize(0.5);
959  ggHFAuAuSysErr->SetMarkerColor(kBlack);
960  ggHFAuAuSysErr->SetLineColor(kBlack);
961  ggHFAuAuSysErr->SetLineWidth(2);
962  ggHFAuAuSysErr->SetFillColor(kGray+1);
963 
964  ggrRAA_ppg077SysErr->SetFillColor(0);
965  ggrRAA_ppg077SysErr->SetMarkerStyle(6);
966  ggrRAA_ppg077SysErr->SetMarkerSize(0.5);
967  ggrRAA_ppg077SysErr->SetMarkerColor(kBlack);
968  ggrRAA_ppg077SysErr->SetLineColor(kBlack);
969  ggrRAA_ppg077SysErr->SetLineWidth(2);
970  ggrRAA_ppg077SysErr->SetFillColor(kGray+1);
971 
972  ggrRAA_ppg077->SetLineWidth(2);
973  ggrRAA_ppg077->SetMarkerStyle(20);
974  ggrRAA_ppg077->SetMarkerSize(1.0);
975  ggrRAA_ppg077->SetMarkerColor(kBlack);
976 
977  ggRAASysErr_charm->SetFillColor(kGreen+2);
978  ggRAASysErr_charm->SetFillStyle(3001);
979  ggRAA_charm->SetLineWidth(3);
980  ggRAA_charm->SetLineColor(kGreen+2);
981  ggRAASysErr_bottom->SetFillColor(kBlue);
982  ggRAASysErr_bottom->SetFillStyle(3002);
983  ggRAA_bottom->SetLineWidth(3);
984  ggRAA_bottom->SetLineColor(kBlue);
985 
986 // ggFracFONLL_b->Write();
987  ggFrac_c->Write();
988  ggFrac_b->Write();
989  ggRAA_charm->Write();
990  ggRAA_bottom->Write();
991  ggRAASysErr_charm->Write();
992  ggRAASysErr_bottom->Write();
993  ggCharmAuAu->Write();
994  ggBottomAuAu->Write();
995  ggHFAuAu->Write();
996  ggCharmAuAuSysErr->Write();
997  ggBottomAuAuSysErr->Write();
998  ggHFAuAuSysErr->Write();
999  ggrRAA_ppg077SysErr->Write();
1000  ggrRAA_ppg077->Write();
1001 
1002  TF1 ffnorm_bottom = TF1(*fnorm_bottom);
1003  TF1 ffnorm_charm = TF1(*fnorm_charm);
1004  TF1 ffnorm_hf = TF1(*fnorm_hf);
1005  ffnorm_bottom.SetName("ffnorm_bottom");
1006  ffnorm_charm.SetName("ffnorm_charm");
1007  ffnorm_hf.SetName("ffnorm_hf");
1008  ffnorm_bottom.Write();
1009  ffnorm_charm.Write();
1010  ffnorm_hf.Write();
1011  TF1 ffnorm_bottom_nosupp = TF1(*fnorm_bottom_nosupp);
1012  TF1 ffnorm_charm_nosupp = TF1(*fnorm_charm_nosupp);
1013  TF1 ffnorm_hf_nosupp = TF1(*fnorm_hf_nosupp);
1014  ffnorm_bottom_nosupp.SetName("ffnorm_bottom_nosupp");
1015  ffnorm_charm_nosupp.SetName("ffnorm_charm_nosupp");
1016  ffnorm_hf_nosupp.SetName("ffnorm_hf_nosupp");
1017  ffnorm_bottom_nosupp.Write();
1018  ffnorm_charm_nosupp.Write();
1019  ffnorm_hf_nosupp.Write();
1020 
1021 
1022  TF1 mymybotfrac = TF1(*mybotfrac);
1023  mymybotfrac.SetName("mymybotfrac");
1024  mymybotfrac.SetLineColor(kGreen+2);
1025  mymybotfrac.Write();
1026 
1027  TF1 ppythia_botfrac = TF1(*pythia_botfrac);
1028  ppythia_botfrac.SetName("ppythia_botfrac");
1029  ppythia_botfrac.SetLineColor(kGreen+2);
1030  ppythia_botfrac.Write();
1031 
1032 
1033  TGraphAsymmErrors* ttcharm = new TGraphAsymmErrors(*tcharm);
1034  TGraphAsymmErrors* ttbottom = new TGraphAsymmErrors(*tbottom);
1035  ttcharm->SetFillColor(kGreen+2);
1036  ttcharm->SetFillStyle(3001);
1037  ttcharm->SetLineWidth(1);
1038  ttcharm->SetName("ttcharm");
1039  ttcharm->Write();
1040  ttcharm->SetLineColor(kGreen+2);
1041  ttbottom->SetFillColor(kBlue);
1042  ttbottom->SetFillStyle(3002);
1043  ttbottom->SetLineWidth(1);
1044  ttbottom->SetLineColor(kBlue);
1045  ttbottom->SetName("ttbottom");
1046  ttbottom->Write();
1047 
1048  TGraphAsymmErrors *ggfonll = new TGraphAsymmErrors(*gfonll);
1049  ggfonll->SetLineColor(kBlack);
1050  ggfonll->SetFillColor(kGray);
1051  ggfonll->SetLineWidth(3.0);
1052  ggfonll->SetFillStyle(3001);
1053  ggfonll->SetName("ggfonll");
1054  ggfonll->Write();
1055  TGraph *ggfonllup = new TGraph(*gfonllup);
1056  ggfonllup->SetLineColor(kBlue+2);
1057  ggfonllup->SetLineWidth(2.0);
1058  ggfonllup->SetName("ggfonllup");
1059  ggfonllup->Write();
1060 
1061  TGraphAsymmErrors* bbfrac_cent1 = new TGraphAsymmErrors(*bfrac_cent1);
1062  bbfrac_cent1->SetLineColor(kOrange+8);
1063  bbfrac_cent1->SetFillColor(kOrange+8);
1064  bbfrac_cent1->SetLineWidth(3.0);
1065  bbfrac_cent1->SetFillStyle(3002);
1066  bbfrac_cent1->SetName("bbfrac_cent1");
1067  bbfrac_cent1->Write();
1068 
1069  TGraphErrors *ttRaadata = new TGraphErrors(*tRaadata);
1070  TGraphErrors *ttRaadata_sys = new TGraphErrors(*tRaadata_sys);
1071  ttRaadata_sys->SetMarkerStyle(21);
1072  ttRaadata_sys->SetMarkerSize(1.0);
1073  ttRaadata_sys->SetMarkerColor(kGray+2);
1074  ttRaadata_sys->SetLineColor(kGray+2);
1075  ttRaadata_sys->SetFillStyle(0);
1076  ttRaadata_sys->SetName("ttRaadata_sys");
1077  ttRaadata_sys->Write();
1078  ttRaadata_sys->SetFillStyle(0);
1079  ttRaadata->SetMarkerStyle(21);
1080  ttRaadata->SetMarkerSize(1.0);
1081  ttRaadata->SetMarkerColor(kGray+2);
1082  ttRaadata->SetName("ttRaadata");
1083  ttRaadata->Write();
1084 
1085  TGraph* mymytcharm = new TGraph(*mytcharm);
1086  mymytcharm->SetLineWidth(3);
1087  mymytcharm->SetLineColor(kGreen+2);
1088  mymytcharm->SetName("mymytcharm");
1089  mymytcharm->Write();
1090  TGraph* mymytbottom = new TGraph(*mytbottom);
1091  mymytbottom->SetLineWidth(3);
1092  mymytbottom->SetLineColor(kBlue);
1093  mymytbottom->SetName("mymytbottom");
1094  mymytbottom->Write();
1095 
1096  TGraph* grgr_mysupp_charm = new TGraph(*gr_mysupp_charm);
1097  grgr_mysupp_charm->SetLineWidth(3);
1098  grgr_mysupp_charm->SetLineColor(kGreen+2);
1099  grgr_mysupp_charm->SetName("grgr_mysupp_charm");
1100  grgr_mysupp_charm->Write();
1101  TGraph* grgr_mysupp_bottom = new TGraph(*gr_mysupp_bottom);
1102  grgr_mysupp_bottom->SetLineWidth(3);
1103  grgr_mysupp_bottom->SetLineColor(kGreen+2);
1104  grgr_mysupp_bottom->SetName("grgr_mysupp_bottom");
1105  grgr_mysupp_bottom->Write();
1106 
1107  TGraphErrors *pppgraph_stat = new TGraphErrors(Ndatapp,XX,YY,0,Estat);
1108  TGraphErrors *pppgraph_stat_charm = new TGraphErrors(Ndatapp,XX,YYC,0,ECstat);
1109  TGraphErrors *pppgraph_stat_bottom = new TGraphErrors(Ndatapp,XX,YYB,0,EBstat);
1110  TGraphErrors *pppgraph_sys = new TGraphErrors(Ndatapp,XX,YY,0,Esys);
1111  pppgraph_stat->SetMarkerStyle(24);
1112  pppgraph_stat_charm->SetMarkerStyle(24);
1113  pppgraph_stat_charm->SetMarkerColor(kRed);
1114  pppgraph_stat_charm->SetLineColor(kRed);
1115  pppgraph_stat_bottom->SetMarkerStyle(24);
1116  pppgraph_stat_bottom->SetMarkerColor(kBlue);
1117  pppgraph_stat_bottom->SetLineColor(kBlue);
1118  pppgraph_stat->SetName("pppgraph_stat");
1119  pppgraph_stat->Write();
1120  pppgraph_stat_charm->SetName("pppgraph_stat_charm");
1121  pppgraph_stat_charm->Write();
1122  pppgraph_stat_bottom->SetName("pppgraph_stat_bottom");
1123  pppgraph_stat_bottom->Write();
1124  pppgraph_sys->SetName("pppgraph_sys");
1125  pppgraph_sys->Write();
1126 
1127 TGraph* ggr_mybotfrac = new TGraph(*gr_mybotfrac);
1128 ggr_mybotfrac->SetLineWidth(2);
1129 ggr_mybotfrac->SetLineColor(kGreen+2);
1130 ggr_mybotfrac->SetName("ggr_mybotfrac");
1131 ggr_mybotfrac->Write();
1132 
1133 TGraph* ggr_pythia_botfrac = new TGraph(*gr_pythia_botfrac);
1134 ggr_pythia_botfrac->SetLineWidth(2);
1135 ggr_pythia_botfrac->SetLineColor(kGreen+2);
1136 ggr_pythia_botfrac->SetName("ggr_pythia_botfrac");
1137 ggr_pythia_botfrac->Write();
1138 
1139 
1140 fout->Write();
1141 fout->Close();
1142 
1143 
1144  return 0;
1145 }
1146