--- /dev/null
+#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=3;\r
+Double_t ptlims[nPtBins+1]={2.,4.,8.,16.};\r
+Int_t rebin[nPtBins]={4,6,6};\r
+Double_t sigmapt[nPtBins]={ 0.012, 0.016, 0.018 };\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 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=kDplusKpipi,\r
+ TString fileNameb="Mult/AnalysisResults.root",\r
+ TString fileNamec="MultC/AnalysisResults.root",\r
+ TString fileNamed="MultD/AnalysisResults.root",\r
+ TString fileNamee="MultE/AnalysisResults.root",\r
+ const char *CutsType="StdPid",\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]);\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]);\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
+ 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,s);\r
+ hSignal[j]->SetBinError(iBin+1,errs);\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]->GetXaxis()->SetTitle("Mass (GeV/c^{2})");\r
+ hMass[imult]->SetMarkerColor(2*imult);\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]->GetXaxis()->SetTitle("Sigma (GeV/c^{2})");\r
+ hSigma[imult]->SetMarkerColor(2*imult);\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==0) { \r
+ hBackgroundNormSigma[i]->SetMarkerColor(kBlack);\r
+ hBackgroundNormSigma[i]->Draw("PA");\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%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<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("hNtrVsZvtxCorr");\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(3,2);\r
+ for(Int_t i=0; i<nFiles; i++){\r
+ cNtrVsZvtx->cd(i);\r
+ // hNtrackVsVtxZ[i]->Fit("pol4");\r
+ hNtrackVsVtxZ[i]->Draw("colz");\r
+ }\r
+\r
+ TCanvas *cNtrVsZvtxCorr = new TCanvas("cNtrVsZvtxCorr","Ntr Vs Zvtx Corr");\r
+ cNtrVsZvtxCorr->Divide(3,2);\r
+ for(Int_t i=0; i<nFiles; i++){\r
+ cNtrVsZvtxCorr->cd(i);\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(3,2);\r
+ for(Int_t i=0; i<nFiles; i++){\r
+ cZvtx->cd(i);\r
+ firstbin = hNtrAxis->FindBin( multlims[i] );\r
+ lastbin = hNtrAxis->FindBin( multlims[i+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(3,2);\r
+ for(Int_t i=0; i<nMultbins; i++){\r
+ cZvtxCorr->cd(i);\r
+ firstbin = hNtrAxis->FindBin( multlims[i] );\r
+ lastbin = hNtrAxis->FindBin( multlims[i+1] );\r
+ hZvtxCorr[i] = (TH1F*)hNtrackVsVtxZCorr[i]->ProjectionX(Form("hZvtx_%d",i),firstbin,lastbin);\r
+ hZvtxCorr[i]->Draw();\r
+ }\r
+ */\r
+ return kTRUE;\r
+}\r
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff