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