PWGHF/vertexingHF/macros/CalculateAveragePt.C

   1 // MACRO TO COMPUTE THE <pt> from data for D mesons
   2 // works for D0, Ds and D+; D* case should be tested since the side bands definition is less clear
   3 //
   4 // Requires
   5 // a 2D (pt,mass) histogram, with fine pt bins in pt (finer than those used for the spectrum, e.g. 100 or 200 MeV/c)
   6 // a histogram of background yield vs pt (with the pt bins in which the <pt> has to be computed)
   7 // a histogram of the invariant mass sigma (with the pt bins in which the <pt> has to be computed)
   8 // optionally a histogram of signal yield (not fully needed) (with the pt bins in which the <pt> has to be computed)
   9 // optionally a histogram of the mean of the signal peak obtained from the fit (with the pt bins in which the <pt> has to be computed)
  10 // The last 4 histograms are among the output of the standard macro used for fitting the mass distributions
  11 // optionally (but strongly recommended) an histograms with prompt D reco efficiencies (in finer bins than those in which the <pt> is wanted, otherwise they will have no effect).
  12 //            The binning of the efficiency histo should match that of the pt axis of the 2D (pt,mass) plots x rebin (where rebin is used to rebin the data histogram with the method Rebin)
  13 //            e.g. if you have data in 200 MeV/c and efficiency in 500 MeV/c bins than the minimum binning for which data and efficiencis match is 1 GeV/c
  14 //                    -> you have to use rebin=5 (or multiples of 5) to have the macro working (it stops otherwise)
  15 // Options:
  16 // rebin (see above)
  17 // usefitforsubtraction: with this switched on the pt shape of the backgrond is fitted with an exponential and the fit function is used for the background subtraction in the signal region
  18 //            (useful if there are stat fluctuations in the background pt shape distribution, but take care and check the fit quality in the control plots)
  19 // correct for efficiency: if switched on, the pt distribution of the signal (i.e. that in the signal peak after back subtraction) is weighted by the efficiency
  20 //                         it accounts for the pt dependence of the efficiency -> compulsory if a realistic (corrected) <pt> has to be computed
  21 //
  22 // useParGenAccOverLimAcc: uses a parametrization for correcting the efficiencies by GenAcc/LimAcc. It is useful since the GenAcc/LimAcc factor might
  23 //                         not be available in very fine bins, especially at high pt
  24 //
  25 //
  26 // There are standard methods hard-coded for D0 (DoStandardForD0), D+ and Ds as an example (of course you have to change the paths of the files and the name of the histos!!).
  27 //     To run them (e.g. for D0):
  28 //    -1 (my advice): copy the macro in your working directory and change the name of the paths and histo there (it helps to retrieve the results in the future).
  29 //                    Or copy the ingredient files in the working directory,
  30 //     0- set the path of the files
  31 //     1- .L CalculateAveragePt.C++  (compile the macro, not really needed)
  32 //     2- DoStandardForD0(5,kFALSE,kTRUE,kTRUE,0,3)
  33 //    ... it should be easy :)
  34
  35 //
  36 //
  37 //  For problems write to: andrea.rossi@cern.ch
  38 //
  39 ////////////////////////////////////////////////////////////////////////////////////
  40
  41 #include <TH2F.h>
  42 #include <TH1F.h>
  43 #include <TH2D.h>
  44 #include <TH1D.h>
  45 #include <TCanvas.h>
  46 #include <TF1.h>
  47 #include <TFile.h>
  48 #include <TDirectory.h>
  49 #include <TMath.h>
  50 #include <TString.h>
  51 #include <TGraphAsymmErrors.h>
  52 #include <TAxis.h>
  53 #include <TLegend.h>
  54
  55 TH2F *hPtInvMass=0x0;
  56 Int_t nptbins=0;
  57 Double_t *ptbinlimits=0x0;
  58 Double_t *rawsignal=0x0,*rawback=0x0,*sigma=0x0,*meansignal=0x0;
  59 Double_t nsigmaSignal=3.,nsigmaSBstart=5.;
  60 Double_t mesonMass=1.8645;
  61 Int_t nbinsx;//=hPtInvMass->GetNbinsX();
  62 Int_t nbinsy;//=hPtInvMass->GetNbinsY();
  63 Double_t binwidthpt;//=hPtInvMass->GetYaxis()->GetBinWidth(1);// ASSUME ALL BINS ARE WITH THE SAME WIDTH
  64 Double_t binwidthInvMss;//=hPtInvMass->GetXaxis()->GetBinWidth(1);// ASSUME ALL BINS ARE WITH THE SAME WIDTH
  65 Double_t ptmin;//=hPtInvMass->GetYaxis()->GetBinLowEdge(1);
  66 Double_t ptmax;//=hPtInvMass->GetYaxis()->GetBinUpEdge(nbinsy);
  67 Bool_t useFitForSubtraction=kFALSE;
  68 Bool_t useParGenAccOverLimAcc=kFALSE;
  69 Bool_t subtractSB=kTRUE;
  70 Bool_t corrForEff=kFALSE;
  71 void SetSubtractSB(Bool_t subtract){subtractSB=subtract;}
  72 void SetCorrForEff(Bool_t correff){corrForEff=correff;}
  73 void SetUseFitForSubtraction(Bool_t useFit){useFitForSubtraction=useFit;}
  74 void SetNsigmaForSignal(Double_t nsigm){nsigmaSignal=nsigm;}
  75 void SetNsigmaStartSB(Double_t nsigm){nsigmaSBstart=nsigm;
  76 }
  77 void CalculateAveragePt(Int_t rebin=1,Int_t firstbin=0,Int_t lastbin=0);
  78 TF1* ParametricGenAccOverLimAccCorr(){// parametric form of GenAcc/LimAcc correction for D0->Kpi
  79   TF1 *fGenAccOverLimAcc=new TF1("fat7","[2]+([0]*TMath::ATan([1]*x+[3])+[4]*TMath::Log(1.+x/[5]))",0,24);
  80   fGenAccOverLimAcc->SetParameters(2.76153e-01,6.69658e-01,8.45865e-01,-1.87709e+00,2.89845e+03,1.39651e+05);
  81   return fGenAccOverLimAcc;
  82 }
  83
  84 Int_t standrebin[8]={2,2,2,2,2,2,4,5};
  85 TCanvas **cPtDistrNoSubtr;
  86 TF1 **fitfunc;
  87 TH1D *hAvRawYieldSpectrum,*hSignal,*hMeanSignal,*hEfficNum,*hEfficDenum,*hBack,*hSigma;
  88 void SetHistosEfficiency(TH1D *hNum,TH1D *hDenum){
  89   hEfficNum=(TH1D*)hNum;
  90   hEfficDenum=(TH1D*)hDenum;
  91
  92 }
  93 TGraphAsymmErrors *grAvRawYieldSpectrum;
  94 TGraphAsymmErrors *grAvPtVSPtmean;
  95
  96 TGraphAsymmErrors *grBackAvPtVSPtmean;
  97
  98 void SetPtBinLimits(const Int_t npt,Double_t *ptbinlim){
  99   if(nptbins>0){
 100     delete ptbinlimits;ptbinlimits=0x0;
 101   }
 102   nptbins=npt;
 103   ptbinlimits=new Double_t[nptbins+1];
 104   for(Int_t bin=0;bin<=nptbins;bin++){
 105     ptbinlimits[bin]=ptbinlim[bin];
 106   }
 107   cPtDistrNoSubtr=new TCanvas*[nptbins];
 108   fitfunc=new TF1*[nptbins];
 109 }
 110
 111 void SetPtBinLimits(TH1 *histo){
 112
 113   if(nptbins>0){
 114     delete ptbinlimits;ptbinlimits=0x0;
 115   }
 116   TAxis *xax=histo->GetXaxis();
 117   nptbins=xax->GetNbins();
 118   ptbinlimits=new Double_t[nptbins+1];
 119   for(Int_t bin=0;bin<=nptbins;bin++){
 120     ptbinlimits[bin]=xax->GetBinLowEdge(bin+1);
 121   }
 122    cPtDistrNoSubtr=new TCanvas*[nptbins];
 123    fitfunc=new TF1*[nptbins];
 124 }
 125
 126 void SetHistRawSignal(TH1D *hS){hSignal=hS;
 127  rawsignal=&(hSignal->GetArray())[1];
 128  return;}
 129
 130 void SetHistMean(TH1D *hM){hMeanSignal=hM;
 131  meansignal=&(hMeanSignal->GetArray())[1];
 132  return;
 133 }
 134
 135 void SetHistRawBack(TH1D *hB){hBack=hB;
 136  rawback=&(hBack->GetArray())[1];
 137  hAvRawYieldSpectrum=new TH1D(*hBack);
 138  hAvRawYieldSpectrum->SetName("hAvRawYieldDistrTotPeak");
 139  hAvRawYieldSpectrum->SetTitle("Average Raw Yield in bins after subtraction from bin counting");
 140  hAvRawYieldSpectrum->Reset(0);
 141  hAvRawYieldSpectrum->SetLineColor(kRed);
 142  return;
 143 }
 144 void SetHistSigma(TH1D *hSig){hSigma=hSig;
 145  sigma=&(hSigma->GetArray())[1];
 146  return;
 147 }
 148
 149 void SetHistoMassPt(TH2F *h2){
 150   hPtInvMass=new TH2F(*h2);
 151   nbinsx=hPtInvMass->GetNbinsX();
 152   nbinsy=hPtInvMass->GetNbinsY();
 153   binwidthpt=hPtInvMass->GetYaxis()->GetBinWidth(1);// ASSUME ALL BINS ARE WITH THE SAME WIDTH
 154   binwidthInvMss=hPtInvMass->GetXaxis()->GetBinWidth(1);// ASSUME ALL BINS ARE WITH THE SAME WIDTH
 155   ptmin=hPtInvMass->GetYaxis()->GetBinLowEdge(1);
 156   ptmax=hPtInvMass->GetYaxis()->GetBinUpEdge(nbinsy);
 157
 158 }
 159
 160 TH1D* CheckBinningAndMerge(TH1D *hA,TH1D *hB,Double_t precision=0.001,Double_t minX=-9999.,Double_t maxX=-9999.){// Checks that histo hB binning is finer and fits hA binning
 161   // and return an Histo with hB bins integrated to match than hA bins
 162   // TO BE ADDED A CHECK ON THE MIN AND MAX VALUES OF THE AXES IN THE 2 PLOTS
 163   // NO PARTICULAR TREATMENT OF THE ERROR IN THE INTEGRATION (SHOULD BE OK BUT FOR EFFIENCIES)
 164
 165   TH1D *hWork=new TH1D(*hA);
 166   hWork->SetName(Form("%sNewBins",hB->GetName()));
 167   hWork->SetTitle(hB->GetTitle());
 168
 169   Double_t lebinA,uebinA;
 170   Int_t binBleA=0,binBueA=0;
 171   Int_t binBleALast=0,binBueALast=0;
 172   Int_t maxbinA=hA->GetNbinsX();
 173   Int_t minbinA=1;
 174   if(maxX!=-9999.){
 175     maxbinA=hA->FindBin(maxX);
 176   }
 177   if(minX!=-9999.){
 178     minbinA=hA->FindBin(minX);
 179   }
 180
 181   for(Int_t binA=minbinA;binA<=maxbinA;binA++){
 182     uebinA=hA->GetXaxis()->GetBinUpEdge(binA);
 183     lebinA=hA->GetXaxis()->GetBinLowEdge(binA);
 184     // Check that Low Edge of binA is close to the low or upper edge of a bin in hB
 185     // within precition
 186     binBleA=hB->FindBin(lebinA);
 187
 188     if(TMath::Abs(hB->GetXaxis()->GetBinLowEdge(binBleA)-lebinA)>precision&&TMath::Abs(hB->GetXaxis()->GetBinUpEdge(binBleA)-lebinA)>precision){
 189       printf("The bins of the efficiency histo do not match those of the pt spectrum histo (edges 1):\n Cannot correct for efficiensies \n");
 190       delete hWork;
 191       return 0x0;
 192     }
 193     if(TMath::Abs(hB->GetXaxis()->GetBinLowEdge(binBleA)-lebinA)>TMath::Abs(hB->GetXaxis()->GetBinUpEdge(binBleA)-lebinA))binBleA++;
 194     if(binBleA<=binBleALast){
 195       printf("Problems with rebinning Hist B, non compatibility with different bin of Hist A (1)\n");
 196       delete hWork;
 197       return 0x0;
 198     }
 199     if(hB->GetBinWidth(binBleA)-hA->GetBinWidth(binA)>precision){
 200       printf("The bins of the efficiency histo do not match those of the pt spectrum histo (bin width 1):\n Cannot correct for efficiensies \n");
 201       delete hWork;
 202       return 0x0;
 203     }
 204     // Check that Upper Edge of binA is close to the low or upper edge of a bin in hB
 205     // within precition
 206     binBueA=hB->FindBin(uebinA);
 207     if(TMath::Abs(hB->GetXaxis()->GetBinLowEdge(binBueA)-uebinA)>precision&&TMath::Abs(hB->GetXaxis()->GetBinUpEdge(binBueA)-uebinA)>precision){
 208       printf("The bins of the efficiency histo do not match those of the pt spectrum histo (edges 2) (bin %d Vs %d):\n Cannot correct for efficiensies \n",binA,binBueA);
 209       printf("%f Vs %f or %f \n",uebinA,hB->GetXaxis()->GetBinLowEdge(binBueA),hB->GetXaxis()->GetBinUpEdge(binBueA));
 210       delete hWork;
 211       return 0x0;
 212     }
 213     if(TMath::Abs(hB->GetXaxis()->GetBinLowEdge(binBueA)-uebinA)<TMath::Abs(hB->GetXaxis()->GetBinUpEdge(binBueA)-uebinA))binBueA--;
 214     if(binBueA<=binBueALast){
 215       printf("Problems with rebinning Hist B, non compatibility with different bin of Hist A (2) \n");
 216       delete hWork;
 217       return 0x0;
 218     }
 219     if(hB->GetBinWidth(binBueA)-hA->GetBinWidth(binA)>precision){
 220       printf("The bins of the efficiency histo do not match those of the pt spectrum histo (bin width 2):\n Cannot correct for efficiensies \n");
 221       delete hWork;
 222       return 0x0;
 223     }
 224     if(binBleA>binBueA){
 225       printf("Mathcing failure (probably a bug in the code \n");
 226       delete hWork;
 227       return 0x0;
 228     }
 229
 230     hWork->SetBinContent(binA,hB->Integral(binBleA,binBueA));
 231   }
 232
 233   return hWork;
 234
 235 }
 236
 237 Bool_t CorrectForEfficiency(TH1D *hPtHisto){
 238   TH1D *hWorkNum=0x0,*hWorkDenum=0x0;
 239   //  hWorkNum=CheckBinningAndMerge(hPtHisto,hEfficNum,0.01,1.,23.9);
 240   //  printf("ptmin and pt max: %f, %f \n",ptbinlimits[0],ptbinlimits[nptbins]);
 241   hWorkNum=CheckBinningAndMerge(hPtHisto,hEfficNum,0.01,ptbinlimits[0]*1.0001,ptbinlimits[nptbins]*0.9999);
 242   hWorkNum->Sumw2();
 243
 244   if(!hWorkNum)return kFALSE;
 245   hWorkDenum=CheckBinningAndMerge(hPtHisto,hEfficDenum,0.01,ptbinlimits[0]*1.0001,ptbinlimits[nptbins]*0.9999);
 246   //  hWorkDenum=CheckBinningAndMerge(hPtHisto,hEfficDenum,0.01,1.,23.9);;
 247   hWorkDenum->Sumw2();
 248   if(!hWorkDenum){
 249     delete hWorkNum;
 250     return kFALSE;
 251   }
 252   TH1D *hWorkNumCl=(TH1D*)hWorkDenum->Clone("hWorkNumCl");
 253   hWorkNumCl->Divide(hWorkNum,hWorkDenum,1,1,"B");
 254   if(useParGenAccOverLimAcc){
 255     TF1 *fGenAccoOverLimAcc=ParametricGenAccOverLimAccCorr();
 256     for(Int_t j=1;j<=hWorkNumCl->GetNbinsX();j++){
 257       printf(" EFFF for %f: %f * %f \n",hWorkNumCl->GetBinCenter(j), hWorkNumCl->GetBinContent(j),fGenAccoOverLimAcc->Eval(hWorkNumCl->GetBinCenter(j)));
 258       hWorkNumCl->SetBinContent(j,hWorkNumCl->GetBinContent(j)*fGenAccoOverLimAcc->Eval(hWorkNumCl->GetBinCenter(j)));
 259       hWorkNumCl->SetBinError(j,hWorkNumCl->GetBinError(j)*fGenAccoOverLimAcc->Eval(hWorkNumCl->GetBinCenter(j)));
 260
 261     }
 262   }
 263   hPtHisto->Divide(hWorkNumCl);
 264
 265   delete hWorkNum;
 266   delete hWorkNumCl;
 267   delete hWorkDenum;
 268   return kTRUE;
 269 }
 270
 271 TH1D *HistoPtShapeFromData(Int_t ptbin,Int_t rebin=1.){
 272   if(!hPtInvMass){
 273     printf("No histo set \n");
 274     return 0x0;
 275   }
 276   Double_t mean=mesonMass;
 277   if(meansignal!=0x0){
 278     mean=meansignal[ptbin];
 279     printf("USING MEAN FROM MASS FIT: %f\n",mean);
 280   }
 281   TString namehist="hPtDistrSB_Bin",funcname;
 282   namehist+=ptbin;
 283
 284   TString cname="cPtDistrNoSubBin";
 285   cname+=ptbin;
 286   cPtDistrNoSubtr[ptbin]=new TCanvas(cname.Data(),"CPtDistrNoSub",700,700);
 287
 288   printf("After Setting Canva \n");
 289   Int_t binleftpt=-1,binrightpt=-1;
 290   TH1D *hPtDistrSB=new TH1D(namehist.Data(),"Side band Pt distribution",nbinsy,ptmin,ptmax);
 291   namehist="hPtDistrTotPeak_Bin";
 292   namehist+=ptbin;
 293   printf("After Setting SB hist \n");
 294   TH1D *hPtDistrTotPeak=new TH1D(namehist.Data(),"Total Pt distribution in signal mass reagion",nbinsy,ptmin,ptmax);
 295   for(Int_t jbi=1;jbi<hPtDistrTotPeak->GetNbinsX();jbi++){
 296     hPtDistrTotPeak->SetBinContent(jbi,0);
 297     hPtDistrTotPeak->SetBinError(jbi,0);
 298     hPtDistrSB->SetBinContent(jbi,0);
 299     hPtDistrSB->SetBinError(jbi,0);
 300   }
 301
 302
 303
 304   printf("Before loop on xbins \n");
 305   for(Int_t j=1;j<nbinsx;j++){
 306     if((hPtInvMass->GetXaxis()->GetBinLowEdge(j)<=mean-nsigmaSBstart*sigma[ptbin])||(hPtInvMass->GetXaxis()->GetBinUpEdge(j)>=mean+nsigmaSBstart*sigma[ptbin])){
 307       for(Int_t k=1;k<nbinsy;k++){
 308         if(hPtInvMass->GetYaxis()->GetBinLowEdge(k)>=ptbinlimits[ptbin]&&hPtInvMass->GetYaxis()->GetBinUpEdge(k)<=ptbinlimits[ptbin+1]){
 309           hPtDistrSB->SetBinContent(k,hPtDistrSB->GetBinContent(k)+hPtInvMass->GetBinContent(j,k));
 310           hPtDistrSB->SetBinError(k,TMath::Sqrt(hPtDistrSB->GetBinError(k)*hPtDistrSB->GetBinError(k)+hPtInvMass->GetBinError(j,k)*hPtInvMass->GetBinError(j,k)));
 311         }
 312       }
 313     }
 314     else if((hPtInvMass->GetXaxis()->GetBinLowEdge(j)>=mean-nsigmaSignal*sigma[ptbin])&&(hPtInvMass->GetXaxis()->GetBinUpEdge(j)<=mean+nsigmaSignal*sigma[ptbin])){
 315       for(Int_t k=1;k<nbinsy;k++){
 316         if(hPtInvMass->GetYaxis()->GetBinLowEdge(k)>=ptbinlimits[ptbin]&&hPtInvMass->GetYaxis()->GetBinUpEdge(k)<=ptbinlimits[ptbin+1]){
 317           if(binleftpt==-1)binleftpt=k;
 318           hPtDistrTotPeak->SetBinContent(k,hPtDistrTotPeak->GetBinContent(k)+hPtInvMass->GetBinContent(j,k));
 319           hPtDistrTotPeak->SetBinError(k,TMath::Sqrt(hPtDistrTotPeak->GetBinError(k)*hPtDistrTotPeak->GetBinError(k)+hPtInvMass->GetBinError(j,k)*hPtInvMass->GetBinError(j,k)));
 320           binrightpt=k;
 321         }
 322       }
 323     }
 324   }
 325   hPtDistrSB->Sumw2();
 326   hPtDistrTotPeak->Sumw2();
 327   if(rebin!=1){//Rebin and recalculate needed info
 328     hPtDistrSB->Rebin(rebin);
 329     hPtDistrTotPeak->Rebin(rebin);
 330     binleftpt=-1;
 331     binleftpt=-1;
 332     for(Int_t k=1;k<=hPtDistrTotPeak->GetNbinsX();k++){
 333       if(hPtDistrTotPeak->GetBinLowEdge(k)>=ptbinlimits[ptbin]&&hPtDistrTotPeak->GetBinLowEdge(k+1)<=ptbinlimits[ptbin+1]){
 334         if(binleftpt==-1)binleftpt=k;
 335         binrightpt=k;
 336       }
 337     }
 338     binwidthpt=hPtDistrTotPeak->GetBinWidth(1);
 339   }
 340
 341   Double_t avptbin=0.,errA=0.,errB=0.,errC=0.,errisq=0.;
 342   // CALCULATE BACKGROUND AVERAGE PT
 343   for(Int_t bbin=1;bbin<=hPtDistrSB->GetNbinsX();bbin++){
 344     avptbin+=hPtDistrSB->GetBinCenter(bbin)*hPtDistrSB->GetBinContent(bbin)*hPtDistrSB->GetBinWidth(bbin);
 345     errisq=hPtDistrSB->GetBinError(bbin)*hPtDistrSB->GetBinError(bbin);
 346     errA+=hPtDistrSB->GetBinCenter(bbin)*hPtDistrSB->GetBinCenter(bbin)*errisq;
 347     errB+=hPtDistrSB->GetBinCenter(bbin)*errisq;
 348     errC+=errisq;
 349     // printf("ERROR VALUES: %f, %f \n",errA,errB);
 350   }
 351   avptbin/=hPtDistrSB->Integral("width");
 352   errisq=TMath::Sqrt(errA-2.*avptbin*errB+avptbin*avptbin*errC)/hPtDistrSB->Integral("width");
 353   if(!grBackAvPtVSPtmean){
 354     grBackAvPtVSPtmean=new TGraphAsymmErrors();
 355     grBackAvPtVSPtmean->SetName("grBackAvPtVSPtmean");
 356   }
 357   Int_t grbin=grBackAvPtVSPtmean->GetN();
 358   grBackAvPtVSPtmean->SetPoint(grbin,0.5*(ptbinlimits[ptbin+1]+ptbinlimits[ptbin]),avptbin);
 359   grBackAvPtVSPtmean->SetPointError(grbin,0.5*(ptbinlimits[ptbin+1]-ptbinlimits[ptbin]),0.5*(ptbinlimits[ptbin+1]-ptbinlimits[ptbin]),errisq,errisq);
 360
 361   TH1D *hPtSignalFromSubtraction=new TH1D(*hPtDistrTotPeak);
 362   namehist="hPtSignalFromSubtraction_Bin";
 363   namehist+=ptbin;
 364   hPtSignalFromSubtraction->SetName(namehist.Data());
 365   hPtSignalFromSubtraction->SetTitle("Signal Pt distr from subtraction");
 366
 367
 368   if(subtractSB){
 369     if(!useFitForSubtraction)hPtSignalFromSubtraction->Add(hPtDistrSB,-rawback[ptbin]*(nsigmaSignal/3.)/hPtDistrSB->Integral());// ROUGH LINEAR SCALING OF BACKGOUND!!
 370     else{
 371       funcname="expofit";
 372       funcname+=ptbin;
 373       fitfunc[ptbin]=new TF1(funcname.Data(),"expo",ptbinlimits[ptbin],ptbinlimits[ptbin+1]);
 374       TCanvas *cTempt=new TCanvas();
 375       cTempt->cd();
 376       hPtDistrSB->Fit(fitfunc[ptbin],"RLEV","",ptbinlimits[ptbin],ptbinlimits[ptbin+1]);
 377       delete cTempt;
 378       for(Int_t jb=binleftpt;jb<=binrightpt;jb++){
 379         printf("Err before sub: %f \n",hPtSignalFromSubtraction->GetBinError(jb));
 380         hPtSignalFromSubtraction->SetBinContent(jb,hPtSignalFromSubtraction->GetBinContent(jb)-fitfunc[ptbin]->Eval(hPtSignalFromSubtraction->GetBinCenter(jb))/fitfunc[ptbin]->Integral(ptbinlimits[ptbin],ptbinlimits[ptbin+1])*binwidthpt*rawback[ptbin]*(nsigmaSignal/3.));// ROUGH LINEAR SCALING OF BACKGOUND!!
 381         printf("Err after sub: %f \n\n",hPtSignalFromSubtraction->GetBinError(jb));
 382       }
 383     }
 384   }
 385
 386
 387   cPtDistrNoSubtr[ptbin]->Divide(1,2);
 388   cPtDistrNoSubtr[ptbin]->cd(1);
 389   hPtDistrSB->Draw();
 390   cPtDistrNoSubtr[ptbin]->cd(2);
 391   if(corrForEff){
 392     if(!CorrectForEfficiency(hPtSignalFromSubtraction)){
 393       printf("Correction for efficieny failed \n");
 394     }
 395   }
 396   hPtDistrTotPeak->Draw();
 397   hPtSignalFromSubtraction->Draw("sames");
 398
 399
 400   return hPtSignalFromSubtraction;
 401
 402 }
 403
 404 void CalculateAveragePt(TH2* hMassPt,TH1D *hB,TH1D *hSigm,TH1D *hEffNum=0x0,TH1D *hEffDenum=0x0,TH1D *hS=0x0,TH1D *hMean=0x0,Int_t rebin=1,Int_t firstbin=3,Int_t lastbin=8){
 405
 406   SetPtBinLimits(hB);
 407   SetHistoMassPt((TH2F*)hMassPt);
 408   SetHistRawBack(hB);
 409   SetHistSigma(hSigm);
 410   if(((!hEffNum)&&(hEffDenum))||((hEffNum)&&(!hEffDenum))){
 411     printf("Two histos are needed for the efficiency: numerator and denumerator (for rebinning) \n");
 412   }
 413   else if(hEffNum){
 414     SetHistosEfficiency(hEffNum,hEffDenum);
 415     //    SetHistEffNum(hEffNum);
 416     // SetHistEffDenum(hEffDenum);
 417   }
 418   if(hS)SetHistRawSignal(hS);
 419   if(hMean)SetHistMean(hMean);
 420   printf("Histos set \n");
 421   CalculateAveragePt(rebin,firstbin,lastbin);
 422 }
 423
 424 void CalculateAveragePt(Int_t rebin,Int_t firstbin,Int_t lastbin){
 425
 426   TH1D *hPtSpectra;
 427   if(rebin>0){
 428     hPtSpectra=new TH1D("hPtDistrTotPeak","Total Pt distribution in signal mass reagion",nbinsy,ptmin,ptmax);
 429     if(rebin!=1)hPtSpectra->Rebin(rebin);
 430   }
 431   else if(rebin==-1){
 432     printf("Standard rebinning not implemented yet  \n");return;
 433   }
 434   else {
 435     printf("Wrong rebin numnber \n");return;
 436   }
 437
 438   grAvRawYieldSpectrum=new TGraphAsymmErrors();
 439   grAvPtVSPtmean=new TGraphAsymmErrors();
 440   grBackAvPtVSPtmean=new TGraphAsymmErrors();
 441   grBackAvPtVSPtmean->SetName("grBackAvPtVSPtmean");
 442
 443
 444   Double_t avptbin=0.,errA,errB,errC,errisq;
 445   for(Int_t bin=firstbin;bin<=lastbin;bin++){
 446     //    printf("Before Getting Subtracted Histo \n");
 447     avptbin=0.;
 448     errA=0.;
 449     errB=0.;
 450     errC=0.;
 451     errisq=0.;
 452
 453     TH1D *h=HistoPtShapeFromData(bin,rebin);
 454     hPtSpectra->Add(h);
 455     hAvRawYieldSpectrum->SetBinContent(bin+1,h->Integral("width")/hAvRawYieldSpectrum->GetBinWidth(bin+1));
 456     for(Int_t bbin=1;bbin<=h->GetNbinsX();bbin++){
 457       avptbin+=h->GetBinCenter(bbin)*h->GetBinContent(bbin)*h->GetBinWidth(bbin);
 458       errisq=h->GetBinError(bbin)*h->GetBinError(bbin);
 459       errA+=h->GetBinCenter(bbin)*h->GetBinCenter(bbin)*errisq;
 460       errB+=h->GetBinCenter(bbin)*errisq;
 461       errC+=errisq;
 462       // printf("ERROR VALUES: %f, %f \n",errA,errB);
 463     }
 464     avptbin/=h->Integral("width");
 465     errisq=TMath::Sqrt(errA-2.*avptbin*errB+avptbin*avptbin*errC)/h->Integral("width");
 466
 467     grAvRawYieldSpectrum->SetPoint(bin-firstbin,avptbin,h->Integral("width")/hAvRawYieldSpectrum->GetBinWidth(bin+1));
 468     grAvRawYieldSpectrum->SetPointError(bin-firstbin,avptbin-ptbinlimits[bin],ptbinlimits[bin+1]-avptbin,errisq/2.,errisq/2.);
 469
 470     grAvPtVSPtmean->SetPoint(bin-firstbin,0.5*(ptbinlimits[bin+1]+ptbinlimits[bin]),avptbin);
 471     grAvPtVSPtmean->SetPointError(bin-firstbin,0.5*(ptbinlimits[bin+1]-ptbinlimits[bin]),0.5*(ptbinlimits[bin+1]-ptbinlimits[bin]),errisq,errisq);
 472
 473
 474     hSignal->SetBinContent(bin+1,hSignal->GetBinContent(bin+1)*hPtSpectra->GetBinWidth(1)/hSignal->GetBinWidth(bin+1));
 475     hSignal->SetBinError(bin+1,hSignal->GetBinError(bin+1)*hPtSpectra->GetBinWidth(1)/hSignal->GetBinWidth(bin+1));
 476     //    hPtSpectra->Add(HistoPtShapeFromData(bin,rebin),1./efficiency[bin]/nev*factorTheoryCompare*primaryCharmCorrect[bin]/binwidthpt);//(ptbinlimits[bin+1]-ptbinlimits[bin])
 477     cPtDistrNoSubtr[bin]->Draw();
 478   }
 479
 480    TString nameout="ptCorrection";
 481    TString namegrRawAvPt="grAvRawVsAvPt",namegrAvPtPtmean="grAvPtMeanPt";
 482    if(useFitForSubtraction){
 483      nameout.Append("FitSB");
 484      namegrRawAvPt.Append("FitSB");
 485      namegrAvPtPtmean.Append("FitSB");
 486    }
 487    if(!subtractSB){
 488      namegrRawAvPt.Append("NoSBsub.root");
 489      namegrAvPtPtmean.Append("NoSBsub.root");
 490      nameout.Append("NoSBsubtract.root");
 491    }
 492    else nameout.Append(".root");
 493
 494
 495    grAvRawYieldSpectrum->SetName(namegrRawAvPt.Data());
 496    grAvRawYieldSpectrum->SetMarkerStyle(21);
 497    grAvRawYieldSpectrum->SetMarkerSize(1.2);
 498    grAvRawYieldSpectrum->SetMarkerColor(kBlue);
 499    grAvRawYieldSpectrum->SetLineColor(kBlue);
 500
 501    grAvPtVSPtmean->SetName(namegrAvPtPtmean.Data());
 502    grAvPtVSPtmean->SetMarkerStyle(21);
 503    grAvPtVSPtmean->SetMarkerSize(1.2);
 504    grAvPtVSPtmean->SetMarkerColor(kBlue);
 505    grAvPtVSPtmean->SetLineColor(kBlue);
 506
 507
 508
 509    //  TFile *ftheory=TFile::Open("/Users/administrator/ALICE/CHARM/ppData_2010/PtSelInBin/predictions_D0.root");
 510    // TH1D *histoThCompareCentral=(TH1D*)ftheory->Get("hD0Kpipred_central");
 511    // histoThCompareCentral->Rebin(5);
 512
 513
 514
 515
 516   TCanvas *cPtSpectra=new TCanvas("cPtSpectra","cPtSpectra",700,700);
 517   cPtSpectra->cd();
 518   hPtSpectra->Draw();
 519   //hAvRawYieldSpectrum->Draw("same");
 520   hSignal->SetLineColor(kGreen);
 521   hSignal->Draw("same");
 522   grAvRawYieldSpectrum->Draw("p");
 523
 524   TCanvas *cAvPtVSPtmean=new TCanvas("cAvPtVSPtmean","cAvPtVSPtmean",700,700);
 525   cAvPtVSPtmean->cd();
 526   grAvPtVSPtmean->Draw("ap");
 527   grAvPtVSPtmean->GetYaxis()->SetTitle("<p_{t}> (GeV/c)");
 528   grAvPtVSPtmean->GetXaxis()->SetTitle("bin centre p_{t} (GeV/c)");
 529
 530   grBackAvPtVSPtmean->SetMarkerColor(kRed);
 531   grBackAvPtVSPtmean->SetMarkerStyle(24);
 532   grBackAvPtVSPtmean->SetLineColor(kRed);
 533   grBackAvPtVSPtmean->GetYaxis()->SetTitle("<p_{t}> (GeV/c)");
 534   grBackAvPtVSPtmean->GetXaxis()->SetTitle("bin centre p_{t} (GeV/c)");
 535   grBackAvPtVSPtmean->Draw("p");
 536
 537
 538   TLegend *leg=new TLegend(0.7,0.5,0.95,0.2,"","NDC");
 539   leg->AddEntry(grAvPtVSPtmean,"signal","lp");
 540   leg->AddEntry(grBackAvPtVSPtmean,"background","lp");
 541   leg->Draw();
 542
 543   // PRINTING RESULTS
 544   Double_t *ygrSignal,*ygrBack,*eygrSignal,*eygrBack;
 545   ygrSignal=grAvPtVSPtmean->GetY();
 546   eygrSignal=grAvPtVSPtmean->GetEYhigh();
 547   ygrBack=grBackAvPtVSPtmean->GetY();
 548   eygrBack=grBackAvPtVSPtmean->GetEYhigh();
 549   printf("Av Pt for signal: \n");
 550   for(Int_t in=0;in<grAvPtVSPtmean->GetN();in++){
 551     cout<<ygrSignal[in]<<endl;
 552   }
 553   printf("Error on Av Pt for signal: \n");
 554   for(Int_t in=0;in<grAvPtVSPtmean->GetN();in++){
 555     cout<<eygrSignal[in]<<endl;
 556   }
 557
 558  printf("Av Pt for back: \n");
 559   for(Int_t in=0;in<grAvPtVSPtmean->GetN();in++){
 560     cout<<ygrBack[in]<<endl;
 561   }
 562   printf("Error on Av Pt for back: \n");
 563   for(Int_t in=0;in<grAvPtVSPtmean->GetN();in++){
 564     cout<<eygrBack[in]<<endl;
 565   }
 566
 567
 568
 569
 570   nameout="ptCorrection";
 571   if(useFitForSubtraction)nameout.Append("FitSB");
 572   if(!subtractSB)nameout.Append("NoSBsubtract.root");
 573   else nameout.Append(".root");
 574   TFile *fout=new TFile(nameout.Data(),"RECREATE");
 575   fout->cd();
 576   grAvPtVSPtmean->Write();
 577   grAvRawYieldSpectrum->Write();
 578   grBackAvPtVSPtmean->Write();
 579   cPtSpectra->Write();
 580   for(Int_t bin=firstbin;bin<=lastbin;bin++){
 581     cPtDistrNoSubtr[bin]->Write();
 582   }
 583
 584
 585
 586   fout->Close();
 587
 588   //  histoThCompareCentral->Draw("sames");
 589 }
 590
 591
 592 void DoStandardForD0(Int_t rebin,Bool_t usefit,Bool_t corrforeff=kTRUE,Bool_t useParGenAccLimacc=kTRUE,Int_t firstbin=0,Int_t lastbin=3){// firstbin and lastbin are w.r.t. the signal histo plot, starting from 0
 593
 594   //N.B.  ADAPT THE PART BELOW WITH YOUR PATHS AND HISTO NAMES
 595   TFile *fData=TFile::Open("/Users/administrator/ALICE/CHARM/PbPBdata_2011/TestTrain/2013June4TrainData92_MC61/Data/Standard/RAAvsNPart/FinalMassPlots_v2/RawYieldBoth_tight.root");
 596   TFile *fCF=TFile::Open("/Users/administrator/ALICE/CHARM/PbPBdata_2011/2013Jun8MCptWeightDataMoreSplit/MergeWithMyCode/Standard/EffPrompt/fileEff_D0_from_c.root");//EffForHFspectrumLHC10f6NOTd4Corr.root");
 597   TFile *fDataList=TFile::Open("/Users/administrator/ALICE/CHARM/PbPBdata_2011/TestTrain/2013June4TrainData92_MC61/Data/AnalysisResults.root");
 598   TDirectory *fdir=(TDirectory*)fDataList->Get("PWG3_D2H_d0D0");
 599   TList *lslist=(TList*)fdir->Get("clistTGHCsign_d0D0");
 600   TH2F *hMassPt=(TH2F*)lslist->FindObject("hInvMassPtTGHCsign");
 601   TH1D *hSig=(TH1D*)fData->Get("hSigma");
 602   TH1D *hBackground=(TH1D*)fData->Get("hBackground");
 603   TH1D *hSign=(TH1D*)fData->Get("hSignal");
 604   TH1D *hMean=(TH1D*)fData->Get("hMass");
 605
 606   mesonMass=1.86484;
 607
 608   TH1D *hEffNum=(TH1D*)fCF->Get("hRecoPIDpt");
 609   TH1D *hEffDeNum=(TH1D*)fCF->Get("hMCAccpt");
 610
 611   /////////////////////////////////////////////////////////////////
 612   //   FOR DEBUGGING: CHECK EFFICIENCY PLOTS REBINNING
 613   //   TCanvas *cEff=new TCanvas();
 614   //   cEff->cd();
 615   //   TH1D *hRatio=new TH1D(*(TH1D*)hEffNum);
 616   //   hRatio->Divide(hEffDeNum);
 617
 618   //   TH1D *hPt=(TH1D*)(hMassPt->ProjectionY());
 619
 620   //   hPt->Rebin(rebin);cout<< hPt->GetBinWidth(1)<<endl;
 621   //   cout<< hEffNum->GetBinWidth(1)<<endl;
 622   //   TH1D *hWorkNum=0x0;
 623   //   hWorkNum=CheckBinningAndMerge(hPt,(TH1D*)hEffNum,0.01,23.9);
 624   //   TH1D *hWorkDenum=0x0;
 625   //   hWorkDenum=CheckBinningAndMerge(hPt,(TH1D*)hEffDeNum,0.01.,23.9);
 626   //   TH1D *hWorkRatio=new TH1D(*hWorkNum);
 627   //   hWorkRatio->Divide(hWorkDenum);
 628
 629
 630   //   hWorkRatio->Draw();
 631   //   hPt->Draw("sames");
 632
 633   //   hRatio->Draw("sames");
 634   /////////////////////////////////////////////////////////////////
 635
 636   SetNsigmaForSignal(3.);
 637   SetNsigmaStartSB(5.);
 638   SetUseFitForSubtraction(usefit);
 639   corrForEff=corrforeff;
 640   useParGenAccOverLimAcc=useParGenAccLimacc;
 641   CalculateAveragePt(hMassPt,hBackground,hSig,hEffNum,hEffDeNum,hSign,hMean,rebin,firstbin,lastbin);
 642
 643
 644
 645
 646 }
 647
 648
 649 void DoStandardForDs(Int_t rebin,Bool_t usefit,Bool_t corrforeff=kTRUE,Int_t firstbin=0,Int_t lastbin=3){// firstbin and lastbin are w.r.t. the signal histo plot, starting from 0
 650
 651   //N.B.  ADAPT THE PART BELOW WITH YOUR PATHS AND HISTO NAMES
 652 TFile *fData=TFile::Open("RawYieldBoth.root");
 653 TFile *fCF=TFile::Open("fileEff_Ds_CommonFramework_from_c_Enriched.root");
 654
 655 TH1D *hSig=(TH1D*)fData->Get("hSigma");
 656 TH1D *hMean=(TH1D*)fData->Get("hMass");
 657 TH1D *hBackground=(TH1D*)fData->Get("hBackground");
 658 TH1D *hSign=(TH1D*)fData->Get("hSignal");
 659
 660 TH1D *hEffNum=(TH1D*)fCF->Get("hRecoPIDpt");
 661 TH1D *hEffDeNum=(TH1D*)fCF->Get("hMCAccpt");
 662  if(!hEffNum)return;
 663  if(!hEffDeNum)return;
 664  TFile *file=new TFile("AnalysisResults.root");
 665  TDirectory *dirFile=(TDirectory*)file->Get("PWG3_D2H_InvMassDs");
 666  TList *cOutput = (TList*)dirFile->Get("coutputDs");
 667  TH2F *hMassPt=(TH2F*)cOutput->FindObject("hPtVsMass");
 668  //TH2F *hMassPt=(TH2F*)fData->Get("hPtVsMass");
 669
 670  mesonMass=1.96847;
 671
 672  SetUseFitForSubtraction(usefit);
 673  corrForEff=corrforeff;
 674  CalculateAveragePt(hMassPt,hBackground,hSig,hEffNum,hEffDeNum,hSign,hMean,rebin,firstbin,lastbin);
 675
 676
 677 }
 678 void DoStandardForDplus(Int_t rebin,Bool_t usefit,Bool_t corrforeff=kTRUE,Int_t firstbin=0,Int_t lastbin=2){
 679 // firstbin and lastbin are w.r.t. the signal histo plot, starting from 0
 680
 681   //N.B.  ADAPT THE PART BELOW WITH YOUR PATHS AND HISTO NAMES
 682
 683   TFile *fData=TFile::Open("/Users/administrator/ALICE/CHARM/ppData_2010/PtSelInBin/testWithEff/Dplus/RawYieldBoth.root");
 684   TFile *fCF=TFile::Open("/Users/administrator/ALICE/CHARM/ppData_2010/PtSelInBin/testWithEff/Dplus/CFEfficiency.root");
 685   TH2F *hMassPt=(TH2F*)fData->Get("hPtVsMassTC");
 686   TH1D *hSig=(TH1D*)fData->Get("hSigma");
 687   TH1D *hMean=(TH1D*)fData->Get("hMean");
 688   TH1D *hBackground=(TH1D*)fData->Get("hBackground");
 689   TH1D *hSign=(TH1D*)fData->Get("hSignal");
 690
 691   TH1D *hEffNum=(TH1D*)fCF->Get("RecoPID");
 692   TH1D *hEffDeNum=(TH1D*)fCF->Get("GenAcc");
 693
 694   mesonMass=1.86962;
 695
 696   /////////////////////////////////////////////////////////////////
 697   //   FOR DEBUGGING: CHECK EFFICIENCY PLOTS REBINNING
 698   //   TCanvas *cEff=new TCanvas();
 699   //   cEff->cd();
 700   //   TH1D *hRatio=new TH1D(*(TH1D*)hEffNum);
 701   //   hRatio->Divide(hEffDeNum);
 702
 703   //   TH1D *hPt=(TH1D*)(hMassPt->ProjectionY());
 704
 705   //   hPt->Rebin(rebin);cout<< hPt->GetBinWidth(1)<<endl;
 706   //   cout<< hEffNum->GetBinWidth(1)<<endl;
 707   //   TH1D *hWorkNum=0x0;
 708   //   hWorkNum=CheckBinningAndMerge(hPt,(TH1D*)hEffNum,0.01,23.9);
 709   //   TH1D *hWorkDenum=0x0;
 710   //   hWorkDenum=CheckBinningAndMerge(hPt,(TH1D*)hEffDeNum,0.01.,23.9);
 711   //   TH1D *hWorkRatio=new TH1D(*hWorkNum);
 712   //   hWorkRatio->Divide(hWorkDenum);
 713   //   hWorkRatio->Draw();
 714   //   hPt->Draw("sames");
 715   //   hRatio->Draw("sames");
 716   /////////////////////////////////////////////////////////////////
 717
 718   SetUseFitForSubtraction(usefit);
 719   corrForEff=corrforeff;
 720   CalculateAveragePt(hMassPt,hBackground,hSig,hEffNum,hEffDeNum,hSign,hMean,rebin,firstbin,lastbin);
 721
 722
 723
 724
 725 }
 726
 727
 728
 729
 730 TH1D *SmearEffHisto(TH1D *hInput,TString name="hEffNum",Double_t maxPt=40.,Double_t step=0.1){
 731   Int_t nbins=(Int_t)(maxPt/step);
 732   Int_t maxBinA=hInput->GetNbinsX();
 733   Int_t binBl,binBu;
 734   Double_t binlea,binuea,contInMin,contInMax;
 735
 736   TH1D *h=new TH1D(name.Data(),name.Data(),nbins,0.,maxPt);
 737   for(Int_t bin=1;bin<=maxBinA;bin++){
 738     contInMin=hInput->GetBinContent(bin);
 739     if(bin==maxBinA) contInMax=contInMin*1.2;
 740     else contInMax=hInput->GetBinContent(bin+1);
 741     binlea=hInput->GetBinLowEdge(bin);
 742     binuea=hInput->GetXaxis()->GetBinUpEdge(bin);
 743     binBl=h->FindBin(binlea);
 744     binBu=h->FindBin(binuea);
 745     for(Int_t bb=binBl;bb<binBu;bb++){
 746       h->SetBinContent(bb,contInMin+(bb-binBl)*(contInMax-contInMin)/(binBu-binBl));
 747     }
 748   }
 749   return h;
 750
 751 }
 752
 753
 754 void AverageD0DplusResults(TString fileD0="/Users/administrator/ALICE/CHARM/ppData_2010/2011_Jul_05/data/LHC10bcdeAOD057/AvPt/MassRegSel3PkMore5SBEffCorrMeanFit/ptCorrectionFitSB.root",TString fileDplus="/Users/administrator/ALICE/CHARM/ppData_2010/2011_Jul_05/Dplus/AvPt/2011Jul26Renu/average_ptNew.root"){
 755
 756   Double_t *yD0,*eyD0,*yDplus,*eyDplus,*x,*exLow,*exHigh;
 757   Double_t error=0.,average=0.;
 758
 759   TFile *fD0=TFile::Open(fileD0.Data());
 760   // TCanvas *cD0=(TCanvas*)fD0->Get("cAvPtVSPtmean");
 761   //  TGraphAsymmErrors *grD0=(TGraphAsymmErrors*)cD0->FindObject("grAvPtMeanPtFitSB");
 762   TGraphAsymmErrors *grD0=(TGraphAsymmErrors*)fD0->Get("grAvPtMeanPtFitSB");
 763
 764   grD0->SetName("grAvPtMeanPtFitSB_D0");
 765   yD0=grD0->GetY();
 766   eyD0=grD0->GetEYhigh();
 767
 768
 769   TFile *fDplus=TFile::Open(fileDplus.Data());
 770   TCanvas *cDplus=(TCanvas*)fDplus->Get("cAvPtVSPtmean");
 771   TGraphAsymmErrors *grDplus=(TGraphAsymmErrors*)cDplus->FindObject("grAvPtMeanPtFitSB");
 772   grDplus->SetName("grAvPtMeanPtFitSB_Dplus");
 773   yDplus=grDplus->GetY();
 774   eyDplus=grDplus->GetEYhigh();
 775
 776
 777   Int_t nlmax=grDplus->GetN();
 778   Int_t nlD0=grD0->GetN();
 779   Int_t nlDplus=grDplus->GetN();
 780
 781   if(grD0->GetN()>nlmax){
 782     nlmax=grD0->GetN();
 783     x=grD0->GetX();
 784     exLow=grD0->GetEXlow();
 785     exHigh=grD0->GetEXhigh();
 786   }
 787   else {
 788     x=grDplus->GetX();
 789     exLow=grDplus->GetEXlow();
 790     exHigh=grDplus->GetEXhigh();
 791   }
 792
 793
 794   TGraphAsymmErrors *grAverage=new TGraphAsymmErrors();
 795   grAverage->SetName("grAverageD0Dplus");
 796
 797   for(Int_t j=0;j<nlmax;j++){
 798     if(j<nlD0&&j<nlDplus){
 799       average=(yD0[j]/(eyD0[j]*eyD0[j])+yDplus[j]/(eyDplus[j]*eyDplus[j]))/(1./(eyD0[j]*eyD0[j])+1./(eyDplus[j]*eyDplus[j]));
 800       error=TMath::Sqrt(1./(1./(eyD0[j]*eyD0[j])+1./(eyDplus[j]*eyDplus[j])));
 801     }
 802     else if(j>=nlD0){
 803       average=yDplus[j];
 804       error=eyDplus[j];
 805     }
 806     else if(j>=nlDplus){
 807       average=yD0[j];
 808       error=eyD0[j];
 809     }
 810     grAverage->SetPoint(j,x[j],average);
 811     grAverage->SetPointError(j,exLow[j],exHigh[j],error,error);
 812   }
 813
 814   // Printing Info: D0
 815   printf("Average Value D0: \n");
 816   for(Int_t j=0;j<nlD0;j++){
 817     cout<<yD0[j]<<endl;
 818
 819   }
 820
 821  printf("Average Errors D0: \n");
 822   for(Int_t j=0;j<nlD0;j++){
 823     cout<<eyD0[j]<<endl;
 824
 825   }
 826
 827   // Printing Info: Dplus
 828   printf("Average Value Dplus: \n");
 829   for(Int_t j=0;j<nlDplus;j++){
 830     cout<<yDplus[j]<<endl;
 831
 832   }
 833
 834   printf("Average Errors Dplus: \n");
 835   for(Int_t j=0;j<nlDplus;j++){
 836     cout<<eyDplus[j]<<endl;
 837
 838   }
 839
 840
 841   // Printing Info: Dplus D0 average
 842   yDplus=grAverage->GetY();
 843   eyDplus=grAverage->GetEYhigh();
 844   printf("Average Value D0 + Dplus: \n");
 845   for(Int_t j=0;j<nlDplus;j++){
 846     cout<<yDplus[j]<<endl;
 847
 848   }
 849
 850   printf("Average Errors D0+Dplus: \n");
 851   for(Int_t j=0;j<nlDplus;j++){
 852     cout<<eyDplus[j]<<endl;
 853
 854   }
 855
 856
 857
 858
 859   TCanvas *cCompare=new TCanvas("cCompareD0Dplus","cCompareD0Dplus",700,700);
 860   cCompare->cd();
 861
 862   grDplus->SetMarkerStyle(21);
 863   grDplus->SetMarkerSize(1.2);
 864   grDplus->SetMarkerColor(kBlue);
 865   grDplus->SetLineColor(kBlue);
 866   grDplus->Draw("ap");
 867
 868   grD0->SetMarkerStyle(22);
 869   grD0->SetMarkerSize(1.2);
 870   grD0->SetMarkerColor(kRed);
 871   grD0->SetLineColor(kRed);
 872   grD0->Draw("p");
 873
 874   grAverage->SetMarkerStyle(20);
 875   grAverage->SetMarkerSize(1.2);
 876   grAverage->SetMarkerColor(kBlack);
 877   grAverage->SetLineColor(kBlack);
 878   grAverage->Draw("p");
 879
 880
 881
 882 }