PWGHF/vertexingHF/macros/ReadDvsMultiplicity.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)\r
   2 #include <TInterpreter.h>\r
   3 #include <TString.h>\r
   4 #include <TObjString.h>\r
   5 #include <TObjArray.h>\r
   6 #include <TMath.h>\r
   7 #include <TFile.h>\r
   8 #include <TCanvas.h>\r
   9 #include <TH1.h>\r
  10 #include <TH1F.h>\r
  11 #include <TH2F.h>\r
  12 #include <TH3.h>\r
  13 #include <TH3F.h>\r
  14 #include <TH1D.h>\r
  15 #include <TF1.h>\r
  16 #include <TSystem.h>\r
  17 #include <TStyle.h>\r
  18 #include <TLegend.h>\r
  19 #include <TList.h>\r
  20 #include <TLegendEntry.h>\r
  21 #include <TDatabasePDG.h>\r
  22 #include <TGraph.h>\r
  23 \r
  24 #include "AliAODRecoDecayHF.h"\r
  25 #include "AliRDHFCuts.h"\r
  26 #include "AliRDHFCutsDplustoKpipi.h"\r
  27 #include "AliRDHFCutsD0toKpi.h"\r
  28 #include "AliHFMassFitter.h"\r
  29 #include "AliNormalizationCounter.h"\r
  30 #endif\r
  31 \r
  32 /* $Id$ */ \r
  33 \r
  34 enum {kD0toKpi, kDplusKpipi};\r
  35 enum {kBoth, kParticleOnly, kAntiParticleOnly};\r
  36 enum {kExpo=0, kLinear, kPol2};\r
  37 enum {kGaus=0, kDoubleGaus};\r
  38 enum {kCorr=0, kUnCorr, kNoPid};\r
  39 \r
  40 \r
  41 // Common variables: to be configured by the user\r
  42 const Int_t nPtBins=5;\r
  43 Double_t ptlims[nPtBins+1]={2.,4.,6.,8.,12.,24.};\r
  44 Int_t rebin[nPtBins]={4,4,6,6,8};\r
  45 Double_t sigmapt[nPtBins]={ 0.010, 0.012, 0.016, 0.018, 0.20 };\r
  46 Bool_t fixPeakSigma = kFALSE;\r
  47 //\r
  48 const Int_t nMultbins=6;\r
  49 Double_t multlims[nMultbins+1]={1.,9.,14.,20.,31.,49.,100.};\r
  50 // const Int_t nMultbins=1;\r
  51 // Double_t multlims[nMultbins+1]={0.,500.};\r
  52 //\r
  53 Int_t firstUsedBin[nPtBins]={-1,-1,-1,-1,-1}; // -1 uses all bins, >=1 to set the lower bin to be accepted from original histo\r
  54 //\r
  55 Bool_t isMC=kFALSE;\r
  56 Int_t typeb=kExpo;\r
  57 Int_t types=kGaus;\r
  58 Int_t optPartAntiPart=kBoth;\r
  59 Int_t factor4refl=0;\r
  60 Float_t massRangeForCounting=0.05; // GeV\r
  61 TH2F* hPtMass=0x0;\r
  62 TString suffix="StdPid";\r
  63 //for D0only\r
  64 const Int_t nsamples=2;//3;\r
  65 Int_t nevents[nsamples]={1.18860695e+08 /*LHC10dnewTPCpid*/,9.0374946e+07 /*LHC10b+c*/};\r
  66 \r
  67 // Functions\r
  68 Bool_t LoadDplusHistos(TObjArray* listFiles, TH3F** hPtMassMult, TH2F** hNtrZvtx, TH2F** hNtrZvtxCorr, const char *CutsType, Int_t Option);\r
  69 Bool_t LoadD0toKpiHistos(TObjArray* listFiles, TH3F** hPtMassMult, TH2F** hNtrZvtx, TH2F** hNtrZvtxCorr, AliNormalizationCounter *counter,\r
  70                          const char *CutsType, Int_t Option);\r
  71 Bool_t CheckNtrVsZvtx(TH2F** hNtrZvtx, TH2F** hNtrZvtxCorr, Int_t nFiles);\r
  72 TH1F* RebinHisto(TH1F* hOrig, Int_t reb, Int_t firstUse=-1);\r
  73 \r
  74 \r
  75 void ReadDvsMultiplicity(Int_t analysisType=kD0toKpi,\r
  76                          TString fileNameb="AnalysisResults.root",\r
  77                          TString fileNamec="",\r
  78                          TString fileNamed="",\r
  79                          TString fileNamee="",\r
  80                          const char *CutsType="",\r
  81                          Int_t Option=kCorr)\r
  82 {\r
  83   // gSystem->SetIncludePath("-I. -I$ROOTSYS/include -I$ALICE_ROOT -I$ALICE_ROOT/include -I$ALICE_ROOT/ITS -I$ALICE_ROOT/TPC -I$ALICE_ROOT/CONTAINERS -I$ALICE_ROOT/STEER/STEER -I$ALICE_ROOT/STEER/STEERBase -I$ALICE_ROOT/STEER/ESD -I$ALICE_ROOT/STEER/AOD -I$ALICE_ROOT/TRD -I$ALICE_ROOT/macros -I$ALICE_ROOT/ANALYSIS  -I$ALICE_ROOT/OADB -I$ALICE_ROOT/PWGHF -I$ALICE_ROOT/PWGHF/base -I$ALICE_ROOT/PWGHF/vertexingHF -I$ALICE_ROOT/PWG/FLOW/Case -I$ALICE_ROOT/PWG/FLOW/Tasks -g");\r
  84 \r
  85   // gInterpreter->ExecuteMacro("$ALICE_ROOT/PWGHF/vertexingHF/macros/LoadLibraries.C");\r
  86   gStyle->SetOptTitle(1);\r
  87 \r
  88   Int_t nFiles=0;\r
  89   TObjArray* listFiles=new TObjArray();\r
  90   if(fileNameb!="") { listFiles->AddLast(new TObjString(fileNameb.Data())); nFiles++; }\r
  91   if(fileNamec!="") { listFiles->AddLast(new TObjString(fileNamec.Data())); nFiles++; }\r
  92   if(fileNamed!="") { listFiles->AddLast(new TObjString(fileNamed.Data())); nFiles++; }\r
  93   if(fileNamee!="") { listFiles->AddLast(new TObjString(fileNamee.Data())); nFiles++; }\r
  94   if(listFiles->GetEntries()==0){\r
  95     printf("Missing file names in input\n");\r
  96     return;\r
  97   }\r
  98   TH3F** hPtMassMult=new TH3F*[1];\r
  99   hPtMassMult[0]=0x0;\r
 100   TH1F** hmass=new TH1F*[nPtBins*nMultbins];\r
 101   for(Int_t i=0; i<nPtBins*nMultbins; i++) hmass[i]=0x0;\r
 102   TH2F** hNtrZvtx=new TH2F*[4];\r
 103   for(Int_t i=0; i<4; i++) hNtrZvtx[i]=0x0;\r
 104   TH2F** hNtrZvtxCorr=new TH2F*[4];\r
 105   for(Int_t i=0; i<4; i++) hNtrZvtxCorr[i]=0x0;\r
 106   AliNormalizationCounter *counter=0x0;\r
 107   TH1F * hNormalization = new TH1F("hNormalization","Events in the norm counter",nMultbins,multlims);\r
 108 \r
 109   Float_t massD;\r
 110   Bool_t retCode;\r
 111   if(analysisType==kD0toKpi) {\r
 112     massD=1.86484003067016602e+00;//(Float_t)(TDatabasePDG::Instance()->GetParticle(421)->Mass());\r
 113     retCode=LoadD0toKpiHistos(listFiles,hPtMassMult,hNtrZvtx,hNtrZvtxCorr,counter,CutsType,Option);\r
 114   }\r
 115   else if(analysisType==kDplusKpipi) {\r
 116     massD=1.86961996555328369e+00;//(Float_t)(TDatabasePDG::Instance()->GetParticle(411)->Mass());\r
 117     retCode=LoadDplusHistos(listFiles,hPtMassMult,hNtrZvtx,hNtrZvtxCorr,CutsType,Option);\r
 118   }\r
 119   else{\r
 120     printf("Wrong analysis type parameter\n");\r
 121     return;\r
 122   }\r
 123   if(!retCode){\r
 124     printf("ERROR in reading input files\n");\r
 125     return;\r
 126   }\r
 127   //\r
 128   // Check the multiplicity variables\r
 129   //\r
 130   Bool_t isMult = CheckNtrVsZvtx(hNtrZvtx,hNtrZvtxCorr,nFiles);\r
 131   if(!isMult) return;\r
 132   //\r
 133   //\r
 134   printf(" Preparing the mass fits");\r
 135   TH1F* hCntSig1[nMultbins];\r
 136   TH1F* hCntSig2[nMultbins];\r
 137   TH1F* hNDiffCntSig1[nMultbins];\r
 138   TH1F* hNDiffCntSig2[nMultbins];\r
 139   TH1F* hSignal[nMultbins];\r
 140   TH1F* hRelErrSig[nMultbins];\r
 141   TH1F* hInvSignif[nMultbins];\r
 142   TH1F* hBackground[nMultbins];\r
 143   TH1F* hBackgroundNormSigma[nMultbins];\r
 144   TH1F* hSignificance[nMultbins];\r
 145   TH1F* hMass[nMultbins];\r
 146   TH1F* hSigma[nMultbins];\r
 147   for(Int_t i=0; i<nMultbins; i++){\r
 148     hCntSig1[i]=new TH1F(Form("hCntSig1_%d",i),Form("hCntSig1_%d",i),nPtBins,ptlims);\r
 149     hCntSig2[i]=new TH1F(Form("hCntSig2_%d",i),Form("hCntSig2_%d",i),nPtBins,ptlims);\r
 150     hNDiffCntSig1[i]=new TH1F(Form("hNDiffCntSig1_%d",i),Form("hNDiffCntSig1_%d",i),nPtBins,ptlims);\r
 151     hNDiffCntSig2[i]=new TH1F(Form("hNDiffCntSig2_%d",i),Form("hNDiffCntSig2_%d",i),nPtBins,ptlims);\r
 152     hSignal[i]=new TH1F(Form("hSignal_%d",i),Form("hSignal_%d",i),nPtBins,ptlims);\r
 153     hRelErrSig[i]=new TH1F(Form("hRelErrSig_%d",i),Form("hRelErrSig_%d",i),nPtBins,ptlims);\r
 154     hInvSignif[i]=new TH1F(Form("hInvSignif_%d",i),Form("hInvSignif_%d",i),nPtBins,ptlims);\r
 155     hBackground[i]=new TH1F(Form("hBackground_%d",i),Form("hBackground_%d",i),nPtBins,ptlims);\r
 156     hBackgroundNormSigma[i]=new TH1F(Form("hBackgroundNormSigma_%d",i),Form("hBackgroundNormSigma_%d",i),nPtBins,ptlims);\r
 157     hSignificance[i]=new TH1F(Form("hSignificance_%d",i),Form("hSignificance_%d",i),nPtBins,ptlims);\r
 158     hMass[i]=new TH1F(Form("hMass_%d",i),Form("hMass_%d",i),nPtBins,ptlims);\r
 159     hSigma[i]=new TH1F(Form("hSigma_%d",i),Form("hSigma_%d",i),nPtBins,ptlims);\r
 160   }\r
 161   printf(", defined...\n");\r
 162 \r
 163   //  std::cout << " htemp :"<<hPtMassMult[0]<<std::endl;\r
 164   TH1F* hptaxis = (TH1F*)hPtMassMult[0]->ProjectionZ("hptaxis");\r
 165   TH1F* hmassaxis = (TH1F*)hPtMassMult[0]->ProjectionY("hmassaxis");\r
 166   TH1F* hmultaxis = (TH1F*)hPtMassMult[0]->ProjectionX("hmultaxis");\r
 167   Int_t nMassBins=hmassaxis->GetNbinsX();\r
 168   Double_t hmin=hmassaxis->GetBinLowEdge(3);\r
 169   Double_t hmax=hmassaxis->GetBinLowEdge(nMassBins-2) + hmassaxis->GetBinWidth(nMassBins-2);\r
 170   Float_t minBinSum=hmassaxis->FindBin(massD-massRangeForCounting);\r
 171   Float_t maxBinSum=hmassaxis->FindBin(massD+massRangeForCounting);\r
 172   Int_t iPad=1;\r
 173   \r
 174   printf("Now initializing the fit functions\n");\r
 175   TF1* funBckStore1=0x0;\r
 176   TF1* funBckStore2=0x0;\r
 177   TF1* funBckStore3=0x0;\r
 178 \r
 179   Int_t nPtMultbins = nPtBins*nMultbins;\r
 180   AliHFMassFitter** fitter=new AliHFMassFitter*[nPtMultbins];\r
 181   Double_t arrchisquare0[nPtBins];\r
 182   Double_t arrchisquare1[nPtBins];\r
 183   Double_t arrchisquare2[nPtBins];\r
 184   Double_t arrchisquare3[nPtBins];\r
 185   Double_t arrchisquare4[nPtBins];\r
 186   Double_t arrchisquare5[nPtBins];\r
 187   for(Int_t i=0; i<nPtBins; i++){\r
 188     arrchisquare0[i]=0.;\r
 189     arrchisquare1[i]=0.;\r
 190     arrchisquare2[i]=0.;\r
 191     arrchisquare3[i]=0.;\r
 192     arrchisquare4[i]=0.;\r
 193     arrchisquare5[i]=0.;\r
 194   }\r
 195   \r
 196   TCanvas* canvas[nMultbins];\r
 197   Int_t nx=2, ny=2;\r
 198   if(nPtBins>4){\r
 199     nx=3;\r
 200     ny=2;\r
 201   }\r
 202   if(nPtBins>6){\r
 203     nx=4;\r
 204     ny=3;\r
 205   }\r
 206   if(nPtBins>12){\r
 207     nx=5;\r
 208     ny=4;\r
 209   }\r
 210   for(Int_t i=0; i<nMultbins; i++ ){\r
 211     canvas[i] = new TCanvas(Form("c%d",i),Form("summary canvas for mult bin %d",i));\r
 212     canvas[i]->Divide(nx,ny);\r
 213   }\r
 214   TCanvas *myCanvas[nPtMultbins];\r
 215   \r
 216   //\r
 217   //\r
 218   // Loop on multiplicity bins\r
 219   //\r
 220   Int_t massBin=0;\r
 221   Double_t sig,errsig,s,errs,b,errb;\r
 222   for(Int_t j=0; j<nMultbins; j++){\r
 223     //    printf(" Studying multiplicity bin %d\n",j);\r
 224     Int_t multbinlow = hmultaxis->FindBin(multlims[j]);\r
 225     Int_t multbinhigh = hmultaxis->FindBin(multlims[j+1])-1;\r
 226     Float_t val = multbinlow + (multbinhigh-multbinlow)/2.;\r
 227     Int_t hnbin = hNormalization->FindBin(val);\r
 228     Float_t nevents = 0.;\r
 229     if(counter) { nevents = counter->GetNEventsForNorm(multbinlow,multbinhigh); cout << endl<<endl<<" Nevents ("<<multbinlow<<","<<multbinhigh<<") ="<<nevents<<endl<<endl<<endl;}\r
 230     hNormalization->SetBinContent(hnbin,nevents);\r
 231     //\r
 232     // Loop on pt bins\r
 233     //\r
 234     iPad=1;\r
 235     for(Int_t iBin=0; iBin<nPtBins; iBin++){\r
 236       canvas[j]->cd(iPad++);\r
 237       //      printf(" projecting to the mass histogram\n");\r
 238       Int_t ptbinlow = hptaxis->FindBin(ptlims[iBin]);\r
 239       Int_t ptbinhigh = hptaxis->FindBin(ptlims[iBin+1])-1;\r
 240       hmass[massBin] = (TH1F*)hPtMassMult[0]->ProjectionY(Form("hmass_%d_%d",j,iBin),multbinlow,multbinhigh,ptbinlow,ptbinhigh);\r
 241       if(  hmass[massBin]->GetEntries() < 60 ) {\r
 242         massBin++;\r
 243         continue;\r
 244       }\r
 245       Int_t origNbins=hmass[massBin]->GetNbinsX(); \r
 246       //      printf(" rebinning the mass histogram\n");\r
 247       TH1F* hRebinned=RebinHisto(hmass[massBin],rebin[iBin],firstUsedBin[iBin]);\r
 248       hmin=hRebinned->GetBinLowEdge(2);\r
 249       hmax=hRebinned->GetBinLowEdge(hRebinned->GetNbinsX());\r
 250       //      printf(" define the mass fitter and options \n");\r
 251       fitter[massBin] = new AliHFMassFitter(hRebinned,hmin, hmax,1,typeb,types);\r
 252       fitter[massBin]->SetRangeFit(1.65,2.10);\r
 253       Int_t rebinItem = origNbins/fitter[massBin]->GetBinN();\r
 254       fitter[massBin]->SetReflectionSigmaFactor(factor4refl);\r
 255       fitter[massBin]->SetInitialGaussianMean(massD);\r
 256       fitter[massBin]->SetInitialGaussianSigma(sigmapt[iBin]);\r
 257       if(fixPeakSigma) {\r
 258         fitter[massBin]->SetFixGaussianMean(massD);\r
 259         fitter[massBin]->SetFixGaussianSigma(sigmapt[iBin],kTRUE);\r
 260       }\r
 261       Bool_t out=fitter[massBin]->MassFitter(0);\r
 262       if(!out) {\r
 263         fitter[massBin]->GetHistoClone()->Draw();\r
 264         massBin++;\r
 265         continue;\r
 266       }\r
 267       //      printf(" getting the fit parameters\n");\r
 268       Double_t mass=fitter[massBin]->GetMean();\r
 269       Double_t massUnc=fitter[massBin]->GetMeanUncertainty();\r
 270       Double_t sigma=fitter[massBin]->GetSigma();\r
 271       Double_t sigmaUnc=fitter[massBin]->GetSigmaUncertainty();\r
 272       if(j==0) arrchisquare0[iBin]=fitter[massBin]->GetReducedChiSquare();\r
 273       else if(j==1) arrchisquare1[iBin]=fitter[massBin]->GetReducedChiSquare();\r
 274       else if(j==2) arrchisquare2[iBin]=fitter[massBin]->GetReducedChiSquare();\r
 275       else if(j==3) arrchisquare3[iBin]=fitter[massBin]->GetReducedChiSquare();\r
 276       else if(j==4) arrchisquare4[iBin]=fitter[massBin]->GetReducedChiSquare();\r
 277       else if(j==5) arrchisquare5[iBin]=fitter[massBin]->GetReducedChiSquare();\r
 278       TF1* fB1=fitter[massBin]->GetBackgroundFullRangeFunc();\r
 279       TF1* fB2=fitter[massBin]->GetBackgroundRecalcFunc();\r
 280       TF1* fM=fitter[massBin]->GetMassFunc();\r
 281       if(iBin==0 && fB1) funBckStore1=new TF1(*fB1);\r
 282       if(iBin==0 && fB2) funBckStore2=new TF1(*fB2);\r
 283       if(iBin==0 && fM) funBckStore3=new TF1(*fM);\r
 284 \r
 285       fitter[massBin]->DrawHere(gPad);\r
 286       fitter[massBin]->Signal(3,s,errs);\r
 287       fitter[massBin]->Background(3,b,errb);\r
 288       fitter[massBin]->Significance(3,sig,errsig);\r
 289       Double_t ry=fitter[massBin]->GetRawYield();\r
 290       Double_t ery=fitter[massBin]->GetRawYieldError();\r
 291       myCanvas[massBin] = new TCanvas(Form("myCanvas_%d_%d",j,iBin),Form("Invariant mass mult bin %d, pt bin %d",j,iBin));\r
 292       fitter[massBin]->DrawHere(gPad);\r
 293     \r
 294       Float_t cntSig1=0.;\r
 295       Float_t cntSig2=0.;\r
 296       Float_t cntErr=0.;\r
 297       for(Int_t iMB=minBinSum; iMB<=maxBinSum; iMB++){\r
 298         Float_t bkg1=fB1 ? fB1->Eval(hmass[massBin]->GetBinCenter(iMB))/rebinItem : 0;\r
 299         Float_t bkg2=fB2 ? fB2->Eval(hmass[massBin]->GetBinCenter(iMB))/rebinItem : 0;\r
 300         cntSig1+=(hmass[massBin]->GetBinContent(iMB)-bkg1);\r
 301         cntSig2+=(hmass[massBin]->GetBinContent(iMB)-bkg2);\r
 302         cntErr+=(hmass[massBin]->GetBinContent(iMB));\r
 303       }\r
 304       hCntSig1[j]->SetBinContent(iBin+1,cntSig1);\r
 305       hCntSig1[j]->SetBinError(iBin+1,TMath::Sqrt(cntErr));\r
 306       hNDiffCntSig1[j]->SetBinContent(iBin+1,(s-cntSig1)/s);\r
 307       hNDiffCntSig1[j]->SetBinError(iBin+1,TMath::Sqrt(cntErr)/s);\r
 308       hCntSig2[j]->SetBinContent(iBin+1,cntSig2);\r
 309       hNDiffCntSig2[j]->SetBinContent(iBin+1,(s-cntSig2)/s);\r
 310       hNDiffCntSig2[j]->SetBinError(iBin+1,TMath::Sqrt(cntErr)/s);\r
 311       hCntSig2[j]->SetBinError(iBin+1,TMath::Sqrt(cntErr));\r
 312       hSignal[j]->SetBinContent(iBin+1,ry);\r
 313       hSignal[j]->SetBinError(iBin+1,ery);\r
 314       hRelErrSig[j]->SetBinContent(iBin+1,errs/s);\r
 315       hInvSignif[j]->SetBinContent(iBin+1,1/sig);\r
 316       hInvSignif[j]->SetBinError(iBin+1,errsig/(sig*sig));\r
 317       hBackground[j]->SetBinContent(iBin+1,b); //consider sigma\r
 318       hBackground[j]->SetBinError(iBin+1,errb);\r
 319       hBackgroundNormSigma[j]->SetBinContent(iBin+1,b/(3*fitter[iBin]->GetSigma())*(3*0.012)); //consider sigma\r
 320       hBackgroundNormSigma[j]->SetBinError(iBin+1,errb);\r
 321       hSignificance[j]->SetBinContent(iBin+1,sig);\r
 322       hSignificance[j]->SetBinError(iBin+1,errsig);\r
 323       hMass[j]->SetBinContent(iBin+1,mass);\r
 324       hMass[j]->SetBinError(iBin+1,massUnc);\r
 325       hSigma[j]->SetBinContent(iBin+1,sigma);\r
 326       hSigma[j]->SetBinError(iBin+1,sigmaUnc);\r
 327 \r
 328       massBin++;\r
 329       delete hRebinned;\r
 330     }// end loop on pt bins\r
 331 \r
 332     canvas[j]->Update();\r
 333     canvas[j]->SaveAs(Form("hMass%s_%d_%d.eps",CutsType,typeb,j));\r
 334     \r
 335   }// end loop on multiplicity bins\r
 336 \r
 337 \r
 338   TCanvas *cpar=new TCanvas("cpar","Fit params",1200,600);\r
 339   cpar->Divide(2,1);\r
 340   cpar->cd(1);\r
 341   for(Int_t imult=0; imult<nMultbins; imult++) {\r
 342     hMass[imult]->SetMarkerStyle(20);\r
 343     hMass[imult]->GetXaxis()->SetTitle("Pt (GeV/c)");\r
 344     hMass[imult]->GetYaxis()->SetTitle("Mass (GeV/c^{2})");\r
 345     hMass[imult]->SetMarkerColor(2*imult);\r
 346     if(imult==5) hMass[imult]->SetMarkerColor(2*imult-3);\r
 347     if(imult==0) {\r
 348       hMass[imult]->SetMarkerColor(kBlack);\r
 349       hMass[imult]->Draw("PE");\r
 350     }\r
 351     else  hMass[imult]->Draw("PEsame");\r
 352   }\r
 353   cpar->cd(2);\r
 354   for(Int_t imult=0; imult<nMultbins; imult++) {\r
 355     hSigma[imult]->SetMarkerStyle(20);\r
 356     //  hSigma[0]->Draw("PE");\r
 357     hSigma[imult]->GetXaxis()->SetTitle("Pt (GeV/c)");\r
 358     hSigma[imult]->GetYaxis()->SetTitle("Sigma (GeV/c^{2})");\r
 359     hSigma[imult]->SetMarkerColor(2*imult);\r
 360     if(imult==5) hSigma[imult]->SetMarkerColor(2*imult-3);\r
 361     if(imult==0) {\r
 362       hSigma[imult]->SetMarkerColor(kBlack);\r
 363       hSigma[imult]->Draw("PE");\r
 364     }\r
 365     else  hSigma[imult]->Draw("PEsame");\r
 366   }\r
 367   cpar->Update();\r
 368 \r
 369   /*\r
 370   TCanvas** csig;//= new TCanvas*[nMultbins];\r
 371   TCanvas** cDiffS;//=new TCanvas*[nMultbins];\r
 372   for(Int_t i=0; i<nMultbins; i++){\r
 373     csig[i] =new TCanvas(Form("csig_%d",i),Form("Results, multiplicity bin %d",i),1200,600);\r
 374     csig[i]->Divide(3,1);\r
 375     csig[i]->cd(1);\r
 376     hSignal[i]->SetMarkerStyle(20);\r
 377     hSignal[i]->SetMarkerColor(4);\r
 378     hSignal[i]->SetLineColor(4);\r
 379     hSignal[i]->GetXaxis()->SetTitle("Pt (GeV/c)");\r
 380     hSignal[i]->GetYaxis()->SetTitle("Signal");\r
 381     hSignal[i]->Draw("P");\r
 382     hCntSig1[i]->SetMarkerStyle(26);\r
 383     hCntSig1[i]->SetMarkerColor(2);\r
 384     hCntSig1[i]->SetLineColor(2);\r
 385     hCntSig1[i]->Draw("PSAME");\r
 386     hCntSig2[i]->SetMarkerStyle(29);\r
 387     hCntSig2[i]->SetMarkerColor(kGray+1);\r
 388     hCntSig2[i]->SetLineColor(kGray+1);\r
 389     hCntSig2[i]->Draw("PSAME");\r
 390     TLegend* leg=new TLegend(0.4,0.7,0.89,0.89);\r
 391     leg->SetFillColor(0);\r
 392     TLegendEntry* ent=leg->AddEntry(hSignal[i],"From Fit","PL");\r
 393     ent->SetTextColor(hSignal[i]->GetMarkerColor());\r
 394     ent=leg->AddEntry(hCntSig1[i],"From Counting1","PL");\r
 395     ent->SetTextColor(hCntSig1[i]->GetMarkerColor());\r
 396     ent=leg->AddEntry(hCntSig2[i],"From Counting2","PL");\r
 397     ent->SetTextColor(hCntSig2[i]->GetMarkerColor());\r
 398     leg->Draw();\r
 399     csig[i]->cd(2);\r
 400     hBackground[i]->SetMarkerStyle(20);\r
 401     hBackground[i]->Draw("P");\r
 402     hBackground[i]->GetXaxis()->SetTitle("Pt (GeV/c)");\r
 403     hBackground[i]->GetYaxis()->SetTitle("Background");\r
 404     csig[i]->cd(3);\r
 405     hSignificance[i]->SetMarkerStyle(20);\r
 406     hSignificance[i]->Draw("P");\r
 407     hSignificance[i]->GetXaxis()->SetTitle("Pt (GeV/c)");\r
 408     hSignificance[i]->GetYaxis()->SetTitle("Significance");\r
 409     \r
 410     cDiffS[i]=new TCanvas(Form("cDiffS_%d",i),Form("Difference, multiplicity bin %d",i),1200,600);\r
 411     cDiffS[i]->Divide(2,1);\r
 412     cDiffS[i]->cd(1);\r
 413     hRelErrSig[i]->SetMarkerStyle(20); //fullcircle\r
 414     hRelErrSig[i]->SetTitleOffset(1.2);  \r
 415     hRelErrSig[i]->SetTitle("Rel Error from Fit;p_{t} (GeV/c);Signal Relative Error");\r
 416     hRelErrSig[i]->Draw("P");\r
 417     hInvSignif[i]->SetMarkerStyle(21); //fullsquare\r
 418     hInvSignif[i]->SetMarkerColor(kMagenta+1);\r
 419     hInvSignif[i]->SetLineColor(kMagenta+1);\r
 420     hInvSignif[i]->Draw("PSAMES");\r
 421     TLegend* leg2=new TLegend(0.4,0.7,0.89,0.89);\r
 422     leg2->SetFillColor(0);\r
 423     TLegendEntry* ent2=leg2->AddEntry(hRelErrSig[i],"From Fit","P");\r
 424     ent2->SetTextColor(hRelErrSig[i]->GetMarkerColor());\r
 425     ent2=leg2->AddEntry(hInvSignif[i],"1/Significance","PL");\r
 426     ent2->SetTextColor(hInvSignif[i]->GetMarkerColor());\r
 427     leg2->Draw();\r
 428 \r
 429     cDiffS[i]->cd(2);\r
 430     hNDiffCntSig1[i]->SetMarkerStyle(26);\r
 431     hNDiffCntSig1[i]->SetMarkerColor(2);\r
 432     hNDiffCntSig1[i]->SetLineColor(2);\r
 433     hNDiffCntSig1[i]->SetTitle("Cmp Fit-Count;p_{t} (GeV/c);(S_{fit}-S_{count})/S_{fit}");\r
 434     hNDiffCntSig1[i]->Draw("P");\r
 435     hNDiffCntSig2[i]->SetMarkerStyle(29);\r
 436     hNDiffCntSig2[i]->SetMarkerColor(kGray+1);\r
 437     hNDiffCntSig2[i]->SetLineColor(kGray+1);\r
 438     hNDiffCntSig2[i]->Draw("PSAME");\r
 439     TLegend* leg1=new TLegend(0.4,0.7,0.89,0.89);\r
 440     leg1->SetFillColor(0);\r
 441     TLegendEntry* ent1=leg1->AddEntry(hNDiffCntSig1[i],"From Counting1","PL");\r
 442     ent1->SetTextColor(hNDiffCntSig1[i]->GetMarkerColor());\r
 443     ent1=leg1->AddEntry(hNDiffCntSig2[i],"From Counting2","PL");\r
 444     ent1->SetTextColor(hNDiffCntSig2[i]->GetMarkerColor());\r
 445     leg1->Draw();\r
 446   }\r
 447 */\r
 448 \r
 449   TGraph* grReducedChiSquare0=new TGraph(nPtBins,ptlims,arrchisquare0);\r
 450   grReducedChiSquare0->SetName("grReducedChiSquare0");\r
 451   grReducedChiSquare0->SetTitle("Reduced Chi2;p_{t} (GeV/c);#tilde{#chi}^{2}");\r
 452   TGraph* grReducedChiSquare1=new TGraph(nPtBins,ptlims,arrchisquare1);\r
 453   grReducedChiSquare1->SetName("grReducedChiSquare1");\r
 454   grReducedChiSquare1->SetTitle("Reduced Chi2;p_{t} (GeV/c);#tilde{#chi}^{2}");\r
 455   TGraph* grReducedChiSquare2=new TGraph(nPtBins,ptlims,arrchisquare2);\r
 456   grReducedChiSquare2->SetName("grReducedChiSquare2");\r
 457   grReducedChiSquare2->SetTitle("Reduced Chi2;p_{t} (GeV/c);#tilde{#chi}^{2}");\r
 458   TGraph* grReducedChiSquare3=new TGraph(nPtBins,ptlims,arrchisquare3);\r
 459   grReducedChiSquare3->SetName("grReducedChiSquare3");\r
 460   grReducedChiSquare3->SetTitle("Reduced Chi2;p_{t} (GeV/c);#tilde{#chi}^{2}");\r
 461   TCanvas *cChi2=new TCanvas("cChi2","reduced chi square",600,600);\r
 462   cChi2->cd();\r
 463   grReducedChiSquare0->SetMarkerStyle(21);\r
 464   grReducedChiSquare0->Draw("AP");\r
 465   grReducedChiSquare1->SetMarkerStyle(22);\r
 466   grReducedChiSquare1->Draw("Psame");\r
 467   grReducedChiSquare2->SetMarkerStyle(23);\r
 468   grReducedChiSquare2->Draw("Psame");\r
 469   grReducedChiSquare3->SetMarkerStyle(24);\r
 470   grReducedChiSquare3->Draw("Psame");\r
 471 \r
 472   TCanvas* cbkgNormSigma=new TCanvas("cbkgNormSigma","Background normalized to sigma",400,600);\r
 473   cbkgNormSigma->cd();\r
 474   for(Int_t i=0; i<nMultbins; i++){\r
 475     hBackgroundNormSigma[i]->SetMarkerStyle(20);\r
 476     hBackgroundNormSigma[i]->GetXaxis()->SetTitle("Pt (GeV/c)");\r
 477     hBackgroundNormSigma[i]->GetYaxis()->SetTitle("Background #times 3 #times 0.012/ (3 #times #sigma)");\r
 478     hBackgroundNormSigma[i]->SetMarkerColor(2*i);\r
 479     if(i==5) hBackgroundNormSigma[i]->SetMarkerColor(2*i-3);\r
 480     if(i==0) { \r
 481       hBackgroundNormSigma[i]->SetMarkerColor(kBlack);\r
 482       hBackgroundNormSigma[i]->Draw("PE");\r
 483     }\r
 484     else  hBackgroundNormSigma[i]->Draw("Psame");\r
 485   }\r
 486 \r
 487 \r
 488   TString partname="Both";\r
 489   if(optPartAntiPart==kParticleOnly) {\r
 490     if(analysisType==kD0toKpi) partname="D0";\r
 491     if(analysisType==kDplusKpipi) partname="Dplus";\r
 492   }\r
 493   if(optPartAntiPart==kAntiParticleOnly) {\r
 494     if(analysisType==kD0toKpi) partname="D0bar";\r
 495     if(analysisType==kDplusKpipi) partname="Dminus";\r
 496   }\r
 497 \r
 498   TString outfilename = Form("RawYield_Mult_%s_%s",partname.Data(),CutsType);\r
 499   if(fixPeakSigma) outfilename += "_SigmaFixed";\r
 500   if(typeb==0) outfilename += "_Expo.root";\r
 501   else if(typeb==1) outfilename += "_Linear.root";\r
 502   else if(typeb==2) outfilename += "_Pol2.root";\r
 503 \r
 504   TFile* outf=new TFile(outfilename,"recreate");\r
 505   outf->cd(); \r
 506   hNormalization->Write();\r
 507   for(Int_t j=0; j<massBin; j++) hmass[j]->Write();\r
 508   for(Int_t j=0; j<nMultbins; j++){\r
 509     hMass[j]->Write();\r
 510     hSigma[j]->Write();\r
 511     hCntSig1[j]->Write();\r
 512     hCntSig2[j]->Write();\r
 513     hNDiffCntSig1[j]->Write();\r
 514     hNDiffCntSig2[j]->Write();\r
 515     hSignal[j]->Write();\r
 516     hRelErrSig[j]->Write();\r
 517     hInvSignif[j]->Write();\r
 518     hBackground[j]->Write();\r
 519     hBackgroundNormSigma[j]->Write();\r
 520     hSignificance[j]->Write();\r
 521   }\r
 522   grReducedChiSquare0->Write();\r
 523   grReducedChiSquare1->Write();\r
 524   grReducedChiSquare2->Write();\r
 525   grReducedChiSquare3->Write();\r
 526   //  hPtMass->Write();\r
 527   outf->Close();\r
 528   \r
 529 }\r
 530 \r
 531 //_____________________________________________________________________________________________\r
 532 Bool_t LoadD0toKpiHistos(TObjArray* listFiles, TH3F** hPtMassMult, TH2F** hNtrZvtx, TH2F** hNtrZvtxCorr, AliNormalizationCounter *counter, const char *CutsType, Int_t Option)\r
 533 {\r
 534   Int_t nFiles=listFiles->GetEntries();\r
 535   TList **hlist=new TList*[nFiles];\r
 536   TList **hlistNorm=new TList*[nFiles];\r
 537   AliRDHFCutsD0toKpi** dcuts=new AliRDHFCutsD0toKpi*[nFiles];\r
 538 \r
 539   Int_t nReadFiles=0;\r
 540   for(Int_t iFile=0; iFile<nFiles; iFile++){\r
 541     TString fName=((TObjString*)listFiles->At(iFile))->GetString();    \r
 542     TFile *f=TFile::Open(fName.Data());\r
 543     if(!f){\r
 544       printf("ERROR: file %s does not exist\n",fName.Data());\r
 545       continue;\r
 546     }\r
 547     printf("Open File %s\n",f->GetName());\r
 548 \r
 549     TString dirname="PWG3_D2H_DMult_D0";\r
 550     if(optPartAntiPart==kParticleOnly) dirname+="D0";\r
 551     else if(optPartAntiPart==kAntiParticleOnly) dirname+="D0bar";\r
 552     dirname += CutsType;\r
 553     TDirectory *dir = (TDirectory*)f->Get(dirname);\r
 554     if(!dir){\r
 555       printf("ERROR: directory %s not found in %s\n",dirname.Data(),fName.Data());\r
 556       continue;\r
 557     }\r
 558 \r
 559     TString listmassname="coutputD0";\r
 560     if(optPartAntiPart==kParticleOnly) listmassname+="D0";\r
 561     else if(optPartAntiPart==kAntiParticleOnly) listmassname+="D0bar";\r
 562     listmassname += CutsType;\r
 563     printf("List mass name %s\n",listmassname.Data());\r
 564     hlist[nReadFiles]=(TList*)dir->Get(listmassname);\r
 565 \r
 566     TString listnorm="coutputNormD0";\r
 567     if(optPartAntiPart==kParticleOnly) listnorm+="D0";\r
 568     else if(optPartAntiPart==kAntiParticleOnly) listnorm+="D0bar";\r
 569     listnorm += CutsType;\r
 570     printf("List norm name %s\n",listnorm.Data());\r
 571     hlistNorm[nReadFiles]=(TList*)dir->Get(listnorm);\r
 572     //    AliNormalizationCounter *tmpcounter = (AliNormalizationCounter*)hlistNorm[nReadFiles]->FindObject("NormCounterCorrMult");\r
 573     //    counter->Add(tmpcounter);\r
 574     //    delete tmpcounter;\r
 575     //    counter = tmpcounter;\r
 576 \r
 577     TString cutsobjname="coutputCutsD0";\r
 578     if(optPartAntiPart==kParticleOnly) cutsobjname+="D0";\r
 579     else if(optPartAntiPart==kAntiParticleOnly) cutsobjname+="D0bar";\r
 580     cutsobjname += CutsType;\r
 581     printf("Cuts name %s\n",cutsobjname.Data());\r
 582     dcuts[nReadFiles]=(AliRDHFCutsD0toKpi*)dir->Get(cutsobjname);\r
 583     if(!dcuts[nReadFiles]) {\r
 584       printf("ERROR: Cut objects do not match\n");\r
 585       return kFALSE;\r
 586     }\r
 587     /*\r
 588     if(nReadFiles>0){\r
 589       Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]);\r
 590       if(!sameCuts){\r
 591         printf("ERROR: Cut objects do not match\n");\r
 592         return kFALSE;\r
 593       }\r
 594     }\r
 595     */\r
 596     nReadFiles++;\r
 597   }\r
 598   if(nReadFiles<nFiles){\r
 599     printf("WARNING: not all requested files have been found\n");\r
 600     if (nReadFiles==0) {\r
 601       printf("ERROR: Any file/dir found\n");\r
 602       return kFALSE;\r
 603     }\r
 604   }\r
 605   //  printf("Cuts type %s, particle/antipart %d\n",CutsType,optPartAntiPart);\r
 606 \r
 607   /*\r
 608   Int_t nPtBinsCuts=dcuts[0]->GetNPtBins();\r
 609   printf("Number of pt bins for cut object = %d\n",nPtBins);\r
 610   Float_t *ptlimsCuts=dcuts[0]->GetPtBinLimits();\r
 611   ptlimsCuts[nPtBinsCuts]=ptlimsCuts[nPtBinsCuts-1]+4.;\r
 612   */\r
 613 \r
 614   printf("Get the 3D histogram \n");\r
 615   const char *histoname="";\r
 616   if(optPartAntiPart==kParticleOnly) histoname= "hPtVsMassvsMultPart";\r
 617   else if(optPartAntiPart==kAntiParticleOnly) histoname="hPtVsMassvsMultAntiPart";\r
 618   else if(optPartAntiPart==kBoth) histoname="hPtVsMassvsMult";\r
 619   if(Option==kUnCorr) histoname="hPtVsMassvsMultUncorr";\r
 620   if(Option==kNoPid) histoname="hPtVsMassvsMultNoPid";\r
 621 \r
 622   TH3F * htemp;\r
 623   for(Int_t iFile=0; iFile<nReadFiles; iFile++){\r
 624     printf(" Looking for histo histMass %s for file %d\n",histoname,iFile);\r
 625     htemp=(TH3F*)hlist[iFile]->FindObject(Form("%s",histoname));\r
 626     //    cout << " htemp :"<<htemp<<endl;\r
 627     if(!hPtMassMult[0]){\r
 628       hPtMassMult[0]=new TH3F(*htemp);\r
 629     }else{\r
 630       hPtMassMult[0]->Add(htemp);\r
 631     }\r
 632     hNtrZvtx[iFile] = (TH2F*)hlist[iFile]->FindObject("hNtrVsZvtx");\r
 633     hNtrZvtxCorr[iFile] = (TH2F*)hlist[iFile]->FindObject("hNtrCorrVsZvtx");\r
 634   }\r
 635   \r
 636   //  cout<<" hPtMassMult:"<<hPtMassMult[0]<<endl;\r
 637 \r
 638   TString partname="Both";\r
 639   if(optPartAntiPart==kParticleOnly) partname="D0";\r
 640   if(optPartAntiPart==kAntiParticleOnly) partname="D0bar";\r
 641 \r
 642   TString outfilename = Form("RawYield%s_%s",partname.Data(),CutsType);\r
 643   if(fixPeakSigma) outfilename += "_SigmaFixed";\r
 644   if(typeb==0) outfilename += "_Expo.root";\r
 645   else if(typeb==1) outfilename += "_Linear.root";\r
 646   else if(typeb==2) outfilename += "_Pol2.root";\r
 647   TFile* outf=new TFile(outfilename,"recreate");\r
 648   outf->cd();\r
 649   dcuts[0]->Write();\r
 650   outf->Close();\r
 651   return kTRUE;\r
 652 \r
 653 }\r
 654 \r
 655 //_____________________________________________________________________________________________\r
 656 Bool_t LoadDplusHistos(TObjArray* listFiles, TH3F** hPtMassMult, TH2F** hNtrZvtx, TH2F** hNtrZvtxCorr, const char *CutsType, Int_t Option)\r
 657 {\r
 658 Int_t nFiles=listFiles->GetEntries();\r
 659   TList **hlist=new TList*[nFiles];\r
 660   TList **hlistNorm=new TList*[nFiles];\r
 661   AliRDHFCutsDplustoKpipi** dcuts=new AliRDHFCutsDplustoKpipi*[nFiles];\r
 662 \r
 663   Int_t nReadFiles=0;\r
 664   for(Int_t iFile=0; iFile<nFiles; iFile++){\r
 665     TString fName=((TObjString*)listFiles->At(iFile))->GetString();    \r
 666     TFile *f=TFile::Open(fName.Data());\r
 667     if(!f){\r
 668       printf("ERROR: file %s does not exist\n",fName.Data());\r
 669       continue;\r
 670     }\r
 671     printf("Open File %s\n",f->GetName());\r
 672  TDirectory *dir = (TDirectory*)f->Get(Form("PWG3_D2H_DMult_Dplus%s",suffix.Data()));\r
 673     // TDirectory *dir = (TDirectory*)f->Get("PWG3_D2H_DMult");\r
 674     if(!dir){\r
 675       printf("ERROR: directory PWG3_D2H_DMult not found in %s\n",fName.Data());\r
 676       continue;\r
 677     }\r
 678     hlist[nReadFiles]=(TList*)dir->Get(Form("coutputDplus%s",suffix.Data()));\r
 679     TList *listcut = (TList*)dir->Get(Form("coutputCutsDplus%s",suffix.Data()));\r
 680     TList *listNorm = (TList*)dir->Get(Form("coutputNormDplus%s",suffix.Data()));\r
 681     dcuts[nReadFiles]=(AliRDHFCutsDplustoKpipi*)listcut->At(0);\r
 682      if(nReadFiles>0){\r
 683       Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]);\r
 684       if(!sameCuts){\r
 685  \r
 686         printf("ERROR: Cut objects do not match\n");\r
 687         return kFALSE;\r
 688       }\r
 689      }\r
 690 \r
 691 \r
 692 \r
 693  nReadFiles++;\r
 694       }\r
 695   if(nReadFiles<nFiles){\r
 696     printf("WARNING: not all requested files have been found\n");\r
 697     if (nReadFiles==0) {\r
 698       printf("ERROR: Any file/dir found\n");\r
 699       return kFALSE;\r
 700     }\r
 701   }\r
 702 \r
 703 \r
 704  printf("Get the 3D histogram \n");\r
 705   const char *histoname="";\r
 706   if(optPartAntiPart==kParticleOnly) histoname= "hPtVsMassvsMultPart";\r
 707   else if(optPartAntiPart==kAntiParticleOnly) histoname="hPtVsMassvsMultAntiPart";\r
 708   else if(optPartAntiPart==kBoth) histoname="hPtVsMassvsMult";\r
 709   if(Option==kUnCorr) histoname="hPtVsMassvsMultUncorr";\r
 710   if(Option==kNoPid) histoname="hPtVsMassvsMultNoPid";\r
 711 \r
 712   TH3F * htemp;\r
 713   for(Int_t iFile=0; iFile<nReadFiles; iFile++){\r
 714     printf(" Looking for histo histMass %s for file %d\n",histoname,iFile);\r
 715     htemp=(TH3F*)hlist[iFile]->FindObject(Form("%s",histoname));\r
 716     //    cout << " htemp :"<<htemp<<endl;\r
 717     if(!hPtMassMult[0]){\r
 718       hPtMassMult[0]=new TH3F(*htemp);\r
 719     }else{\r
 720       hPtMassMult[0]->Add(htemp);\r
 721     }\r
 722     //  (TH2F*)hNtrZvtx[iFile] = (TH2F*)hlist[iFile]->FindObject("hNtrVsZvtx"); \r
 723     //(TH2F*)hNtrZvtxCorr[iFile] = (TH2F*)hlist[iFile]->FindObject("hNtrVsZvtxCorr");\r
 724   }\r
 725   \r
 726   //  cout<<" hPtMassMult:"<<hPtMassMult[0]<<endl;\r
 727 \r
 728   TString partname="Both";\r
 729   if(optPartAntiPart==kParticleOnly) partname="Dplus";\r
 730   if(optPartAntiPart==kAntiParticleOnly) partname="Dminus";\r
 731 \r
 732   TString outfilename = Form("RawYield%s_%s",partname.Data(),CutsType);\r
 733   if(fixPeakSigma) outfilename += "_SigmaFixed";\r
 734   if(typeb==0) outfilename += "_Expo.root";\r
 735   else if(typeb==1) outfilename += "_Linear.root";\r
 736   else if(typeb==2) outfilename += "_Pol2.root";\r
 737   TFile* outf=new TFile(outfilename,"recreate");\r
 738   outf->cd();\r
 739   dcuts[0]->Write();\r
 740   outf->Close();\r
 741   return kTRUE;\r
 742   \r
 743 }\r
 744 \r
 745 //_____________________________________________________________________________________________\r
 746 TH1F* RebinHisto(TH1F* hOrig, Int_t reb, Int_t firstUse){\r
 747   // Rebin histogram, from bin firstUse to lastUse\r
 748   // Use all bins if firstUse=-1\r
 749   \r
 750   Int_t nBinOrig=hOrig->GetNbinsX();\r
 751   Int_t firstBinOrig=1;\r
 752   Int_t lastBinOrig=nBinOrig;\r
 753   Int_t nBinOrigUsed=nBinOrig;\r
 754   Int_t nBinFinal=nBinOrig/reb;\r
 755   if(firstUse>=1){\r
 756     firstBinOrig=firstUse;\r
 757     nBinFinal=(nBinOrig-firstUse+1)/reb;\r
 758     nBinOrigUsed=nBinFinal*reb;\r
 759     lastBinOrig=firstBinOrig+nBinOrigUsed-1;\r
 760   }else{\r
 761     Int_t exc=nBinOrigUsed%reb;\r
 762     if(exc!=0){\r
 763       nBinOrigUsed-=exc;\r
 764       firstBinOrig+=exc/2;\r
 765       lastBinOrig=firstBinOrig+nBinOrigUsed-1;\r
 766     }\r
 767   }\r
 768   \r
 769   printf("Rebin from %d bins to %d bins -- Used bins=%d in range %d-%d\n",nBinOrig,nBinFinal,nBinOrigUsed,firstBinOrig,lastBinOrig);\r
 770   Float_t lowLim=hOrig->GetXaxis()->GetBinLowEdge(firstBinOrig);\r
 771   Float_t hiLim=hOrig->GetXaxis()->GetBinUpEdge(lastBinOrig);\r
 772   TH1F* hRebin=new TH1F(Form("%s-rebin",hOrig->GetName()),hOrig->GetTitle(),nBinFinal,lowLim,hiLim);\r
 773   Int_t lastSummed=firstBinOrig-1;\r
 774   for(Int_t iBin=1;iBin<=nBinFinal; iBin++){\r
 775     Float_t sum=0.;\r
 776     for(Int_t iOrigBin=0;iOrigBin<reb;iOrigBin++){\r
 777       sum+=hOrig->GetBinContent(lastSummed+1);\r
 778       lastSummed++;\r
 779     }\r
 780     hRebin->SetBinContent(iBin,sum);\r
 781   }\r
 782   return hRebin;\r
 783 }\r
 784 \r
 785 //_____________________________________________________________________________________________\r
 786 Bool_t CheckNtrVsZvtx(TH2F** hNtrackVsVtxZ, TH2F** hNtrackVsVtxZCorr, Int_t nFiles)\r
 787 {\r
 788 \r
 789   TCanvas *cNtrVsZvtx = new TCanvas("cNtrVsZvtx","Ntr Vs Zvtx");\r
 790   cNtrVsZvtx->Divide(2,2);\r
 791   for(Int_t i=0; i<nFiles; i++){\r
 792     cNtrVsZvtx->cd(i+1);\r
 793     //    hNtrackVsVtxZ[i]->Fit("pol4");\r
 794     hNtrackVsVtxZ[i]->Draw("colz");\r
 795     cNtrVsZvtx->Update();\r
 796   }\r
 797 \r
 798   TCanvas *cNtrVsZvtxCorr = new TCanvas("cNtrVsZvtxCorr","Ntr Vs Zvtx Corr");\r
 799   cNtrVsZvtxCorr->Divide(2,2);\r
 800   for(Int_t i=0; i<nFiles; i++){\r
 801     cNtrVsZvtxCorr->cd(i+1);\r
 802     //    hNtrackVsVtxZCorr[i]->Fit("pol4");\r
 803     hNtrackVsVtxZCorr[i]->Draw("colz");\r
 804   }\r
 805 \r
 806   TH1F *hNtrAxis = (TH1F*)hNtrackVsVtxZ[0]->ProjectionY("hNtrAxis");\r
 807   TH1F *hZvtx[nMultbins];\r
 808   Int_t firstbin=0, lastbin=0;\r
 809   TCanvas *cZvtx = new TCanvas("cZvtx","Zvtx projections");\r
 810   cZvtx->Divide(2,2);\r
 811   for(Int_t i=0; i<nFiles; i++){\r
 812     cZvtx->cd(i+1);\r
 813     firstbin = hNtrAxis->FindBin( multlims[i] );\r
 814     lastbin = hNtrAxis->FindBin( multlims[i+1] ) -1;\r
 815     hZvtx[i] = (TH1F*)hNtrackVsVtxZ[i]->ProjectionX(Form("hZvtx_%d",i),firstbin,lastbin);\r
 816     hZvtx[i]->Draw();\r
 817   }\r
 818   TH1F *hZvtxCorr[nMultbins]; \r
 819   TCanvas *cZvtxCorr = new TCanvas("cZvtxCorr","Zvtx projections Corr");\r
 820   cZvtxCorr->Divide(2,2);\r
 821   for(Int_t i=0; i<nFiles; i++){\r
 822     cZvtxCorr->cd(i+1);\r
 823     firstbin = hNtrAxis->FindBin( multlims[i] );\r
 824     lastbin = hNtrAxis->FindBin( multlims[i+1] ) -1;\r
 825     hZvtxCorr[i] = (TH1F*)hNtrackVsVtxZCorr[i]->ProjectionX(Form("hZvtxCorr_%d",i),firstbin,lastbin);\r
 826     hZvtxCorr[i]->Draw();\r
 827   }\r
 828 \r
 829   return kTRUE;\r
 830 }\r