]>
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> | |
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" | |
470a7b55 | 22 | #include "AliRDHFCutsDStartoKpipi.h" |
23851d0b | 23 | #include "AliRDHFCutsDstoKKpi.h" |
cb2336fe | 24 | #include "AliRDHFCutsD0toKpi.h" |
b92c0195 | 25 | #include "AliHFMassFitter.h" |
642de1ee | 26 | #include "AliNormalizationCounter.h" |
b92c0195 | 27 | #endif |
28 | ||
29 | ||
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) | |
cb2336fe | 33 | // D0: C. Bianchin (cbianchi@pd.infn.it) |
b92c0195 | 34 | |
35 | ||
36 | ||
37 | // | |
470a7b55 | 38 | enum {kD0toKpi, kDplusKpipi, kDStarD0pi, kDsKKpi}; |
cb2336fe | 39 | enum {kBoth, kParticleOnly, kAntiParticleOnly}; |
470a7b55 | 40 | enum {kExpo=0, kThrExpo=5, kLinear, kPol2}; |
b92c0195 | 41 | enum {kGaus=0, kDoubleGaus}; |
42 | ||
43 | ||
44 | // Common variables: to be configured by the user | |
0587ab27 | 45 | const Int_t nPtBins=6; |
46 | Double_t ptlims[nPtBins+1]={2.,3.,4.,5.,6.,8.,12.}; | |
1918f6ea | 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 | |
36f9ecb7 | 49 | |
642de1ee | 50 | TString suffix="Loose_SecondSet1236_ForCF08"; |
51 | ||
52 | ||
0587ab27 | 53 | //const Int_t nPtBins=7;//6; |
54 | //Double_t ptlims[nPtBins+1]={1.,2.,3.,4.,5.,6.,8.,12.}; | |
5f1ad8a6 | 55 | //Int_t rebin[nPtBins+1]={8,6,10,10,10,10,10,10}; //for looser cuts |
0587ab27 | 56 | //Int_t rebin[nPtBins+1]={10,10,10,10,10,10,10,10}; //for Chiara's cuts |
5f1ad8a6 | 57 | Int_t typeb=kExpo; |
b92c0195 | 58 | Int_t types=kGaus; |
cb2336fe | 59 | Int_t optPartAntiPart=kBoth; |
b92c0195 | 60 | Int_t factor4refl=0; |
e654a7da | 61 | Double_t minMassForFit=1.7; |
62 | Double_t maxMassForFit=2.1; | |
b92c0195 | 63 | Float_t massRangeForCounting=0.05; // GeV |
0587ab27 | 64 | TH2F* hPtMass=0x0; |
642de1ee | 65 | Double_t nEventsForNorm=0.; |
b92c0195 | 66 | |
cb2336fe | 67 | //for D0only |
36f9ecb7 | 68 | Bool_t cutsappliedondistr=kTRUE; |
5f1ad8a6 | 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*/}; | |
b92c0195 | 76 | // Functions |
77 | ||
78 | Bool_t LoadDplusHistos(TObjArray* listFiles, TH1F** hMass); | |
23851d0b | 79 | Bool_t LoadDsHistos(TObjArray* listFiles, TH1F** hMass); |
b92c0195 | 80 | Bool_t LoadD0toKpiHistos(TObjArray* listFiles, TH1F** hMass); |
470a7b55 | 81 | Bool_t LoadDstarD0piHistos(TObjArray* listFiles, TH1F** hMass); |
82 | ||
1918f6ea | 83 | TH1F* RebinHisto(TH1F* hOrig, Int_t reb, Int_t firstUse=-1); |
b92c0195 | 84 | |
85 | ||
86 | ||
87 | ||
88 | void FitMassSpectra(Int_t analysisType=kDplusKpipi, | |
23851d0b | 89 | TString fileNameb="", |
90 | TString fileNamec="", | |
91 | TString fileNamed="", | |
92 | TString fileNamee="" | |
b92c0195 | 93 | ){ |
94 | // | |
95 | ||
1918f6ea | 96 | gInterpreter->ExecuteMacro("$ALICE_ROOT/PWGHF/vertexingHF/macros/LoadLibraries.C"); |
b92c0195 | 97 | gStyle->SetOptTitle(0); |
98 | ||
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())); | |
23851d0b | 103 | if(fileNamee!="") listFiles->AddLast(new TObjString(fileNamee.Data())); |
b92c0195 | 104 | if(listFiles->GetEntries()==0){ |
105 | printf("Missing file names in input\n"); | |
106 | return; | |
107 | } | |
108 | ||
b92c0195 | 109 | TH1F** hmass=new TH1F*[nPtBins]; |
110 | for(Int_t i=0;i<nPtBins;i++) hmass[i]=0x0; | |
cb2336fe | 111 | |
470a7b55 | 112 | Float_t massD, massD0_fDstar; |
b92c0195 | 113 | Bool_t retCode; |
114 | if(analysisType==kD0toKpi){ | |
115 | retCode=LoadD0toKpiHistos(listFiles,hmass); | |
116 | massD=TDatabasePDG::Instance()->GetParticle(421)->Mass(); | |
117 | } | |
118 | else if(analysisType==kDplusKpipi){ | |
119 | retCode=LoadDplusHistos(listFiles,hmass); | |
120 | massD=TDatabasePDG::Instance()->GetParticle(411)->Mass(); | |
121 | } | |
470a7b55 | 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; | |
127 | } | |
23851d0b | 128 | else if(analysisType==kDsKKpi){ |
129 | retCode=LoadDsHistos(listFiles,hmass); | |
130 | massD=TDatabasePDG::Instance()->GetParticle(431)->Mass(); | |
131 | } | |
b92c0195 | 132 | else{ |
23851d0b | 133 | printf("Wrong analysis type parameter\n"); |
b92c0195 | 134 | return; |
135 | } | |
136 | if(!retCode){ | |
137 | printf("ERROR in reading input files\n"); | |
138 | return; | |
139 | } | |
140 | ||
141 | ||
142 | ||
0587ab27 | 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); | |
155 | ||
b92c0195 | 156 | |
1918f6ea | 157 | Int_t nMassBins=hmass[0]->GetNbinsX(); |
e654a7da | 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)); | |
1918f6ea | 160 | Float_t minBinSum=hmass[0]->FindBin(massD-massRangeForCounting); |
161 | Float_t maxBinSum=hmass[0]->FindBin(massD+massRangeForCounting); | |
b92c0195 | 162 | Int_t iPad=1; |
163 | ||
cb2336fe | 164 | TF1* funBckStore1=0x0; |
165 | TF1* funBckStore2=0x0; | |
166 | TF1* funBckStore3=0x0; | |
b92c0195 | 167 | |
168 | AliHFMassFitter** fitter=new AliHFMassFitter*[nPtBins]; | |
5f1ad8a6 | 169 | Double_t arrchisquare[nPtBins]; |
b92c0195 | 170 | TCanvas* c1= new TCanvas("c1","MassSpectra",1000,800); |
171 | Int_t nx=3, ny=2; | |
172 | if(nPtBins>6){ | |
173 | nx=4; | |
174 | ny=3; | |
175 | } | |
176 | if(nPtBins>12){ | |
177 | nx=5; | |
178 | ny=4; | |
179 | } | |
180 | c1->Divide(nx,ny); | |
181 | ||
182 | Double_t sig,errsig,s,errs,b,errb; | |
183 | for(Int_t iBin=0; iBin<nPtBins; iBin++){ | |
184 | c1->cd(iPad++); | |
53d1574e | 185 | Int_t origNbins=hmass[iBin]->GetNbinsX(); |
1918f6ea | 186 | TH1F* hRebinned=RebinHisto(hmass[iBin],rebin[iBin],firstUsedBin[iBin]); |
e654a7da | 187 | hmin=TMath::Max(minMassForFit,hRebinned->GetBinLowEdge(2)); |
188 | hmax=TMath::Min(maxMassForFit,hRebinned->GetBinLowEdge(hRebinned->GetNbinsX())); | |
1918f6ea | 189 | fitter[iBin]=new AliHFMassFitter(hRebinned,hmin, hmax,1,typeb,types); |
53d1574e | 190 | rebin[iBin]=origNbins/fitter[iBin]->GetBinN(); |
b92c0195 | 191 | fitter[iBin]->SetReflectionSigmaFactor(factor4refl); |
192 | fitter[iBin]->SetInitialGaussianMean(massD); | |
470a7b55 | 193 | if(analysisType==kDStarD0pi) fitter[iBin]->SetInitialGaussianSigma(0.00065); |
b92c0195 | 194 | Bool_t out=fitter[iBin]->MassFitter(0); |
5f1ad8a6 | 195 | if(!out) { |
196 | fitter[iBin]->GetHistoClone()->Draw(); | |
197 | continue; | |
198 | } | |
0587ab27 | 199 | Double_t mass=fitter[iBin]->GetMean(); |
200 | Double_t sigma=fitter[iBin]->GetSigma(); | |
5f1ad8a6 | 201 | arrchisquare[iBin]=fitter[iBin]->GetReducedChiSquare(); |
b92c0195 | 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); | |
208 | ||
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); | |
e654a7da | 213 | Double_t ry=fitter[iBin]->GetRawYield(); |
214 | Double_t ery=fitter[iBin]->GetRawYieldError(); | |
b92c0195 | 215 | Float_t cntSig1=0.; |
216 | Float_t cntSig2=0.; | |
217 | Float_t cntErr=0.; | |
218 | for(Int_t iMB=minBinSum; iMB<=maxBinSum; iMB++){ | |
cb2336fe | 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; | |
b92c0195 | 221 | cntSig1+=(hmass[iBin]->GetBinContent(iMB)-bkg1); |
222 | cntSig2+=(hmass[iBin]->GetBinContent(iMB)-bkg2); | |
223 | cntErr+=(hmass[iBin]->GetBinContent(iMB)); | |
224 | } | |
225 | hCntSig1->SetBinContent(iBin+1,cntSig1); | |
226 | hCntSig1->SetBinError(iBin+1,TMath::Sqrt(cntErr)); | |
36f9ecb7 | 227 | hNDiffCntSig1->SetBinContent(iBin+1,(s-cntSig1)/s); |
228 | hNDiffCntSig1->SetBinError(iBin+1,TMath::Sqrt(cntErr)/s); | |
b92c0195 | 229 | hCntSig2->SetBinContent(iBin+1,cntSig2); |
36f9ecb7 | 230 | hNDiffCntSig2->SetBinContent(iBin+1,(s-cntSig2)/s); |
231 | hNDiffCntSig2->SetBinError(iBin+1,TMath::Sqrt(cntErr)/s); | |
b92c0195 | 232 | hCntSig2->SetBinError(iBin+1,TMath::Sqrt(cntErr)); |
e654a7da | 233 | hSignal->SetBinContent(iBin+1,ry); |
234 | hSignal->SetBinError(iBin+1,ery); | |
36f9ecb7 | 235 | hRelErrSig->SetBinContent(iBin+1,errs/s); |
236 | hInvSignif->SetBinContent(iBin+1,1/sig); | |
237 | hInvSignif->SetBinError(iBin+1,errsig/(sig*sig)); | |
5f1ad8a6 | 238 | hBackground->SetBinContent(iBin+1,b); //consider sigma |
b92c0195 | 239 | hBackground->SetBinError(iBin+1,errb); |
5f1ad8a6 | 240 | hBackgroundNormSigma->SetBinContent(iBin+1,b/(3*fitter[iBin]->GetSigma())*(3*0.012)); //consider sigma |
241 | hBackgroundNormSigma->SetBinError(iBin+1,errb); | |
b92c0195 | 242 | hSignificance->SetBinContent(iBin+1,sig); |
243 | hSignificance->SetBinError(iBin+1,errsig); | |
0587ab27 | 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); | |
b92c0195 | 248 | |
249 | } | |
0587ab27 | 250 | |
5f1ad8a6 | 251 | /* |
b92c0195 | 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"); | |
5f1ad8a6 | 256 | */ |
0587ab27 | 257 | |
258 | TCanvas *cpar=new TCanvas("cpar","Fit params",1200,600); | |
259 | cpar->Divide(2,1); | |
260 | cpar->cd(1); | |
261 | hMass->SetMarkerStyle(20); | |
262 | hMass->Draw("PE"); | |
263 | hMass->GetXaxis()->SetTitle("Pt (GeV/c)"); | |
264 | hMass->GetXaxis()->SetTitle("Mass (GeV/c^{2})"); | |
265 | cpar->cd(2); | |
266 | hSigma->SetMarkerStyle(20); | |
267 | hSigma->Draw("PE"); | |
268 | hSigma->GetXaxis()->SetTitle("Pt (GeV/c)"); | |
269 | hSigma->GetXaxis()->SetTitle("Sigma (GeV/c^{2})"); | |
270 | ||
b92c0195 | 271 | TCanvas* csig=new TCanvas("csig","Results",1200,600); |
272 | csig->Divide(3,1); | |
273 | csig->cd(1); | |
274 | hSignal->SetMarkerStyle(20); | |
275 | hSignal->SetMarkerColor(4); | |
276 | hSignal->SetLineColor(4); | |
277 | hSignal->GetXaxis()->SetTitle("Pt (GeV/c)"); | |
278 | hSignal->GetYaxis()->SetTitle("Signal"); | |
279 | hSignal->Draw("P"); | |
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()); | |
296 | leg->Draw(); | |
297 | csig->cd(2); | |
298 | hBackground->SetMarkerStyle(20); | |
299 | hBackground->Draw("P"); | |
300 | hBackground->GetXaxis()->SetTitle("Pt (GeV/c)"); | |
301 | hBackground->GetYaxis()->SetTitle("Background"); | |
302 | csig->cd(3); | |
303 | hSignificance->SetMarkerStyle(20); | |
304 | hSignificance->Draw("P"); | |
305 | hSignificance->GetXaxis()->SetTitle("Pt (GeV/c)"); | |
306 | hSignificance->GetYaxis()->SetTitle("Significance"); | |
307 | ||
36f9ecb7 | 308 | TCanvas* cDiffS=new TCanvas("cDiffS","Difference",1200,600); |
309 | cDiffS->Divide(2,1); | |
310 | cDiffS->cd(1); | |
311 | hRelErrSig->SetMarkerStyle(20); //fullcircle | |
312 | hRelErrSig->SetTitleOffset(1.2); | |
5f1ad8a6 | 313 | hRelErrSig->SetTitle("Rel Error from Fit;p_{t} (GeV/c);Signal Relative Error"); |
36f9ecb7 | 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()); | |
325 | leg2->Draw(); | |
326 | ||
327 | cDiffS->cd(2); | |
328 | hNDiffCntSig1->SetMarkerStyle(26); | |
329 | hNDiffCntSig1->SetMarkerColor(2); | |
330 | hNDiffCntSig1->SetLineColor(2); | |
5f1ad8a6 | 331 | hNDiffCntSig1->SetTitle("Cmp Fit-Count;p_{t} (GeV/c);(S_{fit}-S_{count})/S_{fit}"); |
36f9ecb7 | 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()); | |
343 | leg1->Draw(); | |
344 | ||
5f1ad8a6 | 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); | |
349 | cChi2->cd(); | |
350 | grReducedChiSquare->SetMarkerStyle(21); | |
351 | grReducedChiSquare->Draw("AP"); | |
352 | ||
353 | TCanvas* cbkgNormSigma=new TCanvas("cbkgNormSigma","Background normalized to sigma",400,600); | |
354 | cbkgNormSigma->cd(); | |
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(); | |
360 | ||
361 | ||
362 | TString partname="Both"; | |
363 | if(optPartAntiPart==kParticleOnly) { | |
364 | if(analysisType==kD0toKpi) partname="D0"; | |
365 | if(analysisType==kDplusKpipi) partname="Dplus"; | |
23851d0b | 366 | if(analysisType==kDsKKpi) partname="Dsplus"; |
5f1ad8a6 | 367 | } |
368 | if(optPartAntiPart==kAntiParticleOnly) { | |
369 | if(analysisType==kD0toKpi) partname="D0bar"; | |
370 | if(analysisType==kDplusKpipi) partname="Dminus"; | |
23851d0b | 371 | if(analysisType==kDsKKpi) partname="Dsminus"; |
5f1ad8a6 | 372 | } |
36f9ecb7 | 373 | |
642de1ee | 374 | printf("Events for norm = %f\n",nEventsForNorm); |
375 | TH1F* hNEvents=new TH1F("hNEvents","",1,0.,1.); | |
376 | hNEvents->SetBinContent(1,nEventsForNorm); | |
377 | ||
5f1ad8a6 | 378 | TFile* outf=new TFile(Form("RawYield%s.root",partname.Data()),"update"); |
b92c0195 | 379 | outf->cd(); |
642de1ee | 380 | hNEvents->Write(); |
0587ab27 | 381 | hMass->Write(); |
382 | hSigma->Write(); | |
b92c0195 | 383 | hCntSig1->Write(); |
384 | hCntSig2->Write(); | |
36f9ecb7 | 385 | hNDiffCntSig1->Write(); |
386 | hNDiffCntSig2->Write(); | |
b92c0195 | 387 | hSignal->Write(); |
36f9ecb7 | 388 | hRelErrSig->Write(); |
389 | hInvSignif->Write(); | |
b92c0195 | 390 | hBackground->Write(); |
5f1ad8a6 | 391 | hBackgroundNormSigma->Write(); |
b92c0195 | 392 | hSignificance->Write(); |
5f1ad8a6 | 393 | grReducedChiSquare->Write(); |
0587ab27 | 394 | hPtMass->Write(); |
b92c0195 | 395 | outf->Close(); |
396 | } | |
397 | ||
398 | ||
399 | Bool_t LoadDplusHistos(TObjArray* listFiles, TH1F** hMass){ | |
400 | ||
401 | Int_t nFiles=listFiles->GetEntries(); | |
402 | TList **hlist=new TList*[nFiles]; | |
403 | AliRDHFCutsDplustoKpipi** dcuts=new AliRDHFCutsDplustoKpipi*[nFiles]; | |
404 | ||
405 | Int_t nReadFiles=0; | |
406 | for(Int_t iFile=0; iFile<nFiles; iFile++){ | |
407 | TString fName=((TObjString*)listFiles->At(iFile))->GetString(); | |
408 | TFile *f=TFile::Open(fName.Data()); | |
409 | if(!f){ | |
410 | printf("ERROR: file %s does not exist\n",fName.Data()); | |
411 | continue; | |
412 | } | |
413 | printf("Open File %s\n",f->GetName()); | |
414 | TDirectory *dir = (TDirectory*)f->Get("PWG3_D2H_InvMassDplus"); | |
415 | if(!dir){ | |
cb2336fe | 416 | printf("ERROR: directory PWG3_D2H_InvMassDplus not found in %s\n",fName.Data()); |
b92c0195 | 417 | continue; |
418 | } | |
642de1ee | 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); | |
b92c0195 | 422 | if(nReadFiles>0){ |
423 | Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]); | |
424 | if(!sameCuts){ | |
425 | printf("ERROR: Cut objects do not match\n"); | |
426 | return kFALSE; | |
427 | } | |
428 | } | |
642de1ee | 429 | AliNormalizationCounter* c=(AliNormalizationCounter*)dir->Get(Form("coutputDplusNorm%s",suffix.Data())); |
430 | printf("Events for normalization = %f\n",c->GetNEventsForNorm()); | |
431 | nEventsForNorm+=c->GetNEventsForNorm(); | |
b92c0195 | 432 | nReadFiles++; |
433 | } | |
434 | if(nReadFiles<nFiles){ | |
435 | printf("WARNING: not all requested files have been found\n"); | |
436 | } | |
437 | ||
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.; | |
442 | ||
443 | Int_t iFinBin=0; | |
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++){ | |
cb2336fe | 450 | TString histoName; |
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()); | |
455 | if(!htemp){ | |
456 | printf("ERROR: Histogram %s not found\n",histoName.Data()); | |
457 | return kFALSE; | |
458 | } | |
b92c0195 | 459 | if(!hMass[iFinBin]){ |
460 | hMass[iFinBin]=new TH1F(*htemp); | |
461 | }else{ | |
462 | hMass[iFinBin]->Add(htemp); | |
463 | } | |
464 | } | |
465 | } | |
466 | } | |
5f1ad8a6 | 467 | TString partname="Both"; |
468 | if(optPartAntiPart==kParticleOnly) partname="Dplus"; | |
469 | if(optPartAntiPart==kAntiParticleOnly) partname="Dminus"; | |
b92c0195 | 470 | |
0587ab27 | 471 | for(Int_t iFile=0; iFile<nReadFiles; iFile++){ |
472 | TH2F* htemp2=(TH2F*)hlist[iFile]->FindObject("hPtVsMassTC"); | |
473 | if(iFile==0){ | |
474 | hPtMass=new TH2F(*htemp2); | |
475 | }else{ | |
476 | hPtMass->Add(htemp2); | |
477 | } | |
478 | } | |
479 | ||
5f1ad8a6 | 480 | TFile* outf=new TFile(Form("RawYield%s.root",partname.Data()),"recreate"); |
b92c0195 | 481 | outf->cd(); |
482 | dcuts[0]->Write(); | |
483 | outf->Close(); | |
484 | ||
485 | return kTRUE; | |
486 | ||
487 | } | |
488 | ||
23851d0b | 489 | |
490 | Bool_t LoadDsHistos(TObjArray* listFiles, TH1F** hMass){ | |
491 | ||
492 | Int_t nFiles=listFiles->GetEntries(); | |
493 | TList **hlist=new TList*[nFiles]; | |
494 | AliRDHFCutsDstoKKpi** dcuts=new AliRDHFCutsDstoKKpi*[nFiles]; | |
495 | ||
496 | Int_t nReadFiles=0; | |
497 | for(Int_t iFile=0; iFile<nFiles; iFile++){ | |
498 | TString fName=((TObjString*)listFiles->At(iFile))->GetString(); | |
499 | TFile *f=TFile::Open(fName.Data()); | |
500 | if(!f){ | |
501 | printf("ERROR: file %s does not exist\n",fName.Data()); | |
502 | continue; | |
503 | } | |
504 | printf("Open File %s\n",f->GetName()); | |
505 | TDirectory *dir = (TDirectory*)f->Get("PWG3_D2H_InvMassDs"); | |
506 | if(!dir){ | |
507 | printf("ERROR: directory PWG3_D2H_InvMassDs not found in %s\n",fName.Data()); | |
508 | continue; | |
509 | } | |
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; | |
514 | if(nReadFiles>0){ | |
515 | Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]); | |
516 | if(!sameCuts){ | |
517 | printf("ERROR: Cut objects do not match\n"); | |
518 | return kFALSE; | |
519 | } | |
520 | } | |
521 | nReadFiles++; | |
522 | } | |
523 | if(nReadFiles<nFiles){ | |
524 | printf("WARNING: not all requested files have been found\n"); | |
525 | } | |
526 | ||
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.; | |
531 | ||
532 | Int_t iFinBin=0; | |
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++){ | |
539 | TString histoName; | |
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); | |
544 | } | |
545 | else if(optPartAntiPart==kAntiParticleOnly){ | |
546 | printf("Particle/Antiparticle not yet enabled for Ds"); | |
547 | histoName.Form("hMassAllPt%dphi",i); | |
548 | } | |
549 | TH1F* htemp=(TH1F*)hlist[iFile]->FindObject(histoName.Data()); | |
550 | if(!htemp){ | |
551 | printf("ERROR: Histogram %s not found\n",histoName.Data()); | |
552 | return kFALSE; | |
553 | } | |
554 | if(!hMass[iFinBin]){ | |
555 | hMass[iFinBin]=new TH1F(*htemp); | |
556 | }else{ | |
557 | hMass[iFinBin]->Add(htemp); | |
558 | } | |
559 | } | |
560 | } | |
561 | } | |
562 | TString partname="Both"; | |
563 | if(optPartAntiPart==kParticleOnly) partname="Both"; | |
564 | if(optPartAntiPart==kAntiParticleOnly) partname="Both"; | |
565 | ||
566 | for(Int_t iFile=0; iFile<nReadFiles; iFile++){ | |
567 | TH2F* htemp2=(TH2F*)hlist[iFile]->FindObject("hPtVsMassPhi"); | |
568 | if(iFile==0){ | |
569 | hPtMass=new TH2F(*htemp2); | |
570 | }else{ | |
571 | hPtMass->Add(htemp2); | |
572 | } | |
573 | } | |
574 | ||
575 | TFile* outf=new TFile(Form("RawYield%s.root",partname.Data()),"recreate"); | |
576 | outf->cd(); | |
577 | dcuts[0]->Write(); | |
578 | outf->Close(); | |
579 | ||
580 | return kTRUE; | |
581 | ||
582 | } | |
583 | ||
584 | ||
b92c0195 | 585 | Bool_t LoadD0toKpiHistos(TObjArray* listFiles, TH1F** hMass){ |
cb2336fe | 586 | // |
587 | Int_t nFiles=listFiles->GetEntries(); | |
588 | TList **hlist=new TList*[nFiles]; | |
589 | AliRDHFCutsD0toKpi** dcuts=new AliRDHFCutsD0toKpi*[nFiles]; | |
590 | ||
591 | Int_t nReadFiles=0; | |
592 | for(Int_t iFile=0; iFile<nFiles; iFile++){ | |
593 | TString fName=((TObjString*)listFiles->At(iFile))->GetString(); | |
594 | TFile *f=TFile::Open(fName.Data()); | |
595 | if(!f){ | |
596 | printf("ERROR: file %s does not exist\n",fName.Data()); | |
597 | continue; | |
598 | } | |
599 | printf("Open File %s\n",f->GetName()); | |
600 | ||
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); | |
606 | if(!dir){ | |
607 | printf("ERROR: directory %s not found in %s\n",dirname.Data(),fName.Data()); | |
608 | continue; | |
609 | } | |
610 | TString listmassname="coutputmassD0Mass"; | |
5f1ad8a6 | 611 | if(optPartAntiPart==kParticleOnly) listmassname+="D0"; |
612 | else if(optPartAntiPart==kAntiParticleOnly) listmassname+="D0bar"; | |
cb2336fe | 613 | if(cutsappliedondistr) listmassname+="C"; |
614 | ||
615 | hlist[nReadFiles]=(TList*)dir->Get(listmassname); | |
616 | ||
617 | TString cutsobjname="cutsD0"; | |
5f1ad8a6 | 618 | if(optPartAntiPart==kParticleOnly) cutsobjname+="D0"; |
619 | else if(optPartAntiPart==kAntiParticleOnly) cutsobjname+="D0bar"; | |
cb2336fe | 620 | if(cutsappliedondistr) cutsobjname+="C"; |
621 | ||
622 | dcuts[nReadFiles]=(AliRDHFCutsD0toKpi*)dir->Get(cutsobjname); | |
623 | if(nReadFiles>0){ | |
624 | Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]); | |
625 | if(!sameCuts){ | |
626 | printf("ERROR: Cut objects do not match\n"); | |
627 | return kFALSE; | |
628 | } | |
629 | } | |
630 | nReadFiles++; | |
631 | } | |
632 | if(nReadFiles<nFiles){ | |
633 | printf("WARNING: not all requested files have been found\n"); | |
36f9ecb7 | 634 | if (nReadFiles==0) { |
635 | printf("ERROR: Any file/dir found\n"); | |
636 | return kFALSE; | |
637 | } | |
cb2336fe | 638 | } |
639 | ||
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.; | |
644 | ||
645 | Int_t iFinBin=0; | |
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)); | |
653 | if(!hMass[iFinBin]){ | |
654 | hMass[iFinBin]=new TH1F(*htemp); | |
655 | }else{ | |
656 | hMass[iFinBin]->Add(htemp); | |
657 | } | |
658 | } | |
659 | } | |
660 | } | |
5f1ad8a6 | 661 | TString partname="Both"; |
662 | if(optPartAntiPart==kParticleOnly) partname="D0"; | |
663 | if(optPartAntiPart==kAntiParticleOnly) partname="D0bar"; | |
cb2336fe | 664 | |
5f1ad8a6 | 665 | TFile* outf=new TFile(Form("RawYield%s.root",partname.Data()),"recreate"); |
cb2336fe | 666 | outf->cd(); |
667 | dcuts[0]->Write(); | |
668 | outf->Close(); | |
669 | return kTRUE; | |
b92c0195 | 670 | } |
5f1ad8a6 | 671 | |
470a7b55 | 672 | Bool_t LoadDstarD0piHistos(TObjArray* listFiles, TH1F** hMass){ |
673 | // | |
674 | Int_t nFiles=listFiles->GetEntries(); | |
675 | TList **hlist=new TList*[nFiles]; | |
676 | AliRDHFCutsDStartoKpipi** dcuts=new AliRDHFCutsDStartoKpipi*[nFiles]; | |
677 | Int_t nReadFiles=0; | |
678 | for(Int_t iFile=0; iFile<nFiles; iFile++){ | |
679 | TString fName=((TObjString*)listFiles->At(iFile))->GetString(); | |
680 | TFile *f=TFile::Open(fName.Data()); | |
681 | if(!f){ | |
682 | printf("ERROR: file %s does not exist\n",fName.Data()); | |
683 | continue; | |
684 | } | |
685 | printf("Open File %s\n",f->GetName()); | |
686 | TString dirname="PWG3_D2H_DStarSpectra"; | |
687 | TDirectory *dir = (TDirectory*)f->Get(dirname); | |
688 | if(!dir){ | |
689 | printf("ERROR: directory %s not found in %s\n",dirname.Data(),fName.Data()); | |
690 | continue; | |
691 | } | |
692 | TString listmassname="DStarPID00"; | |
693 | ||
694 | hlist[nReadFiles]=(TList*)dir->Get(listmassname); | |
695 | TString cutsobjname="DStartoKpipiCuts"; | |
696 | dcuts[nReadFiles]=(AliRDHFCutsDStartoKpipi*)dir->Get(cutsobjname); | |
697 | if(nReadFiles>0){ | |
698 | Bool_t sameCuts=dcuts[nReadFiles]->CompareCuts(dcuts[0]); | |
699 | if(!sameCuts){ | |
700 | printf("ERROR: Cut objects do not match\n"); | |
701 | return kFALSE; | |
702 | } | |
703 | } | |
704 | nReadFiles++; | |
705 | } | |
706 | if(nReadFiles<nFiles){ | |
707 | printf("WARNING: not all requested files have been found\n"); | |
708 | if (nReadFiles==0) { | |
709 | printf("ERROR: Any file/dir found\n"); | |
710 | return kFALSE; | |
711 | } | |
712 | } | |
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.; | |
717 | ||
718 | ||
719 | ||
720 | ||
721 | Int_t iFinBin=0; | |
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)); | |
729 | if(!hMass[iFinBin]){ | |
730 | hMass[iFinBin]=new TH1F(*htemp); | |
731 | }else{ | |
732 | hMass[iFinBin]->Add(htemp); | |
733 | } | |
734 | } | |
735 | } | |
736 | } | |
737 | TString partname="Both"; | |
738 | if(optPartAntiPart==kParticleOnly) partname="DStar"; | |
739 | if(optPartAntiPart==kAntiParticleOnly) partname="DStarbar"; | |
740 | ||
741 | TFile* outf=new TFile(Form("RawYield%s.root",partname.Data()),"recreate"); | |
742 | outf->cd(); | |
743 | dcuts[0]->Write(); | |
744 | outf->Close(); | |
745 | return kTRUE; | |
746 | } | |
747 | ||
5f1ad8a6 | 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 | |
751 | //arguments: | |
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) | |
756 | ||
757 | gStyle->SetOptStat(0); | |
758 | gStyle->SetOptTitle(0); | |
759 | gStyle->SetCanvasColor(0); | |
760 | gStyle->SetFrameFillColor(0); | |
761 | gStyle->SetTitleFillColor(0); | |
762 | gStyle->SetFrameBorderMode(0); | |
763 | ||
764 | if(!filenameYield) filenameYield=new TString[ncmp]; | |
765 | ||
766 | for(Int_t k=0;k<ncmp;k++){ | |
767 | if(!filenameYield) filenameYield[k]="RawYield.root"; | |
768 | filenameYield[k].Prepend(paths[k]); | |
769 | } | |
770 | ||
771 | TCanvas* cSig=new TCanvas("cSig","Results",1200,600); | |
772 | cSig->Divide(3,1); | |
773 | TCanvas* cBkgN=new TCanvas("cBkgN","Background normalized to sigma",400,600); | |
774 | TCanvas* cDiffS=new TCanvas("cDiffS","Difference",1200,600); | |
775 | cDiffS->Divide(2,1); | |
776 | TCanvas *cChi2=new TCanvas("cChi2","reduced chi square",600,600); | |
777 | ||
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); | |
784 | ||
785 | for(Int_t iTypes=0;iTypes<ncmp;iTypes++){ | |
786 | TFile* fin=new TFile(filenameYield[iTypes]); | |
787 | if(!fin){ | |
788 | printf("WARNING: %s not found",filenameYield[iTypes].Data()); | |
789 | continue; | |
790 | } | |
791 | ||
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)); | |
800 | ||
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); | |
809 | ||
810 | TLegendEntry* ent4=leg4->AddEntry(hSignal,Form("%s",legtext[iTypes].Data()),"PL"); | |
811 | ent4->SetTextColor(hSignal->GetMarkerColor()); | |
812 | ||
813 | if(ncmp==nsamples){ | |
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])); | |
819 | } | |
820 | ||
821 | if(iTypes==0){ | |
822 | cSig->cd(1); | |
823 | hSignal->DrawClone("P"); | |
824 | cSig->cd(2); | |
825 | hBackground->DrawClone("P"); | |
826 | cSig->cd(3); | |
827 | hSignificance->DrawClone("P"); | |
828 | cBkgN->cd(); | |
829 | hBackgroundNormSigma->DrawClone("P"); | |
830 | } else{ | |
831 | cSig->cd(1); | |
832 | hSignal->DrawClone("Psames"); | |
833 | cSig->cd(2); | |
834 | hBackground->DrawClone("Psames"); | |
835 | cSig->cd(3); | |
836 | hSignificance->DrawClone("Psames"); | |
837 | cBkgN->cd(); | |
838 | hBackgroundNormSigma->DrawClone("Psames"); | |
839 | } | |
840 | ||
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)); | |
845 | ||
846 | hRelErrSig->SetMarkerColor(iTypes+2); | |
847 | hRelErrSig->SetLineColor(iTypes+2); | |
848 | hInvSignif->SetMarkerColor(iTypes+2); | |
849 | hInvSignif->SetLineColor(iTypes+2); | |
850 | ||
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()); | |
855 | ||
856 | cDiffS->cd(1); | |
857 | if(iTypes==0){ | |
858 | hRelErrSig->DrawClone("P"); | |
859 | hInvSignif->DrawClone(); | |
860 | } else{ | |
861 | hRelErrSig->DrawClone("Psames"); | |
862 | hInvSignif->DrawClone("sames"); | |
863 | } | |
864 | ||
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)); | |
869 | ||
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()); | |
878 | ||
879 | cDiffS->cd(2); | |
880 | if(iTypes==0){ | |
881 | hNDiffCntSig1->DrawClone(); | |
882 | hNDiffCntSig2->DrawClone(); | |
883 | }else{ | |
884 | hNDiffCntSig1->DrawClone("sames"); | |
885 | hNDiffCntSig2->DrawClone("sames"); | |
886 | } | |
887 | ||
888 | TGraph* grReducedChiSquare=(TGraph*)fin->Get("grReducedChiSquare"); | |
889 | grReducedChiSquare->SetName(Form("%s%d",grReducedChiSquare->GetName(),iTypes)); | |
890 | ||
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()); | |
895 | ||
896 | cChi2->cd(); | |
897 | if(iTypes==0) grReducedChiSquare->DrawClone("AP"); | |
898 | else grReducedChiSquare->DrawClone("P"); | |
899 | } | |
900 | ||
901 | cSig->cd(1); | |
902 | leg4->Draw(); | |
903 | ||
904 | cDiffS->cd(1); | |
905 | leg1->Draw(); | |
906 | ||
907 | cDiffS->cd(2); | |
908 | leg2->Draw(); | |
909 | ||
910 | cChi2->cd(); | |
911 | leg3->Draw(); | |
912 | ||
913 | TFile* fout=new TFile("ComparisonRawYield.root","RECREATE"); | |
914 | fout->cd(); | |
915 | cDiffS->Write(); | |
916 | cChi2->Write(); | |
917 | fout->Close(); | |
918 | } | |
1918f6ea | 919 | |
920 | ||
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 | |
924 | ||
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; | |
930 | if(firstUse>=1){ | |
931 | firstBinOrig=firstUse; | |
932 | nBinFinal=(nBinOrig-firstUse+1)/reb; | |
933 | nBinOrigUsed=nBinFinal*reb; | |
934 | lastBinOrig=firstBinOrig+nBinOrigUsed-1; | |
935 | }else{ | |
936 | Int_t exc=nBinOrigUsed%reb; | |
937 | if(exc!=0){ | |
938 | nBinOrigUsed-=exc; | |
939 | firstBinOrig+=exc/2; | |
940 | lastBinOrig=firstBinOrig+nBinOrigUsed-1; | |
941 | } | |
942 | } | |
943 | ||
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++){ | |
950 | Float_t sum=0.; | |
951 | for(Int_t iOrigBin=0;iOrigBin<reb;iOrigBin++){ | |
952 | sum+=hOrig->GetBinContent(lastSummed+1); | |
953 | lastSummed++; | |
954 | } | |
955 | hRebin->SetBinContent(iBin,sum); | |
956 | } | |
957 | return hRebin; | |
958 | } |