-#if !defined(__CINT__) || defined(__MAKECINT__)\r
-#include <TInterpreter.h>\r
-#include <TString.h>\r
-#include <TObjString.h>\r
-#include <TObjArray.h>\r
-#include <TMath.h>\r
-#include <TFile.h>\r
-#include <TCanvas.h>\r
-#include <TH1.h>\r
-#include <TH1F.h>\r
-#include <TH2F.h>\r
-#include <TH3.h>\r
-#include <TH3F.h>\r
-#include <TH1D.h>\r
-#include <TF1.h>\r
-#include <TSystem.h>\r
-#include <TStyle.h>\r
-#include <TLegend.h>\r
-#include <TList.h>\r
-#include <TLegendEntry.h>\r
-#include <TDatabasePDG.h>\r
-#include <TGraph.h>\r
-\r
-#include "AliAODRecoDecayHF.h"\r
-#include "AliRDHFCuts.h"\r
-#include "AliRDHFCutsDplustoKpipi.h"\r
-#include "AliRDHFCutsD0toKpi.h"\r
-#include "AliHFMassFitter.h"\r
-#endif\r
-\r
-/* $Id$ */ \r
-\r
-enum {kD0toKpi, kDplusKpipi};\r
-enum {kBoth, kParticleOnly, kAntiParticleOnly};\r
-enum {kExpo=0, kLinear, kPol2};\r
-enum {kGaus=0, kDoubleGaus};\r
-enum {kCorr=0, kUnCorr, kNoPid};\r
-\r
-\r
-// Common variables: to be configured by the user\r
-const Int_t nPtBins=5;\r
-Double_t ptlims[nPtBins+1]={2.,4.,6.,8.,12.,16.};\r
-Int_t rebin[nPtBins]={4,6,6,6,8};\r
-Double_t sigmapt[nPtBins]={ 0.012, 0.016, 0.016, 0.018, 0.20 };\r
-Bool_t fixPeakSigma = kFALSE;\r
-//\r
-const Int_t nMultbins=6;\r
-Double_t multlims[nMultbins+1]={1.,9.,14.,20.,31.,49.,100.};\r
-//\r
-Int_t firstUsedBin[nPtBins]={-1,-1,-1,-1,-1}; // -1 uses all bins, >=1 to set the lower bin to be accepted from original histo\r
-//\r
-Bool_t isMC=kFALSE;\r
-Int_t typeb=kExpo;\r
-Int_t types=kGaus;\r
-Int_t optPartAntiPart=kBoth;\r
-Int_t factor4refl=0;\r
-Float_t massRangeForCounting=0.05; // GeV\r
-TH2F* hPtMass=0x0;\r
-TString suffix="StdPid";\r
-//for D0only\r
-const Int_t nsamples=2;//3;\r
-Int_t nevents[nsamples]={1.18860695e+08 /*LHC10dnewTPCpid*/,9.0374946e+07 /*LHC10b+c*/};\r
-\r
-// Functions\r
-Bool_t LoadDplusHistos(TObjArray* listFiles, TH3F** hPtMassMult, TH2F** hNtrZvtx, TH2F** hNtrZvtxCorr, const char *CutsType, Int_t Option);\r
-Bool_t LoadD0toKpiHistos(TObjArray* listFiles, TH3F** hPtMassMult, TH2F** hNtrZvtx, TH2F** hNtrZvtxCorr, \r
- const char *CutsType, Int_t Option);\r
-Bool_t CheckNtrVsZvtx(TH2F** hNtrZvtx, TH2F** hNtrZvtxCorr, Int_t nFiles);\r
-TH1F* RebinHisto(TH1F* hOrig, Int_t reb, Int_t firstUse=-1);\r
-\r
-\r
-void ReadDvsMultiplicity(Int_t analysisType=kD0toKpi,\r
- TString fileNameb="AnalysisResults.root",\r
- TString fileNamec="",\r
- TString fileNamed="",\r
- TString fileNamee="",\r
- const char *CutsType="",\r
- Int_t Option=kCorr)\r
-{\r
- // gSystem->SetIncludePath("-I. -I$ROOTSYS/include -I$ALICE_ROOT -I$ALICE_ROOT/include -I$ALICE_ROOT/ITS -I$ALICE_ROOT/TPC -I$ALICE_ROOT/CONTAINERS -I$ALICE_ROOT/STEER/STEER -I$ALICE_ROOT/STEER/STEERBase -I$ALICE_ROOT/STEER/ESD -I$ALICE_ROOT/STEER/AOD -I$ALICE_ROOT/TRD -I$ALICE_ROOT/macros -I$ALICE_ROOT/ANALYSIS -I$ALICE_ROOT/OADB -I$ALICE_ROOT/PWGHF -I$ALICE_ROOT/PWGHF/base -I$ALICE_ROOT/PWGHF/vertexingHF -I$ALICE_ROOT/PWG/FLOW/Case -I$ALICE_ROOT/PWG/FLOW/Tasks -g");\r
-\r
- // gInterpreter->ExecuteMacro("$ALICE_ROOT/PWGHF/vertexingHF/macros/LoadLibraries.C");\r
- gStyle->SetOptTitle(1);\r
-\r
- Int_t nFiles=0;\r
- TObjArray* listFiles=new TObjArray();\r
- if(fileNameb!="") { listFiles->AddLast(new TObjString(fileNameb.Data())); nFiles++; }\r
- if(fileNamec!="") { listFiles->AddLast(new TObjString(fileNamec.Data())); nFiles++; }\r
- if(fileNamed!="") { listFiles->AddLast(new TObjString(fileNamed.Data())); nFiles++; }\r
- if(fileNamee!="") { listFiles->AddLast(new TObjString(fileNamee.Data())); nFiles++; }\r
- if(listFiles->GetEntries()==0){\r
- printf("Missing file names in input\n");\r
- return;\r
- }\r
- TH3F** hPtMassMult=new TH3F*[1];\r
- hPtMassMult[0]=0x0;\r
- TH1F** hmass=new TH1F*[nPtBins*nMultbins];\r
- for(Int_t i=0; i<nPtBins*nMultbins; i++) hmass[i]=0x0;\r
- TH2F** hNtrZvtx=new TH2F*[4];\r
- for(Int_t i=0; i<4; i++) hNtrZvtx[i]=0x0;\r
- TH2F** hNtrZvtxCorr=new TH2F*[4];\r
- for(Int_t i=0; i<4; i++) hNtrZvtxCorr[i]=0x0;\r
-\r
- Float_t massD;\r
- Bool_t retCode;\r
- if(analysisType==kD0toKpi) {\r
- massD=1.86484003067016602e+00;//(Float_t)(TDatabasePDG::Instance()->GetParticle(421)->Mass());\r
- retCode=LoadD0toKpiHistos(listFiles,hPtMassMult,hNtrZvtx,hNtrZvtxCorr,CutsType,Option);\r
- }\r
- else if(analysisType==kDplusKpipi) {\r
- massD=1.86961996555328369e+00;//(Float_t)(TDatabasePDG::Instance()->GetParticle(411)->Mass());\r
- retCode=LoadDplusHistos(listFiles,hPtMassMult,hNtrZvtx,hNtrZvtxCorr,CutsType,Option);\r
- }\r
- else{\r
- printf("Wrong analysis type parameter\n");\r
- return;\r
- }\r
- if(!retCode){\r
- printf("ERROR in reading input files\n");\r
- return;\r
- }\r
- //\r
- // Check the multiplicity variables\r
- //\r
- Bool_t isMult = CheckNtrVsZvtx(hNtrZvtx,hNtrZvtxCorr,nFiles);\r
- if(!isMult) return;\r
- //\r
- //\r
- printf(" Preparing the mass fits");\r
- TH1F* hCntSig1[nMultbins];\r
- TH1F* hCntSig2[nMultbins];\r
- TH1F* hNDiffCntSig1[nMultbins];\r
- TH1F* hNDiffCntSig2[nMultbins];\r
- TH1F* hSignal[nMultbins];\r
- TH1F* hRelErrSig[nMultbins];\r
- TH1F* hInvSignif[nMultbins];\r
- TH1F* hBackground[nMultbins];\r
- TH1F* hBackgroundNormSigma[nMultbins];\r
- TH1F* hSignificance[nMultbins];\r
- TH1F* hMass[nMultbins];\r
- TH1F* hSigma[nMultbins];\r
- for(Int_t i=0; i<nMultbins; i++){\r
- hCntSig1[i]=new TH1F(Form("hCntSig1_%d",i),Form("hCntSig1_%d",i),nPtBins,ptlims);\r
- hCntSig2[i]=new TH1F(Form("hCntSig2_%d",i),Form("hCntSig2_%d",i),nPtBins,ptlims);\r
- hNDiffCntSig1[i]=new TH1F(Form("hNDiffCntSig1_%d",i),Form("hNDiffCntSig1_%d",i),nPtBins,ptlims);\r
- hNDiffCntSig2[i]=new TH1F(Form("hNDiffCntSig2_%d",i),Form("hNDiffCntSig2_%d",i),nPtBins,ptlims);\r
- hSignal[i]=new TH1F(Form("hSignal_%d",i),Form("hSignal_%d",i),nPtBins,ptlims);\r
- hRelErrSig[i]=new TH1F(Form("hRelErrSig_%d",i),Form("hRelErrSig_%d",i),nPtBins,ptlims);\r
- hInvSignif[i]=new TH1F(Form("hInvSignif_%d",i),Form("hInvSignif_%d",i),nPtBins,ptlims);\r
- hBackground[i]=new TH1F(Form("hBackground_%d",i),Form("hBackground_%d",i),nPtBins,ptlims);\r
- hBackgroundNormSigma[i]=new TH1F(Form("hBackgroundNormSigma_%d",i),Form("hBackgroundNormSigma_%d",i),nPtBins,ptlims);\r
- hSignificance[i]=new TH1F(Form("hSignificance_%d",i),Form("hSignificance_%d",i),nPtBins,ptlims);\r
- hMass[i]=new TH1F(Form("hMass_%d",i),Form("hMass_%d",i),nPtBins,ptlims);\r
- hSigma[i]=new TH1F(Form("hSigma_%d",i),Form("hSigma_%d",i),nPtBins,ptlims);\r
- }\r
- printf(", defined...\n");\r
-\r
- // std::cout << " htemp :"<<hPtMassMult[0]<<std::endl;\r
- TH1F* hptaxis = (TH1F*)hPtMassMult[0]->ProjectionZ("hptaxis");\r
- TH1F* hmassaxis = (TH1F*)hPtMassMult[0]->ProjectionY("hmassaxis");\r
- TH1F* hmultaxis = (TH1F*)hPtMassMult[0]->ProjectionX("hmultaxis");\r
- Int_t nMassBins=hmassaxis->GetNbinsX();\r
- Double_t hmin=hmassaxis->GetBinLowEdge(3);\r
- Double_t hmax=hmassaxis->GetBinLowEdge(nMassBins-2) + hmassaxis->GetBinWidth(nMassBins-2);\r
- Float_t minBinSum=hmassaxis->FindBin(massD-massRangeForCounting);\r
- Float_t maxBinSum=hmassaxis->FindBin(massD+massRangeForCounting);\r
- Int_t iPad=1;\r
- \r
- printf("Now initializing the fit functions\n");\r
- TF1* funBckStore1=0x0;\r
- TF1* funBckStore2=0x0;\r
- TF1* funBckStore3=0x0;\r
-\r
- Int_t nPtMultbins = nPtBins*nMultbins;\r
- AliHFMassFitter** fitter=new AliHFMassFitter*[nPtMultbins];\r
- Double_t arrchisquare0[nPtBins];\r
- Double_t arrchisquare1[nPtBins];\r
- Double_t arrchisquare2[nPtBins];\r
- Double_t arrchisquare3[nPtBins];\r
- Double_t arrchisquare4[nPtBins];\r
- Double_t arrchisquare5[nPtBins];\r
- for(Int_t i=0; i<nPtBins; i++){\r
- arrchisquare0[i]=0.;\r
- arrchisquare1[i]=0.;\r
- arrchisquare2[i]=0.;\r
- arrchisquare3[i]=0.;\r
- arrchisquare4[i]=0.;\r
- arrchisquare5[i]=0.;\r
- }\r
- \r
- TCanvas* canvas[nMultbins];\r
- Int_t nx=2, ny=2;\r
- if(nPtBins>4){\r
- nx=3;\r
- ny=2;\r
- }\r
- if(nPtBins>6){\r
- nx=4;\r
- ny=3;\r
- }\r
- if(nPtBins>12){\r
- nx=5;\r
- ny=4;\r
- }\r
- for(Int_t i=0; i<nMultbins; i++ ){\r
- canvas[i] = new TCanvas(Form("c%d",i),Form("summary canvas for mult bin %d",i));\r
- canvas[i]->Divide(nx,ny);\r
- }\r
- TCanvas *myCanvas[nPtMultbins];\r
- \r
- //\r
- //\r
- // Loop on multiplicity bins\r
- //\r
- Int_t massBin=0;\r
- Double_t sig,errsig,s,errs,b,errb;\r
- for(Int_t j=0; j<nMultbins; j++){\r
- // printf(" Studying multiplicity bin %d\n",j);\r
- Int_t multbinlow = hmultaxis->FindBin(multlims[j]);\r
- Int_t multbinhigh = hmultaxis->FindBin(multlims[j+1])-1;\r
- //\r
- // Loop on pt bins\r
- //\r
- iPad=1;\r
- for(Int_t iBin=0; iBin<nPtBins; iBin++){\r
- canvas[j]->cd(iPad++);\r
- // printf(" projecting to the mass histogram\n");\r
- Int_t ptbinlow = hptaxis->FindBin(ptlims[iBin]);\r
- Int_t ptbinhigh = hptaxis->FindBin(ptlims[iBin+1])-1;\r
- hmass[massBin] = (TH1F*)hPtMassMult[0]->ProjectionY(Form("hmass_%d_%d",j,iBin),multbinlow,multbinhigh,ptbinlow,ptbinhigh);\r
- if( hmass[massBin]->GetEntries() < 60 ) {\r
- massBin++;\r
- continue;\r
- }\r
- Int_t origNbins=hmass[massBin]->GetNbinsX(); \r
- // printf(" rebinning the mass histogram\n");\r
- TH1F* hRebinned=RebinHisto(hmass[massBin],rebin[iBin],firstUsedBin[iBin]);\r
- hmin=hRebinned->GetBinLowEdge(2);\r
- hmax=hRebinned->GetBinLowEdge(hRebinned->GetNbinsX());\r
- // printf(" define the mass fitter and options \n");\r
- fitter[massBin]=new AliHFMassFitter(hRebinned,hmin, hmax,1,typeb,types);\r
- fitter[massBin]->SetRangeFit(1.65,2.10);\r
- rebin[massBin]=origNbins/fitter[iBin]->GetBinN();\r
- fitter[massBin]->SetReflectionSigmaFactor(factor4refl);\r
- fitter[massBin]->SetInitialGaussianMean(massD);\r
- fitter[massBin]->SetInitialGaussianSigma(sigmapt[iBin]);\r
- if(fixPeakSigma) {\r
- fitter[massBin]->SetFixGaussianMean(massD);\r
- fitter[massBin]->SetFixGaussianSigma(sigmapt[iBin],kTRUE);\r
- }\r
- Bool_t out=fitter[massBin]->MassFitter(0);\r
- if(!out) {\r
- fitter[massBin]->GetHistoClone()->Draw();\r
- massBin++;\r
- continue;\r
- }\r
- // printf(" getting the fit parameters\n");\r
- Double_t mass=fitter[massBin]->GetMean();\r
- Double_t sigma=fitter[massBin]->GetSigma();\r
- if(j==0) arrchisquare0[iBin]=fitter[massBin]->GetReducedChiSquare();\r
- else if(j==1) arrchisquare1[iBin]=fitter[massBin]->GetReducedChiSquare();\r
- else if(j==2) arrchisquare2[iBin]=fitter[massBin]->GetReducedChiSquare();\r
- else if(j==3) arrchisquare3[iBin]=fitter[massBin]->GetReducedChiSquare();\r
- else if(j==4) arrchisquare4[iBin]=fitter[massBin]->GetReducedChiSquare();\r
- else if(j==5) arrchisquare5[iBin]=fitter[massBin]->GetReducedChiSquare();\r
- TF1* fB1=fitter[massBin]->GetBackgroundFullRangeFunc();\r
- TF1* fB2=fitter[massBin]->GetBackgroundRecalcFunc();\r
- TF1* fM=fitter[massBin]->GetMassFunc();\r
- if(iBin==0 && fB1) funBckStore1=new TF1(*fB1);\r
- if(iBin==0 && fB2) funBckStore2=new TF1(*fB2);\r
- if(iBin==0 && fM) funBckStore3=new TF1(*fM);\r
-\r
- fitter[massBin]->DrawHere(gPad);\r
- fitter[massBin]->Signal(3,s,errs);\r
- fitter[massBin]->Background(3,b,errb);\r
- fitter[massBin]->Significance(3,sig,errsig);\r
- Double_t ry=fitter[iBin]->GetRawYield();\r
- Double_t ery=fitter[iBin]->GetRawYieldError();\r
- myCanvas[massBin] = new TCanvas(Form("myCanvas_%d_%d",j,iBin),Form("Invariant mass mult bin %d, pt bin %d",j,iBin));\r
- fitter[massBin]->DrawHere(gPad);\r
- \r
- Float_t cntSig1=0.;\r
- Float_t cntSig2=0.;\r
- Float_t cntErr=0.;\r
- for(Int_t iMB=minBinSum; iMB<=maxBinSum; iMB++){\r
- Float_t bkg1=fB1 ? fB1->Eval(hmass[massBin]->GetBinCenter(iMB))/rebin[iBin] : 0;\r
- Float_t bkg2=fB2 ? fB2->Eval(hmass[massBin]->GetBinCenter(iMB))/rebin[iBin] : 0;\r
- cntSig1+=(hmass[massBin]->GetBinContent(iMB)-bkg1);\r
- cntSig2+=(hmass[massBin]->GetBinContent(iMB)-bkg2);\r
- cntErr+=(hmass[massBin]->GetBinContent(iMB));\r
- }\r
- hCntSig1[j]->SetBinContent(iBin+1,cntSig1);\r
- hCntSig1[j]->SetBinError(iBin+1,TMath::Sqrt(cntErr));\r
- hNDiffCntSig1[j]->SetBinContent(iBin+1,(s-cntSig1)/s);\r
- hNDiffCntSig1[j]->SetBinError(iBin+1,TMath::Sqrt(cntErr)/s);\r
- hCntSig2[j]->SetBinContent(iBin+1,cntSig2);\r
- hNDiffCntSig2[j]->SetBinContent(iBin+1,(s-cntSig2)/s);\r
- hNDiffCntSig2[j]->SetBinError(iBin+1,TMath::Sqrt(cntErr)/s);\r
- hCntSig2[j]->SetBinError(iBin+1,TMath::Sqrt(cntErr));\r
- hSignal[j]->SetBinContent(iBin+1,ry);\r
- hSignal[j]->SetBinError(iBin+1,ery);\r
- hRelErrSig[j]->SetBinContent(iBin+1,errs/s);\r
- hInvSignif[j]->SetBinContent(iBin+1,1/sig);\r
- hInvSignif[j]->SetBinError(iBin+1,errsig/(sig*sig));\r
- hBackground[j]->SetBinContent(iBin+1,b); //consider sigma\r
- hBackground[j]->SetBinError(iBin+1,errb);\r
- hBackgroundNormSigma[j]->SetBinContent(iBin+1,b/(3*fitter[iBin]->GetSigma())*(3*0.012)); //consider sigma\r
- hBackgroundNormSigma[j]->SetBinError(iBin+1,errb);\r
- hSignificance[j]->SetBinContent(iBin+1,sig);\r
- hSignificance[j]->SetBinError(iBin+1,errsig);\r
- hMass[j]->SetBinContent(iBin+1,mass);\r
- hMass[j]->SetBinError(iBin+1,0.0001);\r
- hSigma[j]->SetBinContent(iBin+1,sigma);\r
- hSigma[j]->SetBinError(iBin+1,0.0001);\r
-\r
- massBin++;\r
- }// end loop on pt bins\r
-\r
- canvas[j]->Update();\r
- canvas[j]->SaveAs(Form("hMass%s_%d_%d.eps",CutsType,typeb,j));\r
- \r
- }// end loop on multiplicity bins\r
-\r
-\r
- TCanvas *cpar=new TCanvas("cpar","Fit params",1200,600);\r
- cpar->Divide(2,1);\r
- cpar->cd(1);\r
- for(Int_t imult=0; imult<nMultbins; imult++) {\r
- hMass[imult]->SetMarkerStyle(20);\r
- hMass[imult]->GetXaxis()->SetTitle("Pt (GeV/c)");\r
- hMass[imult]->GetYaxis()->SetTitle("Mass (GeV/c^{2})");\r
- hMass[imult]->SetMarkerColor(2*imult);\r
- if(imult==5) hMass[imult]->SetMarkerColor(2*imult-3);\r
- if(imult==0) {\r
- hMass[imult]->SetMarkerColor(kBlack);\r
- hMass[imult]->Draw("PE");\r
- }\r
- else hMass[imult]->Draw("PEsame");\r
- }\r
- cpar->cd(2);\r
- for(Int_t imult=0; imult<nMultbins; imult++) {\r
- hSigma[imult]->SetMarkerStyle(20);\r
- // hSigma[0]->Draw("PE");\r
- hSigma[imult]->GetXaxis()->SetTitle("Pt (GeV/c)");\r
- hSigma[imult]->GetYaxis()->SetTitle("Sigma (GeV/c^{2})");\r
- hSigma[imult]->SetMarkerColor(2*imult);\r
- if(imult==5) hSigma[imult]->SetMarkerColor(2*imult-3);\r
- if(imult==0) {\r
- hSigma[imult]->SetMarkerColor(kBlack);\r
- hSigma[imult]->Draw("PE");\r
- }\r
- else hSigma[imult]->Draw("PEsame");\r
- }\r
- cpar->Update();\r
-\r
- /*\r
- TCanvas** csig;//= new TCanvas*[nMultbins];\r
- TCanvas** cDiffS;//=new TCanvas*[nMultbins];\r
- for(Int_t i=0; i<nMultbins; i++){\r
- csig[i] =new TCanvas(Form("csig_%d",i),Form("Results, multiplicity bin %d",i),1200,600);\r
- csig[i]->Divide(3,1);\r
- csig[i]->cd(1);\r
- hSignal[i]->SetMarkerStyle(20);\r
- hSignal[i]->SetMarkerColor(4);\r
- hSignal[i]->SetLineColor(4);\r
- hSignal[i]->GetXaxis()->SetTitle("Pt (GeV/c)");\r
- hSignal[i]->GetYaxis()->SetTitle("Signal");\r
- hSignal[i]->Draw("P");\r
- hCntSig1[i]->SetMarkerStyle(26);\r
- hCntSig1[i]->SetMarkerColor(2);\r
- hCntSig1[i]->SetLineColor(2);\r
- hCntSig1[i]->Draw("PSAME");\r
- hCntSig2[i]->SetMarkerStyle(29);\r
- hCntSig2[i]->SetMarkerColor(kGray+1);\r
- hCntSig2[i]->SetLineColor(kGray+1);\r
- hCntSig2[i]->Draw("PSAME");\r
- TLegend* leg=new TLegend(0.4,0.7,0.89,0.89);\r
- leg->SetFillColor(0);\r
- TLegendEntry* ent=leg->AddEntry(hSignal[i],"From Fit","PL");\r
- ent->SetTextColor(hSignal[i]->GetMarkerColor());\r
- ent=leg->AddEntry(hCntSig1[i],"From Counting1","PL");\r
- ent->SetTextColor(hCntSig1[i]->GetMarkerColor());\r
- ent=leg->AddEntry(hCntSig2[i],"From Counting2","PL");\r
- ent->SetTextColor(hCntSig2[i]->GetMarkerColor());\r
- leg->Draw();\r
- csig[i]->cd(2);\r
- hBackground[i]->SetMarkerStyle(20);\r
- hBackground[i]->Draw("P");\r
- hBackground[i]->GetXaxis()->SetTitle("Pt (GeV/c)");\r
- hBackground[i]->GetYaxis()->SetTitle("Background");\r
- csig[i]->cd(3);\r
- hSignificance[i]->SetMarkerStyle(20);\r
- hSignificance[i]->Draw("P");\r
- hSignificance[i]->GetXaxis()->SetTitle("Pt (GeV/c)");\r
- hSignificance[i]->GetYaxis()->SetTitle("Significance");\r
- \r
- cDiffS[i]=new TCanvas(Form("cDiffS_%d",i),Form("Difference, multiplicity bin %d",i),1200,600);\r
- cDiffS[i]->Divide(2,1);\r
- cDiffS[i]->cd(1);\r
- hRelErrSig[i]->SetMarkerStyle(20); //fullcircle\r
- hRelErrSig[i]->SetTitleOffset(1.2); \r
- hRelErrSig[i]->SetTitle("Rel Error from Fit;p_{t} (GeV/c);Signal Relative Error");\r
- hRelErrSig[i]->Draw("P");\r
- hInvSignif[i]->SetMarkerStyle(21); //fullsquare\r
- hInvSignif[i]->SetMarkerColor(kMagenta+1);\r
- hInvSignif[i]->SetLineColor(kMagenta+1);\r
- hInvSignif[i]->Draw("PSAMES");\r
- TLegend* leg2=new TLegend(0.4,0.7,0.89,0.89);\r
- leg2->SetFillColor(0);\r
- TLegendEntry* ent2=leg2->AddEntry(hRelErrSig[i],"From Fit","P");\r
- ent2->SetTextColor(hRelErrSig[i]->GetMarkerColor());\r
- ent2=leg2->AddEntry(hInvSignif[i],"1/Significance","PL");\r
- ent2->SetTextColor(hInvSignif[i]->GetMarkerColor());\r
- leg2->Draw();\r
-\r
- cDiffS[i]->cd(2);\r
- hNDiffCntSig1[i]->SetMarkerStyle(26);\r
- hNDiffCntSig1[i]->SetMarkerColor(2);\r
- hNDiffCntSig1[i]->SetLineColor(2);\r
- hNDiffCntSig1[i]->SetTitle("Cmp Fit-Count;p_{t} (GeV/c);(S_{fit}-S_{count})/S_{fit}");\r
- hNDiffCntSig1[i]->Draw("P");\r
- hNDiffCntSig2[i]->SetMarkerStyle(29);\r
- hNDiffCntSig2[i]->SetMarkerColor(kGray+1);\r
- hNDiffCntSig2[i]->SetLineColor(kGray+1);\r
- hNDiffCntSig2[i]->Draw("PSAME");\r
- TLegend* leg1=new TLegend(0.4,0.7,0.89,0.89);\r
- leg1->SetFillColor(0);\r
- TLegendEntry* ent1=leg1->AddEntry(hNDiffCntSig1[i],"From Counting1","PL");\r
- ent1->SetTextColor(hNDiffCntSig1[i]->GetMarkerColor());\r
- ent1=leg1->AddEntry(hNDiffCntSig2[i],"From Counting2","PL");\r
- ent1->SetTextColor(hNDiffCntSig2[i]->GetMarkerColor());\r
- leg1->Draw();\r
- }\r
-*/\r
-\r
- TGraph* grReducedChiSquare0=new TGraph(nPtBins,ptlims,arrchisquare0);\r
- grReducedChiSquare0->SetName("grReducedChiSquare0");\r
- grReducedChiSquare0->SetTitle("Reduced Chi2;p_{t} (GeV/c);#tilde{#chi}^{2}");\r
- TGraph* grReducedChiSquare1=new TGraph(nPtBins,ptlims,arrchisquare1);\r
- grReducedChiSquare1->SetName("grReducedChiSquare1");\r
- grReducedChiSquare1->SetTitle("Reduced Chi2;p_{t} (GeV/c);#tilde{#chi}^{2}");\r
- TGraph* grReducedChiSquare2=new TGraph(nPtBins,ptlims,arrchisquare2);\r
- grReducedChiSquare2->SetName("grReducedChiSquare2");\r
- grReducedChiSquare2->SetTitle("Reduced Chi2;p_{t} (GeV/c);#tilde{#chi}^{2}");\r
- TGraph* grReducedChiSquare3=new TGraph(nPtBins,ptlims,arrchisquare3);\r
- grReducedChiSquare3->SetName("grReducedChiSquare3");\r
- grReducedChiSquare3->SetTitle("Reduced Chi2;p_{t} (GeV/c);#tilde{#chi}^{2}");\r
- TCanvas *cChi2=new TCanvas("cChi2","reduced chi square",600,600);\r
- cChi2->cd();\r
- grReducedChiSquare0->SetMarkerStyle(21);\r
- grReducedChiSquare0->Draw("AP");\r
- grReducedChiSquare1->SetMarkerStyle(22);\r
- grReducedChiSquare1->Draw("Psame");\r
- grReducedChiSquare2->SetMarkerStyle(23);\r
- grReducedChiSquare2->Draw("Psame");\r
- grReducedChiSquare3->SetMarkerStyle(24);\r
- grReducedChiSquare3->Draw("Psame");\r
-\r
- TCanvas* cbkgNormSigma=new TCanvas("cbkgNormSigma","Background normalized to sigma",400,600);\r
- cbkgNormSigma->cd();\r
- for(Int_t i=0; i<nMultbins; i++){\r
- hBackgroundNormSigma[i]->SetMarkerStyle(20);\r
- hBackgroundNormSigma[i]->GetXaxis()->SetTitle("Pt (GeV/c)");\r
- hBackgroundNormSigma[i]->GetYaxis()->SetTitle("Background #times 3 #times 0.012/ (3 #times #sigma)");\r
- hBackgroundNormSigma[i]->SetMarkerColor(2*i);\r
- if(i==5) hBackgroundNormSigma[i]->SetMarkerColor(2*i-3);\r
- if(i==0) { \r
- hBackgroundNormSigma[i]->SetMarkerColor(kBlack);\r
- hBackgroundNormSigma[i]->Draw("PE");\r
- }\r
- else hBackgroundNormSigma[i]->Draw("Psame");\r
- }\r
-\r
-\r
- TString partname="Both";\r
- if(optPartAntiPart==kParticleOnly) {\r
- if(analysisType==kD0toKpi) partname="D0";\r
- if(analysisType==kDplusKpipi) partname="Dplus";\r
- }\r
- if(optPartAntiPart==kAntiParticleOnly) {\r
- if(analysisType==kD0toKpi) partname="D0bar";\r
- if(analysisType==kDplusKpipi) partname="Dminus";\r
- }\r
-\r
- TString outfilename = Form("RawYield_MultInt_%s_%s",partname.Data(),CutsType);\r
- if(fixPeakSigma) outfilename += "_SigmaFixed";\r
- if(typeb==0) outfilename += "_Expo.root";\r
- else if(typeb==1) outfilename += "_Linear.root";\r
- else if(typeb==2) outfilename += "_Pol2.root";\r
-\r
- TFile* outf=new TFile(outfilename,"update");\r
- outf->cd(); \r
- for(Int_t j=0; j<massBin; j++) hmass[j]->Write();\r
- for(Int_t j=0; j<nMultbins; j++){\r
- hMass[j]->Write();\r
- hSigma[j]->Write();\r
- hCntSig1[j]->Write();\r
- hCntSig2[j]->Write();\r
- hNDiffCntSig1[j]->Write();\r
- hNDiffCntSig2[j]->Write();\r
- hSignal[j]->Write();\r
- hRelErrSig[j]->Write();\r
- hInvSignif[j]->Write();\r
- hBackground[j]->Write();\r
- hBackgroundNormSigma[j]->Write();\r
- hSignificance[j]->Write();\r
- }\r
- grReducedChiSquare0->Write();\r
- grReducedChiSquare1->Write();\r
- grReducedChiSquare2->Write();\r
- grReducedChiSquare3->Write();\r
- // hPtMass->Write();\r
- outf->Close();\r
- \r
-}\r
-\r
-//_____________________________________________________________________________________________\r
-Bool_t LoadD0toKpiHistos(TObjArray* listFiles, TH3F** hPtMassMult, TH2F** hNtrZvtx, TH2F** hNtrZvtxCorr, const char *CutsType, Int_t Option)\r
-{\r
- Int_t nFiles=listFiles->GetEntries();\r
- TList **hlist=new TList*[nFiles];\r
- TList **hlistNorm=new TList*[nFiles];\r
- AliRDHFCutsD0toKpi** dcuts=new AliRDHFCutsD0toKpi*[nFiles];\r
-\r
- Int_t nReadFiles=0;\r
- for(Int_t iFile=0; iFile<nFiles; iFile++){\r
- TString fName=((TObjString*)listFiles->At(iFile))->GetString(); \r
- TFile *f=TFile::Open(fName.Data());\r
- if(!f){\r
- printf("ERROR: file %s does not exist\n",fName.Data());\r
- continue;\r
- }\r
- printf("Open File %s\n",f->GetName());\r
-\r
- TString dirname="PWG3_D2H_DMult_D0";\r
- if(optPartAntiPart==kParticleOnly) dirname+="D0";\r
- else if(optPartAntiPart==kAntiParticleOnly) dirname+="D0bar";\r
- dirname += CutsType;\r
- TDirectory *dir = (TDirectory*)f->Get(dirname);\r
- if(!dir){\r
- printf("ERROR: directory %s not found in %s\n",dirname.Data(),fName.Data());\r
- continue;\r
- }\r
-\r
- TString listmassname="coutputD0";\r
- if(optPartAntiPart==kParticleOnly) listmassname+="D0";\r
- else if(optPartAntiPart==kAntiParticleOnly) listmassname+="D0bar";\r
- listmassname += CutsType;\r
- printf("List mass name %s\n",listmassname.Data());\r
- hlist[nReadFiles]=(TList*)dir->Get(listmassname);\r
-\r
- TString listnorm="coutputNormD0";\r
- if(optPartAntiPart==kParticleOnly) listnorm+="D0";\r
- else if(optPartAntiPart==kAntiParticleOnly) listnorm+="D0bar";\r
- listnorm += CutsType;\r
- printf("List norm name %s\n",listnorm.Data());\r
- hlistNorm[nReadFiles]=(TList*)dir->Get(listnorm);\r
-\r
- TString cutsobjname="coutputCutsD0";\r
- if(optPartAntiPart==kParticleOnly) cutsobjname+="D0";\r
- else if(optPartAntiPart==kAntiParticleOnly) cutsobjname+="D0bar";\r
- cutsobjname += CutsType;\r
- printf("Cuts name %s\n",cutsobjname.Data());\r
- dcuts[nReadFiles]=(AliRDHFCutsD0toKpi*)dir->Get(cutsobjname);\r
- if(!dcuts[nReadFiles]) {\r
- printf("ERROR: Cut objects do not match\n");\r
- return kFALSE;\r
- }\r
- /*\r
- if(nReadFiles>0){\r
- Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]);\r
- if(!sameCuts){\r
- printf("ERROR: Cut objects do not match\n");\r
- return kFALSE;\r
- }\r
- }\r
- */\r
- nReadFiles++;\r
- }\r
- if(nReadFiles<nFiles){\r
- printf("WARNING: not all requested files have been found\n");\r
- if (nReadFiles==0) {\r
- printf("ERROR: Any file/dir found\n");\r
- return kFALSE;\r
- }\r
- }\r
- // printf("Cuts type %s, particle/antipart %d\n",CutsType,optPartAntiPart);\r
-\r
- /*\r
- Int_t nPtBinsCuts=dcuts[0]->GetNPtBins();\r
- printf("Number of pt bins for cut object = %d\n",nPtBins);\r
- Float_t *ptlimsCuts=dcuts[0]->GetPtBinLimits();\r
- ptlimsCuts[nPtBinsCuts]=ptlimsCuts[nPtBinsCuts-1]+4.;\r
- */\r
-\r
- printf("Get the 3D histogram \n");\r
- const char *histoname="";\r
- if(optPartAntiPart==kParticleOnly) histoname= "hPtVsMassvsMultPart";\r
- else if(optPartAntiPart==kAntiParticleOnly) histoname="hPtVsMassvsMultAntiPart";\r
- else if(optPartAntiPart==kBoth) histoname="hPtVsMassvsMult";\r
- if(Option==kUnCorr) histoname="hPtVsMassvsMultUncorr";\r
- if(Option==kNoPid) histoname="hPtVsMassvsMultNoPid";\r
-\r
- TH3F * htemp;\r
- for(Int_t iFile=0; iFile<nReadFiles; iFile++){\r
- printf(" Looking for histo histMass %s for file %d\n",histoname,iFile);\r
- htemp=(TH3F*)hlist[iFile]->FindObject(Form("%s",histoname));\r
- // cout << " htemp :"<<htemp<<endl;\r
- if(!hPtMassMult[0]){\r
- hPtMassMult[0]=new TH3F(*htemp);\r
- }else{\r
- hPtMassMult[0]->Add(htemp);\r
- }\r
- hNtrZvtx[iFile] = (TH2F*)hlist[iFile]->FindObject("hNtrVsZvtx");\r
- hNtrZvtxCorr[iFile] = (TH2F*)hlist[iFile]->FindObject("hNtrCorrVsZvtx");\r
- }\r
- \r
- // cout<<" hPtMassMult:"<<hPtMassMult[0]<<endl;\r
-\r
- TString partname="Both";\r
- if(optPartAntiPart==kParticleOnly) partname="D0";\r
- if(optPartAntiPart==kAntiParticleOnly) partname="D0bar";\r
-\r
- TString outfilename = Form("RawYield%s_%s",partname.Data(),CutsType);\r
- if(fixPeakSigma) outfilename += "_SigmaFixed";\r
- if(typeb==0) outfilename += "_Expo.root";\r
- else if(typeb==1) outfilename += "_Linear.root";\r
- else if(typeb==2) outfilename += "_Pol2.root";\r
- TFile* outf=new TFile(outfilename,"recreate");\r
- outf->cd();\r
- dcuts[0]->Write();\r
- outf->Close();\r
- return kTRUE;\r
-\r
-}\r
-\r
-//_____________________________________________________________________________________________\r
-Bool_t LoadDplusHistos(TObjArray* listFiles, TH3F** hPtMassMult, TH2F** hNtrZvtx, TH2F** hNtrZvtxCorr, const char *CutsType, Int_t Option)\r
-{\r
-Int_t nFiles=listFiles->GetEntries();\r
- TList **hlist=new TList*[nFiles];\r
- TList **hlistNorm=new TList*[nFiles];\r
- AliRDHFCutsDplustoKpipi** dcuts=new AliRDHFCutsDplustoKpipi*[nFiles];\r
-\r
- Int_t nReadFiles=0;\r
- for(Int_t iFile=0; iFile<nFiles; iFile++){\r
- TString fName=((TObjString*)listFiles->At(iFile))->GetString(); \r
- TFile *f=TFile::Open(fName.Data());\r
- if(!f){\r
- printf("ERROR: file %s does not exist\n",fName.Data());\r
- continue;\r
- }\r
- printf("Open File %s\n",f->GetName());\r
- TDirectory *dir = (TDirectory*)f->Get(Form("PWG3_D2H_DMult_Dplus%s",suffix.Data()));\r
- // TDirectory *dir = (TDirectory*)f->Get("PWG3_D2H_DMult");\r
- if(!dir){\r
- printf("ERROR: directory PWG3_D2H_DMult not found in %s\n",fName.Data());\r
- continue;\r
- }\r
- hlist[nReadFiles]=(TList*)dir->Get(Form("coutputDplus%s",suffix.Data()));\r
- TList *listcut = (TList*)dir->Get(Form("coutputCutsDplus%s",suffix.Data()));\r
- TList *listNorm = (TList*)dir->Get(Form("coutputNormDplus%s",suffix.Data()));\r
- dcuts[nReadFiles]=(AliRDHFCutsDplustoKpipi*)listcut->At(0);\r
- if(nReadFiles>0){\r
- Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]);\r
- if(!sameCuts){\r
- \r
- printf("ERROR: Cut objects do not match\n");\r
- return kFALSE;\r
- }\r
- }\r
-\r
-\r
-\r
- nReadFiles++;\r
- }\r
- if(nReadFiles<nFiles){\r
- printf("WARNING: not all requested files have been found\n");\r
- if (nReadFiles==0) {\r
- printf("ERROR: Any file/dir found\n");\r
- return kFALSE;\r
- }\r
- }\r
-\r
-\r
- printf("Get the 3D histogram \n");\r
- const char *histoname="";\r
- if(optPartAntiPart==kParticleOnly) histoname= "hPtVsMassvsMultPart";\r
- else if(optPartAntiPart==kAntiParticleOnly) histoname="hPtVsMassvsMultAntiPart";\r
- else if(optPartAntiPart==kBoth) histoname="hPtVsMassvsMult";\r
- if(Option==kUnCorr) histoname="hPtVsMassvsMultUncorr";\r
- if(Option==kNoPid) histoname="hPtVsMassvsMultNoPid";\r
-\r
- TH3F * htemp;\r
- for(Int_t iFile=0; iFile<nReadFiles; iFile++){\r
- printf(" Looking for histo histMass %s for file %d\n",histoname,iFile);\r
- htemp=(TH3F*)hlist[iFile]->FindObject(Form("%s",histoname));\r
- // cout << " htemp :"<<htemp<<endl;\r
- if(!hPtMassMult[0]){\r
- hPtMassMult[0]=new TH3F(*htemp);\r
- }else{\r
- hPtMassMult[0]->Add(htemp);\r
- }\r
- // (TH2F*)hNtrZvtx[iFile] = (TH2F*)hlist[iFile]->FindObject("hNtrVsZvtx"); \r
- //(TH2F*)hNtrZvtxCorr[iFile] = (TH2F*)hlist[iFile]->FindObject("hNtrVsZvtxCorr");\r
- }\r
- \r
- // cout<<" hPtMassMult:"<<hPtMassMult[0]<<endl;\r
-\r
- TString partname="Both";\r
- if(optPartAntiPart==kParticleOnly) partname="Dplus";\r
- if(optPartAntiPart==kAntiParticleOnly) partname="Dminus";\r
-\r
- TString outfilename = Form("RawYield%s_%s",partname.Data(),CutsType);\r
- if(fixPeakSigma) outfilename += "_SigmaFixed";\r
- if(typeb==0) outfilename += "_Expo.root";\r
- else if(typeb==1) outfilename += "_Linear.root";\r
- else if(typeb==2) outfilename += "_Pol2.root";\r
- TFile* outf=new TFile(outfilename,"recreate");\r
- outf->cd();\r
- dcuts[0]->Write();\r
- outf->Close();\r
- return kTRUE;\r
- \r
- }\r
-\r
-//_____________________________________________________________________________________________\r
-TH1F* RebinHisto(TH1F* hOrig, Int_t reb, Int_t firstUse){\r
- // Rebin histogram, from bin firstUse to lastUse\r
- // Use all bins if firstUse=-1\r
- \r
- Int_t nBinOrig=hOrig->GetNbinsX();\r
- Int_t firstBinOrig=1;\r
- Int_t lastBinOrig=nBinOrig;\r
- Int_t nBinOrigUsed=nBinOrig;\r
- Int_t nBinFinal=nBinOrig/reb;\r
- if(firstUse>=1){\r
- firstBinOrig=firstUse;\r
- nBinFinal=(nBinOrig-firstUse+1)/reb;\r
- nBinOrigUsed=nBinFinal*reb;\r
- lastBinOrig=firstBinOrig+nBinOrigUsed-1;\r
- }else{\r
- Int_t exc=nBinOrigUsed%reb;\r
- if(exc!=0){\r
- nBinOrigUsed-=exc;\r
- firstBinOrig+=exc/2;\r
- lastBinOrig=firstBinOrig+nBinOrigUsed-1;\r
- }\r
- }\r
- \r
- printf("Rebin from %d bins to %d bins -- Used bins=%d in range %d-%d\n",nBinOrig,nBinFinal,nBinOrigUsed,firstBinOrig,lastBinOrig);\r
- Float_t lowLim=hOrig->GetXaxis()->GetBinLowEdge(firstBinOrig);\r
- Float_t hiLim=hOrig->GetXaxis()->GetBinUpEdge(lastBinOrig);\r
- TH1F* hRebin=new TH1F(Form("%s-rebin",hOrig->GetName()),hOrig->GetTitle(),nBinFinal,lowLim,hiLim);\r
- Int_t lastSummed=firstBinOrig-1;\r
- for(Int_t iBin=1;iBin<=nBinFinal; iBin++){\r
- Float_t sum=0.;\r
- for(Int_t iOrigBin=0;iOrigBin<reb;iOrigBin++){\r
- sum+=hOrig->GetBinContent(lastSummed+1);\r
- lastSummed++;\r
- }\r
- hRebin->SetBinContent(iBin,sum);\r
- }\r
- return hRebin;\r
-}\r
-\r
-//_____________________________________________________________________________________________\r
-Bool_t CheckNtrVsZvtx(TH2F** hNtrackVsVtxZ, TH2F** hNtrackVsVtxZCorr, Int_t nFiles)\r
-{\r
-\r
- TCanvas *cNtrVsZvtx = new TCanvas("cNtrVsZvtx","Ntr Vs Zvtx");\r
- cNtrVsZvtx->Divide(2,2);\r
- for(Int_t i=0; i<nFiles; i++){\r
- cNtrVsZvtx->cd(i+1);\r
- // hNtrackVsVtxZ[i]->Fit("pol4");\r
- hNtrackVsVtxZ[i]->Draw("colz");\r
- cNtrVsZvtx->Update();\r
- }\r
-\r
- TCanvas *cNtrVsZvtxCorr = new TCanvas("cNtrVsZvtxCorr","Ntr Vs Zvtx Corr");\r
- cNtrVsZvtxCorr->Divide(2,2);\r
- for(Int_t i=0; i<nFiles; i++){\r
- cNtrVsZvtxCorr->cd(i+1);\r
- // hNtrackVsVtxZCorr[i]->Fit("pol4");\r
- hNtrackVsVtxZCorr[i]->Draw("colz");\r
- }\r
-\r
- TH1F *hNtrAxis = (TH1F*)hNtrackVsVtxZ[0]->ProjectionY("hNtrAxis");\r
- TH1F *hZvtx[nMultbins];\r
- Int_t firstbin=0, lastbin=0;\r
- TCanvas *cZvtx = new TCanvas("cZvtx","Zvtx projections");\r
- cZvtx->Divide(2,2);\r
- for(Int_t i=0; i<nFiles; i++){\r
- cZvtx->cd(i+1);\r
- firstbin = hNtrAxis->FindBin( multlims[i] );\r
- lastbin = hNtrAxis->FindBin( multlims[i+1] ) -1;\r
- hZvtx[i] = (TH1F*)hNtrackVsVtxZ[i]->ProjectionX(Form("hZvtx_%d",i),firstbin,lastbin);\r
- hZvtx[i]->Draw();\r
- }\r
- TH1F *hZvtxCorr[nMultbins]; \r
- TCanvas *cZvtxCorr = new TCanvas("cZvtxCorr","Zvtx projections Corr");\r
- cZvtxCorr->Divide(2,2);\r
- for(Int_t i=0; i<nFiles; i++){\r
- cZvtxCorr->cd(i+1);\r
- firstbin = hNtrAxis->FindBin( multlims[i] );\r
- lastbin = hNtrAxis->FindBin( multlims[i+1] ) -1;\r
- hZvtxCorr[i] = (TH1F*)hNtrackVsVtxZCorr[i]->ProjectionX(Form("hZvtxCorr_%d",i),firstbin,lastbin);\r
- hZvtxCorr[i]->Draw();\r
- }\r
-\r
- return kTRUE;\r
-}\r
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <TInterpreter.h>
+#include <TString.h>
+#include <TObjString.h>
+#include <TObjArray.h>
+#include <TMath.h>
+#include <TFile.h>
+#include <TCanvas.h>
+#include <TH1.h>
+#include <TH1F.h>
+#include <TH2F.h>
+#include <TH3.h>
+#include <TH3F.h>
+#include <TH1D.h>
+#include <TF1.h>
+#include <TSystem.h>
+#include <TStyle.h>
+#include <TLegend.h>
+#include <TList.h>
+#include <TLegendEntry.h>
+#include <TDatabasePDG.h>
+#include <TGraph.h>
+
+#include "AliAODRecoDecayHF.h"
+#include "AliRDHFCuts.h"
+#include "AliRDHFCutsDplustoKpipi.h"
+#include "AliRDHFCutsD0toKpi.h"
+#include "AliHFMassFitter.h"
+#include "AliNormalizationCounter.h"
+#endif
+
+/* $Id$ */
+
+enum {kD0toKpi, kDplusKpipi};
+enum {kBoth, kParticleOnly, kAntiParticleOnly};
+enum {kExpo=0, kLinear, kPol2};
+enum {kGaus=0, kDoubleGaus};
+enum {kCorr=0, kUnCorr, kNoPid};
+
+
+// Common variables: to be configured by the user
+// const Int_t nPtBins=6;
+// Double_t ptlims[nPtBins+1]={1., 2.,4.,6.,8.,12.,24.};
+// Int_t rebin[nPtBins]={4,4,6,6,8,8};
+// Double_t sigmapt[nPtBins]={ 0.008, 0.010, 0.012, 0.016, 0.018, 0.020 };
+const Int_t nPtBins=5;
+Double_t ptlims[nPtBins+1]={1.,2.,4.,8.,12.,24.};
+Int_t rebin[nPtBins]={4,4,6,8,8};
+Double_t sigmapt[nPtBins]={ 0.008, 0.014, 0.019, 0.027, 0.033 };
+Bool_t fixPeakSigma = kTRUE;
+//
+const Int_t nMultbins=7;
+Double_t multlims[nMultbins+1]={1.,9.,14.,20.,31.,50.,81.,100.};
+// const Int_t nMultbins=1;
+// Double_t multlims[nMultbins+1]={0.,500.};
+//
+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
+//
+Bool_t isMC=kFALSE;
+Int_t typeb=kExpo;
+Int_t types=kGaus;
+Int_t optPartAntiPart=kBoth;
+Int_t factor4refl=0;
+Float_t massRangeForCounting=0.05; // GeV --> it is 3 sigmapt[binpt]
+Float_t nSigmaRangeForCounting=3.0; // 3 sigmapt[binpt]
+TH2F* hPtMass=0x0;
+TString suffix="StdPid";
+
+
+// Functions
+Bool_t LoadDplusHistos(TObjArray* listFiles, TH3F** hPtMassMult, TH2F** hNtrZvtx, TH2F** hNtrZvtxCorr, const char *CutsType, Int_t Option);
+Bool_t LoadD0toKpiHistos(TObjArray* listFiles, TH3F** hPtMassMult, TH2F** hNtrZvtx, TH2F** hNtrZvtxCorr, AliNormalizationCounter *counter, const char *CutsType, Int_t Option);
+Bool_t CheckNtrVsZvtx(TH2F** hNtrZvtx, TH2F** hNtrZvtxCorr, Int_t nFiles);
+TH1F* RebinHisto(TH1F* hOrig, Int_t reb, Int_t firstUse=-1);
+
+
+void ReadDvsMultiplicity(Int_t analysisType=kD0toKpi,
+ TString fileNameb="AnalysisResults.root",
+ TString fileNamec="",
+ TString fileNamed="",
+ TString fileNamee="",
+ const char *CutsType="",
+ Int_t Option=kCorr)
+{
+ // 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");
+
+ // gInterpreter->ExecuteMacro("$ALICE_ROOT/PWGHF/vertexingHF/macros/LoadLibraries.C");
+ gStyle->SetOptTitle(1);
+
+ TString ntrkname="Ntracklets";
+ // for(int j=0; j<=nMultbins; j++) multlims[j] *= (68./8.8);
+ // ntrkname="Nvzero";
+
+ Int_t nFiles=0;
+ TObjArray* listFiles=new TObjArray();
+ if(fileNameb!="") { listFiles->AddLast(new TObjString(fileNameb.Data())); nFiles++; }
+ if(fileNamec!="") { listFiles->AddLast(new TObjString(fileNamec.Data())); nFiles++; }
+ if(fileNamed!="") { listFiles->AddLast(new TObjString(fileNamed.Data())); nFiles++; }
+ if(fileNamee!="") { listFiles->AddLast(new TObjString(fileNamee.Data())); nFiles++; }
+ if(listFiles->GetEntries()==0){
+ printf("Missing file names in input\n");
+ return;
+ }
+ TH3F** hPtMassMult=new TH3F*[1];
+ hPtMassMult[0]=0x0;
+ TH1F** hmass=new TH1F*[nPtBins*nMultbins];
+ for(Int_t i=0; i<nPtBins*nMultbins; i++) hmass[i]=0x0;
+ TH2F** hNtrZvtx=new TH2F*[4];
+ for(Int_t i=0; i<4; i++) hNtrZvtx[i]=0x0;
+ TH2F** hNtrZvtxCorr=new TH2F*[4];
+ for(Int_t i=0; i<4; i++) hNtrZvtxCorr[i]=0x0;
+ AliNormalizationCounter *counter=0x0;
+ TH1F * hNormalization = new TH1F("hNormalization","Events in the norm counter",nMultbins,multlims);
+
+ Float_t massD;
+ Bool_t retCode;
+ if(analysisType==kD0toKpi) {
+ massD=1.86484003067016602e+00;//(Float_t)(TDatabasePDG::Instance()->GetParticle(421)->Mass());
+ retCode=LoadD0toKpiHistos(listFiles,hPtMassMult,hNtrZvtx,hNtrZvtxCorr,counter,CutsType,Option);
+ }
+ else if(analysisType==kDplusKpipi) {
+ massD=1.86961996555328369e+00;//(Float_t)(TDatabasePDG::Instance()->GetParticle(411)->Mass());
+ retCode=LoadDplusHistos(listFiles,hPtMassMult,hNtrZvtx,hNtrZvtxCorr,CutsType,Option);
+ }
+ else{
+ printf("Wrong analysis type parameter\n");
+ return;
+ }
+ if(!retCode){
+ printf("ERROR in reading input files\n");
+ return;
+ }
+ //
+ // Check the multiplicity variables
+ //
+ Bool_t isMult = CheckNtrVsZvtx(hNtrZvtx,hNtrZvtxCorr,nFiles);
+ if(!isMult) return;
+ //
+ //
+ printf(" Preparing the mass fits");
+ TH1F* hCntSig1[nMultbins];
+ TH1F* hCntSig2[nMultbins];
+ TH1F* hNDiffCntSig1[nMultbins];
+ TH1F* hNDiffCntSig2[nMultbins];
+ TH1F* hSignal[nMultbins];
+ TH1F* hRelErrSig[nMultbins];
+ TH1F* hInvSignif[nMultbins];
+ TH1F* hBackground[nMultbins];
+ TH1F* hBackgroundNormSigma[nMultbins];
+ TH1F* hSignificance[nMultbins];
+ TH1F* hMass[nMultbins];
+ TH1F* hSigma[nMultbins];
+ for(Int_t i=0; i<nMultbins; i++){
+ hCntSig1[i]=new TH1F(Form("hCntSig1_%d",i),Form("hCntSig1_%d",i),nPtBins,ptlims);
+ hCntSig2[i]=new TH1F(Form("hCntSig2_%d",i),Form("hCntSig2_%d",i),nPtBins,ptlims);
+ hNDiffCntSig1[i]=new TH1F(Form("hNDiffCntSig1_%d",i),Form("hNDiffCntSig1_%d",i),nPtBins,ptlims);
+ hNDiffCntSig2[i]=new TH1F(Form("hNDiffCntSig2_%d",i),Form("hNDiffCntSig2_%d",i),nPtBins,ptlims);
+ hSignal[i]=new TH1F(Form("hSignal_%d",i),Form("hSignal_%d",i),nPtBins,ptlims);
+ hRelErrSig[i]=new TH1F(Form("hRelErrSig_%d",i),Form("hRelErrSig_%d",i),nPtBins,ptlims);
+ hInvSignif[i]=new TH1F(Form("hInvSignif_%d",i),Form("hInvSignif_%d",i),nPtBins,ptlims);
+ hBackground[i]=new TH1F(Form("hBackground_%d",i),Form("hBackground_%d",i),nPtBins,ptlims);
+ hBackgroundNormSigma[i]=new TH1F(Form("hBackgroundNormSigma_%d",i),Form("hBackgroundNormSigma_%d",i),nPtBins,ptlims);
+ hSignificance[i]=new TH1F(Form("hSignificance_%d",i),Form("hSignificance_%d",i),nPtBins,ptlims);
+ hMass[i]=new TH1F(Form("hMass_%d",i),Form("hMass_%d",i),nPtBins,ptlims);
+ hSigma[i]=new TH1F(Form("hSigma_%d",i),Form("hSigma_%d",i),nPtBins,ptlims);
+ }
+ printf(", defined...\n");
+
+ // std::cout << " htemp :"<<hPtMassMult[0]<<std::endl;
+ TH1F* hptaxis = (TH1F*)hPtMassMult[0]->ProjectionZ("hptaxis");
+ TH1F* hmassaxis = (TH1F*)hPtMassMult[0]->ProjectionY("hmassaxis");
+ TH1F* hmultaxis = (TH1F*)hPtMassMult[0]->ProjectionX("hmultaxis");
+ Int_t nMassBins=hmassaxis->GetNbinsX();
+ Double_t hmin=hmassaxis->GetBinLowEdge(3);
+ Double_t hmax=hmassaxis->GetBinLowEdge(nMassBins-2) + hmassaxis->GetBinWidth(nMassBins-2);
+ Int_t iPad=1;
+
+ printf("Now initializing the fit functions\n");
+ // TF1* funBckStore1=0x0;
+ // TF1* funBckStore2=0x0;
+ // TF1* funBckStore3=0x0;
+
+ Int_t nPtMultbins = nPtBins*nMultbins;
+ AliHFMassFitter** fitter=new AliHFMassFitter*[nPtMultbins];
+ Double_t arrchisquare0[nPtBins];
+ Double_t arrchisquare1[nPtBins];
+ Double_t arrchisquare2[nPtBins];
+ Double_t arrchisquare3[nPtBins];
+ Double_t arrchisquare4[nPtBins];
+ Double_t arrchisquare5[nPtBins];
+ for(Int_t i=0; i<nPtBins; i++){
+ arrchisquare0[i]=0.;
+ arrchisquare1[i]=0.;
+ arrchisquare2[i]=0.;
+ arrchisquare3[i]=0.;
+ arrchisquare4[i]=0.;
+ arrchisquare5[i]=0.;
+ }
+
+ TCanvas* canvas[nMultbins];
+ Int_t nx=2, ny=2;
+ if(nPtBins>4){
+ nx=3;
+ ny=2;
+ }
+ if(nPtBins>6){
+ nx=4;
+ ny=3;
+ }
+ if(nPtBins>12){
+ nx=5;
+ ny=4;
+ }
+ for(Int_t i=0; i<nMultbins; i++ ){
+ canvas[i] = new TCanvas(Form("c%d",i),Form("summary canvas for mult bin %d",i));
+ canvas[i]->Divide(nx,ny);
+ }
+ TCanvas *myCanvas[nPtMultbins];
+
+ //
+ //
+ // Loop on multiplicity bins
+ //
+ Int_t massBin=0;
+ for(Int_t j=0; j<nMultbins; j++){
+ Double_t sig,errsig,s,errs,b,errb;
+ // printf(" Studying multiplicity bin %d\n",j);
+ Int_t multbinlow = hmultaxis->FindBin(multlims[j]);
+ Int_t multbinhigh = hmultaxis->FindBin(multlims[j+1])-1;
+ Float_t val = multbinlow + (multbinhigh-multbinlow)/2.;
+ Int_t hnbin = hNormalization->FindBin(val);
+ Float_t nevents = 0.;
+ if(counter) { nevents = counter->GetNEventsForNorm(multbinlow,multbinhigh); std::cout << std::endl<<std::endl<<" Nevents ("<<multbinlow<<","<<multbinhigh<<") ="<<nevents<<std::endl<<std::endl<<std::endl;}
+ hNormalization->SetBinContent(hnbin,nevents);
+ //
+ // Loop on pt bins
+ //
+ iPad=1;
+ for(Int_t iBin=0; iBin<nPtBins; iBin++){
+ canvas[j]->cd(iPad++);
+ // printf(" projecting to the mass histogram\n");
+ Int_t ptbinlow = hptaxis->FindBin(ptlims[iBin]);
+ Int_t ptbinhigh = hptaxis->FindBin(ptlims[iBin+1])-1;
+ hmass[massBin] = (TH1F*)hPtMassMult[0]->ProjectionY(Form("hmass_%d_%d",j,iBin),multbinlow,multbinhigh,ptbinlow,ptbinhigh);
+ 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]) );
+ // 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;
+ if( hmass[massBin]->GetEntries() < 60 ) {
+ massBin++;
+ continue;
+ }
+ Int_t origNbins=hmass[massBin]->GetNbinsX();
+ // printf(" rebinning the mass histogram\n");
+ TH1F* hRebinned=RebinHisto(hmass[massBin],rebin[iBin],firstUsedBin[iBin]);
+ hmin=hRebinned->GetBinLowEdge(2);
+ hmax=hRebinned->GetBinLowEdge(hRebinned->GetNbinsX());
+ // printf(" define the mass fitter and options \n");
+ fitter[massBin] = new AliHFMassFitter(hRebinned,hmin, hmax,1,typeb,types);
+ fitter[massBin]->SetRangeFit(1.65,2.10);
+ Int_t rebinItem = origNbins/fitter[massBin]->GetBinN();
+ fitter[massBin]->SetReflectionSigmaFactor(factor4refl);
+ fitter[massBin]->SetInitialGaussianMean(massD);
+ fitter[massBin]->SetInitialGaussianSigma(sigmapt[iBin]);
+ if(fixPeakSigma) {
+ fitter[massBin]->SetFixGaussianMean(massD);
+ fitter[massBin]->SetFixGaussianSigma(sigmapt[iBin],kTRUE);
+ }
+ Bool_t out=fitter[massBin]->MassFitter(0);
+ if(!out) {
+ fitter[massBin]->GetHistoClone()->Draw();
+ massBin++;
+ continue;
+ }
+ // printf(" getting the fit parameters\n");
+ Double_t mass=fitter[massBin]->GetMean();
+ Double_t massUnc=fitter[massBin]->GetMeanUncertainty();
+ Double_t sigma=fitter[massBin]->GetSigma();
+ Double_t sigmaUnc=fitter[massBin]->GetSigmaUncertainty();
+ if(j==0) arrchisquare0[iBin]=fitter[massBin]->GetReducedChiSquare();
+ else if(j==1) arrchisquare1[iBin]=fitter[massBin]->GetReducedChiSquare();
+ else if(j==2) arrchisquare2[iBin]=fitter[massBin]->GetReducedChiSquare();
+ else if(j==3) arrchisquare3[iBin]=fitter[massBin]->GetReducedChiSquare();
+ else if(j==4) arrchisquare4[iBin]=fitter[massBin]->GetReducedChiSquare();
+ else if(j==5) arrchisquare5[iBin]=fitter[massBin]->GetReducedChiSquare();
+ TF1* fB1=fitter[massBin]->GetBackgroundFullRangeFunc();
+ TF1* fB2=fitter[massBin]->GetBackgroundRecalcFunc();
+ // TF1* fM=fitter[massBin]->GetMassFunc();
+ // if(iBin==0 && fB1) funBckStore1=new TF1(*fB1);
+ // if(iBin==0 && fB2) funBckStore2=new TF1(*fB2);
+ // if(iBin==0 && fM) funBckStore3=new TF1(*fM);
+
+ fitter[massBin]->DrawHere(gPad);
+ fitter[massBin]->Signal(3,s,errs);
+ fitter[massBin]->Background(3,b,errb);
+ fitter[massBin]->Significance(3,sig,errsig);
+ Double_t ry=fitter[massBin]->GetRawYield();
+ Double_t ery=fitter[massBin]->GetRawYieldError();
+ myCanvas[massBin] = new TCanvas(Form("myCanvas_%d_%d",j,iBin),Form("Invariant mass mult bin %d, pt bin %d",j,iBin));
+ fitter[massBin]->DrawHere(gPad);
+
+ Float_t cntSig1=0.;
+ Float_t cntSig2=0.;
+ Float_t cntErr=0.;
+ massRangeForCounting = nSigmaRangeForCounting*sigmapt[iBin];
+ // std::cout << " pt bin "<< iBin << " mass range = "<< massRangeForCounting<<std::endl;
+ Float_t minBinSum=hmassaxis->FindBin(massD-massRangeForCounting);
+ Float_t maxBinSum=hmassaxis->FindBin(massD+massRangeForCounting);
+ for(Int_t iMB=minBinSum; iMB<=maxBinSum; iMB++){
+ Float_t bkg1=fB1 ? fB1->Eval(hmass[massBin]->GetBinCenter(iMB))/rebinItem : 0;
+ Float_t bkg2=fB2 ? fB2->Eval(hmass[massBin]->GetBinCenter(iMB))/rebinItem : 0;
+ cntSig1+=(hmass[massBin]->GetBinContent(iMB)-bkg1);
+ cntSig2+=(hmass[massBin]->GetBinContent(iMB)-bkg2);
+ cntErr+=(hmass[massBin]->GetBinContent(iMB));
+ }
+ hCntSig1[j]->SetBinContent(iBin+1,cntSig1);
+ hCntSig1[j]->SetBinError(iBin+1,TMath::Sqrt(cntErr));
+ hNDiffCntSig1[j]->SetBinContent(iBin+1,(s-cntSig1)/s);
+ hNDiffCntSig1[j]->SetBinError(iBin+1,TMath::Sqrt(cntErr)/s);
+ hCntSig2[j]->SetBinContent(iBin+1,cntSig2);
+ hNDiffCntSig2[j]->SetBinContent(iBin+1,(s-cntSig2)/s);
+ hNDiffCntSig2[j]->SetBinError(iBin+1,TMath::Sqrt(cntErr)/s);
+ hCntSig2[j]->SetBinError(iBin+1,TMath::Sqrt(cntErr));
+ hSignal[j]->SetBinContent(iBin+1,ry);
+ hSignal[j]->SetBinError(iBin+1,ery);
+ hRelErrSig[j]->SetBinContent(iBin+1,errs/s);
+ hInvSignif[j]->SetBinContent(iBin+1,1/sig);
+ hInvSignif[j]->SetBinError(iBin+1,errsig/(sig*sig));
+ hBackground[j]->SetBinContent(iBin+1,b); //consider sigma
+ hBackground[j]->SetBinError(iBin+1,errb);
+ hBackgroundNormSigma[j]->SetBinContent(iBin+1,b/(3*fitter[massBin]->GetSigma())*(3*0.012)); //consider sigma
+ hBackgroundNormSigma[j]->SetBinError(iBin+1,errb);
+ hSignificance[j]->SetBinContent(iBin+1,sig);
+ hSignificance[j]->SetBinError(iBin+1,errsig);
+ hMass[j]->SetBinContent(iBin+1,mass);
+ hMass[j]->SetBinError(iBin+1,massUnc);
+ hSigma[j]->SetBinContent(iBin+1,sigma);
+ hSigma[j]->SetBinError(iBin+1,sigmaUnc);
+
+ massBin++;
+ delete hRebinned;
+ }// end loop on pt bins
+
+ canvas[j]->Update();
+ canvas[j]->SaveAs(Form("hMass%s_%d_%d.eps",CutsType,typeb,j));
+ // canvas[j]->SaveAs(Form("hMass%s_%d_%d_MultInt.eps",CutsType,typeb,j));
+
+ }// end loop on multiplicity bins
+
+
+ TCanvas *cpar=new TCanvas("cpar","Fit params",1200,600);
+ cpar->Divide(2,1);
+ cpar->cd(1);
+ for(Int_t imult=0; imult<nMultbins; imult++) {
+ hMass[imult]->SetMarkerStyle(20);
+ hMass[imult]->GetXaxis()->SetTitle("Pt (GeV/c)");
+ hMass[imult]->GetYaxis()->SetTitle("Mass (GeV/c^{2})");
+ hMass[imult]->SetMarkerColor(2*imult);
+ if(imult==5) hMass[imult]->SetMarkerColor(2*imult-3);
+ if(imult==0) {
+ hMass[imult]->SetMarkerColor(kBlack);
+ hMass[imult]->Draw("PE");
+ }
+ else hMass[imult]->Draw("PEsame");
+ }
+ cpar->cd(2);
+ for(Int_t imult=0; imult<nMultbins; imult++) {
+ hSigma[imult]->SetMarkerStyle(20);
+ // hSigma[0]->Draw("PE");
+ hSigma[imult]->GetXaxis()->SetTitle("Pt (GeV/c)");
+ hSigma[imult]->GetYaxis()->SetTitle("Sigma (GeV/c^{2})");
+ hSigma[imult]->SetMarkerColor(2*imult);
+ if(imult==5) hSigma[imult]->SetMarkerColor(2*imult-3);
+ if(imult==0) {
+ hSigma[imult]->SetMarkerColor(kBlack);
+ hSigma[imult]->Draw("PE");
+ }
+ else hSigma[imult]->Draw("PEsame");
+ }
+ cpar->Update();
+
+ /*
+ TCanvas** csig;//= new TCanvas*[nMultbins];
+ TCanvas** cDiffS;//=new TCanvas*[nMultbins];
+ for(Int_t i=0; i<nMultbins; i++){
+ csig[i] =new TCanvas(Form("csig_%d",i),Form("Results, multiplicity bin %d",i),1200,600);
+ csig[i]->Divide(3,1);
+ csig[i]->cd(1);
+ hSignal[i]->SetMarkerStyle(20);
+ hSignal[i]->SetMarkerColor(4);
+ hSignal[i]->SetLineColor(4);
+ hSignal[i]->GetXaxis()->SetTitle("Pt (GeV/c)");
+ hSignal[i]->GetYaxis()->SetTitle("Signal");
+ hSignal[i]->Draw("P");
+ hCntSig1[i]->SetMarkerStyle(26);
+ hCntSig1[i]->SetMarkerColor(2);
+ hCntSig1[i]->SetLineColor(2);
+ hCntSig1[i]->Draw("PSAME");
+ hCntSig2[i]->SetMarkerStyle(29);
+ hCntSig2[i]->SetMarkerColor(kGray+1);
+ hCntSig2[i]->SetLineColor(kGray+1);
+ hCntSig2[i]->Draw("PSAME");
+ TLegend* leg=new TLegend(0.4,0.7,0.89,0.89);
+ leg->SetFillColor(0);
+ TLegendEntry* ent=leg->AddEntry(hSignal[i],"From Fit","PL");
+ ent->SetTextColor(hSignal[i]->GetMarkerColor());
+ ent=leg->AddEntry(hCntSig1[i],"From Counting1","PL");
+ ent->SetTextColor(hCntSig1[i]->GetMarkerColor());
+ ent=leg->AddEntry(hCntSig2[i],"From Counting2","PL");
+ ent->SetTextColor(hCntSig2[i]->GetMarkerColor());
+ leg->Draw();
+ csig[i]->cd(2);
+ hBackground[i]->SetMarkerStyle(20);
+ hBackground[i]->Draw("P");
+ hBackground[i]->GetXaxis()->SetTitle("Pt (GeV/c)");
+ hBackground[i]->GetYaxis()->SetTitle("Background");
+ csig[i]->cd(3);
+ hSignificance[i]->SetMarkerStyle(20);
+ hSignificance[i]->Draw("P");
+ hSignificance[i]->GetXaxis()->SetTitle("Pt (GeV/c)");
+ hSignificance[i]->GetYaxis()->SetTitle("Significance");
+
+ cDiffS[i]=new TCanvas(Form("cDiffS_%d",i),Form("Difference, multiplicity bin %d",i),1200,600);
+ cDiffS[i]->Divide(2,1);
+ cDiffS[i]->cd(1);
+ hRelErrSig[i]->SetMarkerStyle(20); //fullcircle
+ hRelErrSig[i]->SetTitleOffset(1.2);
+ hRelErrSig[i]->SetTitle("Rel Error from Fit;p_{t} (GeV/c);Signal Relative Error");
+ hRelErrSig[i]->Draw("P");
+ hInvSignif[i]->SetMarkerStyle(21); //fullsquare
+ hInvSignif[i]->SetMarkerColor(kMagenta+1);
+ hInvSignif[i]->SetLineColor(kMagenta+1);
+ hInvSignif[i]->Draw("PSAMES");
+ TLegend* leg2=new TLegend(0.4,0.7,0.89,0.89);
+ leg2->SetFillColor(0);
+ TLegendEntry* ent2=leg2->AddEntry(hRelErrSig[i],"From Fit","P");
+ ent2->SetTextColor(hRelErrSig[i]->GetMarkerColor());
+ ent2=leg2->AddEntry(hInvSignif[i],"1/Significance","PL");
+ ent2->SetTextColor(hInvSignif[i]->GetMarkerColor());
+ leg2->Draw();
+
+ cDiffS[i]->cd(2);
+ hNDiffCntSig1[i]->SetMarkerStyle(26);
+ hNDiffCntSig1[i]->SetMarkerColor(2);
+ hNDiffCntSig1[i]->SetLineColor(2);
+ hNDiffCntSig1[i]->SetTitle("Cmp Fit-Count;p_{t} (GeV/c);(S_{fit}-S_{count})/S_{fit}");
+ hNDiffCntSig1[i]->Draw("P");
+ hNDiffCntSig2[i]->SetMarkerStyle(29);
+ hNDiffCntSig2[i]->SetMarkerColor(kGray+1);
+ hNDiffCntSig2[i]->SetLineColor(kGray+1);
+ hNDiffCntSig2[i]->Draw("PSAME");
+ TLegend* leg1=new TLegend(0.4,0.7,0.89,0.89);
+ leg1->SetFillColor(0);
+ TLegendEntry* ent1=leg1->AddEntry(hNDiffCntSig1[i],"From Counting1","PL");
+ ent1->SetTextColor(hNDiffCntSig1[i]->GetMarkerColor());
+ ent1=leg1->AddEntry(hNDiffCntSig2[i],"From Counting2","PL");
+ ent1->SetTextColor(hNDiffCntSig2[i]->GetMarkerColor());
+ leg1->Draw();
+ }
+*/
+
+ TGraph* grReducedChiSquare0=new TGraph(nPtBins,ptlims,arrchisquare0);
+ grReducedChiSquare0->SetName("grReducedChiSquare0");
+ grReducedChiSquare0->SetTitle("Reduced Chi2;p_{t} (GeV/c);#tilde{#chi}^{2}");
+ TGraph* grReducedChiSquare1=new TGraph(nPtBins,ptlims,arrchisquare1);
+ grReducedChiSquare1->SetName("grReducedChiSquare1");
+ grReducedChiSquare1->SetTitle("Reduced Chi2;p_{t} (GeV/c);#tilde{#chi}^{2}");
+ TGraph* grReducedChiSquare2=new TGraph(nPtBins,ptlims,arrchisquare2);
+ grReducedChiSquare2->SetName("grReducedChiSquare2");
+ grReducedChiSquare2->SetTitle("Reduced Chi2;p_{t} (GeV/c);#tilde{#chi}^{2}");
+ TGraph* grReducedChiSquare3=new TGraph(nPtBins,ptlims,arrchisquare3);
+ grReducedChiSquare3->SetName("grReducedChiSquare3");
+ grReducedChiSquare3->SetTitle("Reduced Chi2;p_{t} (GeV/c);#tilde{#chi}^{2}");
+ TCanvas *cChi2=new TCanvas("cChi2","reduced chi square",600,600);
+ cChi2->cd();
+ grReducedChiSquare0->SetMarkerStyle(21);
+ grReducedChiSquare0->Draw("AP");
+ grReducedChiSquare1->SetMarkerStyle(22);
+ grReducedChiSquare1->Draw("Psame");
+ grReducedChiSquare2->SetMarkerStyle(23);
+ grReducedChiSquare2->Draw("Psame");
+ grReducedChiSquare3->SetMarkerStyle(24);
+ grReducedChiSquare3->Draw("Psame");
+
+ TCanvas* cbkgNormSigma=new TCanvas("cbkgNormSigma","Background normalized to sigma",400,600);
+ cbkgNormSigma->cd();
+ for(Int_t i=0; i<nMultbins; i++){
+ hBackgroundNormSigma[i]->SetMarkerStyle(20);
+ hBackgroundNormSigma[i]->GetXaxis()->SetTitle("Pt (GeV/c)");
+ hBackgroundNormSigma[i]->GetYaxis()->SetTitle("Background #times 3 #times 0.012/ (3 #times #sigma)");
+ hBackgroundNormSigma[i]->SetMarkerColor(2*i);
+ if(i==5) hBackgroundNormSigma[i]->SetMarkerColor(2*i-3);
+ if(i==0) {
+ hBackgroundNormSigma[i]->SetMarkerColor(kBlack);
+ hBackgroundNormSigma[i]->Draw("PE");
+ }
+ else hBackgroundNormSigma[i]->Draw("Psame");
+ }
+
+
+ TString partname="Both";
+ if(optPartAntiPart==kParticleOnly) {
+ if(analysisType==kD0toKpi) partname="D0";
+ if(analysisType==kDplusKpipi) partname="Dplus";
+ }
+ if(optPartAntiPart==kAntiParticleOnly) {
+ if(analysisType==kD0toKpi) partname="D0bar";
+ if(analysisType==kDplusKpipi) partname="Dminus";
+ }
+
+ TString outfilename = Form("RawYield_Mult_%s_%s",partname.Data(),CutsType);
+ // outfilename += "_MultInt";
+ if(fixPeakSigma) outfilename += "_SigmaFixed";
+ outfilename += Form("_BCin%1.1fSigma",nSigmaRangeForCounting);
+ if(typeb==0) outfilename += "_Expo.root";
+ else if(typeb==1) outfilename += "_Linear.root";
+ else if(typeb==2) outfilename += "_Pol2.root";
+
+ TFile* outf=new TFile(outfilename,"recreate");
+ outf->cd();
+ hNormalization->Write();
+ for(Int_t j=0; j<massBin; j++) hmass[j]->Write();
+ for(Int_t j=0; j<nMultbins; j++){
+ hMass[j]->Write();
+ hSigma[j]->Write();
+ hCntSig1[j]->Write();
+ hCntSig2[j]->Write();
+ hNDiffCntSig1[j]->Write();
+ hNDiffCntSig2[j]->Write();
+ hSignal[j]->Write();
+ hRelErrSig[j]->Write();
+ hInvSignif[j]->Write();
+ hBackground[j]->Write();
+ hBackgroundNormSigma[j]->Write();
+ hSignificance[j]->Write();
+ }
+ grReducedChiSquare0->Write();
+ grReducedChiSquare1->Write();
+ grReducedChiSquare2->Write();
+ grReducedChiSquare3->Write();
+ // hPtMass->Write();
+ outf->Close();
+
+}
+
+//_____________________________________________________________________________________________
+Bool_t LoadD0toKpiHistos(TObjArray* listFiles, TH3F** hPtMassMult, TH2F** hNtrZvtx, TH2F** hNtrZvtxCorr, AliNormalizationCounter *counter, const char *CutsType, Int_t Option)
+{
+ Int_t nFiles=listFiles->GetEntries();
+ TList **hlist=new TList*[nFiles];
+ TList **hlistNorm=new TList*[nFiles];
+ AliRDHFCutsD0toKpi** dcuts=new AliRDHFCutsD0toKpi*[nFiles];
+
+ Int_t nReadFiles=0;
+ for(Int_t iFile=0; iFile<nFiles; iFile++){
+ TString fName=((TObjString*)listFiles->At(iFile))->GetString();
+ TFile *f=TFile::Open(fName.Data());
+ if(!f){
+ printf("ERROR: file %s does not exist\n",fName.Data());
+ continue;
+ }
+ printf("Open File %s\n",f->GetName());
+
+ TString dirname="PWG3_D2H_DMult_D0";
+ if(optPartAntiPart==kParticleOnly) dirname+="D0";
+ else if(optPartAntiPart==kAntiParticleOnly) dirname+="D0bar";
+ dirname += CutsType;
+ TDirectory *dir = (TDirectory*)f->Get(dirname);
+ if(!dir){
+ printf("ERROR: directory %s not found in %s\n",dirname.Data(),fName.Data());
+ continue;
+ }
+
+ TString listmassname="coutputD0";
+ if(optPartAntiPart==kParticleOnly) listmassname+="D0";
+ else if(optPartAntiPart==kAntiParticleOnly) listmassname+="D0bar";
+ listmassname += CutsType;
+ printf("List mass name %s\n",listmassname.Data());
+ hlist[nReadFiles]=(TList*)dir->Get(listmassname);
+
+ TString listnorm="coutputNormD0";
+ if(optPartAntiPart==kParticleOnly) listnorm+="D0";
+ else if(optPartAntiPart==kAntiParticleOnly) listnorm+="D0bar";
+ listnorm += CutsType;
+ printf("List norm name %s\n",listnorm.Data());
+ hlistNorm[nReadFiles]=(TList*)dir->Get(listnorm);
+ // AliNormalizationCounter *tmpcounter = (AliNormalizationCounter*)hlistNorm[nReadFiles]->FindObject("NormCounterCorrMult");
+ // counter->Add(tmpcounter);
+ // delete tmpcounter;
+ // counter = tmpcounter;
+
+ TString cutsobjname="coutputCutsD0";
+ if(optPartAntiPart==kParticleOnly) cutsobjname+="D0";
+ else if(optPartAntiPart==kAntiParticleOnly) cutsobjname+="D0bar";
+ cutsobjname += CutsType;
+ printf("Cuts name %s\n",cutsobjname.Data());
+ dcuts[nReadFiles]=(AliRDHFCutsD0toKpi*)dir->Get(cutsobjname);
+ if(!dcuts[nReadFiles]) {
+ printf("ERROR: Cut objects do not match\n");
+ return kFALSE;
+ }
+ /*
+ if(nReadFiles>0){
+ Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]);
+ if(!sameCuts){
+ printf("ERROR: Cut objects do not match\n");
+ return kFALSE;
+ }
+ }
+ */
+ nReadFiles++;
+ }
+ if(nReadFiles<nFiles){
+ printf("WARNING: not all requested files have been found\n");
+ if (nReadFiles==0) {
+ printf("ERROR: Any file/dir found\n");
+ return kFALSE;
+ }
+ }
+ // printf("Cuts type %s, particle/antipart %d\n",CutsType,optPartAntiPart);
+
+ /*
+ Int_t nPtBinsCuts=dcuts[0]->GetNPtBins();
+ printf("Number of pt bins for cut object = %d\n",nPtBins);
+ Float_t *ptlimsCuts=dcuts[0]->GetPtBinLimits();
+ ptlimsCuts[nPtBinsCuts]=ptlimsCuts[nPtBinsCuts-1]+4.;
+ */
+
+ printf("Get the 3D histogram \n");
+ const char *histoname="";
+ if(optPartAntiPart==kParticleOnly) histoname= "hPtVsMassvsMultPart";
+ else if(optPartAntiPart==kAntiParticleOnly) histoname="hPtVsMassvsMultAntiPart";
+ else if(optPartAntiPart==kBoth) histoname="hPtVsMassvsMult";
+ if(Option==kUnCorr) histoname="hPtVsMassvsMultUncorr";
+ if(Option==kNoPid) histoname="hPtVsMassvsMultNoPid";
+
+ TH3F * htemp;
+ for(Int_t iFile=0; iFile<nReadFiles; iFile++){
+ printf(" Looking for histo histMass %s for file %d\n",histoname,iFile);
+ htemp=(TH3F*)hlist[iFile]->FindObject(Form("%s",histoname));
+ // cout << " htemp :"<<htemp<<endl;
+ if(!hPtMassMult[0]){
+ hPtMassMult[0]=new TH3F(*htemp);
+ }else{
+ hPtMassMult[0]->Add(htemp);
+ }
+ hNtrZvtx[iFile] = (TH2F*)hlist[iFile]->FindObject("hNtrVsZvtx");
+ hNtrZvtxCorr[iFile] = (TH2F*)hlist[iFile]->FindObject("hNtrCorrVsZvtx");
+ }
+
+ // std::cout<<" hPtMassMult:"<<hPtMassMult[0]<<std::endl;
+
+ TString partname="Both";
+ if(optPartAntiPart==kParticleOnly) partname="D0";
+ if(optPartAntiPart==kAntiParticleOnly) partname="D0bar";
+
+ TString outfilename = Form("RawYield%s_%s",partname.Data(),CutsType);
+ if(fixPeakSigma) outfilename += "_SigmaFixed";
+ if(typeb==0) outfilename += "_Expo.root";
+ else if(typeb==1) outfilename += "_Linear.root";
+ else if(typeb==2) outfilename += "_Pol2.root";
+ TFile* outf=new TFile(outfilename,"recreate");
+ outf->cd();
+ dcuts[0]->Write();
+ outf->Close();
+ return kTRUE;
+
+}
+
+//_____________________________________________________________________________________________
+Bool_t LoadDplusHistos(TObjArray* listFiles, TH3F** hPtMassMult, TH2F** hNtrZvtx, TH2F** hNtrZvtxCorr, const char *CutsType, Int_t Option)
+{
+Int_t nFiles=listFiles->GetEntries();
+ TList **hlist=new TList*[nFiles];
+ TList **hlistNorm=new TList*[nFiles];
+ AliRDHFCutsDplustoKpipi** dcuts=new AliRDHFCutsDplustoKpipi*[nFiles];
+
+ Int_t nReadFiles=0;
+ for(Int_t iFile=0; iFile<nFiles; iFile++){
+ TString fName=((TObjString*)listFiles->At(iFile))->GetString();
+ TFile *f=TFile::Open(fName.Data());
+ if(!f){
+ printf("ERROR: file %s does not exist\n",fName.Data());
+ continue;
+ }
+ printf("Open File %s\n",f->GetName());
+ TDirectory *dir = (TDirectory*)f->Get(Form("PWG3_D2H_DMult_Dplus%s",suffix.Data()));
+ // TDirectory *dir = (TDirectory*)f->Get("PWG3_D2H_DMult");
+ if(!dir){
+ printf("ERROR: directory PWG3_D2H_DMult not found in %s\n",fName.Data());
+ continue;
+ }
+ hlist[nReadFiles]=(TList*)dir->Get(Form("coutputDplus%s",suffix.Data()));
+ TList *listcut = (TList*)dir->Get(Form("coutputCutsDplus%s",suffix.Data()));
+ TList *listNorm = (TList*)dir->Get(Form("coutputNormDplus%s",suffix.Data()));
+ dcuts[nReadFiles]=(AliRDHFCutsDplustoKpipi*)listcut->At(0);
+ if(nReadFiles>0){
+ Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]);
+ if(!sameCuts){
+
+ printf("ERROR: Cut objects do not match\n");
+ return kFALSE;
+ }
+ }
+
+
+
+ nReadFiles++;
+ }
+ if(nReadFiles<nFiles){
+ printf("WARNING: not all requested files have been found\n");
+ if (nReadFiles==0) {
+ printf("ERROR: Any file/dir found\n");
+ return kFALSE;
+ }
+ }
+
+
+ printf("Get the 3D histogram \n");
+ const char *histoname="";
+ if(optPartAntiPart==kParticleOnly) histoname= "hPtVsMassvsMultPart";
+ else if(optPartAntiPart==kAntiParticleOnly) histoname="hPtVsMassvsMultAntiPart";
+ else if(optPartAntiPart==kBoth) histoname="hPtVsMassvsMult";
+ if(Option==kUnCorr) histoname="hPtVsMassvsMultUncorr";
+ if(Option==kNoPid) histoname="hPtVsMassvsMultNoPid";
+
+ TH3F * htemp;
+ for(Int_t iFile=0; iFile<nReadFiles; iFile++){
+ printf(" Looking for histo histMass %s for file %d\n",histoname,iFile);
+ htemp=(TH3F*)hlist[iFile]->FindObject(Form("%s",histoname));
+ // cout << " htemp :"<<htemp<<endl;
+ if(!hPtMassMult[0]){
+ hPtMassMult[0]=new TH3F(*htemp);
+ }else{
+ hPtMassMult[0]->Add(htemp);
+ }
+ // (TH2F*)hNtrZvtx[iFile] = (TH2F*)hlist[iFile]->FindObject("hNtrVsZvtx");
+ //(TH2F*)hNtrZvtxCorr[iFile] = (TH2F*)hlist[iFile]->FindObject("hNtrVsZvtxCorr");
+ }
+
+ // cout<<" hPtMassMult:"<<hPtMassMult[0]<<endl;
+
+ TString partname="Both";
+ if(optPartAntiPart==kParticleOnly) partname="Dplus";
+ if(optPartAntiPart==kAntiParticleOnly) partname="Dminus";
+
+ TString outfilename = Form("RawYield%s_%s",partname.Data(),CutsType);
+ if(fixPeakSigma) outfilename += "_SigmaFixed";
+ if(typeb==0) outfilename += "_Expo.root";
+ else if(typeb==1) outfilename += "_Linear.root";
+ else if(typeb==2) outfilename += "_Pol2.root";
+ TFile* outf=new TFile(outfilename,"recreate");
+ outf->cd();
+ dcuts[0]->Write();
+ outf->Close();
+ return kTRUE;
+
+}
+
+//_____________________________________________________________________________________________
+TH1F* RebinHisto(TH1F* hOrig, Int_t reb, Int_t firstUse){
+ // Rebin histogram, from bin firstUse to lastUse
+ // Use all bins if firstUse=-1
+
+ Int_t nBinOrig=hOrig->GetNbinsX();
+ Int_t firstBinOrig=1;
+ Int_t lastBinOrig=nBinOrig;
+ Int_t nBinOrigUsed=nBinOrig;
+ Int_t nBinFinal=nBinOrig/reb;
+ if(firstUse>=1){
+ firstBinOrig=firstUse;
+ nBinFinal=(nBinOrig-firstUse+1)/reb;
+ nBinOrigUsed=nBinFinal*reb;
+ lastBinOrig=firstBinOrig+nBinOrigUsed-1;
+ }else{
+ Int_t exc=nBinOrigUsed%reb;
+ if(exc!=0){
+ nBinOrigUsed-=exc;
+ firstBinOrig+=exc/2;
+ lastBinOrig=firstBinOrig+nBinOrigUsed-1;
+ }
+ }
+
+ printf("Rebin from %d bins to %d bins -- Used bins=%d in range %d-%d\n",nBinOrig,nBinFinal,nBinOrigUsed,firstBinOrig,lastBinOrig);
+ Float_t lowLim=hOrig->GetXaxis()->GetBinLowEdge(firstBinOrig);
+ Float_t hiLim=hOrig->GetXaxis()->GetBinUpEdge(lastBinOrig);
+ TH1F* hRebin=new TH1F(Form("%s-rebin",hOrig->GetName()),hOrig->GetTitle(),nBinFinal,lowLim,hiLim);
+ Int_t lastSummed=firstBinOrig-1;
+ for(Int_t iBin=1;iBin<=nBinFinal; iBin++){
+ Float_t sum=0.;
+ for(Int_t iOrigBin=0;iOrigBin<reb;iOrigBin++){
+ sum+=hOrig->GetBinContent(lastSummed+1);
+ lastSummed++;
+ }
+ hRebin->SetBinContent(iBin,sum);
+ }
+ return hRebin;
+}
+
+//_____________________________________________________________________________________________
+Bool_t CheckNtrVsZvtx(TH2F** hNtrackVsVtxZ, TH2F** hNtrackVsVtxZCorr, Int_t nFiles)
+{
+
+ TCanvas *cNtrVsZvtx = new TCanvas("cNtrVsZvtx","Ntr Vs Zvtx");
+ cNtrVsZvtx->Divide(2,2);
+ for(Int_t i=0; i<nFiles; i++){
+ cNtrVsZvtx->cd(i+1);
+ // hNtrackVsVtxZ[i]->Fit("pol4");
+ hNtrackVsVtxZ[i]->Draw("colz");
+ cNtrVsZvtx->Update();
+ }
+
+ TCanvas *cNtrVsZvtxCorr = new TCanvas("cNtrVsZvtxCorr","Ntr Vs Zvtx Corr");
+ cNtrVsZvtxCorr->Divide(2,2);
+ for(Int_t i=0; i<nFiles; i++){
+ cNtrVsZvtxCorr->cd(i+1);
+ // hNtrackVsVtxZCorr[i]->Fit("pol4");
+ hNtrackVsVtxZCorr[i]->Draw("colz");
+ }
+
+ TH1F *hNtrAxis = (TH1F*)hNtrackVsVtxZ[0]->ProjectionY("hNtrAxis");
+ TH1F *hZvtx[nMultbins];
+ Int_t firstbin=0, lastbin=0;
+ TCanvas *cZvtx = new TCanvas("cZvtx","Zvtx projections");
+ cZvtx->Divide(2,2);
+ for(Int_t i=0; i<nFiles; i++){
+ cZvtx->cd(i+1);
+ firstbin = hNtrAxis->FindBin( multlims[i] );
+ lastbin = hNtrAxis->FindBin( multlims[i+1] ) -1;
+ hZvtx[i] = (TH1F*)hNtrackVsVtxZ[i]->ProjectionX(Form("hZvtx_%d",i),firstbin,lastbin);
+ hZvtx[i]->Draw();
+ }
+ TH1F *hZvtxCorr[nMultbins];
+ TCanvas *cZvtxCorr = new TCanvas("cZvtxCorr","Zvtx projections Corr");
+ cZvtxCorr->Divide(2,2);
+ for(Int_t i=0; i<nFiles; i++){
+ cZvtxCorr->cd(i+1);
+ firstbin = hNtrAxis->FindBin( multlims[i] );
+ lastbin = hNtrAxis->FindBin( multlims[i+1] ) -1;
+ hZvtxCorr[i] = (TH1F*)hNtrackVsVtxZCorr[i]->ProjectionX(Form("hZvtxCorr_%d",i),firstbin,lastbin);
+ hZvtxCorr[i]->Draw();
+ }
+
+ return kTRUE;
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff