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