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