PWG3/hfe/macros/MakeD2eSpectra.C

   1 #include <TF1.h>
   2 #include <TFormula.h>
   3 void MakeD2eSpectra(){
   4
   5   // load libs
   6   gSystem->Load("libANALYSIS.so");
   7   gSystem->Load("libANALYSISalice.so");
   8   gSystem->Load("libCORRFW.so");
   9   gSystem->Load("libPWG0base.so");
  10   gSystem->Load("/alidata10/alice_u/minjung/032.heavy/hfe/bg4conv/util/hfe/libPWG3hfeDevel.so");
  11
  12   setGeneralStyle();
  13
  14   // read input files
  15   //char filename0[]="/alidata10/alice_u/minjung/train/GSI/V006.MC_pp/2011-03-08_2238.5742/LHC10f6a/HFEtask.root"; // pwg3 talk(linear binning)
  16   //char filename1[]="../TableOldnorm/HFPtSpectrum_D0Kpi_rebinnedth_combinedFD_021210_DecayBRSubtracted.root";     // pwg3 talk(use of 71.4mb)
  17   //char filename2[]="../TableOldnorm/HFPtSpectrum_DplusKpipi_combinedFD_021210_DecayBRSubtracted.root";           // pwg3 talk(use of 71.4mb)
  18   //char filename0[]="/alidata10/alice_u/minjung/train/GSI/V006.MC_pp/2011-03-15_2121.5954/LHC10f6a/HFEtask.root";   // preliminary plot (log binning)
  19   char filename0[]="./HFEtaskLHC10f6awExtendedDpt.root";   // QM plot
  20   char filename1[]="../TableOldnorm/HFPtSpectrum_D0Kpi_combinedFD_rebinnedth_150311_newsigma_woDecayBR.root";      // preliminary plot(use of 62.1mb)
  21   char filename2[]="../TableOldnorm/HFPtSpectrum_DplusKpipi_combinedFD_150311_newsigma_woDecayBR.root";            // preliminary plot(use of 62.1mb)
  22
  23   TFile *_file[3];
  24   _file[0] = TFile::Open(filename0);
  25   _file[1] = TFile::Open(filename1);
  26   _file[2] = TFile::Open(filename2);
  27
  28   TList *tl = (TList *)_file[0]->Get("TRD_HFE_QA2")->FindObject("MCqa");
  29
  30   //count # of events
  31   TList *tl_result = (TList *)_file[0]->Get("TRD_HFE_Results2");
  32   AliCFContainer *containermc = (AliCFContainer *) tl_result->FindObject("eventContainer");
  33   AliCFDataGrid *mcGrid1 = new AliCFDataGrid("mcgrid1","",*containermc,AliHFEcuts::kEventStepGenerated);
  34   AliCFDataGrid *mcGrid2 = new AliCFDataGrid("mcgrid2","",*containermc,AliHFEcuts::kEventStepRecNoCut);
  35   printf("mc # of entries kEventStepGenerated= %lf kEventStepRecNoCut= %lf\n",mcGrid1->GetEntries(),mcGrid2->GetEntries());
  36
  37   Int_t nEvt = (Int_t) mcGrid1->GetEntries();
  38   Double_t nevtnorm = 1./float(nEvt); //normalize with total number of event
  39
  40
  41   // cross section for the old D meson data
  42   //Double_t Xsectrion4e = 71400000000.;
  43   //Double_t Xsectrion4D = 71400.; // consider the data point is represented as umb
  44   Double_t Xsectrion4e = 62100000000.;
  45   Double_t Xsectrion4D = 62100.; // consider the data point is represented as umb
  46
  47
  48   TString kDType[7];
  49   kDType[0]="421";
  50   kDType[1]="411";
  51   kDType[2]="431";
  52   kDType[3]="4122";
  53   kDType[4]="4132";
  54   kDType[5]="4232";
  55   kDType[6]="4332";
  56
  57   TString kDType2[4];
  58   kDType2[0]="";
  59   kDType2[1]="D0";
  60   kDType2[2]="Dp";
  61   kDType2[3]="Drest";
  62
  63   TString hname;
  64
  65   // 2D histo having D, e pt correlation
  66   TH2D *hDeCorr[4]; // 4 groups.  "", "D0", "Dp", "Drest"
  67   for(int i=0; i<4; i++){
  68     hname = "mcqa_barrel_PtCorr"+kDType2[i]+"_c";
  69     hDeCorr[i]=(TH2D*)tl->FindObject(hname);
  70   }
  71
  72   // D mesons' pt, y 2D histograms from PYTHIA
  73   TH2D *hDPYTHIA[7]; // 7 particles. 421:D0, 411:D+, 431:Ds, 4122, 4132, 4232, 4332
  74   for(int i=0; i<7; i++){
  75     hname = "Dmeson"+kDType[i];
  76     hDPYTHIA[i]=(TH2D*)tl->FindObject(hname);
  77   }
  78
  79   // PYTHIA pt distribution in |y|<0.5 of 7 D mesons
  80   TH1D *hDPYTHIApt[7];
  81   for(int i=0; i<7; i++){
  82     hname= "hD"+kDType[i]+"PYTHIApt";
  83     hDPYTHIApt[i]=(TH1D*)hDPYTHIA[i]->ProjectionX(hname,5,5); // |rapidity|<0.5
  84
  85     CorrectFromTheWidth(hDPYTHIApt[i]); // bin width correction
  86     hDPYTHIApt[i]->Scale(nevtnorm);     // # of event normalization
  87     hDPYTHIApt[i]->Scale(0.5);          // (particle + antiparticle)/2
  88     hDPYTHIApt[i]->Scale(Xsectrion4D);  // convert to cross section
  89     hDPYTHIApt[i]->SetMarkerStyle(20+i);
  90     hDPYTHIApt[i]->GetYaxis()->SetTitle("d#sigma/dp_{t} |_{|y|<0.5} [#mu b/GeV/c]");
  91   }
  92
  93   // D measured (be careful with binnings)
  94   TH1D *hD0MeasuredStat= (TH1D *)_file[1]->Get("histoSigmaCorr");                        // data & statistic errors, start from bin number 4
  95   TGraphAsymmErrors *gD0measuredSys= (TGraphAsymmErrors *)_file[1]->Get("gSigmaCorr");   // data & systematic errors, start from bin number 4
  96
  97   TH1D *hDpMeasuredStat = (TH1D *)_file[2]->Get("histoSigmaCorr");                       // data & statistic errors, start from bin number 1
  98   TGraphAsymmErrors *gDpMeasuredSys = (TGraphAsymmErrors *)_file[2]->Get("gSigmaCorr");  // data & statistic errors, start from bin number 1
  99
 100   // drawing -------------------------------------------------------
 101   cPtDmeson = new TCanvas("cPtDmeson","pT of D mesons",0,0,800,500);
 102   cPtDmeson->Divide(2,1);
 103
 104   cPtDmeson->cd(1);
 105   hDPYTHIApt[0]->Draw(""); //0: D0
 106   hD0MeasuredStat->Draw("same");
 107   TLegend *legD0 = new TLegend(0.60,0.75,0.75,0.85);
 108   legD0->AddEntry(hDPYTHIApt[0],"D^{0} PYTHIA","p");
 109   legD0->AddEntry(hD0MeasuredStat,"D^{0} Measured","p");
 110   setLegendStyle(*legD0,0);
 111   legD0->Draw("same");
 112
 113   cPtDmeson->cd(2);
 114   hDPYTHIApt[1]->Draw(""); //1: Dp
 115   hDpMeasuredStat->Draw("same");
 116   TLegend *legDp = new TLegend(0.60,0.75,0.75,0.85);
 117   legDp->AddEntry(hDPYTHIApt[1],"D^{+} PYTHIA","p");
 118   legDp->AddEntry(hDpMeasuredStat,"D^{+} Measured","p");
 119   setLegendStyle(*legDp,0);
 120   legDp->Draw("same");
 121   //----------------------------------------------------------------
 122
 123
 124   Double_t ptval[6];
 125   Double_t XsecD0[6], XsecD0ErrStatFrac[6], XsecD0ErrSysFracLow[6], XsecD0ErrSysFracHigh[6];
 126   Double_t XsecDp[6], XsecDpErrStatFrac[6], XsecDpErrSysFracLow[6], XsecDpErrSysFracHigh[6];
 127   Double_t iM2PD0[6], iM2PDp[6];
 128
 129   Double_t xbins[12]={0.,0.5,1.,2.,3.,4.,5.,6.,8.,12.,16.,20.};
 130   TH1D *hM2PD0= new TH1D("hM2PD0","D0: Measured/Pythia",11,xbins);
 131   TH1D *hM2PDp= new TH1D("hM2PDp","Dp: Measured/Pythia",11,xbins);
 132
 133   for(int i=0; i<6; i++){ // 0 -> first measured point, there are 6 points for pp for the moment
 134      // D0
 135      iM2PD0[i]=0.;
 136      ptval[i]=hD0MeasuredStat->GetBinCenter(i+4); // +4: starting from the first measured point
 137      XsecD0[i]=hD0MeasuredStat->GetBinContent(i+4);
 138      if(hDPYTHIApt[0]->GetBinContent(i+4)>0) iM2PD0[i]=hD0MeasuredStat->GetBinContent(i+4)/hDPYTHIApt[0]->GetBinContent(i+4);
 139      if(hD0MeasuredStat->GetBinContent(i+4)>0){
 140        XsecD0ErrStatFrac[i]    = hD0MeasuredStat->GetBinError(i+4)/hD0MeasuredStat->GetBinContent(i+4);
 141        XsecD0ErrSysFracLow[i]  = gD0measuredSys->GetErrorYlow(i+4)/hD0MeasuredStat->GetBinContent(i+4);
 142        XsecD0ErrSysFracHigh[i] = gD0measuredSys->GetErrorYhigh(i+4)/hD0MeasuredStat->GetBinContent(i+4);
 143      }
 144
 145      // D+
 146      iM2PDp[i]=0.;
 147      XsecDp[i]=hDpMeasuredStat->GetBinContent(i+1); // be careful with the fact that D+ measured starting from bin 1
 148      if(hDPYTHIApt[1]->GetBinContent(i+4)>0) iM2PDp[i]=hDpMeasuredStat->GetBinContent(i+1)/hDPYTHIApt[1]->GetBinContent(i+4); //D+
 149      if(hDpMeasuredStat->GetBinContent(i+1)>0) {
 150        XsecDpErrStatFrac[i]    = hDpMeasuredStat->GetBinError(i+1)/hDpMeasuredStat->GetBinContent(i+1);
 151        XsecDpErrSysFracLow[i]  = gDpMeasuredSys->GetErrorYlow(i+1)/hDpMeasuredStat->GetBinContent(i+1);
 152        XsecDpErrSysFracHigh[i] = gDpMeasuredSys->GetErrorYhigh(i+1)/hDpMeasuredStat->GetBinContent(i+1);
 153      }
 154
 155      hM2PD0->Fill(ptval[i],iM2PD0[i]);
 156      hM2PDp->Fill(ptval[i],iM2PDp[i]);
 157   }
 158   // drawing -------------------------------------------------------
 159   cM2PDmeson = new TCanvas("cM2PDmeson","Measured/Pythia",0,0,800,500);
 160   cM2PDmeson->Divide(2,1);
 161   cM2PDmeson->cd(1);
 162   hM2PD0->SetMarkerStyle(20);
 163   hM2PD0->SetXTitle("p_{t} [GeV/c]");
 164   hM2PD0->SetYTitle("ratio");
 165   hM2PD0->Draw("p");
 166   TLegend *legM2PD0 = new TLegend(0.45,0.75,0.75,0.85);
 167   legM2PD0->AddEntry(hM2PD0,"D0: Measured/PYTHIA","p");
 168   setLegendStyle(*legM2PD0,0);
 169   legM2PD0->Draw("same");
 170
 171   cM2PDmeson->cd(2);
 172   hM2PDp->SetMarkerStyle(20);
 173   hM2PDp->SetXTitle("p_{t} [GeV/c]");
 174   hM2PDp->SetYTitle("ratio");
 175   hM2PDp->Draw("p");
 176   TLegend *legM2PDp = new TLegend(0.45,0.75,0.75,0.85);
 177   legM2PDp->AddEntry(hM2PDp,"D+: Measured/PYTHIA","p");
 178   setLegendStyle(*legM2PDp,0);
 179   legM2PDp->Draw("same");
 180   //----------------------------------------------------------------
 181
 182
 183   // electrons pt spectra from given D meson pt bins (before any scaling)
 184   TString kDPtbin[14];
 185   kDPtbin[0]="Dpt0";
 186   kDPtbin[1]="Dpt05";
 187   kDPtbin[2]="Dpt1";
 188   kDPtbin[3]="Dpt2";
 189   kDPtbin[4]="Dpt3";
 190   kDPtbin[5]="Dpt4";
 191   kDPtbin[6]="Dpt5";
 192   kDPtbin[7]="Dpt6";
 193   kDPtbin[8]="Dpt8";
 194   kDPtbin[9]="Dpt12";
 195   kDPtbin[10]="Dpt16";
 196   kDPtbin[11]="Dpt20";
 197   kDPtbin[12]="Dpt30";
 198   kDPtbin[13]="Dpt40";
 199
 200   Int_t kbinl[14]={1, 6,11,21,31,41,51,61, 81,121,161,201,301,401};
 201   Int_t kbinh[14]={5,10,20,30,40,50,60,80,120,160,200,300,400,500};
 202
 203   TH1D *hD2e[4][14];
 204   TH1D *hD2eNorm[4][14];
 205   TH1D *hD2eNormWm[2][14];
 206   Int_t nRebin = 1;
 207   for(int iDtype=0; iDtype<2; iDtype++){ // 4 groups.  0:"", 1:"D0", 2:"Dp", 3:"Drest"
 208     for(int iDptbin=0; iDptbin<14; iDptbin++){
 209       hname = "h"+kDType2[iDtype+1]+kDPtbin[iDptbin];
 210       hD2e[iDtype][iDptbin] = (TH1D*)hDeCorr[iDtype+1]->ProjectionY(hname,kbinl[iDptbin],kbinh[iDptbin]);
 211       hD2e[iDtype][iDptbin]->Rebin(nRebin);
 212       CorrectFromTheWidth(hD2e[iDtype][iDptbin]);
 213       hname = "hNorm"+kDType2[iDtype]+kDPtbin[iDptbin];
 214       hD2eNorm[iDtype][iDptbin] = (TH1D*)hD2e[iDtype][iDptbin]->Clone(hname);
 215       hD2eNormWm[iDtype][iDptbin] = (TH1D*)hD2e[iDtype][iDptbin]->Clone(hname);
 216       if(iDtype==0) hD2eNorm[2][iDptbin] = (TH1D*)hD2e[iDtype][iDptbin]->Clone(hname); // for Ds, use D0 electron spectra
 217       if(iDtype==0) hD2eNorm[3][iDptbin] = (TH1D*)hD2e[iDtype][iDptbin]->Clone(hname); // for Lc, use D0 electron spectra
 218     }
 219   }
 220
 221   hD0pt0to2=(TH1D*)hD2e[0][0]->Clone("hD0pt0to2");
 222   hD0pt2to12=(TH1D*)hD2e[0][3]->Clone("hD0pt2to12");
 223   hD0pt12to50=(TH1D*)hD2e[0][9]->Clone("hD0pt12to50");
 224   hD0pt0to2->Add(hD2e[0][1]);
 225   hD0pt0to2->Add(hD2e[0][2]);
 226   hD0pt2to12->Add(hD2e[0][4]);
 227   hD0pt2to12->Add(hD2e[0][5]);
 228   hD0pt2to12->Add(hD2e[0][6]);
 229   hD0pt2to12->Add(hD2e[0][7]);
 230   hD0pt2to12->Add(hD2e[0][8]);
 231   hD0pt12to50->Add(hD2e[0][10]);
 232   hD0pt12to50->Add(hD2e[0][11]);
 233   hD0pt12to50->Add(hD2e[0][12]);
 234   hD0pt12to50->Add(hD2e[0][13]);
 235   hD0ptall=(TH1D*)hD0pt0to2->Clone("hD0ptall");
 236   hD0ptall->Add(hD0pt2to12);
 237   hD0ptall->Add(hD0pt12to50);
 238
 239   hD0lowptcontrib=(TH1D*)hD0pt0to2->Clone("hD0lowptcontrib");
 240   hD0lowptcontrib->Divide(hD0pt2to12);
 241
 242   hD0pt0to2->Divide(hD0ptall);
 243   hD0pt2to12->Divide(hD0ptall);
 244   hD0pt12to50->Divide(hD0ptall);
 245
 246   new TCanvas;
 247   hD0pt0to2->SetMarkerStyle(24);
 248   hD0pt2to12->SetMarkerStyle(24);
 249   hD0pt12to50->SetMarkerStyle(24);
 250   hD0pt0to2->SetMarkerColor(1);
 251   hD0pt0to2->SetLineColor(1);
 252   hD0pt2to12->SetMarkerColor(4);
 253   hD0pt2to12->SetLineColor(4);
 254   hD0pt12to50->SetMarkerColor(2);
 255   hD0pt12to50->SetLineColor(2);
 256   hD0pt0to2->Draw();
 257   hD0pt2to12->Draw("same");
 258   hD0pt12to50->Draw("same");
 259   gPad->SetLogy();
 260   gPad->SetGrid();
 261
 262   // drawing -------------------------------------------------------
 263   cPtD2e = new TCanvas("cPtD2e","pT of e from certain D pt bin",0,0,1000,700);
 264   cPtD2e->Divide(2,1);
 265   Int_t colorcode=1;
 266   TLegend *legD2e[4];
 267   for(int iDtype=0; iDtype<2; iDtype++){
 268     legD2e[iDtype] = new TLegend(0.75,0.45,0.88,0.88);
 269     legD2e[iDtype]->AddEntry(hD2e[0][0],kDType2[iDtype+1],"");
 270     for(int iDptbin=0; iDptbin<14; iDptbin++){
 271       cPtD2e->cd(iDtype+1);
 272       colorcode=iDptbin+1;
 273       if(colorcode==10) colorcode=38;
 274       hD2e[iDtype][iDptbin]->SetLineColor(colorcode);
 275       hD2e[iDtype][iDptbin]->SetMarkerColor(colorcode);
 276       hD2e[iDtype][iDptbin]->SetMarkerStyle(20);
 277       if(iDptbin==0) hD2e[iDtype][iDptbin]->Draw();
 278       else hD2e[iDtype][iDptbin]->Draw("samep");
 279       gPad->SetLogy();
 280       gPad->SetLogx();
 281       legD2e[iDtype]->AddEntry(hD2e[iDtype][iDptbin],kDPtbin[iDptbin],"p");
 282       setLegendStyle(*legD2e[iDtype],0);
 283     }
 284     legD2e[iDtype]->Draw("same");
 285   }
 286   //----------------------------------------------------------------
 287
 288   // electrons pt spectra from given D meson pt bins "with scaling based on measured X section"
 289   Double_t accfactorD[4][6] = { {1.73018, 1.72687, 1.76262, 1.94291, 1.82095, 1.69734,}, // acceptance factor of D0->e
 290                                 {1.77726, 1.72664, 1.79331, 1.73913, 1.71084, 1.74054,}, // acceptance factor of Dp->e
 291                                 {1.77726, 1.72664, 1.79331, 1.73913, 1.71084, 1.74054,}, // acceptance factor of Ds->e(copy of Dp)
 292                                 {1.77726, 1.72664, 1.79331, 1.73913, 1.71084, 1.74054} }; // acceptance factor of Lc->e(copy of Dp)
 293   Double_t ptbwdth[6] = {1., 1., 1., 1., 2., 4.}; // pt bin width for binwidth correction
 294
 295   // normalize PYTHIA electron spectra with the measured D meson cross section * branching ratio
 296   Double_t brachratioD[4];
 297   brachratioD[0] = 0.0653; // D0->e branching ratio
 298   brachratioD[1] = 0.16;   // Dp->e branching ratio
 299   brachratioD[2] = 0.08;   // Ds->e branching ratio
 300   brachratioD[3] = 0.045;  // Lc->e branching ratio
 301   //Double_t DiffMthdScaleD0= 1./1.07;  // considering difference between weighting method and X section method for D0 : new
 302   //Double_t DiffMthdScaleDp= 1./1.19;  // considering difference between weighting method and X section method for Dp : new
 303   Double_t DiffMthdScaleD0= 1./1.12;  // considering difference between weighting method and X section method for D0
 304   Double_t DiffMthdScaleDp= 1./1.25;  // considering difference between weighting method and X section method for Dp
 305
 306
 307   Double_t norm, xsec;
 308   Double_t xsecD[4][6];
 309   for(int iptb=0; iptb<6; iptb++){
 310     xsecD[0][iptb] = XsecD0[iptb]; //D0
 311     xsecD[1][iptb] = XsecDp[iptb]; //Dp
 312     xsecD[2][iptb] = (XsecD0[iptb]+XsecDp[iptb])*2.37/(8.49+2.65+5.07); //Ds: use of ZEUS ratio
 313     xsecD[3][iptb] = (XsecD0[iptb]+XsecDp[iptb])*3.59/(8.49+2.65+5.07); //Lc: use of ZEUS ratio
 314   }
 315
 316   for(int iDtype=0; iDtype<4; iDtype++){
 317     for(int iptb=0; iptb<6; iptb++){
 318       xsec = xsecD[iDtype][iptb]*1000000.; // change scale from ub to pb
 319       norm = hD2eNorm[iDtype][iptb+3]->Integral("width")/(xsec*brachratioD[iDtype]*ptbwdth[iptb]); //D0
 320       norm = accfactorD[iDtype][iptb]/norm;
 321       hD2eNorm[iDtype][iptb+3]->Scale(norm);
 322     }
 323   }
 324
 325   // not measured pt data points for D0
 326   hD2eNorm[0][9]->Scale(nevtnorm);
 327   hD2eNorm[0][9]->Scale(0.44);  // measured/pythia for this pt bin
 328   hD2eNorm[0][9]->Scale(0.5);
 329   hD2eNorm[0][9]->Scale(DiffMthdScaleD0);
 330   hD2eNorm[0][9]->Scale(Xsectrion4e);
 331   hD2eNorm[0][10]->Scale(nevtnorm);
 332   hD2eNorm[0][10]->Scale(0.44); // measured/pythia for this pt bin
 333   hD2eNorm[0][10]->Scale(0.5);
 334   hD2eNorm[0][10]->Scale(DiffMthdScaleD0);
 335   hD2eNorm[0][10]->Scale(Xsectrion4e);
 336   hD2eNorm[0][11]->Scale(nevtnorm);
 337   hD2eNorm[0][11]->Scale(0.44); // measured/pythia for this pt bin
 338   hD2eNorm[0][11]->Scale(0.5);
 339   hD2eNorm[0][11]->Scale(DiffMthdScaleD0);
 340   hD2eNorm[0][11]->Scale(Xsectrion4e);
 341   hD2eNorm[0][12]->Scale(nevtnorm);
 342   hD2eNorm[0][12]->Scale(0.44); // measured/pythia for this pt bin
 343   hD2eNorm[0][12]->Scale(0.5);
 344   hD2eNorm[0][12]->Scale(DiffMthdScaleD0);
 345   hD2eNorm[0][12]->Scale(Xsectrion4e);
 346   hD2eNorm[0][13]->Scale(nevtnorm);
 347   hD2eNorm[0][13]->Scale(0.44); // measured/pythia for this pt bin
 348   hD2eNorm[0][13]->Scale(0.5);
 349   hD2eNorm[0][13]->Scale(DiffMthdScaleD0);
 350   hD2eNorm[0][13]->Scale(Xsectrion4e);
 351
 352
 353   // not measured pt data points for D+
 354   hD2eNorm[1][9]->Scale(nevtnorm);
 355   hD2eNorm[1][9]->Scale(0.41);  // measured/pythia for this pt bin
 356   hD2eNorm[1][9]->Scale(0.5);
 357   hD2eNorm[1][9]->Scale(DiffMthdScaleDp);
 358   hD2eNorm[1][9]->Scale(Xsectrion4e);
 359   hD2eNorm[1][10]->Scale(nevtnorm);
 360   hD2eNorm[1][10]->Scale(0.41); // measured/pythia for this pt bin
 361   hD2eNorm[1][10]->Scale(0.5);
 362   hD2eNorm[1][10]->Scale(DiffMthdScaleDp);
 363   hD2eNorm[1][10]->Scale(Xsectrion4e);
 364   hD2eNorm[1][11]->Scale(nevtnorm);
 365   hD2eNorm[1][11]->Scale(0.41); // measured/pythia for this pt bin
 366   hD2eNorm[1][11]->Scale(0.5);
 367   hD2eNorm[1][11]->Scale(DiffMthdScaleDp);
 368   hD2eNorm[1][11]->Scale(Xsectrion4e);
 369   hD2eNorm[1][12]->Scale(nevtnorm);
 370   hD2eNorm[1][12]->Scale(0.41); // measured/pythia for this pt bin
 371   hD2eNorm[1][12]->Scale(0.5);
 372   hD2eNorm[1][12]->Scale(DiffMthdScaleDp);
 373   hD2eNorm[1][12]->Scale(Xsectrion4e);
 374   hD2eNorm[1][13]->Scale(nevtnorm);
 375   hD2eNorm[1][13]->Scale(0.41); // measured/pythia for this pt bin
 376   hD2eNorm[1][13]->Scale(0.5);
 377   hD2eNorm[1][13]->Scale(DiffMthdScaleDp);
 378   hD2eNorm[1][13]->Scale(Xsectrion4e);
 379
 380
 381   // estimation for Ds and Lc based on D0, D+ at high pt
 382   Double_t fragRatio=2.37/(8.49+2.65+5.07);
 383   Double_t brachRatio[2];
 384   Double_t brachRatio[0]=2*brachratioD[2]/(brachratioD[0]+brachratioD[1]);
 385   Double_t brachRatio[1]=2*brachratioD[3]/(brachratioD[0]+brachratioD[1]);
 386   for(int iptb=9; iptb<14; iptb++){
 387     hD2eNorm[2][iptb]->Add(hD2eNorm[0][iptb]);
 388     hD2eNorm[2][iptb]->Add(hD2eNorm[1][iptb]);
 389     hD2eNorm[2][iptb]->Scale(fragRatio);
 390     hD2eNorm[2][iptb]->Scale(brachRatio[0]);
 391     hD2eNorm[3][iptb]->Add(hD2eNorm[0][iptb]);
 392     hD2eNorm[3][iptb]->Add(hD2eNorm[1][iptb]);
 393     hD2eNorm[3][iptb]->Scale(fragRatio);
 394     hD2eNorm[3][iptb]->Scale(brachRatio[1]);
 395   }
 396
 397
 398   int kptmaxb = 11; // D0, D+
 399   TH1D *hD2eNormSum[4];
 400   for(int iDtype=0; iDtype<4; iDtype++){
 401     hname="hD2eNormSum" + iDtype;
 402     hD2eNormSum[iDtype]=(TH1D*)hD2eNorm[iDtype][3]->Clone(hname);
 403 //    if(iDtype==2 || iDtype==3) kptmaxb=6; //Ds, Lc    //mmjj
 404     for(int iptb=1; iptb<kptmaxb; iptb++){
 405       hD2eNormSum[iDtype]->Add(hD2eNorm[iDtype][iptb+3]);
 406     }
 407   }
 408
 409   for(int i=1; i<hD2eNormSum[0]->FindBin(2.0); i++){ // consider low pt contribution
 410      for(int iDtype=0; iDtype<4; iDtype++){
 411        Double_t icontrib =hD0lowptcontrib->GetBinContent(i)*hD2eNormSum[iDtype]->GetBinContent(i);
 412        Double_t ipt =hD2eNormSum[iDtype]->GetBinCenter(i);
 413        hD2eNormSum[iDtype]->Fill(ipt, icontrib);
 414      }
 415   }
 416
 417
 418   TH1D *hD2eNormSumStat[4]; // copy of hD2eNormSum[4] to get mc statistical errors
 419   for(int iDtype=0; iDtype<4; iDtype++){
 420     hname="hD2eNormSumStat" + iDtype;
 421     hD2eNormSumStat[iDtype]=(TH1D*)hD2eNormSum[iDtype]->Clone(hname);
 422   }
 423
 424   // sum
 425   hD2eNormTotalStat=(TH1D*)hD2eNormSumStat[0]->Clone("hD2eNormTotalStat");
 426   hD2eNormTotalStat->Add(hD2eNormSumStat[1]);
 427   hD2eNormTotalStat->Add(hD2eNormSumStat[2]);
 428   hD2eNormTotalStat->Add(hD2eNormSumStat[3]);
 429
 430   // drawing, here the error bars are only for MC statistics ---------------------------
 431   cPtD2eSum = new TCanvas("cPtD2eSum","pT of e from D",0,0,1000,600);
 432   cPtD2eSum->Divide(2,1);
 433   cPtD2eSum->cd(1);
 434   hD2eNormSumStat[0]->GetYaxis()->SetTitle("BRxd#sigma/dp_{t} | |#eta|<0.9  [pb/GeV/c]");
 435   hD2eNormSumStat[0]->SetMarkerStyle(20);
 436   hD2eNormSumStat[0]->Draw();
 437   hD2eNormSumStat[1]->SetMarkerColor(2);
 438   hD2eNormSumStat[1]->SetMarkerStyle(20);
 439   hD2eNormSumStat[1]->Draw("same");
 440   hD2eNormSumStat[2]->SetMarkerColor(4);
 441   hD2eNormSumStat[2]->SetMarkerStyle(20);
 442   hD2eNormSumStat[2]->Draw("same");
 443   hD2eNormSumStat[3]->SetMarkerColor(3);
 444   hD2eNormSumStat[3]->SetMarkerStyle(20);
 445   hD2eNormSumStat[3]->Draw("same");
 446   gPad->SetLogy();
 447   TLegend *legsum = new TLegend(0.25,0.70,0.75,0.85);
 448   legsum->AddEntry(hD2eNormSumStat[0],"D^{0}#rightarrow e 2.0<mother p_{t}<50.0 GeV/c","p");
 449   legsum->AddEntry(hD2eNormSumStat[1],"D^{+}#rightarrow e 2.0<mother p_{t}<50.0 GeV/c","p");
 450   legsum->AddEntry(hD2eNormSumStat[2],"D_{s}#rightarrow e 2.0<mother p_{t}<50.0 GeV/c","p");
 451   legsum->AddEntry(hD2eNormSumStat[3],"#Lambda_{c}#rightarrow e 2.0<mother p_{t}<50.0 GeV/c","p");
 452   setLegendStyle(*legsum,0);  legsum->Draw("same");
 453
 454   cPtD2eSum->cd(2);
 455   hD2eNormTotalStat->SetMarkerStyle(20);
 456   hD2eNormTotalStat->GetYaxis()->SetTitle("BRxd#sigma/dp_{t} | |#eta|<0.9  [pb/GeV/c]");
 457   hD2eNormTotalStat->Draw();
 458   TLegend *legtotal = new TLegend(0.30,0.75,0.75,0.85);
 459   legtotal->AddEntry(hD2eNormTotalStat,"D#rightarrow e 2.0<mother p_{t}<50.0 GeV/c","p");
 460   setLegendStyle(*legtotal,0);
 461   legtotal->Draw("same");
 462   gPad->SetLogy();
 463   //----------------------------------------------------------------
 464
 465
 466   // calculating propagated error of D meson's statistic and systematic errors
 467   TGraphAsymmErrors *gD2eNormSum[4]; // store propagated error from D systematic errors
 468   for(iDtype=0; iDtype<4; iDtype++){
 469     gD2eNormSum[iDtype]= new TGraphAsymmErrors(hD2eNormSum[iDtype]);
 470   }
 471
 472   Double_t statErr[11]; //kptmaxb
 473   Double_t sysErrLow[11]; //kptmaxb
 474   Double_t sysErrHigh[11]; //kptmaxb
 475   Double_t statErrSum[4]; // D0, D+, Ds. Lc
 476   Double_t sysErrLowSum[4]; // D0, D+, Ds. Lc
 477   Double_t sysErrHighSum[4]; // D0, D+, Ds. Lc
 478   Int_t ib4err = 0;
 479   for(int i=0; i<hD2eNormSum[0]->GetNbinsX(); i++){
 480      for(int iDtype=0; iDtype<4; iDtype++){
 481        statErrSum[iDtype]=0.;
 482        sysErrLowSum[iDtype]=0.;
 483        sysErrHighSum[iDtype]=0.;
 484      }
 485
 486      for(int iptb=0; iptb<11; iptb++){
 487        // calculate propagated error from D statistical errors
 488        if(iptb<6) ib4err = iptb;
 489        else ib4err = 5; // last two bins from pythia, so take the error from the last measured bin
 490        // statistical error of D propagated to systematic error of elec
 491        statErr[iptb]=hD2eNorm[0][iptb+3]->GetBinContent(i+1)*XsecD0ErrStatFrac[ib4err]; //D0
 492        if(!(iptb<6)) statErr[iptb] = statErr[iptb]*1.5; //1.5 put conservative error
 493        statErrSum[0] += statErr[iptb]*statErr[iptb];
 494        statErr[iptb]=hD2eNorm[1][iptb+3]->GetBinContent(i+1)*XsecDpErrStatFrac[ib4err]; //Dp
 495        if(!(iptb<6)) statErr[iptb] = statErr[iptb]*1.5;
 496        statErrSum[1] += statErr[iptb]*statErr[iptb];
 497        statErr[iptb]=hD2eNorm[2][iptb+3]->GetBinContent(i+1)*XsecDpErrStatFrac[ib4err]*1.5; //Ds
 498        if(!(iptb<6)) statErr[iptb] = statErr[iptb]*1.5;
 499        statErrSum[2] += statErr[iptb]*statErr[iptb];
 500        statErr[iptb]=hD2eNorm[3][iptb+3]->GetBinContent(i+1)*XsecDpErrStatFrac[ib4err]*1.5; //Lc
 501        if(!(iptb<6)) statErr[iptb] = statErr[iptb]*1.5;
 502        statErrSum[3] += statErr[iptb]*statErr[iptb];
 503
 504        // systematic error of D propagated to systematic error of elec
 505        sysErrLow[iptb]=hD2eNorm[0][iptb+3]->GetBinContent(i+1)*XsecD0ErrSysFracLow[ib4err]; //D0
 506        if(!(iptb<6)) sysErrLow[iptb] = sysErrLow[iptb]*1.5;
 507        sysErrLowSum[0] += sysErrLow[iptb]*sysErrLow[iptb];
 508        sysErrLow[iptb]=hD2eNorm[1][iptb+3]->GetBinContent(i+1)*XsecDpErrSysFracLow[ib4err]; //Dp
 509        if(!(iptb<6)) sysErrLow[iptb] = sysErrLow[iptb]*1.5;
 510        sysErrLowSum[1] += sysErrLow[iptb]*sysErrLow[iptb];
 511        sysErrLow[iptb]=hD2eNorm[2][iptb+3]->GetBinContent(i+1)*XsecDpErrSysFracLow[ib4err]*1.5; //Ds
 512        if(!(iptb<6)) sysErrLow[iptb] = sysErrLow[iptb]*1.5;
 513        sysErrLowSum[2] += sysErrLow[iptb]*sysErrLow[iptb];
 514        sysErrLow[iptb]=hD2eNorm[3][iptb+3]->GetBinContent(i+1)*XsecDpErrSysFracLow[ib4err]*1.5; //Lc
 515        if(!(iptb<6)) sysErrLow[iptb] = sysErrLow[iptb]*1.5;
 516        sysErrLowSum[3] += sysErrLow[iptb]*sysErrLow[iptb];
 517
 518        sysErrHigh[iptb]=hD2eNorm[0][iptb+3]->GetBinContent(i+1)*XsecD0ErrSysFracHigh[ib4err]; //D0
 519        if(!(iptb<6)) sysErrHigh[iptb] = sysErrHigh[iptb]*1.5;
 520        sysErrHighSum[0] += sysErrHigh[iptb]*sysErrHigh[iptb];
 521        sysErrHigh[iptb]=hD2eNorm[1][iptb+3]->GetBinContent(i+1)*XsecDpErrSysFracHigh[ib4err]; //Dp
 522        if(!(iptb<6)) sysErrHigh[iptb] = sysErrHigh[iptb]*1.5;
 523        sysErrHighSum[1] += sysErrHigh[iptb]*sysErrHigh[iptb];
 524        sysErrHigh[iptb]=hD2eNorm[2][iptb+3]->GetBinContent(i+1)*XsecDpErrSysFracHigh[ib4err]*1.5; //Ds
 525        if(!(iptb<6)) sysErrHigh[iptb] = sysErrHigh[iptb]*1.5;
 526        sysErrHighSum[2] += sysErrHigh[iptb]*sysErrHigh[iptb];
 527        sysErrHigh[iptb]=hD2eNorm[3][iptb+3]->GetBinContent(i+1)*XsecDpErrSysFracHigh[ib4err]*1.5; //Lc
 528        if(!(iptb<6)) sysErrHigh[iptb] = sysErrHigh[iptb]*1.5;
 529        sysErrHighSum[3] += sysErrHigh[iptb]*sysErrHigh[iptb];
 530      }
 531
 532      for(iDtype=0; iDtype<4; iDtype++){
 533        hD2eNormSum[iDtype]->SetBinError(i+1,TMath::Sqrt(statErrSum[iDtype]));
 534        gD2eNormSum[iDtype]->SetPointError(i,0,0,TMath::Sqrt(sysErrLowSum[iDtype]),TMath::Sqrt(sysErrHighSum[iDtype]));
 535      }
 536   }
 537
 538   // sum -----------------------------------------------------------
 539   hD2eNormTotal=(TH1D*)hD2eNormSum[0]->Clone("hD2eNormTotal");
 540   hD2eNormTotal->Add(hD2eNormSum[1]);
 541   hD2eNormTotal->Add(hD2eNormSum[2]);
 542   hD2eNormTotal->Add(hD2eNormSum[3]);
 543
 544   // sum up errors -------------------------------------------------
 545   TGraphAsymmErrors *gD2eNormTotal = new TGraphAsymmErrors(hD2eNormTotal);
 546   TGraphAsymmErrors *gD2eNormTotalwMCstat = new TGraphAsymmErrors(hD2eNormTotal);
 547   gD2eNormTotal->SetName("gD2eSys");
 548   gD2eNormTotal->SetTitle("electron spectra based on D measurement with systematic errors");
 549   gD2eNormTotalwMCstat->SetName("gD2eSysStat");
 550   gD2eNormTotalwMCstat->SetTitle("electron spectra based on D measurement with systematic and statistic errors");
 551
 552   Double_t errLow[4], errHigh[4], errDStat[4], errMCStat[4];
 553   Double_t errLowSum, errHighSum, errDStatSum, errMCStatSum;
 554   for(int i=0; i<hD2eNormTotal->GetNbinsX(); i++){
 555     errLowSum = 0.;
 556     errHighSum = 0.;
 557     errDStatSum = 0.;
 558     errMCStatSum = 0.;
 559
 560     for(int iDtype=0; iDtype<2; iDtype++){
 561       errLow[iDtype] = gD2eNormSum[iDtype]->GetErrorYlow(i); // error from D systematic error
 562       errLowSum += errLow[iDtype];
 563
 564       errHigh[iDtype]= gD2eNormSum[iDtype]->GetErrorYhigh(i); // error from D systematic error
 565       errHighSum += errHigh[iDtype];
 566
 567       errDStat[iDtype] = hD2eNormSum[iDtype]->GetBinError(i+1); //error from D statistical error
 568       errDStatSum += errDStat[iDtype];
 569
 570       errMCStat[iDtype] = hD2eNormSumStat[iDtype]->GetBinError(i+1); //error from MC statistical error
 571       errMCStatSum += errMCStat[iDtype];
 572     }
 573     if(i<hD2eNormTotal->FindBin(2.0)){
 574       errDStatSum = (hD0lowptcontrib->GetBinContent(i+1)+1.)*errDStatSum; //put conservative error based on the yield fraction of low pt contribution
 575     }
 576
 577     gD2eNormTotal->SetPointError(i,0,0,TMath::Sqrt(errLowSum*errLowSum + errDStatSum*errDStatSum),TMath::Sqrt(errHighSum*errHighSum + errDStatSum*errDStatSum));
 578     gD2eNormTotalwMCstat->SetPointError(i,0,0,TMath::Sqrt(errLowSum*errLowSum + errDStatSum*errDStatSum + errMCStatSum*errMCStatSum),TMath::Sqrt(errHighSum*errHighSum + errDStatSum*errDStatSum + errMCStatSum*errMCStatSum));
 579
 580   }
 581
 582   // drawing ----------------------------------------------------------
 583   cPtD2eSumwErr = new TCanvas("cPtD2eSumwErr","pT of e from D",0,0,1000,600);
 584   cPtD2eSumwErr->cd(1);
 585   gD2eNormTotalwMCstat->SetLineColor(2);
 586   gD2eNormTotalwMCstat->SetMarkerColor(2);
 587   gD2eNormTotalwMCstat->SetMarkerStyle(20);
 588   gD2eNormTotalwMCstat->Draw("ACP");
 589   gD2eNormTotalwMCstat->GetXaxis()->SetTitle("p_{t} [GeV/c]");
 590   gD2eNormTotalwMCstat->GetYaxis()->SetTitle("BRxd#sigma/dp_{t} | |#eta|<0.9  [pb/GeV/c]");
 591   gD2eNormTotal->SetMarkerStyle(20);
 592   gD2eNormTotal->SetLineStyle(2);
 593   gD2eNormTotal->GetYaxis()->SetTitle("BRxd#sigma/dp_{t} | |#eta|<0.9  [pb/GeV/c]");
 594   gD2eNormTotal->Draw("Psame");
 595   TLegend *legtotalerr = new TLegend(0.50,0.75,0.75,0.85);
 596   legtotalerr->AddEntry(gD2eNormTotal,"D#rightarrow e 2.0<mother p_{t}<50.0 GeV/c","lp");
 597   legtotalerr->AddEntry(gD2eNormTotalwMCstat,"(include MC stat error)","l");
 598   setLegendStyle(*legtotalerr,0);
 599   legtotalerr->Draw("same");
 600   gPad->SetLogy();
 601   // ------------------------------------------------------------------
 602
 603   // weighting method
 604   Double_t iM2PD[6];
 605   for(int iDtype=0; iDtype<2; iDtype++){
 606     for(int iptb=0; iptb<6; iptb++){
 607       if(iDtype==0) iM2PD[iptb]=iM2PD0[iptb];
 608       if(iDtype==1) iM2PD[iptb]=iM2PDp[iptb];
 609       hD2eNormWm[iDtype][iptb+3]->Scale(nevtnorm);
 610       hD2eNormWm[iDtype][iptb+3]->Scale(iM2PD[iptb]);  // measured/pythia for this pt bin
 611       hD2eNormWm[iDtype][iptb+3]->Scale(0.5);
 612       hD2eNormWm[iDtype][iptb+3]->Scale(Xsectrion4e);
 613     }
 614   }
 615
 616   // sum up to 12 GeV/c for D0 and D+ to compare 2 different methods
 617   TH1D *hD2eNormSumM1[2];
 618   TH1D *hD2eNormSumM2[2];
 619   for(int iDtype=0; iDtype<2; iDtype++){
 620     hname="hD2eNormSumM1" + iDtype;
 621     hD2eNormSumM1[iDtype]=(TH1D*)hD2eNorm[iDtype][3]->Clone(hname);
 622     hname="hD2eNormSumM2" + iDtype;
 623     hD2eNormSumM2[iDtype]=(TH1D*)hD2eNormWm[iDtype][3]->Clone(hname);
 624     for(int iptb=1; iptb<6; iptb++){
 625       hD2eNormSumM1[iDtype]->Add(hD2eNorm[iDtype][iptb+3]);
 626       hD2eNormSumM2[iDtype]->Add(hD2eNormWm[iDtype][iptb+3]);
 627     }
 628   }
 629   // drawing -------------------------------------------------------
 630   c2Methods = new TCanvas("c2Methods","e spectra from 2 different methods",0,0,1000,700);
 631   c2Methods->Divide(3,1);
 632   TH1D* hRatio2mtd[2];
 633   TLegend *legM[3];
 634   for(int i=0; i<2; i++){
 635     c2Methods->cd(i+1);
 636     legM[i]= new TLegend(0.25,0.70,0.75,0.85);
 637     hD2eNormSumM1[i]->SetMarkerStyle(20);
 638     hD2eNormSumM1[i]->Draw();
 639     hD2eNormSumM1[i]->GetYaxis()->SetTitle("BRxd#sigma/dp_{t} | |#eta|<0.9  [pb/GeV/c]");
 640     hD2eNormSumM2[i]->SetMarkerColor(2);
 641     hD2eNormSumM2[i]->SetMarkerStyle(20);
 642     hD2eNormSumM2[i]->Draw("same");
 643     hname="hRatio2mtdD" + i;
 644     hRatio2mtd[i]=(TH1D*)hD2eNormSumM2[i]->Clone(hname);
 645     hRatio2mtd[i]->Divide(hD2eNormSumM1[i]);
 646     gPad->SetLogy();
 647   }
 648   legM[2]= new TLegend(0.25,0.70,0.75,0.85);
 649   legM[0]->AddEntry(hD2eNormSumM1[0],"2.0<mother p_{t}<12.0 GeV/c","");
 650   legM[0]->AddEntry(hD2eNormSumM1[0],"D^{0}#rightarrow e(use of Xsection)","p");
 651   legM[0]->AddEntry(hD2eNormSumM2[0],"D^{0}#rightarrow e(use of weighting)","p");
 652   legM[1]->AddEntry(hD2eNormSumM1[1],"2.0<mother p_{t}<12.0 GeV/c","");
 653   legM[1]->AddEntry(hD2eNormSumM1[1],"D^{+}#rightarrow e(use of Xsection)","p");
 654   legM[1]->AddEntry(hD2eNormSumM2[1],"D^{+}#rightarrow e(use of weighting)","p");
 655
 656   legM[2]->AddEntry(hRatio2mtd[0],"2.0<mother p_{t}<12.0 GeV/c","");
 657   legM[2]->AddEntry(hRatio2mtd[0],"D^{0}#rightarrow e: (use of weighting)/(use of Xsection)","p");
 658   legM[2]->AddEntry(hRatio2mtd[1],"D^{+}#rightarrow e: (use of weighting)/(use of Xsection)","p");
 659   c2Methods->cd(1);
 660   legM[0]->Draw("same");
 661   c2Methods->cd(2);
 662   legM[1]->Draw("same");
 663
 664   setLegendStyle(*legM[0],0);
 665   setLegendStyle(*legM[1],0);
 666   setLegendStyle(*legM[2],0);
 667
 668   c2Methods->cd(3);
 669   hRatio2mtd[0]->Draw();
 670   hRatio2mtd[1]->SetLineColor(1);
 671   hRatio2mtd[1]->SetMarkerColor(1);
 672   hRatio2mtd[1]->Draw("same");
 673   legM[2]->Draw("same");
 674   //----------------------------------------------------------------
 675
 676
 677   // produce output file
 678   TFile *f = new TFile("Testoutput.root","recreate");
 679   f->cd();
 680   hD2eNormTotal->Write();
 681   gD2eNormTotal->Write();
 682   gD2eNormTotalwMCstat->Write();
 683   f->Close();
 684
 685
 686   // compare with FONLL ----------------------------------------------
 687   h4FONLL= new TH1F("h4FONLL","h4FONLL",60,0.25,30.25);
 688   ifstream stream1("./e_from_D-FONLL-abseta.lt.0.9.txt");
 689
 690   double b;
 691   if(!stream1)
 692     cout << "While opening a file an error is encountered" << endl;
 693
 694   int i=0;
 695   int ipt2=0;
 696   int k=0;
 697   Double_t val[8][60];
 698   Double_t exl[60];
 699   Double_t exh[60];
 700   Double_t eyl[60];
 701   Double_t eyh[60];
 702   while(!stream1.eof())
 703      {
 704          stream1 >> b;
 705
 706          i=k%8;
 707          ipt2=int(double(k)/8.);
 708          val[i][ipt2] = b;
 709          k++;
 710          if(ipt2==59 && i==7) break;
 711      }
 712
 713   for(int i=0; i<60; i++){
 714     h4FONLL->Fill(val[0][i],val[1][i]);
 715   }
 716
 717   for(int m=0; m<60; m++){
 718     exl[m]=0.0;
 719     exh[m]=0.0;
 720     eyl[m]=val[1][m]-val[2][m];
 721     eyh[m]=val[3][m]-val[1][m];
 722   }
 723
 724   TGraphAsymmErrors* gr=new  TGraphAsymmErrors(60,val[0],val[1],exl,exh,eyl,eyh);
 725   gr->SetMarkerStyle(20);
 726   gr->SetMarkerColor(1);
 727   gr->SetMarkerSize(0.5);
 728   gr->SetLineColor(3);
 729   gr->SetFillColor(5);
 730   gr->GetXaxis()->SetTitle("p_{t} [GeV/c]");
 731   gr->GetYaxis()->SetTitle("BRxd#sigma/dp_{t} | |#eta|<0.9  [pb/GeV/c]");
 732
 733   //gr->Draw("ZACPe3");
 734   //gr->Draw("ZCe3");
 735   //gr->Draw("Psame");
 736   cD2eFONLL = new TCanvas("cD2eFONLL","pT of e from D and FONLL",0,0,1000,600);
 737   cD2eFONLL->cd();
 738   gr->Draw("ACP");
 739   gD2eNormTotalwMCstat->Draw("psame");
 740   TLegend *legwFON = new TLegend(0.50,0.75,0.75,0.85);
 741   legwFON->AddEntry(gr,"D#rightarrow e FONLL","p");
 742   legwFON->AddEntry(gD2eNormTotalwMCstat,"D#rightarrow e 2.0<mother p_{t}<50.0 GeV/c","p");
 743   setLegendStyle(*legwFON,0);
 744   legwFON->Draw("same");
 745   gPad->SetLogy();
 746
 747 }
 748
 749 //--------------------------
 750 void setDataStyle(TH1D &h, Int_t lc, Int_t lw, Int_t ls){
 751
 752   h.SetLineColor(lc);
 753   h.SetLineWidth(lw);
 754   h.SetLineStyle(ls);
 755   h.SetMarkerColor(lc);
 756
 757 }
 758
 759 //--------------------------
 760 void setLegendStyle(TLegend &legend, Int_t bs){
 761
 762   legend.SetBorderSize(bs);
 763   legend.SetFillColor(0);
 764   legend.SetTextFont(132);
 765   legend.SetTextSize(0.04);
 766   legend.SetMargin(0.2);
 767
 768 }
 769
 770 //--------------------------
 771 void setPadStyle(Int_t lvl, TPad *pad){
 772
 773   pad->SetLogy();
 774   if(lvl>0) gPad->SetGridy();
 775   if(lvl>1) gPad->SetGridx();
 776
 777 }
 778
 779 //--------------------------
 780 void setGeneralStyle(){
 781
 782   gStyle->SetPalette(1);
 783   gStyle->SetFrameBorderMode(0);
 784   gStyle->SetFrameFillColor(0);
 785   gStyle->SetPadBorderSize(0);
 786   gStyle->SetPadBorderMode(0);
 787   gStyle->SetCanvasColor(0);
 788   gStyle->SetCanvasBorderSize(10);
 789   gStyle->SetOptStat(0);
 790   gStyle->SetOptTitle(0);
 791   gStyle->SetTitleFillColor(10);
 792   gStyle->SetTitleBorderSize(0);
 793   gStyle->SetTitleSize(0.04,"X");
 794   gStyle->SetTitleSize(0.04,"Y");
 795   gStyle->SetTitleFont(132,"X");
 796   gStyle->SetTitleFont(132,"Y");
 797   gStyle->SetTitleXOffset(0.9);
 798   gStyle->SetTitleYOffset(0.9);
 799   gStyle->SetLabelFont(132,"X");
 800   gStyle->SetLabelFont(132,"Y");
 801   gStyle->SetLabelSize(0.04,"X");
 802   gStyle->SetLabelSize(0.04,"Y");
 803   gStyle->SetTitleSize(0.05,"X");
 804   gStyle->SetTitleSize(0.05,"Y");
 805   gStyle->SetLineWidth(2);
 806
 807 }
 808
 809 void CorrectFromTheWidth(TH1D *h1) const {
 810   //
 811   // Correct from the width of the bins --> dN/dp_{T} (GeV/c)^{-1}
 812   //
 813
 814   TAxis *axis = h1->GetXaxis();
 815   Int_t nbinX = h1->GetNbinsX();
 816
 817   for(Int_t i = 1; i <= nbinX; i++) {
 818
 819     Double_t width = axis->GetBinWidth(i);
 820     Double_t content = h1->GetBinContent(i);
 821     Double_t error = h1->GetBinError(i);
 822     h1->SetBinContent(i,content/width);
 823     h1->SetBinError(i,error/width);
 824     Double_t pt = h1->GetBinCenter(i);
 825   }
 826
 827 }