1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 #include <TInterpreter.h>
4 #include <TObjString.h>
15 #include <TLegendEntry.h>
16 #include <TDatabasePDG.h>
19 #include "AliAODRecoDecayHF.h"
20 #include "AliRDHFCuts.h"
21 #include "AliRDHFCutsDplustoKpipi.h"
22 #include "AliRDHFCutsDStartoKpipi.h"
23 #include "AliRDHFCutsDstoKKpi.h"
24 #include "AliRDHFCutsD0toKpi.h"
25 #include "AliHFMassFitter.h"
26 #include "AliNormalizationCounter.h"
30 // MACRO to perform fits to D meson invariant mass spectra
31 // and store raw yields and cut object into a root output file
32 // Origin: F. Prino (prino@to.infn.it)
33 // D0: C. Bianchin (cbianchi@pd.infn.it)
38 enum {kD0toKpi, kDplusKpipi, kDStarD0pi, kDsKKpi};
39 enum {kBoth, kParticleOnly, kAntiParticleOnly};
40 enum {kExpo=0, kThrExpo=5, kLinear, kPol2};
41 enum {kGaus=0, kDoubleGaus};
44 // Common variables: to be configured by the user
45 const Int_t nPtBins=6;
46 Double_t ptlims[nPtBins+1]={2.,3.,4.,5.,6.,8.,12.};
47 Int_t rebin[nPtBins]={2,4,4,4,4,4};
48 Int_t firstUsedBin[nPtBins]={-1,-1,-1,-1,-1,-1}; // -1 uses all bins, >=1 to set the lower bin to be accepted from original histo
50 TString suffix="Loose_SecondSet1236_ForCF08";
53 //const Int_t nPtBins=7;//6;
54 //Double_t ptlims[nPtBins+1]={1.,2.,3.,4.,5.,6.,8.,12.};
55 //Int_t rebin[nPtBins+1]={8,6,10,10,10,10,10,10}; //for looser cuts
56 //Int_t rebin[nPtBins+1]={10,10,10,10,10,10,10,10}; //for Chiara's cuts
59 Int_t optPartAntiPart=kBoth;
61 Double_t minMassForFit=1.7;
62 Double_t maxMassForFit=2.1;
63 Float_t massRangeForCounting=0.05; // GeV
65 Double_t nEventsForNorm=0.;
68 Bool_t cutsappliedondistr=kTRUE;
69 const Int_t nsamples=2;//3;
70 //Int_t nevents[nsamples]={8.5635859e+07 /*LHC10b+cpart*/,8.9700624e+07/*LHC10dpart*/};
71 //Int_t nevents[nsamples]={9.0374946e+07 /*LHC10b+c*/,9.7593785e+07/*LHC10d*/};
72 //Int_t nevents[nsamples]={1.1777812e+08 /*LHC10dnewTPCpid*/,9.7593785e+07/*LHC10d*/};
73 //Int_t nevents[nsamples]={1.1777812e+08 /*LHC10dnewTPCpid*/,9.7593785e+07/*LHC10d*/,9.0374946e+07 /*LHC10b+c*/};
74 Int_t nevents[nsamples]={1.18860695e+08 /*LHC10dnewTPCpid*/,9.0374946e+07 /*LHC10b+c*/};
75 //Int_t nevents[nsamples]={1. /*LHC10dnewTPCpid*/,1 /*LHC10dnewTPCpidrescale*/};
78 Bool_t LoadDplusHistos(TObjArray* listFiles, TH1F** hMass);
79 Bool_t LoadDsHistos(TObjArray* listFiles, TH1F** hMass);
80 Bool_t LoadD0toKpiHistos(TObjArray* listFiles, TH1F** hMass);
81 Bool_t LoadDstarD0piHistos(TObjArray* listFiles, TH1F** hMass);
83 TH1F* RebinHisto(TH1F* hOrig, Int_t reb, Int_t firstUse=-1);
88 void FitMassSpectra(Int_t analysisType=kDplusKpipi,
96 gInterpreter->ExecuteMacro("$ALICE_ROOT/PWGHF/vertexingHF/macros/LoadLibraries.C");
97 gStyle->SetOptTitle(0);
99 TObjArray* listFiles=new TObjArray();
100 if(fileNameb!="") listFiles->AddLast(new TObjString(fileNameb.Data()));
101 if(fileNamec!="") listFiles->AddLast(new TObjString(fileNamec.Data()));
102 if(fileNamed!="") listFiles->AddLast(new TObjString(fileNamed.Data()));
103 if(fileNamee!="") listFiles->AddLast(new TObjString(fileNamee.Data()));
104 if(listFiles->GetEntries()==0){
105 printf("Missing file names in input\n");
109 TH1F** hmass=new TH1F*[nPtBins];
110 for(Int_t i=0;i<nPtBins;i++) hmass[i]=0x0;
112 Float_t massD, massD0_fDstar;
114 if(analysisType==kD0toKpi){
115 retCode=LoadD0toKpiHistos(listFiles,hmass);
116 massD=TDatabasePDG::Instance()->GetParticle(421)->Mass();
118 else if(analysisType==kDplusKpipi){
119 retCode=LoadDplusHistos(listFiles,hmass);
120 massD=TDatabasePDG::Instance()->GetParticle(411)->Mass();
122 else if(analysisType==kDStarD0pi){
123 retCode=LoadDstarD0piHistos(listFiles,hmass);
124 massD=TDatabasePDG::Instance()->GetParticle(413)->Mass();
125 massD0_fDstar=TDatabasePDG::Instance()->GetParticle(421)->Mass();
126 massD =massD-massD0_fDstar;
128 else if(analysisType==kDsKKpi){
129 retCode=LoadDsHistos(listFiles,hmass);
130 massD=TDatabasePDG::Instance()->GetParticle(431)->Mass();
133 printf("Wrong analysis type parameter\n");
137 printf("ERROR in reading input files\n");
143 TH1D* hCntSig1=new TH1D("hCntSig1","hCntSig1",nPtBins,ptlims);
144 TH1D* hCntSig2=new TH1D("hCntSig2","hCntSig2",nPtBins,ptlims);
145 TH1D* hNDiffCntSig1=new TH1D("hNDiffCntSig1","hNDiffCntSig1",nPtBins,ptlims);
146 TH1D* hNDiffCntSig2=new TH1D("hNDiffCntSig2","hNDiffCntSig2",nPtBins,ptlims);
147 TH1D* hSignal=new TH1D("hSignal","hSignal",nPtBins,ptlims);
148 TH1D* hRelErrSig=new TH1D("hRelErrSig","hRelErrSig",nPtBins,ptlims);
149 TH1D* hInvSignif=new TH1D("hInvSignif","hInvSignif",nPtBins,ptlims);
150 TH1D* hBackground=new TH1D("hBackground","hBackground",nPtBins,ptlims);
151 TH1D* hBackgroundNormSigma=new TH1D("hBackgroundNormSigma","hBackgroundNormSigma",nPtBins,ptlims);
152 TH1D* hSignificance=new TH1D("hSignificance","hSignificance",nPtBins,ptlims);
153 TH1D* hMass=new TH1D("hMass","hMass",nPtBins,ptlims);
154 TH1D* hSigma=new TH1D("hSigma","hSigma",nPtBins,ptlims);
157 Int_t nMassBins=hmass[0]->GetNbinsX();
158 Double_t hmin=TMath::Max(minMassForFit,hmass[0]->GetBinLowEdge(2));
159 Double_t hmax=TMath::Min(maxMassForFit,hmass[0]->GetBinLowEdge(nMassBins-2)+hmass[0]->GetBinWidth(nMassBins-2));
160 Float_t minBinSum=hmass[0]->FindBin(massD-massRangeForCounting);
161 Float_t maxBinSum=hmass[0]->FindBin(massD+massRangeForCounting);
164 TF1* funBckStore1=0x0;
165 TF1* funBckStore2=0x0;
166 TF1* funBckStore3=0x0;
168 AliHFMassFitter** fitter=new AliHFMassFitter*[nPtBins];
169 Double_t arrchisquare[nPtBins];
170 TCanvas* c1= new TCanvas("c1","MassSpectra",1000,800);
182 Double_t sig,errsig,s,errs,b,errb;
183 for(Int_t iBin=0; iBin<nPtBins; iBin++){
185 Int_t origNbins=hmass[iBin]->GetNbinsX();
186 TH1F* hRebinned=RebinHisto(hmass[iBin],rebin[iBin],firstUsedBin[iBin]);
187 hmin=TMath::Max(minMassForFit,hRebinned->GetBinLowEdge(2));
188 hmax=TMath::Min(maxMassForFit,hRebinned->GetBinLowEdge(hRebinned->GetNbinsX()));
189 fitter[iBin]=new AliHFMassFitter(hRebinned,hmin, hmax,1,typeb,types);
190 rebin[iBin]=origNbins/fitter[iBin]->GetBinN();
191 fitter[iBin]->SetReflectionSigmaFactor(factor4refl);
192 fitter[iBin]->SetInitialGaussianMean(massD);
193 if(analysisType==kDStarD0pi) fitter[iBin]->SetInitialGaussianSigma(0.00065);
194 Bool_t out=fitter[iBin]->MassFitter(0);
196 fitter[iBin]->GetHistoClone()->Draw();
199 Double_t mass=fitter[iBin]->GetMean();
200 Double_t sigma=fitter[iBin]->GetSigma();
201 arrchisquare[iBin]=fitter[iBin]->GetReducedChiSquare();
202 TF1* fB1=fitter[iBin]->GetBackgroundFullRangeFunc();
203 TF1* fB2=fitter[iBin]->GetBackgroundRecalcFunc();
204 TF1* fM=fitter[iBin]->GetMassFunc();
205 if(iBin==0 && fB1) funBckStore1=new TF1(*fB1);
206 if(iBin==0 && fB2) funBckStore2=new TF1(*fB2);
207 if(iBin==0 && fM) funBckStore3=new TF1(*fM);
209 fitter[iBin]->DrawHere(gPad);
210 fitter[iBin]->Signal(3,s,errs);
211 fitter[iBin]->Background(3,b,errb);
212 fitter[iBin]->Significance(3,sig,errsig);
213 Double_t ry=fitter[iBin]->GetRawYield();
214 Double_t ery=fitter[iBin]->GetRawYieldError();
218 for(Int_t iMB=minBinSum; iMB<=maxBinSum; iMB++){
219 Float_t bkg1=fB1 ? fB1->Eval(hmass[iBin]->GetBinCenter(iMB))/rebin[iBin] : 0;
220 Float_t bkg2=fB2 ? fB2->Eval(hmass[iBin]->GetBinCenter(iMB))/rebin[iBin] : 0;
221 cntSig1+=(hmass[iBin]->GetBinContent(iMB)-bkg1);
222 cntSig2+=(hmass[iBin]->GetBinContent(iMB)-bkg2);
223 cntErr+=(hmass[iBin]->GetBinContent(iMB));
225 hCntSig1->SetBinContent(iBin+1,cntSig1);
226 hCntSig1->SetBinError(iBin+1,TMath::Sqrt(cntErr));
227 hNDiffCntSig1->SetBinContent(iBin+1,(s-cntSig1)/s);
228 hNDiffCntSig1->SetBinError(iBin+1,TMath::Sqrt(cntErr)/s);
229 hCntSig2->SetBinContent(iBin+1,cntSig2);
230 hNDiffCntSig2->SetBinContent(iBin+1,(s-cntSig2)/s);
231 hNDiffCntSig2->SetBinError(iBin+1,TMath::Sqrt(cntErr)/s);
232 hCntSig2->SetBinError(iBin+1,TMath::Sqrt(cntErr));
233 hSignal->SetBinContent(iBin+1,ry);
234 hSignal->SetBinError(iBin+1,ery);
235 hRelErrSig->SetBinContent(iBin+1,errs/s);
236 hInvSignif->SetBinContent(iBin+1,1/sig);
237 hInvSignif->SetBinError(iBin+1,errsig/(sig*sig));
238 hBackground->SetBinContent(iBin+1,b); //consider sigma
239 hBackground->SetBinError(iBin+1,errb);
240 hBackgroundNormSigma->SetBinContent(iBin+1,b/(3*fitter[iBin]->GetSigma())*(3*0.012)); //consider sigma
241 hBackgroundNormSigma->SetBinError(iBin+1,errb);
242 hSignificance->SetBinContent(iBin+1,sig);
243 hSignificance->SetBinError(iBin+1,errsig);
244 hMass->SetBinContent(iBin+1,mass);
245 hMass->SetBinError(iBin+1,0.0001);
246 hSigma->SetBinContent(iBin+1,sigma);
247 hSigma->SetBinError(iBin+1,0.0001);
252 c1->cd(1); // is some cases the fitting function of 1st bin get lost
253 funBckStore1->Draw("same");
254 funBckStore2->Draw("same");
255 funBckStore3->Draw("same");
258 TCanvas *cpar=new TCanvas("cpar","Fit params",1200,600);
261 hMass->SetMarkerStyle(20);
263 hMass->GetXaxis()->SetTitle("Pt (GeV/c)");
264 hMass->GetXaxis()->SetTitle("Mass (GeV/c^{2})");
266 hSigma->SetMarkerStyle(20);
268 hSigma->GetXaxis()->SetTitle("Pt (GeV/c)");
269 hSigma->GetXaxis()->SetTitle("Sigma (GeV/c^{2})");
271 TCanvas* csig=new TCanvas("csig","Results",1200,600);
274 hSignal->SetMarkerStyle(20);
275 hSignal->SetMarkerColor(4);
276 hSignal->SetLineColor(4);
277 hSignal->GetXaxis()->SetTitle("Pt (GeV/c)");
278 hSignal->GetYaxis()->SetTitle("Signal");
280 hCntSig1->SetMarkerStyle(26);
281 hCntSig1->SetMarkerColor(2);
282 hCntSig1->SetLineColor(2);
283 hCntSig1->Draw("PSAME");
284 hCntSig2->SetMarkerStyle(29);
285 hCntSig2->SetMarkerColor(kGray+1);
286 hCntSig2->SetLineColor(kGray+1);
287 hCntSig2->Draw("PSAME");
288 TLegend* leg=new TLegend(0.4,0.7,0.89,0.89);
289 leg->SetFillColor(0);
290 TLegendEntry* ent=leg->AddEntry(hSignal,"From Fit","PL");
291 ent->SetTextColor(hSignal->GetMarkerColor());
292 ent=leg->AddEntry(hCntSig1,"From Counting1","PL");
293 ent->SetTextColor(hCntSig1->GetMarkerColor());
294 ent=leg->AddEntry(hCntSig2,"From Counting2","PL");
295 ent->SetTextColor(hCntSig2->GetMarkerColor());
298 hBackground->SetMarkerStyle(20);
299 hBackground->Draw("P");
300 hBackground->GetXaxis()->SetTitle("Pt (GeV/c)");
301 hBackground->GetYaxis()->SetTitle("Background");
303 hSignificance->SetMarkerStyle(20);
304 hSignificance->Draw("P");
305 hSignificance->GetXaxis()->SetTitle("Pt (GeV/c)");
306 hSignificance->GetYaxis()->SetTitle("Significance");
308 TCanvas* cDiffS=new TCanvas("cDiffS","Difference",1200,600);
311 hRelErrSig->SetMarkerStyle(20); //fullcircle
312 hRelErrSig->SetTitleOffset(1.2);
313 hRelErrSig->SetTitle("Rel Error from Fit;p_{t} (GeV/c);Signal Relative Error");
314 hRelErrSig->Draw("P");
315 hInvSignif->SetMarkerStyle(21); //fullsquare
316 hInvSignif->SetMarkerColor(kMagenta+1);
317 hInvSignif->SetLineColor(kMagenta+1);
318 hInvSignif->Draw("PSAMES");
319 TLegend* leg2=new TLegend(0.4,0.7,0.89,0.89);
320 leg2->SetFillColor(0);
321 TLegendEntry* ent2=leg2->AddEntry(hRelErrSig,"From Fit","P");
322 ent2->SetTextColor(hRelErrSig->GetMarkerColor());
323 ent2=leg2->AddEntry(hInvSignif,"1/Significance","PL");
324 ent2->SetTextColor(hInvSignif->GetMarkerColor());
328 hNDiffCntSig1->SetMarkerStyle(26);
329 hNDiffCntSig1->SetMarkerColor(2);
330 hNDiffCntSig1->SetLineColor(2);
331 hNDiffCntSig1->SetTitle("Cmp Fit-Count;p_{t} (GeV/c);(S_{fit}-S_{count})/S_{fit}");
332 hNDiffCntSig1->Draw("P");
333 hNDiffCntSig2->SetMarkerStyle(29);
334 hNDiffCntSig2->SetMarkerColor(kGray+1);
335 hNDiffCntSig2->SetLineColor(kGray+1);
336 hNDiffCntSig2->Draw("PSAME");
337 TLegend* leg1=new TLegend(0.4,0.7,0.89,0.89);
338 leg1->SetFillColor(0);
339 TLegendEntry* ent1=leg1->AddEntry(hNDiffCntSig1,"From Counting1","PL");
340 ent1->SetTextColor(hNDiffCntSig1->GetMarkerColor());
341 ent1=leg1->AddEntry(hNDiffCntSig2,"From Counting2","PL");
342 ent1->SetTextColor(hNDiffCntSig2->GetMarkerColor());
345 TGraph* grReducedChiSquare=new TGraph(nPtBins,ptlims,arrchisquare);
346 grReducedChiSquare->SetName("grReducedChiSquare");
347 grReducedChiSquare->SetTitle("Reduced Chi2;p_{t} (GeV/c);#tilde{#chi}^{2}");
348 TCanvas *cChi2=new TCanvas("cChi2","reduced chi square",600,600);
350 grReducedChiSquare->SetMarkerStyle(21);
351 grReducedChiSquare->Draw("AP");
353 TCanvas* cbkgNormSigma=new TCanvas("cbkgNormSigma","Background normalized to sigma",400,600);
355 hBackgroundNormSigma->SetMarkerStyle(20);
356 hBackgroundNormSigma->Draw("P");
357 hBackgroundNormSigma->GetXaxis()->SetTitle("Pt (GeV/c)");
358 hBackgroundNormSigma->GetYaxis()->SetTitle("Background #times 3 #times 0.012/ (3 #times #sigma)");
359 hBackgroundNormSigma->Draw();
362 TString partname="Both";
363 if(optPartAntiPart==kParticleOnly) {
364 if(analysisType==kD0toKpi) partname="D0";
365 if(analysisType==kDplusKpipi) partname="Dplus";
366 if(analysisType==kDsKKpi) partname="Dsplus";
368 if(optPartAntiPart==kAntiParticleOnly) {
369 if(analysisType==kD0toKpi) partname="D0bar";
370 if(analysisType==kDplusKpipi) partname="Dminus";
371 if(analysisType==kDsKKpi) partname="Dsminus";
374 printf("Events for norm = %f\n",nEventsForNorm);
375 TH1F* hNEvents=new TH1F("hNEvents","",1,0.,1.);
376 hNEvents->SetBinContent(1,nEventsForNorm);
378 TFile* outf=new TFile(Form("RawYield%s.root",partname.Data()),"update");
385 hNDiffCntSig1->Write();
386 hNDiffCntSig2->Write();
390 hBackground->Write();
391 hBackgroundNormSigma->Write();
392 hSignificance->Write();
393 grReducedChiSquare->Write();
399 Bool_t LoadDplusHistos(TObjArray* listFiles, TH1F** hMass){
401 Int_t nFiles=listFiles->GetEntries();
402 TList **hlist=new TList*[nFiles];
403 AliRDHFCutsDplustoKpipi** dcuts=new AliRDHFCutsDplustoKpipi*[nFiles];
406 for(Int_t iFile=0; iFile<nFiles; iFile++){
407 TString fName=((TObjString*)listFiles->At(iFile))->GetString();
408 TFile *f=TFile::Open(fName.Data());
410 printf("ERROR: file %s does not exist\n",fName.Data());
413 printf("Open File %s\n",f->GetName());
414 TDirectory *dir = (TDirectory*)f->Get("PWG3_D2H_InvMassDplus");
416 printf("ERROR: directory PWG3_D2H_InvMassDplus not found in %s\n",fName.Data());
419 hlist[nReadFiles]=(TList*)dir->Get(Form("coutputDplus%s",suffix.Data()));
420 TList *listcut = (TList*)dir->Get(Form("coutputDplusCuts%s",suffix.Data()));
421 dcuts[nReadFiles]=(AliRDHFCutsDplustoKpipi*)listcut->At(0);
423 Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]);
425 printf("ERROR: Cut objects do not match\n");
429 AliNormalizationCounter* c=(AliNormalizationCounter*)dir->Get(Form("coutputDplusNorm%s",suffix.Data()));
430 printf("Events for normalization = %f\n",c->GetNEventsForNorm());
431 nEventsForNorm+=c->GetNEventsForNorm();
434 if(nReadFiles<nFiles){
435 printf("WARNING: not all requested files have been found\n");
438 Int_t nPtBinsCuts=dcuts[0]->GetNPtBins();
439 printf("Number of pt bins for cut object = %d\n",nPtBins);
440 Float_t *ptlimsCuts=dcuts[0]->GetPtBinLimits();
441 ptlimsCuts[nPtBinsCuts]=ptlimsCuts[nPtBinsCuts-1]+4.;
444 for(Int_t i=0;i<nPtBinsCuts;i++){
445 if(ptlimsCuts[i]>=ptlims[iFinBin+1]) iFinBin+=1;
446 if(iFinBin>nPtBins) break;
447 if(ptlimsCuts[i]>=ptlims[iFinBin] &&
448 ptlimsCuts[i+1]<=ptlims[iFinBin+1]){
449 for(Int_t iFile=0; iFile<nReadFiles; iFile++){
451 if(optPartAntiPart==kBoth) histoName.Form("hMassPt%dTC",i);
452 else if(optPartAntiPart==kParticleOnly) histoName.Form("hMassPt%dTCPlus",i);
453 else if(optPartAntiPart==kAntiParticleOnly) histoName.Form("hMassPt%dTCMinus",i);
454 TH1F* htemp=(TH1F*)hlist[iFile]->FindObject(histoName.Data());
456 printf("ERROR: Histogram %s not found\n",histoName.Data());
460 hMass[iFinBin]=new TH1F(*htemp);
462 hMass[iFinBin]->Add(htemp);
467 TString partname="Both";
468 if(optPartAntiPart==kParticleOnly) partname="Dplus";
469 if(optPartAntiPart==kAntiParticleOnly) partname="Dminus";
471 for(Int_t iFile=0; iFile<nReadFiles; iFile++){
472 TH2F* htemp2=(TH2F*)hlist[iFile]->FindObject("hPtVsMassTC");
474 hPtMass=new TH2F(*htemp2);
476 hPtMass->Add(htemp2);
480 TFile* outf=new TFile(Form("RawYield%s.root",partname.Data()),"recreate");
490 Bool_t LoadDsHistos(TObjArray* listFiles, TH1F** hMass){
492 Int_t nFiles=listFiles->GetEntries();
493 TList **hlist=new TList*[nFiles];
494 AliRDHFCutsDstoKKpi** dcuts=new AliRDHFCutsDstoKKpi*[nFiles];
497 for(Int_t iFile=0; iFile<nFiles; iFile++){
498 TString fName=((TObjString*)listFiles->At(iFile))->GetString();
499 TFile *f=TFile::Open(fName.Data());
501 printf("ERROR: file %s does not exist\n",fName.Data());
504 printf("Open File %s\n",f->GetName());
505 TDirectory *dir = (TDirectory*)f->Get("PWG3_D2H_InvMassDs");
507 printf("ERROR: directory PWG3_D2H_InvMassDs not found in %s\n",fName.Data());
510 hlist[nReadFiles]=(TList*)dir->Get("coutputDs");
511 TList *listcut = (TList*)dir->Get("coutputDsCuts");
512 dcuts[nReadFiles]=(AliRDHFCutsDstoKKpi*)listcut->At(0);
513 cout<< dcuts[nReadFiles]<<endl;
515 Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]);
517 printf("ERROR: Cut objects do not match\n");
523 if(nReadFiles<nFiles){
524 printf("WARNING: not all requested files have been found\n");
527 Int_t nPtBinsCuts=dcuts[0]->GetNPtBins();
528 printf("Number of pt bins for cut object = %d\n",nPtBins);
529 Float_t *ptlimsCuts=dcuts[0]->GetPtBinLimits();
530 ptlimsCuts[nPtBinsCuts]=ptlimsCuts[nPtBinsCuts-1]+4.;
533 for(Int_t i=0;i<nPtBinsCuts;i++){
534 if(ptlimsCuts[i]>=ptlims[iFinBin+1]) iFinBin+=1;
535 if(iFinBin>nPtBins) break;
536 if(ptlimsCuts[i]>=ptlims[iFinBin] &&
537 ptlimsCuts[i+1]<=ptlims[iFinBin+1]){
538 for(Int_t iFile=0; iFile<nReadFiles; iFile++){
540 if(optPartAntiPart==kBoth) histoName.Form("hMassAllPt%dphi",i);
541 else if(optPartAntiPart==kParticleOnly){
542 printf("Particle/Antiparticle not yet enabled for Ds");
543 histoName.Form("hMassAllPt%dphi",i);
545 else if(optPartAntiPart==kAntiParticleOnly){
546 printf("Particle/Antiparticle not yet enabled for Ds");
547 histoName.Form("hMassAllPt%dphi",i);
549 TH1F* htemp=(TH1F*)hlist[iFile]->FindObject(histoName.Data());
551 printf("ERROR: Histogram %s not found\n",histoName.Data());
555 hMass[iFinBin]=new TH1F(*htemp);
557 hMass[iFinBin]->Add(htemp);
562 TString partname="Both";
563 if(optPartAntiPart==kParticleOnly) partname="Both";
564 if(optPartAntiPart==kAntiParticleOnly) partname="Both";
566 for(Int_t iFile=0; iFile<nReadFiles; iFile++){
567 TH2F* htemp2=(TH2F*)hlist[iFile]->FindObject("hPtVsMassPhi");
569 hPtMass=new TH2F(*htemp2);
571 hPtMass->Add(htemp2);
575 TFile* outf=new TFile(Form("RawYield%s.root",partname.Data()),"recreate");
585 Bool_t LoadD0toKpiHistos(TObjArray* listFiles, TH1F** hMass){
587 Int_t nFiles=listFiles->GetEntries();
588 TList **hlist=new TList*[nFiles];
589 AliRDHFCutsD0toKpi** dcuts=new AliRDHFCutsD0toKpi*[nFiles];
592 for(Int_t iFile=0; iFile<nFiles; iFile++){
593 TString fName=((TObjString*)listFiles->At(iFile))->GetString();
594 TFile *f=TFile::Open(fName.Data());
596 printf("ERROR: file %s does not exist\n",fName.Data());
599 printf("Open File %s\n",f->GetName());
601 TString dirname="PWG3_D2H_D0InvMass";
602 if(optPartAntiPart==kParticleOnly) dirname+="D0";
603 else if(optPartAntiPart==kAntiParticleOnly) dirname+="D0bar";
604 if(cutsappliedondistr) dirname+="C";
605 TDirectory *dir = (TDirectory*)f->Get(dirname);
607 printf("ERROR: directory %s not found in %s\n",dirname.Data(),fName.Data());
610 TString listmassname="coutputmassD0Mass";
611 if(optPartAntiPart==kParticleOnly) listmassname+="D0";
612 else if(optPartAntiPart==kAntiParticleOnly) listmassname+="D0bar";
613 if(cutsappliedondistr) listmassname+="C";
615 hlist[nReadFiles]=(TList*)dir->Get(listmassname);
617 TString cutsobjname="cutsD0";
618 if(optPartAntiPart==kParticleOnly) cutsobjname+="D0";
619 else if(optPartAntiPart==kAntiParticleOnly) cutsobjname+="D0bar";
620 if(cutsappliedondistr) cutsobjname+="C";
622 dcuts[nReadFiles]=(AliRDHFCutsD0toKpi*)dir->Get(cutsobjname);
624 Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]);
626 printf("ERROR: Cut objects do not match\n");
632 if(nReadFiles<nFiles){
633 printf("WARNING: not all requested files have been found\n");
635 printf("ERROR: Any file/dir found\n");
640 Int_t nPtBinsCuts=dcuts[0]->GetNPtBins();
641 printf("Number of pt bins for cut object = %d\n",nPtBins);
642 Float_t *ptlimsCuts=dcuts[0]->GetPtBinLimits();
643 ptlimsCuts[nPtBinsCuts]=ptlimsCuts[nPtBinsCuts-1]+4.;
646 for(Int_t i=0;i<nPtBinsCuts;i++){
647 if(ptlimsCuts[i]>=ptlims[iFinBin+1]) iFinBin+=1;
648 if(iFinBin>nPtBins) break;
649 if(ptlimsCuts[i]>=ptlims[iFinBin] &&
650 ptlimsCuts[i+1]<=ptlims[iFinBin+1]){
651 for(Int_t iFile=0; iFile<nReadFiles; iFile++){
652 TH1F* htemp=(TH1F*)hlist[iFile]->FindObject(Form("histMass_%d",i));
654 hMass[iFinBin]=new TH1F(*htemp);
656 hMass[iFinBin]->Add(htemp);
661 TString partname="Both";
662 if(optPartAntiPart==kParticleOnly) partname="D0";
663 if(optPartAntiPart==kAntiParticleOnly) partname="D0bar";
665 TFile* outf=new TFile(Form("RawYield%s.root",partname.Data()),"recreate");
672 Bool_t LoadDstarD0piHistos(TObjArray* listFiles, TH1F** hMass){
674 Int_t nFiles=listFiles->GetEntries();
675 TList **hlist=new TList*[nFiles];
676 AliRDHFCutsDStartoKpipi** dcuts=new AliRDHFCutsDStartoKpipi*[nFiles];
678 for(Int_t iFile=0; iFile<nFiles; iFile++){
679 TString fName=((TObjString*)listFiles->At(iFile))->GetString();
680 TFile *f=TFile::Open(fName.Data());
682 printf("ERROR: file %s does not exist\n",fName.Data());
685 printf("Open File %s\n",f->GetName());
686 TString dirname="PWG3_D2H_DStarSpectra";
687 TDirectory *dir = (TDirectory*)f->Get(dirname);
689 printf("ERROR: directory %s not found in %s\n",dirname.Data(),fName.Data());
692 TString listmassname="DStarPID00";
694 hlist[nReadFiles]=(TList*)dir->Get(listmassname);
695 TString cutsobjname="DStartoKpipiCuts";
696 dcuts[nReadFiles]=(AliRDHFCutsDStartoKpipi*)dir->Get(cutsobjname);
698 Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]);
700 printf("ERROR: Cut objects do not match\n");
706 if(nReadFiles<nFiles){
707 printf("WARNING: not all requested files have been found\n");
709 printf("ERROR: Any file/dir found\n");
713 Int_t nPtBinsCuts=dcuts[0]->GetNPtBins();
714 printf("Number of pt bins for cut object = %d\n",nPtBins);
715 Float_t *ptlimsCuts=dcuts[0]->GetPtBinLimits();
716 ptlimsCuts[nPtBinsCuts]=ptlimsCuts[nPtBinsCuts-1]+4.;
722 for(Int_t i=0;i<nPtBinsCuts;i++){
723 if(ptlimsCuts[i]>=ptlims[iFinBin+1]) iFinBin+=1;
724 if(iFinBin>nPtBins) break;
725 if(ptlimsCuts[i]>=ptlims[iFinBin] &&
726 ptlimsCuts[i+1]<=ptlims[iFinBin+1]){
727 for(Int_t iFile=0; iFile<nReadFiles; iFile++){
728 TH1F* htemp=(TH1F*)hlist[iFile]->FindObject(Form("histDeltaMass_%d",i));
730 hMass[iFinBin]=new TH1F(*htemp);
732 hMass[iFinBin]->Add(htemp);
737 TString partname="Both";
738 if(optPartAntiPart==kParticleOnly) partname="DStar";
739 if(optPartAntiPart==kAntiParticleOnly) partname="DStarbar";
741 TFile* outf=new TFile(Form("RawYield%s.root",partname.Data()),"recreate");
748 void CompareFitTypes(TString* paths, TString* legtext,Int_t ncmp=3,TString* filenameYield=0x0){
749 //read ncmp RawYield.roots and draw them together
750 //set the global variable nevents before running
752 // - paths= vector of ncmp dimension with the paths of the file RawYield.root
753 // - legtext= vector of ncmp dimension with the label for the legend
754 // - ncmp= number of files to compare (default is 3)
755 // - filenameYield= optional ncmp-dimensional array with the difference between the names of the files to be compared (useful if the 2 files are in the same directory but have different names)
757 gStyle->SetOptStat(0);
758 gStyle->SetOptTitle(0);
759 gStyle->SetCanvasColor(0);
760 gStyle->SetFrameFillColor(0);
761 gStyle->SetTitleFillColor(0);
762 gStyle->SetFrameBorderMode(0);
764 if(!filenameYield) filenameYield=new TString[ncmp];
766 for(Int_t k=0;k<ncmp;k++){
767 if(!filenameYield) filenameYield[k]="RawYield.root";
768 filenameYield[k].Prepend(paths[k]);
771 TCanvas* cSig=new TCanvas("cSig","Results",1200,600);
773 TCanvas* cBkgN=new TCanvas("cBkgN","Background normalized to sigma",400,600);
774 TCanvas* cDiffS=new TCanvas("cDiffS","Difference",1200,600);
776 TCanvas *cChi2=new TCanvas("cChi2","reduced chi square",600,600);
778 TLegend* leg1=new TLegend(0.4,0.7,0.89,0.89);
779 leg1->SetFillColor(0);
780 TLegend* leg2=(TLegend*)leg1->Clone();
781 TLegend* leg3=(TLegend*)leg1->Clone();
782 TLegend* leg4=new TLegend(0.4,0.6,0.8,0.89);
783 leg4->SetFillColor(0);
785 for(Int_t iTypes=0;iTypes<ncmp;iTypes++){
786 TFile* fin=new TFile(filenameYield[iTypes]);
788 printf("WARNING: %s not found",filenameYield[iTypes].Data());
792 TH1F* hSignal=(TH1F*)fin->Get("hSignal");
793 TH1F* hBackground=(TH1F*)fin->Get("hBackground");
794 TH1F* hBackgroundNormSigma=(TH1F*)fin->Get("hBackgroundNormSigma");
795 TH1F* hSignificance=(TH1F*)fin->Get("hSignificance");
796 hSignal->SetName(Form("%s%d",hSignal->GetName(),iTypes));
797 hBackground->SetName(Form("%s%d",hBackground->GetName(),iTypes));
798 hBackgroundNormSigma->SetName(Form("%s%d",hBackgroundNormSigma->GetName(),iTypes));
799 hSignificance->SetName(Form("%s%d",hSignificance->GetName(),iTypes));
801 hSignal->SetMarkerColor(iTypes+2);
802 hSignal->SetLineColor(iTypes+2);
803 hBackground->SetMarkerColor(iTypes+2);
804 hBackground->SetLineColor(iTypes+2);
805 hBackgroundNormSigma->SetMarkerColor(iTypes+2);
806 hBackgroundNormSigma->SetLineColor(iTypes+2);
807 hSignificance->SetMarkerColor(iTypes+2);
808 hSignificance->SetLineColor(iTypes+2);
810 TLegendEntry* ent4=leg4->AddEntry(hSignal,Form("%s",legtext[iTypes].Data()),"PL");
811 ent4->SetTextColor(hSignal->GetMarkerColor());
814 printf("Info: Normalizing signal, background and significance to the number of events\n");
815 hSignal->Scale(1./nevents[iTypes]);
816 hBackground->Scale(1./nevents[iTypes]);
817 hBackgroundNormSigma->Scale(1./nevents[iTypes]);
818 hSignificance->Scale(1./TMath::Sqrt(nevents[iTypes]));
823 hSignal->DrawClone("P");
825 hBackground->DrawClone("P");
827 hSignificance->DrawClone("P");
829 hBackgroundNormSigma->DrawClone("P");
832 hSignal->DrawClone("Psames");
834 hBackground->DrawClone("Psames");
836 hSignificance->DrawClone("Psames");
838 hBackgroundNormSigma->DrawClone("Psames");
841 TH1F* hRelErrSig=(TH1F*)fin->Get("hRelErrSig");
842 TH1F* hInvSignif=(TH1F*)fin->Get("hInvSignif");
843 hRelErrSig->SetName(Form("%s%d",hRelErrSig->GetName(),iTypes));
844 hInvSignif->SetName(Form("%s%d",hInvSignif->GetName(),iTypes));
846 hRelErrSig->SetMarkerColor(iTypes+2);
847 hRelErrSig->SetLineColor(iTypes+2);
848 hInvSignif->SetMarkerColor(iTypes+2);
849 hInvSignif->SetLineColor(iTypes+2);
851 TLegendEntry* ent1=leg1->AddEntry(hRelErrSig,Form("From Fit (%s)",legtext[iTypes].Data()),"P");
852 ent1->SetTextColor(hRelErrSig->GetMarkerColor());
853 ent1=leg1->AddEntry(hInvSignif,Form("1/Significance (%s)",legtext[iTypes].Data()),"PL");
854 ent1->SetTextColor(hInvSignif->GetMarkerColor());
858 hRelErrSig->DrawClone("P");
859 hInvSignif->DrawClone();
861 hRelErrSig->DrawClone("Psames");
862 hInvSignif->DrawClone("sames");
865 TH1F* hNDiffCntSig1=(TH1F*)fin->Get("hNDiffCntSig1");
866 TH1F* hNDiffCntSig2=(TH1F*)fin->Get("hNDiffCntSig2");
867 hNDiffCntSig1->SetName(Form("%s%d",hNDiffCntSig1->GetName(),iTypes));
868 hNDiffCntSig2->SetName(Form("%s%d",hNDiffCntSig2->GetName(),iTypes));
870 hNDiffCntSig1->SetMarkerColor(iTypes+2);
871 hNDiffCntSig1->SetLineColor(iTypes+2);
872 hNDiffCntSig2->SetMarkerColor(iTypes+2);
873 hNDiffCntSig2->SetLineColor(iTypes+2);
874 TLegendEntry* ent2=leg2->AddEntry(hNDiffCntSig1,Form("From Counting1 (%s)",legtext[iTypes].Data()),"PL");
875 ent2->SetTextColor(hNDiffCntSig1->GetMarkerColor());
876 ent2=leg2->AddEntry(hNDiffCntSig2,Form("From Counting2 (%s)",legtext[iTypes].Data()),"PL");
877 ent2->SetTextColor(hNDiffCntSig2->GetMarkerColor());
881 hNDiffCntSig1->DrawClone();
882 hNDiffCntSig2->DrawClone();
884 hNDiffCntSig1->DrawClone("sames");
885 hNDiffCntSig2->DrawClone("sames");
888 TGraph* grReducedChiSquare=(TGraph*)fin->Get("grReducedChiSquare");
889 grReducedChiSquare->SetName(Form("%s%d",grReducedChiSquare->GetName(),iTypes));
891 grReducedChiSquare->SetMarkerColor(iTypes+2);
892 grReducedChiSquare->SetLineColor(iTypes+2);
893 TLegendEntry* ent3=leg3->AddEntry(grReducedChiSquare,Form("%s",legtext[iTypes].Data()),"PL");
894 ent3->SetTextColor(grReducedChiSquare->GetMarkerColor());
897 if(iTypes==0) grReducedChiSquare->DrawClone("AP");
898 else grReducedChiSquare->DrawClone("P");
913 TFile* fout=new TFile("ComparisonRawYield.root","RECREATE");
921 TH1F* RebinHisto(TH1F* hOrig, Int_t reb, Int_t firstUse){
922 // Rebin histogram, from bin firstUse to lastUse
923 // Use all bins if firstUse=-1
925 Int_t nBinOrig=hOrig->GetNbinsX();
926 Int_t firstBinOrig=1;
927 Int_t lastBinOrig=nBinOrig;
928 Int_t nBinOrigUsed=nBinOrig;
929 Int_t nBinFinal=nBinOrig/reb;
931 firstBinOrig=firstUse;
932 nBinFinal=(nBinOrig-firstUse+1)/reb;
933 nBinOrigUsed=nBinFinal*reb;
934 lastBinOrig=firstBinOrig+nBinOrigUsed-1;
936 Int_t exc=nBinOrigUsed%reb;
940 lastBinOrig=firstBinOrig+nBinOrigUsed-1;
944 printf("Rebin from %d bins to %d bins -- Used bins=%d in range %d-%d\n",nBinOrig,nBinFinal,nBinOrigUsed,firstBinOrig,lastBinOrig);
945 Float_t lowLim=hOrig->GetXaxis()->GetBinLowEdge(firstBinOrig);
946 Float_t hiLim=hOrig->GetXaxis()->GetBinUpEdge(lastBinOrig);
947 TH1F* hRebin=new TH1F(Form("%s-rebin",hOrig->GetName()),hOrig->GetTitle(),nBinFinal,lowLim,hiLim);
948 Int_t lastSummed=firstBinOrig-1;
949 for(Int_t iBin=1;iBin<=nBinFinal; iBin++){
951 for(Int_t iOrigBin=0;iOrigBin<reb;iOrigBin++){
952 sum+=hOrig->GetBinContent(lastSummed+1);
955 hRebin->SetBinContent(iBin,sum);