Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:14:51 

0001 
0002 #include "TFile.h"
0003 #include "TF1.h"
0004 #include "TH1.h"
0005 #include "TH2.h"
0006 #include "TGraph.h"
0007 #include "TChain.h"
0008 #include "TRandom3.h"
0009 #include "TGraphAsymmErrors.h"
0010 #include "TGraphErrors.h"
0011 
0012 #include <fstream>
0013 #include <iostream>
0014 //#include <iomanip>
0015 //#include <vector>
0016 
0017 using namespace std;
0018 
0019 int ccbb(double eideff, string ofname);
0020 
0021 TH1D* hcharm_pt;
0022 TH1D* hbottom_pt;
0023 TH1D* hhf_pt;
0024 TH1D* hcharm_pt_nosupp;
0025 TH1D* hbottom_pt_nosupp;
0026 TH1D* hhf_pt_nosupp;
0027 
0028 double norm_charm(double* x, double* par) {
0029   double pt = x[0];
0030   double norm = par[0];
0031   int bin = hcharm_pt->GetXaxis()->FindBin(pt);
0032   double value = hcharm_pt->GetBinContent(bin);
0033   return norm*value;
0034 }
0035 double norm_charm_nosupp(double* x, double* par) {
0036   double pt = x[0];
0037   double norm = par[0];
0038   int bin = hcharm_pt_nosupp->GetXaxis()->FindBin(pt);
0039   double value = hcharm_pt_nosupp->GetBinContent(bin);
0040   return norm*value;
0041 }
0042 
0043 double norm_bottom(double* x, double* par) {
0044   double pt = x[0];
0045   double norm = par[0];
0046   int bin = hbottom_pt->GetXaxis()->FindBin(pt);
0047   double value = hbottom_pt->GetBinContent(bin);
0048   return norm*value;
0049 }
0050 double norm_bottom_nosupp(double* x, double* par) {
0051   double pt = x[0];
0052   double norm = par[0];
0053   int bin = hbottom_pt_nosupp->GetXaxis()->FindBin(pt);
0054   double value = hbottom_pt_nosupp->GetBinContent(bin);
0055   return norm*value;
0056 }
0057 
0058 double norm_hf(double* x, double* par) {
0059   double pt = x[0];
0060   double norm = par[0];
0061   int bin = hhf_pt->GetXaxis()->FindBin(pt);
0062   double value = hhf_pt->GetBinContent(bin);
0063   return norm*value;
0064 }
0065 double norm_hf_nosupp(double* x, double* par) {
0066   double pt = x[0];
0067   double norm = par[0];
0068   int bin = hhf_pt_nosupp->GetXaxis()->FindBin(pt);
0069   double value = hhf_pt_nosupp->GetBinContent(bin);
0070   return norm*value;
0071 }
0072 
0073 
0074 
0075 //====================================================================================
0076 
0077 int main(int argc, char* argv[]) 
0078 {
0079   double eideff = 0.9;
0080   string ofname="ccbb.root";
0081   if(argc==1) cout << argv[0] << " is running with standard parameters..." << endl;
0082   if(argc==2) { eideff = atoi(argv[1]); }
0083   if(argc==3) { eideff = atoi(argv[1]); ofname=argv[2]; }
0084   cout << "eID efficiency = " << eideff << endl;
0085   cout << "output file name: " << ofname << endl;
0086     return ccbb(eideff, ofname);
0087 }
0088 
0089 //------------------------------------------------------------------------------------
0090 
0091 int ccbb(double eideff, string ofname) {
0092 
0093 // 0 = ppg182, min. bias AuAu
0094 // 1 = FONLL, p+p
0095 // 2 = my own function
0096 // 3 = Kazuya, 0-10% central AuAu
0097 int which_botfrac = 3;
0098 
0099 // 0 = Min.Bias from ppg182
0100 // 1 = 0-10% central from Kazuya
0101 // 2 = no suppression
0102 // 3 = my own suppression
0103 int which_supp = 2;
0104 
0105 double ptcut = 2.0;
0106 
0107 double NeventsAuAu = 10.0e+09;
0108 
0109 TRandom* rnd = new TRandom3();
0110 rnd->SetSeed(0);
0111 
0112 const int nbins = 15;
0113 double binlim[nbins+1];
0114 for(int i=0; i<=nbins; i++) {binlim[i]=double(i);}
0115 
0116 double startmass=0.;
0117 double stopmass=20.;
0118 int nchanmass = 400;
0119 double startpt=0.;
0120 double stoppt=10.;
0121 int nchanpt = 100;
0122 double bsizept = (stoppt-startpt)/float(nchanpt);
0123 
0124 char hname[999];
0125 
0126 double erpt[99], erpt2[99]; for(int i=0; i<99; i++) {erpt[i]=0.075; erpt2[i]=0.1;}
0127 std::string tmp0,tmp1,tmp2,tmp3,tmp4,tmp5;
0128 
0129 // last histogram is integrated over pT
0130 TH1D* hcharm[nbins+1];
0131 TH1D* hcharm_ckin3_4[nbins+1];
0132 TH1D* hbottom[nbins+1];
0133 TH1D* hdy[nbins+1];
0134 for(int i=0; i<nbins+1; i++) { sprintf(hname,"hcharm_%d",i);  hcharm[i]  = new TH1D(hname,"",nchanmass,startmass,stopmass); hcharm[i]->Sumw2(); }
0135 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(); }
0136 for(int i=0; i<nbins+1; i++) { sprintf(hname,"hbottom_%d",i); hbottom[i] = new TH1D(hname,"",nchanmass,startmass,stopmass); hbottom[i]->Sumw2(); }
0137 for(int i=0; i<nbins+1; i++) { sprintf(hname,"hdy_%d",i);     hdy[i]     = new TH1D(hname,"",nchanmass,startmass,stopmass); hdy[i]->Sumw2();}
0138 hcharm_pt  = new TH1D("hcharm_pt", "",nchanpt,startpt,stoppt);
0139 hbottom_pt = new TH1D("hbottom_pt","",nchanpt,startpt,stoppt);
0140 hcharm_pt->Sumw2();
0141 hbottom_pt->Sumw2();
0142 hcharm_pt_nosupp  = new TH1D("hcharm_pt_nosupp", "",nchanpt,startpt,stoppt);
0143 hbottom_pt_nosupp = new TH1D("hbottom_pt_nosupp","",nchanpt,startpt,stoppt);
0144 hcharm_pt_nosupp->Sumw2();
0145 hbottom_pt_nosupp->Sumw2();
0146 
0147 float mass,pt,eta,pt1,pt2,eta1,eta2;
0148 float pt_sngl,eta_sngl;
0149 
0150 float fitstart=2.0;
0151 float fitstop=8.0;
0152 TF1* fnorm_charm  = new TF1("fnorm_charm", norm_charm, fitstart,fitstop,1);
0153 TF1* fnorm_bottom = new TF1("fnorm_bottom",norm_bottom,fitstart,fitstop,1);
0154 TF1* fnorm_hf     = new TF1("fnorm_hf"   , norm_hf,    fitstart,fitstop,1);
0155 TF1* fnorm_charm_nosupp  = new TF1("fnorm_charm_nosupp", norm_charm_nosupp, fitstart,fitstop,1);
0156 TF1* fnorm_bottom_nosupp = new TF1("fnorm_bottom_nosupp",norm_bottom_nosupp,fitstart,fitstop,1);
0157 TF1* fnorm_hf_nosupp     = new TF1("fnorm_hf_nosupp"   , norm_hf_nosupp,    fitstart,fitstop,1);
0158 
0159 /*
0160 //--------------------------------------------------------------------------------------
0161 //--------------------------------------------------------------------------------------
0162 TFile *fFONLL=new TFile("FONLL_ratio.root");
0163   TGraph *gYieldFONLL_c = (TGraph*)fFONLL->Get("gRatio_c")->Clone();
0164   TGraph *gYieldFONLL_b = (TGraph*)fFONLL->Get("gRatio_b")->Clone();
0165   int Ndata=gYieldFONLL_c->GetN();
0166   double *Xc=gYieldFONLL_c->GetX();
0167   double *Yc=gYieldFONLL_c->GetY();
0168   double *Xb=gYieldFONLL_b->GetX();
0169   double *Yb=gYieldFONLL_b->GetY();
0170   double XR[999], YR[999], ZR[999];
0171   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]);}
0172   TGraph *gFracFONLL_b  = new TGraph(Ndata,XR,YR);
0173   TGraph *gFracFONLL_c  = new TGraph(Ndata,XR,ZR);
0174 fFONLL->Close();
0175 */
0176 //
0177 // bottom/(charm+bottom) and charm/(charm+bottom) ratios in p+p from FONLL  
0178 // from Darren (ppg182)
0179 //
0180 ifstream fin_fonll;
0181 double pt_fonll[99],fonll[99],fonllup[99],fonll_erup[99],fonll_erdn[99],dtmp0,dtmp1;
0182 fin_fonll.open("FONLL_ratio.txt");
0183 if(!fin_fonll) {cerr << " ERROR: Cannot open input text file. " << endl;}
0184 for(int i=0; i<19; i++) {
0185   fin_fonll >> pt_fonll[i] >> tmp0 >> fonll[i] >> tmp1 >> dtmp0 >> tmp2 >> dtmp1 >> tmp3;
0186   fonllup[i] = dtmp0;
0187   fonll_erup[i] = dtmp0 - fonll[i];
0188   fonll_erdn[i] = fonll[i] - dtmp1;
0189   cout << pt_fonll[i] << " " << fonll[i] << " " << fonll_erup[i] << " " << fonll_erdn[i] << endl;
0190 }
0191 fin_fonll.close();
0192   TGraphAsymmErrors *gfonll  = new TGraphAsymmErrors(19,pt_fonll,fonll,0,0,fonll_erdn,fonll_erup);
0193   TGraph *gfonllup  = new TGraph(19,pt_fonll,fonllup);
0194 
0195 
0196 //----------------------------------------------------------------------
0197 // my bottom/(charm+bottom) fraction just for testing
0198 //----------------------------------------------------------------------
0199 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 
0200 //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
0201 double myx[999];
0202 double myy[999];
0203 double pyy[999];
0204 int myn = 150;
0205 for(int i=0; i<myn; i++) { myx[i] = i*0.1; myy[i] = mybotfrac->Eval(myx[i]); }
0206 //for(int i=0; i<myn; i++) { myx[i] = i*0.1; myy[i] = gfonll->Eval(myx[i])*mybotfrac->Eval(myx[i]); }
0207 TGraph* gr_mybotfrac = new TGraph(myn,myx,myy);
0208 
0209 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
0210 for(int i=0; i<myn; i++) { pyy[i] = pythia_botfrac->Eval(myx[i]); }
0211 TGraph* gr_pythia_botfrac = new TGraph(myn,myx,pyy);
0212 
0213 //--------------------------------------------------------------------------------------
0214 // bottom/(charm+bottom) ratio in Min. Bias Au+Au from ppg182:  
0215 // arXiv:1509.04662 
0216 // A. Adare et al. (PHENIX Collaboration) Phys. Rev. C 93, 034904 (2016)
0217 //--------------------------------------------------------------------------------------
0218 double ratcb[99],ratcberup[99],ratcberdn[99],ratcbpt[99];
0219 for(int i=0; i<15; i++) { ratcbpt[i] = 1.1+0.2*i; }
0220 ratcbpt[15] = 4.25;
0221 ratcbpt[16] = 4.75;
0222 ratcbpt[17] = 5.5;
0223 ratcbpt[18] = 6.5;
0224 ratcbpt[19] = 7.5;
0225 ratcbpt[20] = 8.5;
0226 ratcb[0] = 3.99e-02;
0227 ratcb[1] = 6.15e-02;
0228 ratcb[2] = 9.41e-02;
0229 ratcb[3] = 1.40e-01;
0230 ratcb[4] = 2.01e-01;
0231 ratcb[5] = 2.73e-01;
0232 ratcb[6] = 3.49e-01;
0233 ratcb[7] = 4.22e-01;
0234 ratcb[8] = 4.84e-01;
0235 ratcb[9] = 5.32e-01;
0236 ratcb[10] = 5.64e-01;
0237 ratcb[11] = 5.81e-01;
0238 ratcb[12] = 5.89e-01;
0239 ratcb[13] = 5.87e-01;
0240 ratcb[14] = 5.79e-01;
0241 ratcb[15] = 5.57e-01;
0242 ratcb[16] = 5.13e-01;
0243 ratcb[17] = 4.79e-01;
0244 ratcb[18] = 4.71e-01;
0245 ratcb[19] = 4.61e-01;
0246 ratcb[20] = 4.36e-01;
0247   TGraph *gFrac_b = new TGraph(21,ratcbpt,ratcb);
0248 
0249 double fracc[99];
0250 for(int i=0; i<21; i++) { fracc[i] = 1.0-ratcb[i]; }
0251   TGraph *gFrac_c = new TGraph(21,ratcbpt,fracc);
0252 
0253 
0254 //-------------------------------------------------------------------------------
0255 // bottom fraction from Kazuya Nagashima's talk at QM2017
0256 // PHENIX preliminary for 0-10% central Au+Au 
0257 // http://www.phenix.bnl.gov/phenix/WWW/p/draft/nagasy/QM2017/slide/PHENIX_HF_KazuyaNagashima.pdf
0258 //-------------------------------------------------------------------------------
0259   TGraphAsymmErrors* bfrac_cent1;
0260   TFile *fk = new TFile("kazuya/bfrac_cent1.root");
0261     bfrac_cent1 = (TGraphAsymmErrors*)fk->Get("bfrac_cent1");
0262   fk->Close();
0263 
0264 //-----------------------------------------------------------------------------------------
0265 //  R_AA of HF electrons according to Kazuya Nagashima; 0-10% central; 2004+2014 data
0266 //  Different from ppg077 
0267 //  Probably because it's combined 2004 and 2014 results (ppg077 = 2004 only)
0268 //-----------------------------------------------------------------------------------------
0269   ifstream kfile;
0270   kfile.open("kazuya/csv/ppg077_HFRAA_cent1.csv");
0271   TGraphErrors *tRaadata = new TGraphErrors();
0272   TGraphErrors *tRaadata_sys = new TGraphErrors();
0273   double var1, var2, var3, var4, var5, var6;
0274   for(int i=0;i<28;i++) {
0275     kfile >> var1 >> var2 >> var3 >> var4 >> var5 >> var6;
0276     tRaadata->SetPoint(i,var1,var2);
0277     tRaadata->SetPointError(i,0.,var3);
0278     tRaadata_sys->SetPoint(i,var1,var2);
0279     tRaadata_sys->SetPointError(i,0.1,var5);
0280   }
0281   kfile.close();
0282 
0283 //--------------------------------------------------------------------------------------
0284 // Charm and Bottom R_AA in 0-10% central Au+Au from Kazuya Nagashima's QM2017 presentation  
0285 //--------------------------------------------------------------------------------------
0286   TGraphAsymmErrors* tcharm;
0287   TGraphAsymmErrors* tbottom;
0288   TFile *fk2 = new TFile("kazuya/Raa_cent1.root");
0289     tcharm = (TGraphAsymmErrors*)fk2->Get("tcharm");
0290     tbottom = (TGraphAsymmErrors*)fk2->Get("tbottom");
0291   fk2->Close();
0292 
0293   TGraph* mytcharm = new TGraph();
0294   TGraph* mytbottom = new TGraph();
0295   mytcharm->SetPoint(0,0.0,1.0);
0296   mytcharm->SetPoint(1,1.5,1.0);
0297   mytbottom->SetPoint(0,0.0,1.0);
0298   mytbottom->SetPoint(1,1.5,1.0);
0299   double *tmpx = tcharm->GetX();
0300   double *tmpy = tcharm->GetY();
0301   for(int i=2; i<8; i++) { mytcharm->SetPoint(i,tmpx[i-2],tmpy[i-2]); }
0302   mytcharm->SetPoint(8,15.,tmpy[5]);
0303     tmpx = tbottom->GetX();
0304     tmpy = tbottom->GetY();
0305     for(int i=2; i<8; i++) { mytbottom->SetPoint(i,tmpx[i-2],tmpy[i-2]); }
0306     mytbottom->SetPoint(2,tmpx[0],1.0);
0307     mytbottom->SetPoint(8,15.,tmpy[5]);
0308 
0309 
0310 //--------------------------------------------------------------------------------------
0311 // Charm and Bottom R_AA in Min. Bias Au+Au from ppg182:  
0312 // arXiv:1509.04662 
0313 // A. Adare et al. (PHENIX Collaboration) Phys. Rev. C 93, 034904 (2016)
0314 //--------------------------------------------------------------------------------------
0315 
0316 double pt182[9],raa182charm[9],raa182bottom[9],raa182charm_erdn[9],raa182charm_erup[9],raa182bottom_erdn[9],raa182bottom_erup[9];
0317 double pt182orig[9],raa182charmorig[9],raa182bottomorig[9];
0318 pt182[0]=0.0;
0319 pt182[1]=1.5;
0320 raa182charm[0] = 1.0;
0321 raa182charm[1] = 1.0;
0322 raa182bottom[0] = 1.0;
0323 raa182bottom[1] = 1.0;
0324 
0325 ifstream fin182;
0326 fin182.open("ppg182.txt");
0327 if(!fin182) {cerr << " ERROR: Cannot open input text file. " << endl;}
0328 fin182 >> tmp0 >> tmp1 >> tmp2 >> tmp3;
0329 fin182 >> tmp0 >> tmp1 >> tmp2 >> tmp3;
0330 for(int i=2; i<8; i++) {
0331   fin182 >> pt182[i] >> raa182charm[i] >> raa182charm_erup[i-2] >> raa182charm_erdn[i-2];
0332   pt182orig[i-2]=pt182[i];
0333   raa182charmorig[i-2] = raa182charm[i]; 
0334 }
0335 fin182 >> tmp0 >> tmp1 >> tmp2 >> tmp3;
0336 fin182 >> tmp0 >> tmp1 >> tmp2 >> tmp3;
0337 for(int i=2; i<8; i++) {
0338   fin182 >> pt182[i] >> raa182bottom[i] >> raa182bottom_erup[i-2] >> raa182bottom_erdn[i-2];
0339   raa182bottomorig[i-2] = raa182bottom[i]; 
0340 }
0341 
0342 pt182[8]=15.;
0343 raa182charm[8] = raa182charm[7];
0344 raa182bottom[2] = 1.0;  // measured value is 1.15
0345 raa182bottom[8] = raa182bottom[7];
0346   TGraph *gRAA_charm  = new TGraph(9,pt182,raa182charm);
0347   TGraph *gRAA_bottom = new TGraph(9,pt182,raa182bottom);
0348   TGraphAsymmErrors *gRAASysErr_charm  = new TGraphAsymmErrors(6,pt182orig,raa182charmorig, 0,0,raa182charm_erdn, raa182charm_erup);
0349   TGraphAsymmErrors *gRAASysErr_bottom = new TGraphAsymmErrors(6,pt182orig,raa182bottomorig,0,0,raa182bottom_erdn,raa182bottom_erup);
0350 
0351 
0352 // My own RAA from comparison with data
0353   //TF1* fmyraac =  new TF1("fmyraac","0.303395+18.995424/(66.573738+pow(x,7.237049))",1.0,9.0);  // fit above 3 GeV
0354   //TF1* fmyraab =  new TF1("fmyraab","0.347859+64550.024555/(140970.430941+pow(x,9.780939))",1.0,9.0); // fit above 3 Gev
0355   TF1* fmyraac =  new TF1("fmyraac","0.327075+20.482705/(66.589763+pow(x,7.237551))",1.0,9.0);  // fit above 2 GeV
0356   TF1* fmyraab =  new TF1("fmyraab","0.370667+68699.897905/(140801.861045+pow(x,9.779860))",1.0,9.0); // fit above 2 Gev
0357     double mysuppc[999],mysuppb[999],myxx[999];
0358     int mycount=0;
0359     for(int i=0; i<myn; i++) { if(myx[i]>1.0) { myxx[mycount] = myx[i]; mysuppc[mycount] = fmyraac->Eval(myx[i]); mycount++; } }
0360     TGraph* gr_mysupp_charm = new TGraph(mycount,myxx,mysuppc);
0361     mycount=0;
0362     for(int i=0; i<myn; i++) { if(myx[i]>1.0) { mysuppb[mycount] = fmyraab->Eval(myx[i]); mycount++; } }
0363     TGraph* gr_mysupp_bottom = new TGraph(mycount,myxx,mysuppb);
0364 
0365 
0366 // Choose which suppression to use
0367 TGraph* gCharmSuppression;
0368 TGraph* gBottomSuppression;
0369 if(which_supp==0) {
0370    gCharmSuppression =  new TGraph(*gRAA_charm);
0371    gBottomSuppression = new TGraph(*gRAA_bottom);
0372 } else if(which_supp==1) {
0373    gCharmSuppression =  new TGraph(*mytcharm);
0374    gBottomSuppression = new TGraph(*mytbottom);
0375 } else if(which_supp==2) {
0376    double x1[4] = {0.,5.,10.,15.};
0377    double y1[4] = {1.,1.,1., 1.};
0378    gCharmSuppression =  new TGraph(4,x1,y1);
0379    gBottomSuppression = new TGraph(4,x1,y1);
0380 } else if(which_supp==3) {
0381    gCharmSuppression =  new TGraph(*gr_mysupp_charm);
0382    gBottomSuppression = new TGraph(*gr_mysupp_bottom);
0383 } else {
0384   cerr << "ERROR: Wrong suppression factors !!!" << endl;
0385 }
0386 
0387 // Choose which Bottom Fraction to use
0388 TGraph* gBottomFraction;
0389 if(which_botfrac==0) {        // ppg182, min. bias AuAu
0390  gBottomFraction = new TGraph(*gFrac_b);
0391 } else if(which_botfrac==1) { // FONLL p+p
0392  //gBottomFraction = new TGraph(*gFracFONLL_b);
0393  gBottomFraction = new TGraph(*gfonll);
0394 } else if(which_botfrac==2) { // my own function
0395  gBottomFraction = new TGraph(*gr_mybotfrac);
0396 } else if(which_botfrac==3) { // Kazuya, 0-10% central AuAu
0397  gBottomFraction = new TGraph(*bfrac_cent1);
0398 }
0399 
0400 
0401 
0402 //-----------------------------------------------------------------------------------------
0403 //  R_AA of HF electrons from ppg077; 0-10% central Au+Au; 2004 data
0404 //-----------------------------------------------------------------------------------------
0405 
0406 double ppg077_pt[99],ppg077_raa[99],ppg077_staterup[99],ppg077_staterdn[99],ppg077_syserup[99],ppg077_syserdn[99];
0407 
0408 ifstream finppg077;
0409 finppg077.open("ppg077_0-10.txt");
0410 if(!finppg077) {cerr << " ERROR: Cannot open input text file. " << endl;}
0411 finppg077 >> tmp0 >> tmp1 >> tmp2 >> tmp3 >> tmp4 >> tmp5;
0412 finppg077 >> tmp0 >> tmp1;
0413 finppg077 >> tmp0 >> tmp1 >> tmp2 >> tmp3 >> tmp4 >> tmp5;
0414 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]; }
0415 
0416 TGraphAsymmErrors* grRAA_ppg077 = new TGraphAsymmErrors(28,ppg077_pt,ppg077_raa,0,0,ppg077_staterdn,ppg077_staterup);
0417 TGraphAsymmErrors* grRAA_ppg077SysErr = new TGraphAsymmErrors(28,ppg077_pt,ppg077_raa,erpt2,erpt2,ppg077_syserdn,ppg077_syserup);
0418 
0419 
0420 //------------------------------------------------------------------------------------------
0421 // Invariant yield of HF electrons in Au+Au from ppg066, units are (GeV/c)^-2
0422 // Multiply it by pT, multiply by 2pi. 
0423 // This should be dN/dpT per event per unit of rapidity.
0424 // Then multiply by 2 since we are looking at +-1 unit of rapidity.
0425 // First two errors are statistical (up/down) third and fourth are systematic (up/down)
0426 // 0-10% central
0427 // ppg066 data are in form: 1/2pipT * dN/dpT/dy units are (GeV/c)^{-2}
0428 // these are heavy flavor electrons (c+b)
0429 //------------------------------------------------------------------------------------------
0430 double xau1[99],exau1[99],yau1[99],eyau1up[99],eyau1dn[99],eyau1upsys[99],eyau1dnsys[99];
0431 double ycau1[99],eycau1up[99],eycau1dn[99],eycau1upsys[99],eycau1dnsys[99];
0432 double ybau1[99],eybau1up[99],eybau1dn[99],eybau1upsys[99],eybau1dnsys[99];
0433 double xau2[99],exau2[99],yau2[99],eyau2up[99],eyau2dn[99],eyau2upsys[99],eyau2dnsys[99];
0434 double ycau2[99],eycau2up[99],eycau2dn[99],eycau2upsys[99],eycau2dnsys[99];
0435 double ybau2[99],eybau2up[99],eybau2dn[99],eybau2upsys[99],eybau2dnsys[99];
0436 
0437 ifstream finau1;
0438 finau1.open("ppg066_0-10.txt");
0439 if(!finau1) {cerr << " ERROR: Cannot open input text file. " << endl;}
0440 for(int i=0; i<28; i++) {
0441 finau1 >> xau1[i] >> yau1[i] >> eyau1up[i] >> eyau1dn[i] >> eyau1upsys[i] >> eyau1dnsys[i];
0442 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
0443 eyau1up[i]    = eyau1up[i]*    xau1[i]*2.*3.141592654*2. * NeventsAuAu;
0444 eyau1dn[i]    = eyau1dn[i]*    xau1[i]*2.*3.141592654*2. * NeventsAuAu;
0445 eyau1upsys[i] = eyau1upsys[i]* xau1[i]*2.*3.141592654*2. * NeventsAuAu;
0446 eyau1dnsys[i] = eyau1dnsys[i]* xau1[i]*2.*3.141592654*2. * NeventsAuAu;
0447 
0448 double mybottom = gBottomFraction->Eval(xau1[i]);
0449 double mycharm  = 1.-mybottom;
0450 ycau1[i]       = yau1[i]*       mycharm;  // multiply by c/(c+b) ratio to get charm yield 
0451 eycau1up[i]    = eyau1up[i]*    mycharm;
0452 eycau1dn[i]    = eyau1dn[i]*    mycharm;
0453 eycau1upsys[i] = eyau1upsys[i]* mycharm;
0454 eycau1dnsys[i] = eyau1dnsys[i]* mycharm;
0455 ybau1[i]       = yau1[i]*       mybottom;  // multiply by b/(c+b) ratio to get bottom yield
0456 eybau1up[i]    = eyau1up[i]*    mybottom;
0457 eybau1dn[i]    = eyau1dn[i]*    mybottom;
0458 eybau1upsys[i] = eyau1upsys[i]* mybottom;
0459 eybau1dnsys[i] = eyau1dnsys[i]* mybottom;
0460 
0461 }
0462 finau1.close();
0463 
0464   TGraphAsymmErrors *gHFAuAu      = new TGraphAsymmErrors(28,xau1,yau1,0,0,eyau1dn,eyau1up);
0465   TGraphAsymmErrors *gCharmAuAu   = new TGraphAsymmErrors(28,xau1,ycau1,0,0,eycau1dn,eycau1up);
0466   TGraphAsymmErrors *gBottomAuAu  = new TGraphAsymmErrors(28,xau1,ybau1,0,0,eybau1dn,eybau1up);
0467   TGraphAsymmErrors *gHFAuAuSysErr      = new TGraphAsymmErrors(28,xau1,yau1, erpt,erpt,eyau1dnsys,eyau1upsys);
0468   TGraphAsymmErrors *gCharmAuAuSysErr   = new TGraphAsymmErrors(28,xau1,ycau1,erpt,erpt,eycau1dnsys,eycau1upsys);
0469   TGraphAsymmErrors *gBottomAuAuSysErr  = new TGraphAsymmErrors(28,xau1,ybau1,erpt,erpt,eybau1dnsys,eybau1upsys);
0470 
0471 
0472 //-------------------------------------------------------------------------------------------
0473 // Invariant cross-section of HF electrons from ppg066 (run5 data) in p+p
0474 // Convert to invariant yield by dividing by 42mb (p+p total cross-section)
0475 //-------------------------------------------------------------------------------------------
0476   TGraphErrors *gsys_all4;
0477   TGraphErrors *gxs_all4;
0478   TFile *fpp = new TFile("/phenix/WWW/p/info/ppg/066/material/RUN5pp_best_fit.root");
0479     gxs_all4  = (TGraphErrors*)fpp->Get("gxs_all4");
0480     gsys_all4 = (TGraphErrors*)fpp->Get("gsys_all4");
0481   fpp->Close();
0482   int Ndatapp=gsys_all4->GetN();
0483   double *X1 = gxs_all4->GetX();
0484   double *Y1 = gxs_all4->GetY();
0485   double XX[99], YY[99], Estat[99], Esys[99];
0486   double YYC[99], YYB[99], ECstat[99], EBstat[99], ECsys[99], EBsys[99];
0487   for(int i=0; i<Ndatapp; i++) {
0488     XX[i] = X1[i]; 
0489     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
0490     Estat[i] = gxs_all4->GetErrorY(i)/42.e-03  * XX[i]*2.*3.141592654*2. * NeventsAuAu;
0491     Esys[i]  = gsys_all4->GetErrorY(i)/42.e-03 * XX[i]*2.*3.141592654*2. * NeventsAuAu;
0492     double mybottom = pythia_botfrac->Eval(XX[i]);   // PYTHIA
0493     double mycharm  = 1.-mybottom;
0494     //double mybottom = gfonllup->Eval(XX[i]); // use FONLL upper limit prediction to separate charm and bottom
0495     //double mycharm  = 1.-mybottom;
0496     //double mybottom = gfonll->Eval(XX[i]);  // use FONLL prediction to separate charm and bottom
0497     //double mycharm  = 1.-mybottom;
0498     YYC[i]    = YY[i]    * mycharm;
0499     ECstat[i] = Estat[i] * mycharm;
0500     ECsys[i]  = Esys[i]  * mycharm;
0501     YYB[i]    = YY[i]    * mybottom;
0502     EBstat[i] = Estat[i] * mybottom;
0503     EBsys[i]  = Esys[i]  * mybottom;
0504   }
0505   TGraphErrors *ppgraph_stat  = new TGraphErrors(Ndatapp,XX,YY,0,Estat);
0506   TGraphErrors *ppgraph_stat_charm  = new TGraphErrors(Ndatapp,XX,YYC,0,ECstat);
0507   TGraphErrors *ppgraph_stat_bottom  = new TGraphErrors(Ndatapp,XX,YYB,0,EBstat);
0508   TGraphErrors *ppgraph_sys  = new TGraphErrors(Ndatapp,XX,YY,0,Esys);
0509 
0510 
0511 
0512 
0513 
0514 
0515 
0516 //------------------------------------------------------------------------------------------
0517 // CCbar from PYTHIA
0518 //------------------------------------------------------------------------------------------
0519 
0520 TChain* chain_ccbar = new TChain("ntp2");
0521 
0522 // each of these files has 1.63386e+09 generated event
0523 // cross-section = 0.187161 mb
0524 // 10B 0-10% central Au+Au events before suppression corerespond to 9550B p+p events.
0525 // each file corresponds to 1.63386e+09 * 40/0.187161 = 349.19e+09 Min.Bias p+p events
0526 // need 9550/349.19 = 27.35 files. 
0527 // we can use those 27 files if we believe PYTHIA cross-section
0528 // PYTHIA ccbar cross-section scale up by 3.31994 +- 0.0403619 (fit above 2GeV); (was 2.09),
0529 // so we need 9080 runs (was 5716 runs)
0530 // all files below are for 5716 runs
0531 // 27 files (2700 runs)
0532 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_10.root"); // ckin3_3 is really not ckin3
0533 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_11.root");
0534 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_12.root");
0535 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_13.root");
0536 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_14.root");
0537 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_15.root");
0538 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_16.root");
0539 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_17.root");
0540 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_18.root");
0541 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_19.root");
0542 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_20.root");
0543 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_21.root");
0544 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_22.root");
0545 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_23.root");
0546 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_24.root");
0547 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_25.root");
0548 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_26.root");
0549 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_27.root");
0550 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_28.root");
0551 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_29.root");
0552 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_30.root");
0553 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_31.root");
0554 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_32.root");
0555 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_34.root");
0556 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_35.root");
0557 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_36.root");
0558 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_37.root");
0559 
0560 // 17 more files
0561 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_4.root");
0562 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_5.root");
0563 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_6.root");
0564 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_7.root");
0565 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_8.root");
0566 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_38.root");
0567 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_41.root");
0568 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_42.root");
0569 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_43.root");
0570 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_44.root");
0571 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_48.root");
0572 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_49.root");
0573 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_51.root");
0574 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_53.root");
0575 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_54.root");
0576 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_55.root");
0577 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_56.root");
0578 
0579 // six files below have 560 runs together (5.6 files)
0580 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_9.root");
0581 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_39.root");
0582 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_46.root");
0583 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_47.root");
0584 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_50.root");
0585 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_52.root");
0586 
0587 // two files below have together 388 runs (3.88 files)
0588 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_1.root");
0589 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_3/ccbar_ckin3_3_2.root");
0590 
0591 // ten files below are equivalent to 6 files above (six files = 3.6 files above)
0592 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar/ccbar_0.root");
0593 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar/ccbar_1.root");
0594 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar/ccbar_2.root");
0595 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar/ccbar_3.root");
0596 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar/ccbar_4.root");
0597 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar/ccbar_5.root");
0598 /*
0599 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar/ccbar_6.root");
0600 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar/ccbar_7.root");
0601 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar/ccbar_8.root");
0602 chain_ccbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar/ccbar_9.root");
0603 */
0604 
0605 // Single electron normalization.
0606 // PYTHIA cross-section = 0.18713mb
0607 // measured cs = 544 +/- 39(stat) +/- 142(syst) +/- 200(model) mkb; from ppg085, arXiv:0802.0050, PLB 670 iss. 4-5 p.313 (2009)
0608 // Files 0 through 9 have 1.67471e+07 generated p+p events each 
0609 // Files 10 through 14 have 6.6809e+07 generated events each
0610 // All files correspond to 40/0.18713 * (10*1.67471e+07 + 5*6.6809e+07) = 107.20B Min.Bias p+p events
0611 //   10B 0-10% central Au+Au events (Ncoll = 955) correspond to 9550B p+p events (before suppression)
0612 //   Normalization factor = 9550/107.20 = 89.1
0613 double charm_scale = 89.1;
0614 TChain* chain_ccbar_sngl = new TChain("ntp1");
0615 chain_ccbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_norm/ccbar_norm_0.root");
0616 chain_ccbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_norm/ccbar_norm_1.root");
0617 chain_ccbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_norm/ccbar_norm_2.root");
0618 chain_ccbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_norm/ccbar_norm_3.root");
0619 chain_ccbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_norm/ccbar_norm_4.root");
0620 chain_ccbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_norm/ccbar_norm_5.root");
0621 chain_ccbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_norm/ccbar_norm_6.root");
0622 chain_ccbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_norm/ccbar_norm_7.root");
0623 chain_ccbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_norm/ccbar_norm_8.root");
0624 chain_ccbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_norm/ccbar_norm_9.root");
0625 chain_ccbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_norm/ccbar_norm_10.root");
0626 chain_ccbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_norm/ccbar_norm_11.root");
0627 chain_ccbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_norm/ccbar_norm_12.root");
0628 chain_ccbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_norm/ccbar_norm_13.root");
0629 chain_ccbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_norm/ccbar_norm_14.root");
0630 
0631 cout << "Number of entries in CCbar chain: " << chain_ccbar->GetEntries() << endl;
0632 cout << "Number of entries in CCbar SINGLES chain: " << chain_ccbar_sngl->GetEntries() << endl;
0633 
0634 chain_ccbar->SetBranchAddress("mass",&mass);
0635 chain_ccbar->SetBranchAddress("pt",  &pt);
0636 chain_ccbar->SetBranchAddress("eta", &eta);
0637 chain_ccbar->SetBranchAddress("pt1", &pt1);
0638 chain_ccbar->SetBranchAddress("pt2", &pt2);
0639 chain_ccbar->SetBranchAddress("eta1",&eta1);
0640 chain_ccbar->SetBranchAddress("eta2",&eta2);
0641 
0642 chain_ccbar_sngl->SetBranchAddress("pt",  &pt_sngl);
0643 chain_ccbar_sngl->SetBranchAddress("eta", &eta_sngl);
0644 
0645 // Singles for normalization
0646 for(int j=0; j<chain_ccbar_sngl->GetEntries(); j++){
0647   chain_ccbar_sngl->GetEvent(j);
0648   if(j%100000==0) cout << "entry # " << j << " " << endl;
0649   if(fabs(eta_sngl)<1.0) { hcharm_pt_nosupp->Fill(pt_sngl); }
0650   double random = rnd->Uniform(0.,1.);
0651   double myraa = gCharmSuppression->Eval(pt_sngl);
0652   if(random<=myraa) {  if(fabs(eta_sngl)<1.0) { hcharm_pt->Fill(pt_sngl); } }
0653 }
0654 cout << "Number of entries in charm pT histogram = " << hcharm_pt->GetEntries() << endl;
0655 
0656 hcharm_pt->Scale(1.0/bsizept); // convert to dN/dpT
0657 hcharm_pt->Scale(charm_scale); // scale up to 10B most central AuAu according to cross-section
0658   // Normalize PYTHIA to measured charm
0659   fnorm_charm->SetParameter(0,1.);
0660   gCharmAuAu->Fit(fnorm_charm,"QRN0","",fitstart,fitstop);
0661   double norm_charm = fnorm_charm->GetParameter(0);
0662   cout << "charm normalization = " << norm_charm << endl;
0663 
0664 hcharm_pt_nosupp->Scale(1.0/bsizept); // convert to dN/dpT
0665 hcharm_pt_nosupp->Scale(charm_scale); // scale up to 10B most central AuAu according to cross-section
0666   fnorm_charm_nosupp->SetParameter(0,1.);
0667   ppgraph_stat_charm->Fit(fnorm_charm_nosupp,"QRN0","",fitstart,fitstop);
0668   double norm_charm_nosupp = fnorm_charm_nosupp->GetParameter(0);
0669   cout << "charm normalization (no suppression) = " << norm_charm_nosupp << endl;
0670 
0671 
0672 // invariant mass distributions for CCbar
0673 for(int j=0; j<chain_ccbar->GetEntries()*eideff*eideff; j++){
0674   chain_ccbar->GetEvent(j);
0675   if(j%20000==0) cout << "entry # " << j << endl;
0676   if(pt1>ptcut && pt2>ptcut && fabs(eta1)<1.0 && fabs(eta2)<1.0) {
0677     double random1 = rnd->Uniform(0.,1.);
0678     double myraa1 = gCharmSuppression->Eval(pt1);
0679     double random2 = rnd->Uniform(0.,1.);
0680     double myraa2 = gCharmSuppression->Eval(pt2);
0681     if(random1<=myraa1 && random2<myraa2) {
0682       int mybin = -1;
0683       for(int i=0; i<nbins; i++) { if(pt<binlim[i+1]) { mybin = i; break; } }
0684       if(mybin>=0 && mybin<nbins) { (hcharm[mybin])->Fill(mass); }
0685       (hcharm[nbins])->Fill(mass); // all pT
0686     }
0687   }
0688 }
0689 
0690 //------------------------------------------------------------------------------------------
0691 // CCbar with ckin(3) = 4 from PYTHIA
0692 //------------------------------------------------------------------------------------------
0693 
0694 TChain* chain_ccbar_ckin3_4 = new TChain("ntp2");
0695 // Each "run" has 4.98345e+06 generated events
0696 // Cross-section 0.825489e-03 mb
0697 // ccbar_ckin3_4/ccbar cross-section normalization factor = 241478.0e+09 / 9428.13e+09 = 25.61
0698 // Fit gives scaling factor = 21.6616   +/-  0.539134
0699 // // corerection for acceptance = 1.1824 (we need 1.1824 times more ckin3=4 data because of lower acceptance)
0700 // According to PYTHIA cross-section we need 9550/24147.8 * 1.1824 * 100 = 46.76 runs 
0701 // to get an equivalint of 10B 0-10% central Au+Au events.
0702 // Additional cross-section scaling (from measured HF cross-section) = 3.31994 +- 0.0403619 (was 2.09) 
0703 // This, we need 155.2 runs (was 97.7 runs).
0704 // 
0705 // The file below contains 49 runs, it corresponds to 10B 0-10% AuAu events (before suppression) if we use PYTHIA cross-section.
0706 //chain_ccbar_ckin3_4->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_4/ccbar_ckin3_4_9550B_ppevents.root");
0707 // The file below has 100 runs 
0708 chain_ccbar_ckin3_4->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_4/ccbar_ckin3_4_0.root");
0709 // the file below has 55 runs
0710 chain_ccbar_ckin3_4->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_ccbar_ckin3_4/ccbar_ckin3_4_88.root");
0711 
0712 
0713 cout << "Number of entries in CCbar with ckin(3) = 4 chain: " << chain_ccbar_ckin3_4->GetEntries() << endl;
0714 
0715 chain_ccbar_ckin3_4->SetBranchAddress("mass",&mass);
0716 chain_ccbar_ckin3_4->SetBranchAddress("pt",  &pt);
0717 chain_ccbar_ckin3_4->SetBranchAddress("eta", &eta);
0718 chain_ccbar_ckin3_4->SetBranchAddress("pt1", &pt1);
0719 chain_ccbar_ckin3_4->SetBranchAddress("pt2", &pt2);
0720 chain_ccbar_ckin3_4->SetBranchAddress("eta1",&eta1);
0721 chain_ccbar_ckin3_4->SetBranchAddress("eta2",&eta2);
0722 
0723 // invariant mass distributions for CCbar with ckin(3) = 4
0724 for(int j=0; j<chain_ccbar_ckin3_4->GetEntries()*eideff*eideff; j++){
0725   chain_ccbar_ckin3_4->GetEvent(j);
0726   if(j%20000==0) cout << "entry # " << j << endl;
0727   if(pt1>ptcut && pt2>ptcut && fabs(eta1)<1.0 && fabs(eta2)<1.0) {
0728     double random1 = rnd->Uniform(0.,1.);
0729     double myraa1 = gCharmSuppression->Eval(pt1);
0730     double random2 = rnd->Uniform(0.,1.);
0731     double myraa2 = gCharmSuppression->Eval(pt2);
0732     if(random1<=myraa1 && random2<myraa2) {
0733       int mybin = -1;
0734       for(int i=0; i<nbins; i++) { if(pt<binlim[i+1]) { mybin = i; break; } }
0735       if(mybin>=0 && mybin<nbins) { (hcharm_ckin3_4[mybin])->Fill(mass); }
0736       (hcharm_ckin3_4[nbins])->Fill(mass); // all pT
0737     }
0738   }
0739 }
0740 
0741 //------------------------------------------------------------------------------------------
0742 // BBbar from PYTHIA
0743 //------------------------------------------------------------------------------------------
0744 
0745 TChain* chain_bbbar = new TChain("ntp2");
0746 // PYTHIA cross-section = 0.000734778 mb
0747 // One "run" has 1.189e+06 generated events = 64.734e+09 MinBias p+p events.
0748 // So we need 147.5 runs if we believe pythia cross-section
0749 // But bbbar cross-section in PYTHIA seems to be 3.22271 +- 0.0388839 (was 2.66) times smaller than one measured by PHENIX
0750 // So we need 475.5 runs
0751 // The file below has 148 runs and can be used if we believe PYTHIA cross-section
0752 chain_bbbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_bbbar/bbbar_9550B_ppevents.root"); 
0753 // The file below has 114 runs 
0754 //chain_bbbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_bbbar/bbbar_7350B_ppevents.root"); 
0755 //
0756 // Files 5, 6, 7 contain 100 runs, file 8 contains 93 runs, file 88 contains 28 runs.
0757 // One run has 1.19143e+06 generated events
0758 // all files have 476 runs
0759 chain_bbbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_bbbar/bbbar_5.root"); 
0760 chain_bbbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_bbbar/bbbar_6.root"); 
0761 chain_bbbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_bbbar/bbbar_7.root"); 
0762 chain_bbbar->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_bbbar/bbbar_88.root"); 
0763 
0764 // Single bottom normalization.
0765 // PYTHIA cross-section = 0.00073482mb
0766 // FONLL cs = 1.87 +0.99 -0.67 mkb
0767 // 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]
0768 // measured dSigma/dY at midrapidity = 0.92 +0.34-0.31(stat) +0.39-0.36(sys) mkb [same]
0769 // measured from dielectroins cs = 3.9 +2.5(stat) +3-2(sys) mkb [A. Adare et al., Phys. Lett. B670, 313 (2009)]
0770 // File 0 has 6.17663e+06 generated p+p events
0771 // Corresponds to 40/0.00073482*6.17663e+06 = 336.23B minbias p+p events
0772 //   10B 0-10% central Au+Au events (Ncoll = 955) correspond to 9550B p+p events (before suppression)
0773 //   Bottom normalization factor = 9550/336.23 = 28.4
0774 double bottom_scale = 28.4;
0775 TChain* chain_bbbar_sngl = new TChain("ntp1");
0776 chain_bbbar_sngl->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_bbbar_norm/bbbar_norm_0.root"); 
0777 
0778 cout << "Number of entries in BBbar chain: " << chain_bbbar->GetEntries() << endl;
0779 cout << "Number of entries in BBbar SINGLES chain: " << chain_bbbar_sngl->GetEntries() << endl;
0780 
0781 chain_bbbar->SetBranchAddress("mass",&mass);
0782 chain_bbbar->SetBranchAddress("pt",  &pt);
0783 chain_bbbar->SetBranchAddress("eta", &eta);
0784 chain_bbbar->SetBranchAddress("pt1", &pt1);
0785 chain_bbbar->SetBranchAddress("pt2", &pt2);
0786 chain_bbbar->SetBranchAddress("eta1",&eta1);
0787 chain_bbbar->SetBranchAddress("eta2",&eta2);
0788 
0789 chain_bbbar_sngl->SetBranchAddress("pt",  &pt_sngl);
0790 chain_bbbar_sngl->SetBranchAddress("eta", &eta_sngl);
0791 
0792 // Singles for normalization
0793 for(int j=0; j<chain_bbbar_sngl->GetEntries(); j++){
0794   chain_bbbar_sngl->GetEvent(j);
0795   if(j%100000==0) cout << "entry # " << j << " " << endl;
0796   if(fabs(eta_sngl)<1.0) { hbottom_pt_nosupp->Fill(pt_sngl); }
0797   double random = rnd->Uniform(0.,1.);
0798   double myraa = gBottomSuppression->Eval(pt_sngl);
0799   if(random<=myraa) {  if(fabs(eta_sngl)<1.0) { hbottom_pt->Fill(pt_sngl); } }
0800 }
0801 cout << "Number of entries in bottom pT histogram = " << hbottom_pt->GetEntries() << endl;
0802 
0803 hbottom_pt->Scale(1.0/bsizept);  // convert to dN/dpT
0804 hbottom_pt->Scale(bottom_scale); // scale up to 10B most central AuAu according to cross-section
0805   // normalize PYTHIA to measured bottom
0806   fnorm_bottom->SetParameter(0,1.);
0807   gBottomAuAu->Fit(fnorm_bottom,"QRN0","",fitstart,fitstop);
0808   double norm_bottom = fnorm_bottom->GetParameter(0);
0809   cout << "bottom normalization = " << norm_bottom << endl;
0810 
0811 hbottom_pt_nosupp->Scale(1.0/bsizept);  // convert to dN/dpT
0812 hbottom_pt_nosupp->Scale(bottom_scale); // scale up to 10B most central AuAu according to cross-section
0813   fnorm_bottom_nosupp->SetParameter(0,1.);
0814   ppgraph_stat_bottom->Fit(fnorm_bottom_nosupp,"QRN0","",fitstart,fitstop);
0815   double norm_bottom_nosupp = fnorm_bottom_nosupp->GetParameter(0);
0816   cout << "bottom normalization (no suppression) = " << norm_bottom_nosupp << endl;
0817  
0818   // normalize PYTHIA to measured Heavy Flavor
0819   hhf_pt = (TH1D*)hbottom_pt->Clone("hhf_pt");
0820   hhf_pt->SetName("hhf_pt");
0821   hhf_pt->Add(hcharm_pt);
0822     fnorm_hf->SetParameter(0,1.);
0823     gHFAuAu->Fit(fnorm_hf,"QRN0","",fitstart,fitstop);
0824     double norm_hf = fnorm_hf->GetParameter(0);
0825     cout << "HF normalization = " << norm_hf << endl;
0826 
0827   // HF no suppression
0828   hhf_pt_nosupp = (TH1D*)hbottom_pt_nosupp->Clone("hhf_pt_nosupp");
0829   hhf_pt_nosupp->SetName("hhf_pt_nosupp");
0830   hhf_pt_nosupp->Add(hcharm_pt_nosupp);
0831     fnorm_hf_nosupp->SetParameter(0,1.);
0832     ppgraph_stat->Fit(fnorm_hf_nosupp,"QRN0","",fitstart,fitstop);
0833     double norm_hf_nosupp = fnorm_hf_nosupp->GetParameter(0);
0834     cout << "HF normalization (no suppression) = " << norm_hf_nosupp << endl;
0835 
0836 // invariant mass distributions for BBbar
0837 for(int j=0; j<chain_bbbar->GetEntries()*eideff*eideff; j++){
0838   chain_bbbar->GetEvent(j);
0839   if(j%20000==0) cout << "entry # " << j << endl;
0840   if(pt1>ptcut && pt2>ptcut && fabs(eta1)<1.0 && fabs(eta2)<1.0) {
0841     double random1 = rnd->Uniform(0.,1.);
0842     double myraa1 = gBottomSuppression->Eval(pt1);
0843     double random2 = rnd->Uniform(0.,1.);
0844     double myraa2 = gBottomSuppression->Eval(pt2);
0845     if(random1<=myraa1 && random2<=myraa2) {
0846       int mybin = -1;
0847       for(int i=0; i<nbins; i++) { if(pt<binlim[i+1]) { mybin = i; break; } }
0848       if(mybin>=0 && mybin<nbins) { (hbottom[mybin])->Fill(mass); }
0849       (hbottom[nbins])->Fill(mass);
0850     }
0851   }
0852 }
0853 
0854 
0855 //-----------------------------------------------------------------------
0856 //  DY from PYTHIA
0857 //----------------------------------------------------------------------
0858 TChain* chain_dy = new TChain("ntp2");
0859 // PYTHIA cross-section = 0.000118202 mb
0860 // so, one generated event corresponds to 40/0.000118202 = 338409.5 MinBias p+p events
0861 // one run has 0.465274e+06 generated events
0862 // so, 10B 0-10% central AuAu events (9550B p+p min.bias events) correspond to 60.65 runs; 
0863 // the file below has 61 runs.
0864 // do we need to scale up like charm and bottom?
0865 chain_dy->Add("/phenix/hhj/lebedev/sphenix/phpythia6/ntp_dy/dy_9550B_ppevents.root"); 
0866 
0867 cout << "Number of entries in DY chain: " << chain_dy->GetEntries() << endl;
0868 
0869 chain_dy->SetBranchAddress("mass",&mass);
0870 chain_dy->SetBranchAddress("pt",  &pt);
0871 chain_dy->SetBranchAddress("eta", &eta);
0872 chain_dy->SetBranchAddress("pt1", &pt1);
0873 chain_dy->SetBranchAddress("pt2", &pt2);
0874 chain_dy->SetBranchAddress("eta1",&eta1);
0875 chain_dy->SetBranchAddress("eta2",&eta2);
0876 
0877 // Invariant mass distributions for DY
0878 for(int j=0; j<chain_dy->GetEntries()*eideff*eideff; j++){
0879   chain_dy->GetEvent(j);
0880   if(j%20000==0) cout << "entry # " << j << endl;
0881   if(pt1>ptcut && pt2>ptcut && fabs(eta1)<1.0 && fabs(eta2)<1.0) { 
0882     int mybin = -1;
0883     for(int i=0; i<nbins; i++) { if(pt<binlim[i+1]) { mybin = i; break; } }
0884     if(mybin>=0 && mybin<nbins) { (hdy[mybin])->Fill(mass); }
0885     (hdy[nbins])->Fill(mass);
0886   }
0887 }
0888 
0889 
0890 
0891 //----------------------------------------------------------------------
0892 //  Write out
0893 //----------------------------------------------------------------------
0894 
0895 TFile* fout = new TFile(ofname.c_str(),"RECREATE");
0896 
0897   TH1D* hhcharm[nbins];
0898   TH1D* hhcharm_ckin3_4[nbins];
0899   TH1D* hhbottom[nbins];
0900   TH1D* hhdy[nbins];
0901   for(int i=0; i<nbins+1; i++) {
0902     sprintf(hname,"hhbottom_%d",i); hhbottom[i]  = (TH1D*)(hbottom[i])->Clone(hname);
0903     sprintf(hname,"hhcharm_%d",i);  hhcharm[i]   = (TH1D*)(hcharm[i])->Clone(hname);
0904     sprintf(hname,"hhcharm_ckin3_4_%d",i);  hhcharm_ckin3_4[i]   = (TH1D*)(hcharm_ckin3_4[i])->Clone(hname);
0905     sprintf(hname,"hhdy_%d",i);     hhdy[i]      = (TH1D*)(hdy[i])->Clone(hname);
0906   }
0907 
0908   TH1D* hhbottom_pt  = (TH1D*)hbottom_pt->Clone("hhbottom_pt");
0909   TH1D* hhcharm_pt  = (TH1D*)hcharm_pt->Clone("hhcharm_pt");
0910   TH1D* hhhf_pt = (TH1D*)hhbottom_pt->Clone("hhhf_pt");
0911   hhhf_pt->Add(hhcharm_pt);
0912   hhcharm_pt->SetLineColor(kRed+2);
0913   hhbottom_pt->SetLineColor(kBlue+1);
0914   hhhf_pt->SetLineColor(kBlack);
0915 
0916   TH1D* hhbottom_pt_nosupp  = (TH1D*)hbottom_pt_nosupp->Clone("hhbottom_pt_nosupp");
0917   TH1D* hhcharm_pt_nosupp  = (TH1D*)hcharm_pt_nosupp->Clone("hhcharm_pt_nosupp");
0918   TH1D* hhhf_pt_nosupp = (TH1D*)hhbottom_pt_nosupp->Clone("hhhf_pt_nosupp");
0919   hhhf_pt_nosupp->Add(hhcharm_pt_nosupp);
0920   hhcharm_pt_nosupp->SetLineColor(kRed+2);
0921   hhbottom_pt_nosupp->SetLineColor(kBlue+1);
0922   hhhf_pt_nosupp->SetLineColor(kBlack);
0923 
0924 
0925   //TGraph *ggFracFONLL_b  = new TGraph(Ndata,XR,YR); ggFracFONLL_b->SetName("ggFracFONLL_b");
0926   TGraph *ggFrac_c = new TGraph(21,ratcbpt,fracc); ggFrac_c->SetName("ggFrac_c");
0927   TGraph *ggFrac_b = new TGraph(21,ratcbpt,ratcb); ggFrac_b->SetName("ggFrac_b");
0928   TGraphAsymmErrors *ggCharmAuAu   = new TGraphAsymmErrors(28,xau1,ycau1,0,0,eycau1dn,eycau1up); ggCharmAuAu->SetName("ggCharmAuAu");
0929   TGraphAsymmErrors *ggBottomAuAu  = new TGraphAsymmErrors(28,xau1,ybau1,0,0,eybau1dn,eybau1up); ggBottomAuAu->SetName("ggBottomAuAu");
0930   TGraphAsymmErrors *ggHFAuAu  = new TGraphAsymmErrors(28,xau1,yau1,0,0,eyau1dn,eyau1up); ggHFAuAu->SetName("ggHFAuAu");
0931   TGraphAsymmErrors *ggCharmAuAuSysErr   = new TGraphAsymmErrors(28,xau1,ycau1,erpt,erpt,eycau1dnsys,eycau1upsys); ggCharmAuAuSysErr->SetName("ggCharmAuAuSysErr");
0932   TGraphAsymmErrors *ggBottomAuAuSysErr  = new TGraphAsymmErrors(28,xau1,ybau1,erpt,erpt,eybau1dnsys,eybau1upsys); ggBottomAuAuSysErr->SetName("ggBottomAuAuSysErr");
0933   TGraphAsymmErrors *ggHFAuAuSysErr      = new TGraphAsymmErrors(28,xau1,yau1, erpt,erpt,eyau1dnsys,eyau1upsys); ggHFAuAuSysErr->SetName("ggHFAuAuSysErr");
0934   TGraphAsymmErrors* ggrRAA_ppg077 = new TGraphAsymmErrors(28,ppg077_pt,ppg077_raa,0,0,ppg077_staterdn,ppg077_staterup); ggrRAA_ppg077->SetName("ggrRAA_ppg077");
0935   TGraphAsymmErrors* ggrRAA_ppg077SysErr = new TGraphAsymmErrors(28,ppg077_pt,ppg077_raa,erpt2,erpt2,ppg077_syserdn,ppg077_syserup); ggrRAA_ppg077SysErr->SetName("ggrRAA_ppg077SysErr");
0936 
0937   TGraph *ggRAA_charm  = new TGraph(9,pt182,raa182charm); ggRAA_charm->SetName("ggRAA_charm");
0938   TGraph *ggRAA_bottom = new TGraph(9,pt182,raa182bottom); ggRAA_bottom->SetName("ggRAA_bottom");
0939   TGraphAsymmErrors *ggRAASysErr_charm  = new TGraphAsymmErrors(6,pt182orig,raa182charmorig, 0,0,raa182charm_erdn, raa182charm_erup); ggRAASysErr_charm->SetName("ggRAASysErr_charm");
0940   TGraphAsymmErrors *ggRAASysErr_bottom = new TGraphAsymmErrors(6,pt182orig,raa182bottomorig,0,0,raa182bottom_erdn,raa182bottom_erup); ggRAASysErr_bottom->SetName("ggRAASysErr_bottom");
0941 
0942   ggFrac_b->SetLineColor(kMagenta);
0943   ggFrac_b->SetLineWidth(2);
0944 //  ggFracFONLL_b->SetLineColor(kBlack);
0945 //  ggFracFONLL_b->SetLineWidth(2);
0946 
0947   ggHFAuAu->SetMarkerStyle(20);
0948   ggHFAuAu->SetMarkerColor(kBlack);
0949   ggBottomAuAu->SetMarkerStyle(20);
0950   ggBottomAuAu->SetMarkerColor(kBlue);
0951   ggBottomAuAu->SetLineColor(kBlue);
0952   ggCharmAuAu->SetMarkerStyle(20);
0953   ggCharmAuAu->SetMarkerColor(kRed);
0954   ggCharmAuAu->SetLineColor(kRed);
0955 
0956   ggHFAuAuSysErr->SetFillColor(0);
0957   ggHFAuAuSysErr->SetMarkerStyle(6);
0958   ggHFAuAuSysErr->SetMarkerSize(0.5);
0959   ggHFAuAuSysErr->SetMarkerColor(kBlack);
0960   ggHFAuAuSysErr->SetLineColor(kBlack);
0961   ggHFAuAuSysErr->SetLineWidth(2);
0962   ggHFAuAuSysErr->SetFillColor(kGray+1);
0963 
0964   ggrRAA_ppg077SysErr->SetFillColor(0);
0965   ggrRAA_ppg077SysErr->SetMarkerStyle(6);
0966   ggrRAA_ppg077SysErr->SetMarkerSize(0.5);
0967   ggrRAA_ppg077SysErr->SetMarkerColor(kBlack);
0968   ggrRAA_ppg077SysErr->SetLineColor(kBlack);
0969   ggrRAA_ppg077SysErr->SetLineWidth(2);
0970   ggrRAA_ppg077SysErr->SetFillColor(kGray+1);
0971 
0972   ggrRAA_ppg077->SetLineWidth(2);
0973   ggrRAA_ppg077->SetMarkerStyle(20);
0974   ggrRAA_ppg077->SetMarkerSize(1.0);
0975   ggrRAA_ppg077->SetMarkerColor(kBlack);
0976 
0977   ggRAASysErr_charm->SetFillColor(kGreen+2);
0978   ggRAASysErr_charm->SetFillStyle(3001);
0979   ggRAA_charm->SetLineWidth(3);
0980   ggRAA_charm->SetLineColor(kGreen+2);
0981   ggRAASysErr_bottom->SetFillColor(kBlue);
0982   ggRAASysErr_bottom->SetFillStyle(3002);
0983   ggRAA_bottom->SetLineWidth(3);
0984   ggRAA_bottom->SetLineColor(kBlue);
0985 
0986 //  ggFracFONLL_b->Write();
0987   ggFrac_c->Write();
0988   ggFrac_b->Write();
0989   ggRAA_charm->Write();
0990   ggRAA_bottom->Write();
0991   ggRAASysErr_charm->Write();
0992   ggRAASysErr_bottom->Write();
0993   ggCharmAuAu->Write();
0994   ggBottomAuAu->Write();
0995   ggHFAuAu->Write();
0996   ggCharmAuAuSysErr->Write();
0997   ggBottomAuAuSysErr->Write();
0998   ggHFAuAuSysErr->Write();
0999   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
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!