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