PWGHF/vertexingHF/macros/ReadDvsMultiplicity.C

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