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,-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])-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])-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       Double_t ry=fitter[iBin]->GetRawYield();\r
 278       Double_t ery=fitter[iBin]->GetRawYieldError();\r
 279       myCanvas[massBin] = new TCanvas(Form("myCanvas_%d_%d",j,iBin),Form("Invariant mass mult bin %d, pt bin %d",j,iBin));\r
 280       fitter[massBin]->DrawHere(gPad);\r
 281     \r
 282       Float_t cntSig1=0.;\r
 283       Float_t cntSig2=0.;\r
 284       Float_t cntErr=0.;\r
 285       for(Int_t iMB=minBinSum; iMB<=maxBinSum; iMB++){\r
 286         Float_t bkg1=fB1 ? fB1->Eval(hmass[massBin]->GetBinCenter(iMB))/rebin[iBin] : 0;\r
 287         Float_t bkg2=fB2 ? fB2->Eval(hmass[massBin]->GetBinCenter(iMB))/rebin[iBin] : 0;\r
 288         cntSig1+=(hmass[massBin]->GetBinContent(iMB)-bkg1);\r
 289         cntSig2+=(hmass[massBin]->GetBinContent(iMB)-bkg2);\r
 290         cntErr+=(hmass[massBin]->GetBinContent(iMB));\r
 291       }\r
 292       hCntSig1[j]->SetBinContent(iBin+1,cntSig1);\r
 293       hCntSig1[j]->SetBinError(iBin+1,TMath::Sqrt(cntErr));\r
 294       hNDiffCntSig1[j]->SetBinContent(iBin+1,(s-cntSig1)/s);\r
 295       hNDiffCntSig1[j]->SetBinError(iBin+1,TMath::Sqrt(cntErr)/s);\r
 296       hCntSig2[j]->SetBinContent(iBin+1,cntSig2);\r
 297       hNDiffCntSig2[j]->SetBinContent(iBin+1,(s-cntSig2)/s);\r
 298       hNDiffCntSig2[j]->SetBinError(iBin+1,TMath::Sqrt(cntErr)/s);\r
 299       hCntSig2[j]->SetBinError(iBin+1,TMath::Sqrt(cntErr));\r
 300       hSignal[j]->SetBinContent(iBin+1,ry);\r
 301       hSignal[j]->SetBinError(iBin+1,ery);\r
 302       hRelErrSig[j]->SetBinContent(iBin+1,errs/s);\r
 303       hInvSignif[j]->SetBinContent(iBin+1,1/sig);\r
 304       hInvSignif[j]->SetBinError(iBin+1,errsig/(sig*sig));\r
 305       hBackground[j]->SetBinContent(iBin+1,b); //consider sigma\r
 306       hBackground[j]->SetBinError(iBin+1,errb);\r
 307       hBackgroundNormSigma[j]->SetBinContent(iBin+1,b/(3*fitter[iBin]->GetSigma())*(3*0.012)); //consider sigma\r
 308       hBackgroundNormSigma[j]->SetBinError(iBin+1,errb);\r
 309       hSignificance[j]->SetBinContent(iBin+1,sig);\r
 310       hSignificance[j]->SetBinError(iBin+1,errsig);\r
 311       hMass[j]->SetBinContent(iBin+1,mass);\r
 312       hMass[j]->SetBinError(iBin+1,0.0001);\r
 313       hSigma[j]->SetBinContent(iBin+1,sigma);\r
 314       hSigma[j]->SetBinError(iBin+1,0.0001);\r
 315 \r
 316       massBin++;\r
 317     }// end loop on pt bins\r
 318 \r
 319     canvas[j]->Update();\r
 320     canvas[j]->SaveAs(Form("hMass%s_%d_%d.eps",CutsType,typeb,j));\r
 321     \r
 322   }// end loop on multiplicity bins\r
 323 \r
 324 \r
 325   TCanvas *cpar=new TCanvas("cpar","Fit params",1200,600);\r
 326   cpar->Divide(2,1);\r
 327   cpar->cd(1);\r
 328   for(Int_t imult=0; imult<nMultbins; imult++) {\r
 329     hMass[imult]->SetMarkerStyle(20);\r
 330     hMass[imult]->GetXaxis()->SetTitle("Pt (GeV/c)");\r
 331     hMass[imult]->GetYaxis()->SetTitle("Mass (GeV/c^{2})");\r
 332     hMass[imult]->SetMarkerColor(2*imult);\r
 333     if(imult==5) hMass[imult]->SetMarkerColor(2*imult-3);\r
 334     if(imult==0) {\r
 335       hMass[imult]->SetMarkerColor(kBlack);\r
 336       hMass[imult]->Draw("PE");\r
 337     }\r
 338     else  hMass[imult]->Draw("PEsame");\r
 339   }\r
 340   cpar->cd(2);\r
 341   for(Int_t imult=0; imult<nMultbins; imult++) {\r
 342     hSigma[imult]->SetMarkerStyle(20);\r
 343     //  hSigma[0]->Draw("PE");\r
 344     hSigma[imult]->GetXaxis()->SetTitle("Pt (GeV/c)");\r
 345     hSigma[imult]->GetYaxis()->SetTitle("Sigma (GeV/c^{2})");\r
 346     hSigma[imult]->SetMarkerColor(2*imult);\r
 347     if(imult==5) hSigma[imult]->SetMarkerColor(2*imult-3);\r
 348     if(imult==0) {\r
 349       hSigma[imult]->SetMarkerColor(kBlack);\r
 350       hSigma[imult]->Draw("PE");\r
 351     }\r
 352     else  hSigma[imult]->Draw("PEsame");\r
 353   }\r
 354   cpar->Update();\r
 355 \r
 356   /*\r
 357   TCanvas** csig;//= new TCanvas*[nMultbins];\r
 358   TCanvas** cDiffS;//=new TCanvas*[nMultbins];\r
 359   for(Int_t i=0; i<nMultbins; i++){\r
 360     csig[i] =new TCanvas(Form("csig_%d",i),Form("Results, multiplicity bin %d",i),1200,600);\r
 361     csig[i]->Divide(3,1);\r
 362     csig[i]->cd(1);\r
 363     hSignal[i]->SetMarkerStyle(20);\r
 364     hSignal[i]->SetMarkerColor(4);\r
 365     hSignal[i]->SetLineColor(4);\r
 366     hSignal[i]->GetXaxis()->SetTitle("Pt (GeV/c)");\r
 367     hSignal[i]->GetYaxis()->SetTitle("Signal");\r
 368     hSignal[i]->Draw("P");\r
 369     hCntSig1[i]->SetMarkerStyle(26);\r
 370     hCntSig1[i]->SetMarkerColor(2);\r
 371     hCntSig1[i]->SetLineColor(2);\r
 372     hCntSig1[i]->Draw("PSAME");\r
 373     hCntSig2[i]->SetMarkerStyle(29);\r
 374     hCntSig2[i]->SetMarkerColor(kGray+1);\r
 375     hCntSig2[i]->SetLineColor(kGray+1);\r
 376     hCntSig2[i]->Draw("PSAME");\r
 377     TLegend* leg=new TLegend(0.4,0.7,0.89,0.89);\r
 378     leg->SetFillColor(0);\r
 379     TLegendEntry* ent=leg->AddEntry(hSignal[i],"From Fit","PL");\r
 380     ent->SetTextColor(hSignal[i]->GetMarkerColor());\r
 381     ent=leg->AddEntry(hCntSig1[i],"From Counting1","PL");\r
 382     ent->SetTextColor(hCntSig1[i]->GetMarkerColor());\r
 383     ent=leg->AddEntry(hCntSig2[i],"From Counting2","PL");\r
 384     ent->SetTextColor(hCntSig2[i]->GetMarkerColor());\r
 385     leg->Draw();\r
 386     csig[i]->cd(2);\r
 387     hBackground[i]->SetMarkerStyle(20);\r
 388     hBackground[i]->Draw("P");\r
 389     hBackground[i]->GetXaxis()->SetTitle("Pt (GeV/c)");\r
 390     hBackground[i]->GetYaxis()->SetTitle("Background");\r
 391     csig[i]->cd(3);\r
 392     hSignificance[i]->SetMarkerStyle(20);\r
 393     hSignificance[i]->Draw("P");\r
 394     hSignificance[i]->GetXaxis()->SetTitle("Pt (GeV/c)");\r
 395     hSignificance[i]->GetYaxis()->SetTitle("Significance");\r
 396     \r
 397     cDiffS[i]=new TCanvas(Form("cDiffS_%d",i),Form("Difference, multiplicity bin %d",i),1200,600);\r
 398     cDiffS[i]->Divide(2,1);\r
 399     cDiffS[i]->cd(1);\r
 400     hRelErrSig[i]->SetMarkerStyle(20); //fullcircle\r
 401     hRelErrSig[i]->SetTitleOffset(1.2);  \r
 402     hRelErrSig[i]->SetTitle("Rel Error from Fit;p_{t} (GeV/c);Signal Relative Error");\r
 403     hRelErrSig[i]->Draw("P");\r
 404     hInvSignif[i]->SetMarkerStyle(21); //fullsquare\r
 405     hInvSignif[i]->SetMarkerColor(kMagenta+1);\r
 406     hInvSignif[i]->SetLineColor(kMagenta+1);\r
 407     hInvSignif[i]->Draw("PSAMES");\r
 408     TLegend* leg2=new TLegend(0.4,0.7,0.89,0.89);\r
 409     leg2->SetFillColor(0);\r
 410     TLegendEntry* ent2=leg2->AddEntry(hRelErrSig[i],"From Fit","P");\r
 411     ent2->SetTextColor(hRelErrSig[i]->GetMarkerColor());\r
 412     ent2=leg2->AddEntry(hInvSignif[i],"1/Significance","PL");\r
 413     ent2->SetTextColor(hInvSignif[i]->GetMarkerColor());\r
 414     leg2->Draw();\r
 415 \r
 416     cDiffS[i]->cd(2);\r
 417     hNDiffCntSig1[i]->SetMarkerStyle(26);\r
 418     hNDiffCntSig1[i]->SetMarkerColor(2);\r
 419     hNDiffCntSig1[i]->SetLineColor(2);\r
 420     hNDiffCntSig1[i]->SetTitle("Cmp Fit-Count;p_{t} (GeV/c);(S_{fit}-S_{count})/S_{fit}");\r
 421     hNDiffCntSig1[i]->Draw("P");\r
 422     hNDiffCntSig2[i]->SetMarkerStyle(29);\r
 423     hNDiffCntSig2[i]->SetMarkerColor(kGray+1);\r
 424     hNDiffCntSig2[i]->SetLineColor(kGray+1);\r
 425     hNDiffCntSig2[i]->Draw("PSAME");\r
 426     TLegend* leg1=new TLegend(0.4,0.7,0.89,0.89);\r
 427     leg1->SetFillColor(0);\r
 428     TLegendEntry* ent1=leg1->AddEntry(hNDiffCntSig1[i],"From Counting1","PL");\r
 429     ent1->SetTextColor(hNDiffCntSig1[i]->GetMarkerColor());\r
 430     ent1=leg1->AddEntry(hNDiffCntSig2[i],"From Counting2","PL");\r
 431     ent1->SetTextColor(hNDiffCntSig2[i]->GetMarkerColor());\r
 432     leg1->Draw();\r
 433   }\r
 434 */\r
 435 \r
 436   TGraph* grReducedChiSquare0=new TGraph(nPtBins,ptlims,arrchisquare0);\r
 437   grReducedChiSquare0->SetName("grReducedChiSquare0");\r
 438   grReducedChiSquare0->SetTitle("Reduced Chi2;p_{t} (GeV/c);#tilde{#chi}^{2}");\r
 439   TGraph* grReducedChiSquare1=new TGraph(nPtBins,ptlims,arrchisquare1);\r
 440   grReducedChiSquare1->SetName("grReducedChiSquare1");\r
 441   grReducedChiSquare1->SetTitle("Reduced Chi2;p_{t} (GeV/c);#tilde{#chi}^{2}");\r
 442   TGraph* grReducedChiSquare2=new TGraph(nPtBins,ptlims,arrchisquare2);\r
 443   grReducedChiSquare2->SetName("grReducedChiSquare2");\r
 444   grReducedChiSquare2->SetTitle("Reduced Chi2;p_{t} (GeV/c);#tilde{#chi}^{2}");\r
 445   TGraph* grReducedChiSquare3=new TGraph(nPtBins,ptlims,arrchisquare3);\r
 446   grReducedChiSquare3->SetName("grReducedChiSquare3");\r
 447   grReducedChiSquare3->SetTitle("Reduced Chi2;p_{t} (GeV/c);#tilde{#chi}^{2}");\r
 448   TCanvas *cChi2=new TCanvas("cChi2","reduced chi square",600,600);\r
 449   cChi2->cd();\r
 450   grReducedChiSquare0->SetMarkerStyle(21);\r
 451   grReducedChiSquare0->Draw("AP");\r
 452   grReducedChiSquare1->SetMarkerStyle(22);\r
 453   grReducedChiSquare1->Draw("Psame");\r
 454   grReducedChiSquare2->SetMarkerStyle(23);\r
 455   grReducedChiSquare2->Draw("Psame");\r
 456   grReducedChiSquare3->SetMarkerStyle(24);\r
 457   grReducedChiSquare3->Draw("Psame");\r
 458 \r
 459   TCanvas* cbkgNormSigma=new TCanvas("cbkgNormSigma","Background normalized to sigma",400,600);\r
 460   cbkgNormSigma->cd();\r
 461   for(Int_t i=0; i<nMultbins; i++){\r
 462     hBackgroundNormSigma[i]->SetMarkerStyle(20);\r
 463     hBackgroundNormSigma[i]->GetXaxis()->SetTitle("Pt (GeV/c)");\r
 464     hBackgroundNormSigma[i]->GetYaxis()->SetTitle("Background #times 3 #times 0.012/ (3 #times #sigma)");\r
 465     hBackgroundNormSigma[i]->SetMarkerColor(2*i);\r
 466     if(i==5) hBackgroundNormSigma[i]->SetMarkerColor(2*i-3);\r
 467     if(i==0) { \r
 468       hBackgroundNormSigma[i]->SetMarkerColor(kBlack);\r
 469       hBackgroundNormSigma[i]->Draw("PE");\r
 470     }\r
 471     else  hBackgroundNormSigma[i]->Draw("Psame");\r
 472   }\r
 473 \r
 474 \r
 475   TString partname="Both";\r
 476   if(optPartAntiPart==kParticleOnly) {\r
 477     if(analysisType==kD0toKpi) partname="D0";\r
 478     if(analysisType==kDplusKpipi) partname="Dplus";\r
 479   }\r
 480   if(optPartAntiPart==kAntiParticleOnly) {\r
 481     if(analysisType==kD0toKpi) partname="D0bar";\r
 482     if(analysisType==kDplusKpipi) partname="Dminus";\r
 483   }\r
 484 \r
 485   TString outfilename = Form("RawYield_MultInt_%s_%s",partname.Data(),CutsType);\r
 486   if(fixPeakSigma) outfilename += "_SigmaFixed";\r
 487   if(typeb==0) outfilename += "_Expo.root";\r
 488   else if(typeb==1) outfilename += "_Linear.root";\r
 489   else if(typeb==2) outfilename += "_Pol2.root";\r
 490 \r
 491   TFile* outf=new TFile(outfilename,"update");\r
 492   outf->cd(); \r
 493   for(Int_t j=0; j<massBin; j++) hmass[j]->Write();\r
 494   for(Int_t j=0; j<nMultbins; j++){\r
 495     hMass[j]->Write();\r
 496     hSigma[j]->Write();\r
 497     hCntSig1[j]->Write();\r
 498     hCntSig2[j]->Write();\r
 499     hNDiffCntSig1[j]->Write();\r
 500     hNDiffCntSig2[j]->Write();\r
 501     hSignal[j]->Write();\r
 502     hRelErrSig[j]->Write();\r
 503     hInvSignif[j]->Write();\r
 504     hBackground[j]->Write();\r
 505     hBackgroundNormSigma[j]->Write();\r
 506     hSignificance[j]->Write();\r
 507   }\r
 508   grReducedChiSquare0->Write();\r
 509   grReducedChiSquare1->Write();\r
 510   grReducedChiSquare2->Write();\r
 511   grReducedChiSquare3->Write();\r
 512   //  hPtMass->Write();\r
 513   outf->Close();\r
 514   \r
 515 }\r
 516 \r
 517 //_____________________________________________________________________________________________\r
 518 Bool_t LoadD0toKpiHistos(TObjArray* listFiles, TH3F** hPtMassMult, TH2F** hNtrZvtx, TH2F** hNtrZvtxCorr, const char *CutsType, Int_t Option)\r
 519 {\r
 520   Int_t nFiles=listFiles->GetEntries();\r
 521   TList **hlist=new TList*[nFiles];\r
 522   TList **hlistNorm=new TList*[nFiles];\r
 523   AliRDHFCutsD0toKpi** dcuts=new AliRDHFCutsD0toKpi*[nFiles];\r
 524 \r
 525   Int_t nReadFiles=0;\r
 526   for(Int_t iFile=0; iFile<nFiles; iFile++){\r
 527     TString fName=((TObjString*)listFiles->At(iFile))->GetString();    \r
 528     TFile *f=TFile::Open(fName.Data());\r
 529     if(!f){\r
 530       printf("ERROR: file %s does not exist\n",fName.Data());\r
 531       continue;\r
 532     }\r
 533     printf("Open File %s\n",f->GetName());\r
 534 \r
 535     TString dirname="PWG3_D2H_DMult_D0";\r
 536     if(optPartAntiPart==kParticleOnly) dirname+="D0";\r
 537     else if(optPartAntiPart==kAntiParticleOnly) dirname+="D0bar";\r
 538     dirname += CutsType;\r
 539     TDirectory *dir = (TDirectory*)f->Get(dirname);\r
 540     if(!dir){\r
 541       printf("ERROR: directory %s not found in %s\n",dirname.Data(),fName.Data());\r
 542       continue;\r
 543     }\r
 544 \r
 545     TString listmassname="coutputD0";\r
 546     if(optPartAntiPart==kParticleOnly) listmassname+="D0";\r
 547     else if(optPartAntiPart==kAntiParticleOnly) listmassname+="D0bar";\r
 548     listmassname += CutsType;\r
 549     printf("List mass name %s\n",listmassname.Data());\r
 550     hlist[nReadFiles]=(TList*)dir->Get(listmassname);\r
 551 \r
 552     TString listnorm="coutputNormD0";\r
 553     if(optPartAntiPart==kParticleOnly) listnorm+="D0";\r
 554     else if(optPartAntiPart==kAntiParticleOnly) listnorm+="D0bar";\r
 555     listnorm += CutsType;\r
 556     printf("List norm name %s\n",listnorm.Data());\r
 557     hlistNorm[nReadFiles]=(TList*)dir->Get(listnorm);\r
 558 \r
 559     TString cutsobjname="coutputCutsD0";\r
 560     if(optPartAntiPart==kParticleOnly) cutsobjname+="D0";\r
 561     else if(optPartAntiPart==kAntiParticleOnly) cutsobjname+="D0bar";\r
 562     cutsobjname += CutsType;\r
 563     printf("Cuts name %s\n",cutsobjname.Data());\r
 564     dcuts[nReadFiles]=(AliRDHFCutsD0toKpi*)dir->Get(cutsobjname);\r
 565     if(!dcuts[nReadFiles]) {\r
 566       printf("ERROR: Cut objects do not match\n");\r
 567       return kFALSE;\r
 568     }\r
 569     /*\r
 570     if(nReadFiles>0){\r
 571       Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]);\r
 572       if(!sameCuts){\r
 573         printf("ERROR: Cut objects do not match\n");\r
 574         return kFALSE;\r
 575       }\r
 576     }\r
 577     */\r
 578     nReadFiles++;\r
 579   }\r
 580   if(nReadFiles<nFiles){\r
 581     printf("WARNING: not all requested files have been found\n");\r
 582     if (nReadFiles==0) {\r
 583       printf("ERROR: Any file/dir found\n");\r
 584       return kFALSE;\r
 585     }\r
 586   }\r
 587   //  printf("Cuts type %s, particle/antipart %d\n",CutsType,optPartAntiPart);\r
 588 \r
 589   /*\r
 590   Int_t nPtBinsCuts=dcuts[0]->GetNPtBins();\r
 591   printf("Number of pt bins for cut object = %d\n",nPtBins);\r
 592   Float_t *ptlimsCuts=dcuts[0]->GetPtBinLimits();\r
 593   ptlimsCuts[nPtBinsCuts]=ptlimsCuts[nPtBinsCuts-1]+4.;\r
 594   */\r
 595 \r
 596   printf("Get the 3D histogram \n");\r
 597   const char *histoname="";\r
 598   if(optPartAntiPart==kParticleOnly) histoname= "hPtVsMassvsMultPart";\r
 599   else if(optPartAntiPart==kAntiParticleOnly) histoname="hPtVsMassvsMultAntiPart";\r
 600   else if(optPartAntiPart==kBoth) histoname="hPtVsMassvsMult";\r
 601   if(Option==kUnCorr) histoname="hPtVsMassvsMultUncorr";\r
 602   if(Option==kNoPid) histoname="hPtVsMassvsMultNoPid";\r
 603 \r
 604   TH3F * htemp;\r
 605   for(Int_t iFile=0; iFile<nReadFiles; iFile++){\r
 606     printf(" Looking for histo histMass %s for file %d\n",histoname,iFile);\r
 607     htemp=(TH3F*)hlist[iFile]->FindObject(Form("%s",histoname));\r
 608     //    cout << " htemp :"<<htemp<<endl;\r
 609     if(!hPtMassMult[0]){\r
 610       hPtMassMult[0]=new TH3F(*htemp);\r
 611     }else{\r
 612       hPtMassMult[0]->Add(htemp);\r
 613     }\r
 614     hNtrZvtx[iFile] = (TH2F*)hlist[iFile]->FindObject("hNtrVsZvtx");\r
 615     hNtrZvtxCorr[iFile] = (TH2F*)hlist[iFile]->FindObject("hNtrCorrVsZvtx");\r
 616   }\r
 617   \r
 618   //  cout<<" hPtMassMult:"<<hPtMassMult[0]<<endl;\r
 619 \r
 620   TString partname="Both";\r
 621   if(optPartAntiPart==kParticleOnly) partname="D0";\r
 622   if(optPartAntiPart==kAntiParticleOnly) partname="D0bar";\r
 623 \r
 624   TString outfilename = Form("RawYield%s_%s",partname.Data(),CutsType);\r
 625   if(fixPeakSigma) outfilename += "_SigmaFixed";\r
 626   if(typeb==0) outfilename += "_Expo.root";\r
 627   else if(typeb==1) outfilename += "_Linear.root";\r
 628   else if(typeb==2) outfilename += "_Pol2.root";\r
 629   TFile* outf=new TFile(outfilename,"recreate");\r
 630   outf->cd();\r
 631   dcuts[0]->Write();\r
 632   outf->Close();\r
 633   return kTRUE;\r
 634 \r
 635 }\r
 636 \r
 637 //_____________________________________________________________________________________________\r
 638 Bool_t LoadDplusHistos(TObjArray* listFiles, TH3F** hPtMassMult, TH2F** hNtrZvtx, TH2F** hNtrZvtxCorr, const char *CutsType, Int_t Option)\r
 639 {\r
 640 Int_t nFiles=listFiles->GetEntries();\r
 641   TList **hlist=new TList*[nFiles];\r
 642   TList **hlistNorm=new TList*[nFiles];\r
 643   AliRDHFCutsDplustoKpipi** dcuts=new AliRDHFCutsDplustoKpipi*[nFiles];\r
 644 \r
 645   Int_t nReadFiles=0;\r
 646   for(Int_t iFile=0; iFile<nFiles; iFile++){\r
 647     TString fName=((TObjString*)listFiles->At(iFile))->GetString();    \r
 648     TFile *f=TFile::Open(fName.Data());\r
 649     if(!f){\r
 650       printf("ERROR: file %s does not exist\n",fName.Data());\r
 651       continue;\r
 652     }\r
 653     printf("Open File %s\n",f->GetName());\r
 654  TDirectory *dir = (TDirectory*)f->Get(Form("PWG3_D2H_DMult_Dplus%s",suffix.Data()));\r
 655     // TDirectory *dir = (TDirectory*)f->Get("PWG3_D2H_DMult");\r
 656     if(!dir){\r
 657       printf("ERROR: directory PWG3_D2H_DMult not found in %s\n",fName.Data());\r
 658       continue;\r
 659     }\r
 660     hlist[nReadFiles]=(TList*)dir->Get(Form("coutputDplus%s",suffix.Data()));\r
 661     TList *listcut = (TList*)dir->Get(Form("coutputCutsDplus%s",suffix.Data()));\r
 662     TList *listNorm = (TList*)dir->Get(Form("coutputNormDplus%s",suffix.Data()));\r
 663     dcuts[nReadFiles]=(AliRDHFCutsDplustoKpipi*)listcut->At(0);\r
 664      if(nReadFiles>0){\r
 665       Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]);\r
 666       if(!sameCuts){\r
 667  \r
 668         printf("ERROR: Cut objects do not match\n");\r
 669         return kFALSE;\r
 670       }\r
 671      }\r
 672 \r
 673 \r
 674 \r
 675  nReadFiles++;\r
 676       }\r
 677   if(nReadFiles<nFiles){\r
 678     printf("WARNING: not all requested files have been found\n");\r
 679     if (nReadFiles==0) {\r
 680       printf("ERROR: Any file/dir found\n");\r
 681       return kFALSE;\r
 682     }\r
 683   }\r
 684 \r
 685 \r
 686  printf("Get the 3D histogram \n");\r
 687   const char *histoname="";\r
 688   if(optPartAntiPart==kParticleOnly) histoname= "hPtVsMassvsMultPart";\r
 689   else if(optPartAntiPart==kAntiParticleOnly) histoname="hPtVsMassvsMultAntiPart";\r
 690   else if(optPartAntiPart==kBoth) histoname="hPtVsMassvsMult";\r
 691   if(Option==kUnCorr) histoname="hPtVsMassvsMultUncorr";\r
 692   if(Option==kNoPid) histoname="hPtVsMassvsMultNoPid";\r
 693 \r
 694   TH3F * htemp;\r
 695   for(Int_t iFile=0; iFile<nReadFiles; iFile++){\r
 696     printf(" Looking for histo histMass %s for file %d\n",histoname,iFile);\r
 697     htemp=(TH3F*)hlist[iFile]->FindObject(Form("%s",histoname));\r
 698     //    cout << " htemp :"<<htemp<<endl;\r
 699     if(!hPtMassMult[0]){\r
 700       hPtMassMult[0]=new TH3F(*htemp);\r
 701     }else{\r
 702       hPtMassMult[0]->Add(htemp);\r
 703     }\r
 704     //  (TH2F*)hNtrZvtx[iFile] = (TH2F*)hlist[iFile]->FindObject("hNtrVsZvtx"); \r
 705     //(TH2F*)hNtrZvtxCorr[iFile] = (TH2F*)hlist[iFile]->FindObject("hNtrVsZvtxCorr");\r
 706   }\r
 707   \r
 708   //  cout<<" hPtMassMult:"<<hPtMassMult[0]<<endl;\r
 709 \r
 710   TString partname="Both";\r
 711   if(optPartAntiPart==kParticleOnly) partname="Dplus";\r
 712   if(optPartAntiPart==kAntiParticleOnly) partname="Dminus";\r
 713 \r
 714   TString outfilename = Form("RawYield%s_%s",partname.Data(),CutsType);\r
 715   if(fixPeakSigma) outfilename += "_SigmaFixed";\r
 716   if(typeb==0) outfilename += "_Expo.root";\r
 717   else if(typeb==1) outfilename += "_Linear.root";\r
 718   else if(typeb==2) outfilename += "_Pol2.root";\r
 719   TFile* outf=new TFile(outfilename,"recreate");\r
 720   outf->cd();\r
 721   dcuts[0]->Write();\r
 722   outf->Close();\r
 723   return kTRUE;\r
 724   \r
 725      }\r
 726 \r
 727 //_____________________________________________________________________________________________\r
 728 TH1F* RebinHisto(TH1F* hOrig, Int_t reb, Int_t firstUse){\r
 729   // Rebin histogram, from bin firstUse to lastUse\r
 730   // Use all bins if firstUse=-1\r
 731   \r
 732   Int_t nBinOrig=hOrig->GetNbinsX();\r
 733   Int_t firstBinOrig=1;\r
 734   Int_t lastBinOrig=nBinOrig;\r
 735   Int_t nBinOrigUsed=nBinOrig;\r
 736   Int_t nBinFinal=nBinOrig/reb;\r
 737   if(firstUse>=1){\r
 738     firstBinOrig=firstUse;\r
 739     nBinFinal=(nBinOrig-firstUse+1)/reb;\r
 740     nBinOrigUsed=nBinFinal*reb;\r
 741     lastBinOrig=firstBinOrig+nBinOrigUsed-1;\r
 742   }else{\r
 743     Int_t exc=nBinOrigUsed%reb;\r
 744     if(exc!=0){\r
 745       nBinOrigUsed-=exc;\r
 746       firstBinOrig+=exc/2;\r
 747       lastBinOrig=firstBinOrig+nBinOrigUsed-1;\r
 748     }\r
 749   }\r
 750   \r
 751   printf("Rebin from %d bins to %d bins -- Used bins=%d in range %d-%d\n",nBinOrig,nBinFinal,nBinOrigUsed,firstBinOrig,lastBinOrig);\r
 752   Float_t lowLim=hOrig->GetXaxis()->GetBinLowEdge(firstBinOrig);\r
 753   Float_t hiLim=hOrig->GetXaxis()->GetBinUpEdge(lastBinOrig);\r
 754   TH1F* hRebin=new TH1F(Form("%s-rebin",hOrig->GetName()),hOrig->GetTitle(),nBinFinal,lowLim,hiLim);\r
 755   Int_t lastSummed=firstBinOrig-1;\r
 756   for(Int_t iBin=1;iBin<=nBinFinal; iBin++){\r
 757     Float_t sum=0.;\r
 758     for(Int_t iOrigBin=0;iOrigBin<reb;iOrigBin++){\r
 759       sum+=hOrig->GetBinContent(lastSummed+1);\r
 760       lastSummed++;\r
 761     }\r
 762     hRebin->SetBinContent(iBin,sum);\r
 763   }\r
 764   return hRebin;\r
 765 }\r
 766 \r
 767 //_____________________________________________________________________________________________\r
 768 Bool_t CheckNtrVsZvtx(TH2F** hNtrackVsVtxZ, TH2F** hNtrackVsVtxZCorr, Int_t nFiles)\r
 769 {\r
 770 \r
 771   TCanvas *cNtrVsZvtx = new TCanvas("cNtrVsZvtx","Ntr Vs Zvtx");\r
 772   cNtrVsZvtx->Divide(2,2);\r
 773   for(Int_t i=0; i<nFiles; i++){\r
 774     cNtrVsZvtx->cd(i+1);\r
 775     //    hNtrackVsVtxZ[i]->Fit("pol4");\r
 776     hNtrackVsVtxZ[i]->Draw("colz");\r
 777     cNtrVsZvtx->Update();\r
 778   }\r
 779 \r
 780   TCanvas *cNtrVsZvtxCorr = new TCanvas("cNtrVsZvtxCorr","Ntr Vs Zvtx Corr");\r
 781   cNtrVsZvtxCorr->Divide(2,2);\r
 782   for(Int_t i=0; i<nFiles; i++){\r
 783     cNtrVsZvtxCorr->cd(i+1);\r
 784     //    hNtrackVsVtxZCorr[i]->Fit("pol4");\r
 785     hNtrackVsVtxZCorr[i]->Draw("colz");\r
 786   }\r
 787 \r
 788   TH1F *hNtrAxis = (TH1F*)hNtrackVsVtxZ[0]->ProjectionY("hNtrAxis");\r
 789   TH1F *hZvtx[nMultbins];\r
 790   Int_t firstbin=0, lastbin=0;\r
 791   TCanvas *cZvtx = new TCanvas("cZvtx","Zvtx projections");\r
 792   cZvtx->Divide(2,2);\r
 793   for(Int_t i=0; i<nFiles; i++){\r
 794     cZvtx->cd(i+1);\r
 795     firstbin = hNtrAxis->FindBin( multlims[i] );\r
 796     lastbin = hNtrAxis->FindBin( multlims[i+1] ) -1;\r
 797     hZvtx[i] = (TH1F*)hNtrackVsVtxZ[i]->ProjectionX(Form("hZvtx_%d",i),firstbin,lastbin);\r
 798     hZvtx[i]->Draw();\r
 799   }\r
 800   TH1F *hZvtxCorr[nMultbins]; \r
 801   TCanvas *cZvtxCorr = new TCanvas("cZvtxCorr","Zvtx projections Corr");\r
 802   cZvtxCorr->Divide(2,2);\r
 803   for(Int_t i=0; i<nFiles; i++){\r
 804     cZvtxCorr->cd(i+1);\r
 805     firstbin = hNtrAxis->FindBin( multlims[i] );\r
 806     lastbin = hNtrAxis->FindBin( multlims[i+1] ) -1;\r
 807     hZvtxCorr[i] = (TH1F*)hNtrackVsVtxZCorr[i]->ProjectionX(Form("hZvtxCorr_%d",i),firstbin,lastbin);\r
 808     hZvtxCorr[i]->Draw();\r
 809   }\r
 810 \r
 811   return kTRUE;\r
 812 }\r