PWG3/vertexingHF/macros/FitMassSpectra.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 <TH1F.h>
  10 #include <TH2F.h>
  11 #include <TH1D.h>
  12 #include <TF1.h>
  13 #include <TStyle.h>
  14 #include <TLegend.h>
  15 #include <TLegendEntry.h>
  16 #include <TDatabasePDG.h>
  17 #include <TGraph.h>
  18
  19 #include "AliAODRecoDecayHF.h"
  20 #include "AliRDHFCuts.h"
  21 #include "AliRDHFCutsDplustoKpipi.h"
  22 #include "AliRDHFCutsD0toKpi.h"
  23 #include "AliHFMassFitter.h"
  24 #endif
  25
  26
  27 // MACRO to perform fits to D meson invariant mass spectra
  28 // and store raw yields and cut object into a root output file
  29 // Origin: F. Prino (prino@to.infn.it)
  30 // D0: C. Bianchin (cbianchi@pd.infn.it)
  31
  32
  33
  34 //
  35 enum {kD0toKpi, kDplusKpipi};
  36 enum {kBoth, kParticleOnly, kAntiParticleOnly};
  37 enum {kExpo=0, kLinear, kPol2};
  38 enum {kGaus=0, kDoubleGaus};
  39
  40
  41 // Common variables: to be configured by the user
  42 const Int_t nPtBins=6;
  43 Double_t ptlims[nPtBins+1]={2.,3.,4.,5.,6.,8.,12.};
  44 Int_t rebin[nPtBins+1]={2,4,4,4,4,4,4};
  45
  46 //const Int_t nPtBins=7;//6;
  47 //Double_t ptlims[nPtBins+1]={1.,2.,3.,4.,5.,6.,8.,12.};
  48 //Int_t rebin[nPtBins+1]={8,6,10,10,10,10,10,10}; //for looser cuts
  49 //Int_t rebin[nPtBins+1]={10,10,10,10,10,10,10,10}; //for Chiara's cuts
  50 Int_t typeb=kExpo;
  51 Int_t types=kGaus;
  52 Int_t optPartAntiPart=kBoth;
  53 Int_t factor4refl=0;
  54 Float_t massRangeForCounting=0.05; // GeV
  55 TH2F* hPtMass=0x0;
  56
  57 //for D0only
  58 Bool_t cutsappliedondistr=kTRUE;
  59 const Int_t nsamples=2;//3;
  60 //Int_t nevents[nsamples]={8.5635859e+07 /*LHC10b+cpart*/,8.9700624e+07/*LHC10dpart*/};
  61 //Int_t nevents[nsamples]={9.0374946e+07 /*LHC10b+c*/,9.7593785e+07/*LHC10d*/};
  62 //Int_t nevents[nsamples]={1.1777812e+08 /*LHC10dnewTPCpid*/,9.7593785e+07/*LHC10d*/};
  63 //Int_t nevents[nsamples]={1.1777812e+08 /*LHC10dnewTPCpid*/,9.7593785e+07/*LHC10d*/,9.0374946e+07 /*LHC10b+c*/};
  64 Int_t nevents[nsamples]={1.18860695e+08 /*LHC10dnewTPCpid*/,9.0374946e+07 /*LHC10b+c*/};
  65 //Int_t nevents[nsamples]={1. /*LHC10dnewTPCpid*/,1 /*LHC10dnewTPCpidrescale*/};
  66 // Functions
  67
  68 Bool_t LoadDplusHistos(TObjArray* listFiles, TH1F** hMass);
  69 Bool_t LoadD0toKpiHistos(TObjArray* listFiles, TH1F** hMass);
  70
  71
  72
  73
  74 void FitMassSpectra(Int_t analysisType=kDplusKpipi,
  75                TString fileNameb="",
  76                TString fileNamec="",
  77                TString fileNamed=""
  78                ){
  79   //
  80
  81   gInterpreter->ExecuteMacro("$ALICE_ROOT/PWG3/vertexingHF/macros/LoadLibraries.C");
  82   gStyle->SetOptTitle(0);
  83
  84   TObjArray* listFiles=new TObjArray();
  85   if(fileNameb!="") listFiles->AddLast(new TObjString(fileNameb.Data()));
  86   if(fileNamec!="") listFiles->AddLast(new TObjString(fileNamec.Data()));
  87   if(fileNamed!="") listFiles->AddLast(new TObjString(fileNamed.Data()));
  88   if(listFiles->GetEntries()==0){
  89     printf("Missing file names in input\n");
  90     return;
  91   }
  92
  93   TH1F** hmass=new TH1F*[nPtBins];
  94   for(Int_t i=0;i<nPtBins;i++) hmass[i]=0x0;
  95
  96   Float_t massD;
  97   Bool_t retCode;
  98   if(analysisType==kD0toKpi){
  99     retCode=LoadD0toKpiHistos(listFiles,hmass);
 100     massD=TDatabasePDG::Instance()->GetParticle(421)->Mass();
 101   }
 102   else if(analysisType==kDplusKpipi){
 103     retCode=LoadDplusHistos(listFiles,hmass);
 104     massD=TDatabasePDG::Instance()->GetParticle(411)->Mass();
 105   }
 106   else{
 107     printf("Wronganalysistype parameter\n");
 108     return;
 109   }
 110   if(!retCode){
 111     printf("ERROR in reading input files\n");
 112     return;
 113   }
 114
 115
 116
 117   TH1D* hCntSig1=new TH1D("hCntSig1","hCntSig1",nPtBins,ptlims);
 118   TH1D* hCntSig2=new TH1D("hCntSig2","hCntSig2",nPtBins,ptlims);
 119   TH1D* hNDiffCntSig1=new TH1D("hNDiffCntSig1","hNDiffCntSig1",nPtBins,ptlims);
 120   TH1D* hNDiffCntSig2=new TH1D("hNDiffCntSig2","hNDiffCntSig2",nPtBins,ptlims);
 121   TH1D* hSignal=new TH1D("hSignal","hSignal",nPtBins,ptlims);
 122   TH1D* hRelErrSig=new TH1D("hRelErrSig","hRelErrSig",nPtBins,ptlims);
 123   TH1D* hInvSignif=new TH1D("hInvSignif","hInvSignif",nPtBins,ptlims);
 124   TH1D* hBackground=new TH1D("hBackground","hBackground",nPtBins,ptlims);
 125   TH1D* hBackgroundNormSigma=new TH1D("hBackgroundNormSigma","hBackgroundNormSigma",nPtBins,ptlims);
 126   TH1D* hSignificance=new TH1D("hSignificance","hSignificance",nPtBins,ptlims);
 127   TH1D* hMass=new TH1D("hMass","hMass",nPtBins,ptlims);
 128   TH1D* hSigma=new TH1D("hSigma","hSigma",nPtBins,ptlims);
 129
 130
 131   Int_t nMassBins=hmass[1]->GetNbinsX();
 132   Double_t hmin=hmass[1]->GetBinLowEdge(3);
 133   Double_t hmax=hmass[1]->GetBinLowEdge(nMassBins-2)+hmass[1]->GetBinWidth(nMassBins-2);
 134   Float_t minBinSum=hmass[1]->FindBin(massD-massRangeForCounting);
 135   Float_t maxBinSum=hmass[1]->FindBin(massD+massRangeForCounting);
 136   Int_t iPad=1;
 137
 138   TF1* funBckStore1=0x0;
 139   TF1* funBckStore2=0x0;
 140   TF1* funBckStore3=0x0;
 141
 142   AliHFMassFitter** fitter=new AliHFMassFitter*[nPtBins];
 143   Double_t arrchisquare[nPtBins];
 144   TCanvas* c1= new TCanvas("c1","MassSpectra",1000,800);
 145   Int_t nx=3, ny=2;
 146   if(nPtBins>6){
 147     nx=4;
 148     ny=3;
 149   }
 150   if(nPtBins>12){
 151     nx=5;
 152     ny=4;
 153   }
 154   c1->Divide(nx,ny);
 155
 156   Double_t sig,errsig,s,errs,b,errb;
 157   for(Int_t iBin=0; iBin<nPtBins; iBin++){
 158     c1->cd(iPad++);
 159     Int_t origNbins=hmass[iBin]->GetNbinsX();
 160     fitter[iBin]=new AliHFMassFitter(hmass[iBin],hmin, hmax,rebin[iBin],typeb,types);
 161     rebin[iBin]=origNbins/fitter[iBin]->GetBinN();
 162     fitter[iBin]->SetReflectionSigmaFactor(factor4refl);
 163     fitter[iBin]->SetInitialGaussianMean(massD);
 164     Bool_t out=fitter[iBin]->MassFitter(0);
 165     if(!out) {
 166       fitter[iBin]->GetHistoClone()->Draw();
 167       continue;
 168     }
 169     Double_t mass=fitter[iBin]->GetMean();
 170     Double_t sigma=fitter[iBin]->GetSigma();
 171     arrchisquare[iBin]=fitter[iBin]->GetReducedChiSquare();
 172     TF1* fB1=fitter[iBin]->GetBackgroundFullRangeFunc();
 173     TF1* fB2=fitter[iBin]->GetBackgroundRecalcFunc();
 174     TF1* fM=fitter[iBin]->GetMassFunc();
 175     if(iBin==0 && fB1) funBckStore1=new TF1(*fB1);
 176     if(iBin==0 && fB2) funBckStore2=new TF1(*fB2);
 177     if(iBin==0 && fM) funBckStore3=new TF1(*fM);
 178
 179     fitter[iBin]->DrawHere(gPad);
 180     fitter[iBin]->Signal(3,s,errs);
 181     fitter[iBin]->Background(3,b,errb);
 182     fitter[iBin]->Significance(3,sig,errsig);
 183     Float_t cntSig1=0.;
 184     Float_t cntSig2=0.;
 185     Float_t cntErr=0.;
 186     for(Int_t iMB=minBinSum; iMB<=maxBinSum; iMB++){
 187       Float_t bkg1=fB1 ? fB1->Eval(hmass[iBin]->GetBinCenter(iMB))/rebin[iBin] : 0;
 188       Float_t bkg2=fB2 ? fB2->Eval(hmass[iBin]->GetBinCenter(iMB))/rebin[iBin] : 0;
 189       cntSig1+=(hmass[iBin]->GetBinContent(iMB)-bkg1);
 190       cntSig2+=(hmass[iBin]->GetBinContent(iMB)-bkg2);
 191       cntErr+=(hmass[iBin]->GetBinContent(iMB));
 192     }
 193     hCntSig1->SetBinContent(iBin+1,cntSig1);
 194     hCntSig1->SetBinError(iBin+1,TMath::Sqrt(cntErr));
 195     hNDiffCntSig1->SetBinContent(iBin+1,(s-cntSig1)/s);
 196     hNDiffCntSig1->SetBinError(iBin+1,TMath::Sqrt(cntErr)/s);
 197     hCntSig2->SetBinContent(iBin+1,cntSig2);
 198     hNDiffCntSig2->SetBinContent(iBin+1,(s-cntSig2)/s);
 199     hNDiffCntSig2->SetBinError(iBin+1,TMath::Sqrt(cntErr)/s);
 200     hCntSig2->SetBinError(iBin+1,TMath::Sqrt(cntErr));
 201     hSignal->SetBinContent(iBin+1,s);
 202     hSignal->SetBinError(iBin+1,errs);
 203     hRelErrSig->SetBinContent(iBin+1,errs/s);
 204     hInvSignif->SetBinContent(iBin+1,1/sig);
 205     hInvSignif->SetBinError(iBin+1,errsig/(sig*sig));
 206     hBackground->SetBinContent(iBin+1,b); //consider sigma
 207     hBackground->SetBinError(iBin+1,errb);
 208     hBackgroundNormSigma->SetBinContent(iBin+1,b/(3*fitter[iBin]->GetSigma())*(3*0.012)); //consider sigma
 209     hBackgroundNormSigma->SetBinError(iBin+1,errb);
 210     hSignificance->SetBinContent(iBin+1,sig);
 211     hSignificance->SetBinError(iBin+1,errsig);
 212     hMass->SetBinContent(iBin+1,mass);
 213     hMass->SetBinError(iBin+1,0.0001);
 214     hSigma->SetBinContent(iBin+1,sigma);
 215     hSigma->SetBinError(iBin+1,0.0001);
 216
 217   }
 218
 219   /*
 220   c1->cd(1); // is some cases the fitting function of 1st bin get lost
 221   funBckStore1->Draw("same");
 222   funBckStore2->Draw("same");
 223   funBckStore3->Draw("same");
 224   */
 225
 226   TCanvas *cpar=new TCanvas("cpar","Fit params",1200,600);
 227   cpar->Divide(2,1);
 228   cpar->cd(1);
 229   hMass->SetMarkerStyle(20);
 230   hMass->Draw("PE");
 231   hMass->GetXaxis()->SetTitle("Pt (GeV/c)");
 232   hMass->GetXaxis()->SetTitle("Mass (GeV/c^{2})");
 233   cpar->cd(2);
 234   hSigma->SetMarkerStyle(20);
 235   hSigma->Draw("PE");
 236   hSigma->GetXaxis()->SetTitle("Pt (GeV/c)");
 237   hSigma->GetXaxis()->SetTitle("Sigma (GeV/c^{2})");
 238
 239   TCanvas* csig=new TCanvas("csig","Results",1200,600);
 240   csig->Divide(3,1);
 241   csig->cd(1);
 242   hSignal->SetMarkerStyle(20);
 243   hSignal->SetMarkerColor(4);
 244   hSignal->SetLineColor(4);
 245   hSignal->GetXaxis()->SetTitle("Pt (GeV/c)");
 246   hSignal->GetYaxis()->SetTitle("Signal");
 247   hSignal->Draw("P");
 248   hCntSig1->SetMarkerStyle(26);
 249   hCntSig1->SetMarkerColor(2);
 250   hCntSig1->SetLineColor(2);
 251   hCntSig1->Draw("PSAME");
 252   hCntSig2->SetMarkerStyle(29);
 253   hCntSig2->SetMarkerColor(kGray+1);
 254   hCntSig2->SetLineColor(kGray+1);
 255   hCntSig2->Draw("PSAME");
 256   TLegend* leg=new TLegend(0.4,0.7,0.89,0.89);
 257   leg->SetFillColor(0);
 258   TLegendEntry* ent=leg->AddEntry(hSignal,"From Fit","PL");
 259   ent->SetTextColor(hSignal->GetMarkerColor());
 260   ent=leg->AddEntry(hCntSig1,"From Counting1","PL");
 261   ent->SetTextColor(hCntSig1->GetMarkerColor());
 262   ent=leg->AddEntry(hCntSig2,"From Counting2","PL");
 263   ent->SetTextColor(hCntSig2->GetMarkerColor());
 264   leg->Draw();
 265   csig->cd(2);
 266   hBackground->SetMarkerStyle(20);
 267   hBackground->Draw("P");
 268   hBackground->GetXaxis()->SetTitle("Pt (GeV/c)");
 269   hBackground->GetYaxis()->SetTitle("Background");
 270   csig->cd(3);
 271   hSignificance->SetMarkerStyle(20);
 272   hSignificance->Draw("P");
 273   hSignificance->GetXaxis()->SetTitle("Pt (GeV/c)");
 274   hSignificance->GetYaxis()->SetTitle("Significance");
 275
 276   TCanvas* cDiffS=new TCanvas("cDiffS","Difference",1200,600);
 277   cDiffS->Divide(2,1);
 278   cDiffS->cd(1);
 279   hRelErrSig->SetMarkerStyle(20); //fullcircle
 280   hRelErrSig->SetTitleOffset(1.2);
 281   hRelErrSig->SetTitle("Rel Error from Fit;p_{t} (GeV/c);Signal Relative Error");
 282   hRelErrSig->Draw("P");
 283   hInvSignif->SetMarkerStyle(21); //fullsquare
 284   hInvSignif->SetMarkerColor(kMagenta+1);
 285   hInvSignif->SetLineColor(kMagenta+1);
 286   hInvSignif->Draw("PSAMES");
 287   TLegend* leg2=new TLegend(0.4,0.7,0.89,0.89);
 288   leg2->SetFillColor(0);
 289   TLegendEntry* ent2=leg2->AddEntry(hRelErrSig,"From Fit","P");
 290   ent2->SetTextColor(hRelErrSig->GetMarkerColor());
 291   ent2=leg2->AddEntry(hInvSignif,"1/Significance","PL");
 292   ent2->SetTextColor(hInvSignif->GetMarkerColor());
 293   leg2->Draw();
 294
 295   cDiffS->cd(2);
 296   hNDiffCntSig1->SetMarkerStyle(26);
 297   hNDiffCntSig1->SetMarkerColor(2);
 298   hNDiffCntSig1->SetLineColor(2);
 299   hNDiffCntSig1->SetTitle("Cmp Fit-Count;p_{t} (GeV/c);(S_{fit}-S_{count})/S_{fit}");
 300   hNDiffCntSig1->Draw("P");
 301   hNDiffCntSig2->SetMarkerStyle(29);
 302   hNDiffCntSig2->SetMarkerColor(kGray+1);
 303   hNDiffCntSig2->SetLineColor(kGray+1);
 304   hNDiffCntSig2->Draw("PSAME");
 305   TLegend* leg1=new TLegend(0.4,0.7,0.89,0.89);
 306   leg1->SetFillColor(0);
 307   TLegendEntry* ent1=leg1->AddEntry(hNDiffCntSig1,"From Counting1","PL");
 308   ent1->SetTextColor(hNDiffCntSig1->GetMarkerColor());
 309   ent1=leg1->AddEntry(hNDiffCntSig2,"From Counting2","PL");
 310   ent1->SetTextColor(hNDiffCntSig2->GetMarkerColor());
 311   leg1->Draw();
 312
 313   TGraph* grReducedChiSquare=new TGraph(nPtBins,ptlims,arrchisquare);
 314   grReducedChiSquare->SetName("grReducedChiSquare");
 315   grReducedChiSquare->SetTitle("Reduced Chi2;p_{t} (GeV/c);#tilde{#chi}^{2}");
 316   TCanvas *cChi2=new TCanvas("cChi2","reduced chi square",600,600);
 317   cChi2->cd();
 318   grReducedChiSquare->SetMarkerStyle(21);
 319   grReducedChiSquare->Draw("AP");
 320
 321   TCanvas* cbkgNormSigma=new TCanvas("cbkgNormSigma","Background normalized to sigma",400,600);
 322   cbkgNormSigma->cd();
 323   hBackgroundNormSigma->SetMarkerStyle(20);
 324   hBackgroundNormSigma->Draw("P");
 325   hBackgroundNormSigma->GetXaxis()->SetTitle("Pt (GeV/c)");
 326   hBackgroundNormSigma->GetYaxis()->SetTitle("Background #times 3 #times 0.012/ (3 #times #sigma)");
 327   hBackgroundNormSigma->Draw();
 328
 329
 330   TString partname="Both";
 331   if(optPartAntiPart==kParticleOnly) {
 332     if(analysisType==kD0toKpi) partname="D0";
 333     if(analysisType==kDplusKpipi) partname="Dplus";
 334   }
 335   if(optPartAntiPart==kAntiParticleOnly) {
 336     if(analysisType==kD0toKpi) partname="D0bar";
 337     if(analysisType==kDplusKpipi) partname="Dminus";
 338   }
 339
 340   TFile* outf=new TFile(Form("RawYield%s.root",partname.Data()),"update");
 341   outf->cd();
 342   hMass->Write();
 343   hSigma->Write();
 344   hCntSig1->Write();
 345   hCntSig2->Write();
 346   hNDiffCntSig1->Write();
 347   hNDiffCntSig2->Write();
 348   hSignal->Write();
 349   hRelErrSig->Write();
 350   hInvSignif->Write();
 351   hBackground->Write();
 352   hBackgroundNormSigma->Write();
 353   hSignificance->Write();
 354   grReducedChiSquare->Write();
 355   hPtMass->Write();
 356   outf->Close();
 357 }
 358
 359
 360 Bool_t LoadDplusHistos(TObjArray* listFiles, TH1F** hMass){
 361
 362   Int_t nFiles=listFiles->GetEntries();
 363   TList **hlist=new TList*[nFiles];
 364   AliRDHFCutsDplustoKpipi** dcuts=new AliRDHFCutsDplustoKpipi*[nFiles];
 365
 366   Int_t nReadFiles=0;
 367   for(Int_t iFile=0; iFile<nFiles; iFile++){
 368     TString fName=((TObjString*)listFiles->At(iFile))->GetString();
 369     TFile *f=TFile::Open(fName.Data());
 370     if(!f){
 371       printf("ERROR: file %s does not exist\n",fName.Data());
 372       continue;
 373     }
 374     printf("Open File %s\n",f->GetName());
 375     TDirectory *dir = (TDirectory*)f->Get("PWG3_D2H_InvMassDplus");
 376     if(!dir){
 377       printf("ERROR: directory PWG3_D2H_InvMassDplus not found in %s\n",fName.Data());
 378       continue;
 379     }
 380     hlist[nReadFiles]=(TList*)dir->Get("coutputDplus");
 381     TList *listcut = (TList*)dir->Get("coutputDplusCuts");
 382     dcuts[nReadFiles]=(AliRDHFCutsDplustoKpipi*)listcut->At(1);
 383     if(nReadFiles>0){
 384       Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]);
 385       if(!sameCuts){
 386         printf("ERROR: Cut objects do not match\n");
 387         return kFALSE;
 388       }
 389     }
 390     nReadFiles++;
 391   }
 392   if(nReadFiles<nFiles){
 393     printf("WARNING: not all requested files have been found\n");
 394   }
 395
 396   Int_t nPtBinsCuts=dcuts[0]->GetNPtBins();
 397   printf("Number of pt bins for cut object = %d\n",nPtBins);
 398   Float_t *ptlimsCuts=dcuts[0]->GetPtBinLimits();
 399   ptlimsCuts[nPtBinsCuts]=ptlimsCuts[nPtBinsCuts-1]+4.;
 400
 401   Int_t iFinBin=0;
 402   for(Int_t i=0;i<nPtBinsCuts;i++){
 403     if(ptlimsCuts[i]>=ptlims[iFinBin+1]) iFinBin+=1;
 404     if(iFinBin>nPtBins) break;
 405     if(ptlimsCuts[i]>=ptlims[iFinBin] &&
 406        ptlimsCuts[i+1]<=ptlims[iFinBin+1]){
 407       for(Int_t iFile=0; iFile<nReadFiles; iFile++){
 408         TString histoName;
 409         if(optPartAntiPart==kBoth) histoName.Form("hMassPt%dTC",i);
 410         else if(optPartAntiPart==kParticleOnly) histoName.Form("hMassPt%dTCPlus",i);
 411         else if(optPartAntiPart==kAntiParticleOnly) histoName.Form("hMassPt%dTCMinus",i);
 412         TH1F* htemp=(TH1F*)hlist[iFile]->FindObject(histoName.Data());
 413         if(!htemp){
 414           printf("ERROR: Histogram %s not found\n",histoName.Data());
 415           return kFALSE;
 416         }
 417         if(!hMass[iFinBin]){
 418           hMass[iFinBin]=new TH1F(*htemp);
 419         }else{
 420           hMass[iFinBin]->Add(htemp);
 421         }
 422       }
 423     }
 424   }
 425   TString partname="Both";
 426   if(optPartAntiPart==kParticleOnly) partname="Dplus";
 427   if(optPartAntiPart==kAntiParticleOnly) partname="Dminus";
 428
 429   for(Int_t iFile=0; iFile<nReadFiles; iFile++){
 430     TH2F* htemp2=(TH2F*)hlist[iFile]->FindObject("hPtVsMassTC");
 431     if(iFile==0){
 432       hPtMass=new TH2F(*htemp2);
 433     }else{
 434       hPtMass->Add(htemp2);
 435     }
 436   }
 437
 438   TFile* outf=new TFile(Form("RawYield%s.root",partname.Data()),"recreate");
 439   outf->cd();
 440   dcuts[0]->Write();
 441   outf->Close();
 442
 443   return kTRUE;
 444
 445 }
 446
 447 Bool_t LoadD0toKpiHistos(TObjArray* listFiles, TH1F** hMass){
 448   //
 449   Int_t nFiles=listFiles->GetEntries();
 450   TList **hlist=new TList*[nFiles];
 451   AliRDHFCutsD0toKpi** dcuts=new AliRDHFCutsD0toKpi*[nFiles];
 452
 453   Int_t nReadFiles=0;
 454   for(Int_t iFile=0; iFile<nFiles; iFile++){
 455     TString fName=((TObjString*)listFiles->At(iFile))->GetString();
 456     TFile *f=TFile::Open(fName.Data());
 457     if(!f){
 458       printf("ERROR: file %s does not exist\n",fName.Data());
 459       continue;
 460     }
 461     printf("Open File %s\n",f->GetName());
 462
 463     TString dirname="PWG3_D2H_D0InvMass";
 464     if(optPartAntiPart==kParticleOnly) dirname+="D0";
 465     else if(optPartAntiPart==kAntiParticleOnly) dirname+="D0bar";
 466     if(cutsappliedondistr) dirname+="C";
 467     TDirectory *dir = (TDirectory*)f->Get(dirname);
 468     if(!dir){
 469       printf("ERROR: directory %s not found in %s\n",dirname.Data(),fName.Data());
 470       continue;
 471     }
 472     TString listmassname="coutputmassD0Mass";
 473     if(optPartAntiPart==kParticleOnly) listmassname+="D0";
 474     else if(optPartAntiPart==kAntiParticleOnly) listmassname+="D0bar";
 475     if(cutsappliedondistr) listmassname+="C";
 476
 477     hlist[nReadFiles]=(TList*)dir->Get(listmassname);
 478
 479     TString cutsobjname="cutsD0";
 480     if(optPartAntiPart==kParticleOnly) cutsobjname+="D0";
 481     else if(optPartAntiPart==kAntiParticleOnly) cutsobjname+="D0bar";
 482     if(cutsappliedondistr) cutsobjname+="C";
 483
 484     dcuts[nReadFiles]=(AliRDHFCutsD0toKpi*)dir->Get(cutsobjname);
 485     if(nReadFiles>0){
 486       Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]);
 487       if(!sameCuts){
 488         printf("ERROR: Cut objects do not match\n");
 489         return kFALSE;
 490       }
 491     }
 492     nReadFiles++;
 493   }
 494   if(nReadFiles<nFiles){
 495     printf("WARNING: not all requested files have been found\n");
 496     if (nReadFiles==0) {
 497       printf("ERROR: Any file/dir found\n");
 498       return kFALSE;
 499     }
 500   }
 501
 502   Int_t nPtBinsCuts=dcuts[0]->GetNPtBins();
 503   printf("Number of pt bins for cut object = %d\n",nPtBins);
 504   Float_t *ptlimsCuts=dcuts[0]->GetPtBinLimits();
 505   ptlimsCuts[nPtBinsCuts]=ptlimsCuts[nPtBinsCuts-1]+4.;
 506
 507   Int_t iFinBin=0;
 508   for(Int_t i=0;i<nPtBinsCuts;i++){
 509     if(ptlimsCuts[i]>=ptlims[iFinBin+1]) iFinBin+=1;
 510     if(iFinBin>nPtBins) break;
 511     if(ptlimsCuts[i]>=ptlims[iFinBin] &&
 512        ptlimsCuts[i+1]<=ptlims[iFinBin+1]){
 513       for(Int_t iFile=0; iFile<nReadFiles; iFile++){
 514         TH1F* htemp=(TH1F*)hlist[iFile]->FindObject(Form("histMass_%d",i));
 515         if(!hMass[iFinBin]){
 516           hMass[iFinBin]=new TH1F(*htemp);
 517         }else{
 518           hMass[iFinBin]->Add(htemp);
 519         }
 520       }
 521     }
 522   }
 523   TString partname="Both";
 524   if(optPartAntiPart==kParticleOnly) partname="D0";
 525   if(optPartAntiPart==kAntiParticleOnly) partname="D0bar";
 526
 527   TFile* outf=new TFile(Form("RawYield%s.root",partname.Data()),"recreate");
 528   outf->cd();
 529   dcuts[0]->Write();
 530   outf->Close();
 531   return kTRUE;
 532 }
 533
 534 void CompareFitTypes(TString* paths, TString* legtext,Int_t ncmp=3,TString* filenameYield=0x0){
 535   //read ncmp RawYield.roots and draw them together
 536   //set the global variable nevents before running
 537   //arguments:
 538   // - paths= vector of ncmp dimension with the paths of the file RawYield.root
 539   // - legtext= vector of ncmp dimension with the label for the legend
 540   // - ncmp= number of files to compare (default is 3)
 541   // - filenameYield= optional ncmp-dimensional array with the difference between the names of the files to be compared (useful if the 2 files are in the same directory but have different names)
 542
 543   gStyle->SetOptStat(0);
 544   gStyle->SetOptTitle(0);
 545   gStyle->SetCanvasColor(0);
 546   gStyle->SetFrameFillColor(0);
 547   gStyle->SetTitleFillColor(0);
 548   gStyle->SetFrameBorderMode(0);
 549
 550   if(!filenameYield) filenameYield=new TString[ncmp];
 551
 552   for(Int_t k=0;k<ncmp;k++){
 553     if(!filenameYield) filenameYield[k]="RawYield.root";
 554     filenameYield[k].Prepend(paths[k]);
 555   }
 556
 557   TCanvas* cSig=new TCanvas("cSig","Results",1200,600);
 558   cSig->Divide(3,1);
 559   TCanvas* cBkgN=new TCanvas("cBkgN","Background normalized to sigma",400,600);
 560   TCanvas* cDiffS=new TCanvas("cDiffS","Difference",1200,600);
 561   cDiffS->Divide(2,1);
 562   TCanvas *cChi2=new TCanvas("cChi2","reduced chi square",600,600);
 563
 564   TLegend* leg1=new TLegend(0.4,0.7,0.89,0.89);
 565   leg1->SetFillColor(0);
 566   TLegend* leg2=(TLegend*)leg1->Clone();
 567   TLegend* leg3=(TLegend*)leg1->Clone();
 568   TLegend* leg4=new TLegend(0.4,0.6,0.8,0.89);
 569   leg4->SetFillColor(0);
 570
 571   for(Int_t iTypes=0;iTypes<ncmp;iTypes++){
 572     TFile* fin=new TFile(filenameYield[iTypes]);
 573     if(!fin){
 574       printf("WARNING: %s not found",filenameYield[iTypes].Data());
 575       continue;
 576     }
 577
 578     TH1F* hSignal=(TH1F*)fin->Get("hSignal");
 579     TH1F* hBackground=(TH1F*)fin->Get("hBackground");
 580     TH1F* hBackgroundNormSigma=(TH1F*)fin->Get("hBackgroundNormSigma");
 581     TH1F* hSignificance=(TH1F*)fin->Get("hSignificance");
 582     hSignal->SetName(Form("%s%d",hSignal->GetName(),iTypes));
 583     hBackground->SetName(Form("%s%d",hBackground->GetName(),iTypes));
 584     hBackgroundNormSigma->SetName(Form("%s%d",hBackgroundNormSigma->GetName(),iTypes));
 585     hSignificance->SetName(Form("%s%d",hSignificance->GetName(),iTypes));
 586
 587     hSignal->SetMarkerColor(iTypes+2);
 588     hSignal->SetLineColor(iTypes+2);
 589     hBackground->SetMarkerColor(iTypes+2);
 590     hBackground->SetLineColor(iTypes+2);
 591     hBackgroundNormSigma->SetMarkerColor(iTypes+2);
 592     hBackgroundNormSigma->SetLineColor(iTypes+2);
 593     hSignificance->SetMarkerColor(iTypes+2);
 594     hSignificance->SetLineColor(iTypes+2);
 595
 596     TLegendEntry* ent4=leg4->AddEntry(hSignal,Form("%s",legtext[iTypes].Data()),"PL");
 597     ent4->SetTextColor(hSignal->GetMarkerColor());
 598
 599     if(ncmp==nsamples){
 600       printf("Info: Normalizing signal, background and significance to the number of events\n");
 601       hSignal->Scale(1./nevents[iTypes]);
 602       hBackground->Scale(1./nevents[iTypes]);
 603       hBackgroundNormSigma->Scale(1./nevents[iTypes]);
 604       hSignificance->Scale(1./TMath::Sqrt(nevents[iTypes]));
 605     }
 606
 607     if(iTypes==0){
 608       cSig->cd(1);
 609       hSignal->DrawClone("P");
 610       cSig->cd(2);
 611       hBackground->DrawClone("P");
 612       cSig->cd(3);
 613       hSignificance->DrawClone("P");
 614       cBkgN->cd();
 615       hBackgroundNormSigma->DrawClone("P");
 616     } else{
 617       cSig->cd(1);
 618       hSignal->DrawClone("Psames");
 619       cSig->cd(2);
 620       hBackground->DrawClone("Psames");
 621       cSig->cd(3);
 622       hSignificance->DrawClone("Psames");
 623       cBkgN->cd();
 624       hBackgroundNormSigma->DrawClone("Psames");
 625     }
 626
 627     TH1F* hRelErrSig=(TH1F*)fin->Get("hRelErrSig");
 628     TH1F* hInvSignif=(TH1F*)fin->Get("hInvSignif");
 629     hRelErrSig->SetName(Form("%s%d",hRelErrSig->GetName(),iTypes));
 630     hInvSignif->SetName(Form("%s%d",hInvSignif->GetName(),iTypes));
 631
 632     hRelErrSig->SetMarkerColor(iTypes+2);
 633     hRelErrSig->SetLineColor(iTypes+2);
 634     hInvSignif->SetMarkerColor(iTypes+2);
 635     hInvSignif->SetLineColor(iTypes+2);
 636
 637     TLegendEntry* ent1=leg1->AddEntry(hRelErrSig,Form("From Fit (%s)",legtext[iTypes].Data()),"P");
 638     ent1->SetTextColor(hRelErrSig->GetMarkerColor());
 639     ent1=leg1->AddEntry(hInvSignif,Form("1/Significance (%s)",legtext[iTypes].Data()),"PL");
 640     ent1->SetTextColor(hInvSignif->GetMarkerColor());
 641
 642     cDiffS->cd(1);
 643     if(iTypes==0){
 644       hRelErrSig->DrawClone("P");
 645       hInvSignif->DrawClone();
 646     } else{
 647       hRelErrSig->DrawClone("Psames");
 648       hInvSignif->DrawClone("sames");
 649     }
 650
 651     TH1F* hNDiffCntSig1=(TH1F*)fin->Get("hNDiffCntSig1");
 652     TH1F* hNDiffCntSig2=(TH1F*)fin->Get("hNDiffCntSig2");
 653     hNDiffCntSig1->SetName(Form("%s%d",hNDiffCntSig1->GetName(),iTypes));
 654     hNDiffCntSig2->SetName(Form("%s%d",hNDiffCntSig2->GetName(),iTypes));
 655
 656     hNDiffCntSig1->SetMarkerColor(iTypes+2);
 657     hNDiffCntSig1->SetLineColor(iTypes+2);
 658     hNDiffCntSig2->SetMarkerColor(iTypes+2);
 659     hNDiffCntSig2->SetLineColor(iTypes+2);
 660     TLegendEntry* ent2=leg2->AddEntry(hNDiffCntSig1,Form("From Counting1 (%s)",legtext[iTypes].Data()),"PL");
 661     ent2->SetTextColor(hNDiffCntSig1->GetMarkerColor());
 662     ent2=leg2->AddEntry(hNDiffCntSig2,Form("From Counting2 (%s)",legtext[iTypes].Data()),"PL");
 663     ent2->SetTextColor(hNDiffCntSig2->GetMarkerColor());
 664
 665     cDiffS->cd(2);
 666     if(iTypes==0){
 667       hNDiffCntSig1->DrawClone();
 668       hNDiffCntSig2->DrawClone();
 669     }else{
 670      hNDiffCntSig1->DrawClone("sames");
 671      hNDiffCntSig2->DrawClone("sames");
 672     }
 673
 674     TGraph* grReducedChiSquare=(TGraph*)fin->Get("grReducedChiSquare");
 675     grReducedChiSquare->SetName(Form("%s%d",grReducedChiSquare->GetName(),iTypes));
 676
 677     grReducedChiSquare->SetMarkerColor(iTypes+2);
 678     grReducedChiSquare->SetLineColor(iTypes+2);
 679     TLegendEntry* ent3=leg3->AddEntry(grReducedChiSquare,Form("%s",legtext[iTypes].Data()),"PL");
 680     ent3->SetTextColor(grReducedChiSquare->GetMarkerColor());
 681
 682     cChi2->cd();
 683     if(iTypes==0) grReducedChiSquare->DrawClone("AP");
 684     else grReducedChiSquare->DrawClone("P");
 685   }
 686
 687   cSig->cd(1);
 688   leg4->Draw();
 689
 690   cDiffS->cd(1);
 691   leg1->Draw();
 692
 693   cDiffS->cd(2);
 694   leg2->Draw();
 695
 696   cChi2->cd();
 697   leg3->Draw();
 698
 699   TFile* fout=new TFile("ComparisonRawYield.root","RECREATE");
 700   fout->cd();
 701   cDiffS->Write();
 702   cChi2->Write();
 703   fout->Close();
 704 }