]>
Commit | Line | Data |
---|---|---|
1 | #if !defined(__CINT__) || defined(__MAKECINT__) | |
2 | #include <Riostream.h> | |
3 | #include "TH1D.h" | |
4 | #include "TH1.h" | |
5 | #include "TH2F.h" | |
6 | #include "TNtuple.h" | |
7 | #include "TFile.h" | |
8 | #include "TGraphAsymmErrors.h" | |
9 | #include "TCanvas.h" | |
10 | #include "TROOT.h" | |
11 | #include "TStyle.h" | |
12 | #include "TLegend.h" | |
13 | #include "AliHFSystErr.h" | |
14 | #include "AliHFPtSpectrum.h" | |
15 | #endif | |
16 | ||
17 | /* $Id$ */ | |
18 | ||
19 | // | |
20 | // Macro to use the AliHFPtSpectrum class | |
21 | // to compute the feed-down corrections for heavy-flavor | |
22 | // | |
23 | // Z.Conesa, September 2010 (zconesa@in2p3.fr) | |
24 | // | |
25 | ||
26 | ||
27 | ||
28 | // | |
29 | // Macro execution parameters: | |
30 | // 0) filename with the theoretical predictions (direct & feed-down) | |
31 | // 1) acceptance and reconstruction efficiencies file name (direct & feed-down) | |
32 | // 2) reconstructed spectra file name | |
33 | // 3) output file name | |
34 | // 4) Set the feed-down calculation option flag: knone=none, kfc=fc only, kNb=Nb only | |
35 | // 5-6) Set the luminosity: the number of events analyzed, and the cross-section of the sample [pb] | |
36 | // 7) Set whether the yield is for particle + anti-particles or only one of the 'charges' | |
37 | // 8) Set the centrality class | |
38 | // 9) Flag to decide if there is need to evaluate the dependence on the energy loss | |
39 | // | |
40 | ||
41 | enum centrality{ kpp7, kpp276, k07half, kpPb0100, k010, k1020, k020, k2040, k2030, k3040, k4050, k3050, k5060, k4060, k6080, k4080, k5080, k80100, kpPb020, kpPb2040, kpPb4060, kpPb60100 }; | |
42 | enum centestimator{ kV0M, kV0A, kZNA }; | |
43 | enum BFDSubtrMethod { knone, kfc, kNb }; | |
44 | enum RaavsEP {kPhiIntegrated, kInPlane, kOutOfPlane}; | |
45 | enum rapidity{ kdefault, k08to04, k07to04, k04to01, k01to01, k01to04, k04to07, k04to08 }; | |
46 | ||
47 | void HFPtSpectrum ( const char *mcfilename="FeedDownCorrectionMC.root", | |
48 | const char *efffilename="Efficiencies.root", | |
49 | const char *recofilename="Reconstructed.root", const char *recohistoname="hRawSpectrumD0", | |
50 | const char *outfilename="HFPtSpectrum.root", | |
51 | Int_t fdMethod=kNb, Double_t nevents=1.0, Double_t sigma=1.0, // sigma[pb] | |
52 | Bool_t isParticlePlusAntiParticleYield=true, Int_t cc=kpp7, Bool_t PbPbEloss=false, | |
53 | Int_t ccestimator = kV0M, | |
54 | Int_t isRaavsEP=kPhiIntegrated,const char *epResolfile="", | |
55 | Int_t rapiditySlice=kdefault) { | |
56 | ||
57 | ||
58 | gROOT->Macro("$ALICE_ROOT/PWGHF/vertexingHF/macros/LoadLibraries.C"); | |
59 | ||
60 | // Set if calculation considers asymmetric uncertainties or not | |
61 | Bool_t asym = true; | |
62 | ||
63 | // Set the meson/baryon and decay | |
64 | // (only D0 -> K pi, D+--> K pi pi, D* --> D0 pi, D+s -->KKpi, Lc+ --> pKpi & Lc+ --> pK0S implemented here) | |
65 | Bool_t isD0Kpi = true; | |
66 | Bool_t isDplusKpipi = false; | |
67 | Bool_t isDstarD0pi = false; | |
68 | Bool_t isDsKKpi = false; | |
69 | Bool_t isLctopKpi = false; | |
70 | Bool_t isLcK0Sp = false; | |
71 | ||
72 | Int_t shouldBeOne=0; | |
73 | if(isD0Kpi) shouldBeOne++; | |
74 | if(isDplusKpipi) shouldBeOne++; | |
75 | if(isDstarD0pi) shouldBeOne++; | |
76 | if(isDsKKpi) shouldBeOne++; | |
77 | if(isLctopKpi) shouldBeOne++; | |
78 | if(isLcK0Sp) shouldBeOne++; | |
79 | ||
80 | if (shouldBeOne!=1) { | |
81 | cout << "Sorry, can not deal with more than one correction at the same time"<<endl; | |
82 | return; | |
83 | } | |
84 | ||
85 | Int_t option=3; | |
86 | if (fdMethod==kfc) option=1; | |
87 | else if (fdMethod==kNb) option=2; | |
88 | else if (fdMethod==knone) { option=0; asym=false; } | |
89 | else option=3; | |
90 | ||
91 | if (option>2) { | |
92 | cout<< "Bad calculation option, should be <=2"<<endl; | |
93 | return; | |
94 | } | |
95 | ||
96 | ||
97 | // | |
98 | // Defining the Tab values for the given centrality class | |
99 | // https://twiki.cern.ch/twiki/bin/viewauth/ALICE/CentStudies | |
100 | // | |
101 | Double_t tab = 1., tabUnc = 0.; | |
102 | if( ccestimator == kV0M ) { | |
103 | if ( cc == k07half ) { | |
104 | tab = 24.81; tabUnc = 0.8037; | |
105 | } else if ( cc == k010 ) { | |
106 | tab = 23.48; tabUnc = 0.97; | |
107 | } else if ( cc == k1020 ) { | |
108 | tab = 14.4318; tabUnc = 0.5733; | |
109 | } else if ( cc == k020 ) { | |
110 | tab = 18.93; tabUnc = 0.74; | |
111 | } else if ( cc == k2040 ) { | |
112 | tab = 6.86; tabUnc = 0.28; | |
113 | } else if ( cc == k2030 ) { | |
114 | tab = 8.73769; tabUnc = 0.370219; | |
115 | } else if ( cc == k3040 ) { | |
116 | tab = 5.02755; tabUnc = 0.22099; | |
117 | } else if ( cc == k4050 ) { | |
118 | tab = 2.68327; tabUnc = 0.137073; | |
119 | } else if ( cc == k3050 ) { | |
120 | tab = 3.87011; tabUnc = 0.183847; | |
121 | } else if ( cc == k4060 ) { | |
122 | tab = 2.00; tabUnc= 0.11; | |
123 | } else if ( cc == k4080 ) { | |
124 | tab = 1.20451; tabUnc = 0.071843; | |
125 | } else if ( cc == k5060 ) { | |
126 | tab = 1.32884; tabUnc = 0.0929536; | |
127 | } else if ( cc == k6080 ) { | |
128 | tab = 0.419; tabUnc = 0.033; | |
129 | } else if ( cc == k5080 ) { | |
130 | tab = 0.719; tabUnc = 0.054; | |
131 | } else if ( cc == k80100 ){ | |
132 | tab = 0.0690; tabUnc = 0.0062; | |
133 | } | |
134 | } | |
135 | ||
136 | // pPb Glauber (A. Toia) | |
137 | // https://twiki.cern.ch/twiki/bin/viewauth/ALICE/PACentStudies#Glauber_Calculations_with_sigma | |
138 | if( cc == kpPb0100 ){ | |
139 | tab = 0.098334; tabUnc = 0.0070679; | |
140 | // A=207.2; B=1.; | |
141 | } | |
142 | else if( ccestimator == kV0A ){ | |
143 | if ( cc == kpPb020 ) { | |
144 | tab = 0.183; tabUnc = 0.006245; | |
145 | } else if ( cc == kpPb2040 ) { | |
146 | tab = 0.134; tabUnc = 0.004899; | |
147 | } else if ( cc == kpPb4060 ) { | |
148 | tab = 0.092; tabUnc = 0.004796; | |
149 | } else if ( cc == kpPb60100 ) { | |
150 | tab = 0.041; tabUnc = 0.008832; | |
151 | } | |
152 | } | |
153 | else if( ccestimator == kZNA ){ | |
154 | if ( cc == kpPb020 ) { | |
155 | tab = 0.164; tabUnc = 0.010724; | |
156 | } else if ( cc == kpPb2040 ) { | |
157 | tab = 0.137; tabUnc = 0.005099; | |
158 | } else if ( cc == kpPb4060 ) { | |
159 | tab = 0.1011; tabUnc = 0.006; | |
160 | } else if ( cc == kpPb60100 ) { | |
161 | tab = 0.0459; tabUnc = 0.003162; | |
162 | } | |
163 | } | |
164 | ||
165 | tab *= 1e-9; // to pass from mb^{-1} to pb^{-1} | |
166 | tabUnc *= 1e-9; | |
167 | ||
168 | ||
169 | ||
170 | // | |
171 | // Get the histograms from the files | |
172 | // | |
173 | TH1D *hDirectMCpt=0; // Input MC c-->D spectra | |
174 | TH1D *hFeedDownMCpt=0; // Input MC b-->D spectra | |
175 | TH1D *hDirectMCptMax=0; // Input MC maximum c-->D spectra | |
176 | TH1D *hDirectMCptMin=0; // Input MC minimum c-->D spectra | |
177 | TH1D *hFeedDownMCptMax=0; // Input MC maximum b-->D spectra | |
178 | TH1D *hFeedDownMCptMin=0; // Input MC minimum b-->D spectra | |
179 | // TGraphAsymmErrors *gPrediction=0; // Input MC c-->D spectra | |
180 | TH1D *hDirectEffpt=0; // c-->D Acceptance and efficiency correction | |
181 | TH1D *hFeedDownEffpt=0; // b-->D Acceptance and efficiency correction | |
182 | TH1D *hRECpt=0; // all reconstructed D | |
183 | ||
184 | // | |
185 | // Define/Get the input histograms | |
186 | // | |
187 | Int_t decay=0; | |
188 | TFile * mcfile = new TFile(mcfilename,"read"); | |
189 | if (isD0Kpi){ | |
190 | decay = 1; | |
191 | hDirectMCpt = (TH1D*)mcfile->Get("hD0Kpipred_central"); | |
192 | hFeedDownMCpt = (TH1D*)mcfile->Get("hD0KpifromBpred_central_corr"); | |
193 | hDirectMCptMax = (TH1D*)mcfile->Get("hD0Kpipred_max"); | |
194 | hDirectMCptMin = (TH1D*)mcfile->Get("hD0Kpipred_min"); | |
195 | hFeedDownMCptMax = (TH1D*)mcfile->Get("hD0KpifromBpred_max_corr"); | |
196 | hFeedDownMCptMin = (TH1D*)mcfile->Get("hD0KpifromBpred_min_corr"); | |
197 | // gPrediction = (TGraphAsymmErrors*)mcfile->Get("D0Kpiprediction"); | |
198 | } | |
199 | else if (isDplusKpipi){ | |
200 | decay = 2; | |
201 | hDirectMCpt = (TH1D*)mcfile->Get("hDpluskpipipred_central"); | |
202 | hFeedDownMCpt = (TH1D*)mcfile->Get("hDpluskpipifromBpred_central_corr"); | |
203 | hDirectMCptMax = (TH1D*)mcfile->Get("hDpluskpipipred_max"); | |
204 | hDirectMCptMin = (TH1D*)mcfile->Get("hDpluskpipipred_min"); | |
205 | hFeedDownMCptMax = (TH1D*)mcfile->Get("hDpluskpipifromBpred_max_corr"); | |
206 | hFeedDownMCptMin = (TH1D*)mcfile->Get("hDpluskpipifromBpred_min_corr"); | |
207 | // gPrediction = (TGraphAsymmErrors*)mcfile->Get("Dpluskpipiprediction"); | |
208 | } | |
209 | else if(isDstarD0pi){ | |
210 | decay = 3; | |
211 | hDirectMCpt = (TH1D*)mcfile->Get("hDstarD0pipred_central"); | |
212 | hFeedDownMCpt = (TH1D*)mcfile->Get("hDstarD0pifromBpred_central_corr"); | |
213 | hDirectMCptMax = (TH1D*)mcfile->Get("hDstarD0pipred_max"); | |
214 | hDirectMCptMin = (TH1D*)mcfile->Get("hDstarD0pipred_min"); | |
215 | hFeedDownMCptMax = (TH1D*)mcfile->Get("hDstarD0pifromBpred_max_corr"); | |
216 | hFeedDownMCptMin = (TH1D*)mcfile->Get("hDstarD0pifromBpred_min_corr"); | |
217 | // gPrediction = (TGraphAsymmErrors*)mcfile->Get("DstarD0piprediction"); | |
218 | } | |
219 | else if (isDsKKpi){ | |
220 | decay = 4; | |
221 | hDirectMCpt = (TH1D*)mcfile->Get("hDsPhipitoKkpipred_central"); | |
222 | hFeedDownMCpt = (TH1D*)mcfile->Get("hDsPhipitoKkpifromBpred_central_corr"); | |
223 | hDirectMCptMax = (TH1D*)mcfile->Get("hDsPhipitoKkpipred_max"); | |
224 | hDirectMCptMin = (TH1D*)mcfile->Get("hDsPhipitoKkpipred_min"); | |
225 | hFeedDownMCptMax = (TH1D*)mcfile->Get("hDsPhipitoKkpifromBpred_max_corr"); | |
226 | hFeedDownMCptMin = (TH1D*)mcfile->Get("hDsPhipitoKkpifromBpred_min_corr"); | |
227 | } | |
228 | else if (isLctopKpi){ | |
229 | decay = 5; | |
230 | hDirectMCpt = (TH1D*)mcfile->Get("hLcpkpipred_central"); | |
231 | hFeedDownMCpt = (TH1D*)mcfile->Get("hLcpkpifromBpred_central_corr"); | |
232 | hDirectMCptMax = (TH1D*)mcfile->Get("hLcpkpipred_max"); | |
233 | hDirectMCptMin = (TH1D*)mcfile->Get("hLcpkpipred_min"); | |
234 | hFeedDownMCptMax = (TH1D*)mcfile->Get("hLcpkpifromBpred_max_corr"); | |
235 | hFeedDownMCptMin = (TH1D*)mcfile->Get("hLcpkpifromBpred_min_corr"); | |
236 | } | |
237 | else if (isLcK0Sp){ | |
238 | decay = 6; | |
239 | hDirectMCpt = (TH1D*)mcfile->Get("hLcK0sppred_central"); | |
240 | hFeedDownMCpt = (TH1D*)mcfile->Get("hLcK0spfromBpred_central_corr"); | |
241 | hDirectMCptMax = (TH1D*)mcfile->Get("hLcK0sppred_max"); | |
242 | hDirectMCptMin = (TH1D*)mcfile->Get("hLcK0sppred_min"); | |
243 | hFeedDownMCptMax = (TH1D*)mcfile->Get("hLcK0spfromBpred_max_corr"); | |
244 | hFeedDownMCptMin = (TH1D*)mcfile->Get("hLcK0spfromBpred_min_corr"); | |
245 | } | |
246 | // | |
247 | hDirectMCpt->SetNameTitle("hDirectMCpt","direct MC spectra"); | |
248 | hFeedDownMCpt->SetNameTitle("hFeedDownMCpt","feed-down MC spectra"); | |
249 | hDirectMCptMax->SetNameTitle("hDirectMCptMax","max direct MC spectra"); | |
250 | hDirectMCptMin->SetNameTitle("hDirectMCptMin","min direct MC spectra"); | |
251 | hFeedDownMCptMax->SetNameTitle("hFeedDownMCptMax","max feed-down MC spectra"); | |
252 | hFeedDownMCptMin->SetNameTitle("hFeedDownMCptMin","min feed-down MC spectra"); | |
253 | // | |
254 | // Scale FONLL inputs if we do the analysis in y-slices | |
255 | // | |
256 | if(rapiditySlice!=kdefault){ | |
257 | Double_t scaleFONLL = 1.0; | |
258 | switch(rapiditySlice) { | |
259 | case k08to04: scaleFONLL = (0.093+0.280)/1.0; break; | |
260 | case k07to04: scaleFONLL = 0.280/1.0; break; | |
261 | case k04to01: scaleFONLL = 0.284/1.0; break; | |
262 | case k01to01: scaleFONLL = 0.191/1.0; break; | |
263 | case k01to04: scaleFONLL = 0.288/1.0; break; | |
264 | case k04to07: scaleFONLL = 0.288/1.0; break; | |
265 | case k04to08: scaleFONLL = (0.288+0.096)/1.0; break; | |
266 | } | |
267 | hDirectMCpt->Scale(scaleFONLL); | |
268 | hDirectMCptMax->Scale(scaleFONLL); | |
269 | hDirectMCptMin->Scale(scaleFONLL); | |
270 | switch(rapiditySlice) { | |
271 | case k08to04: scaleFONLL = (0.089+0.274)/1.0; break; | |
272 | case k07to04: scaleFONLL = 0.274/1.0; break; | |
273 | case k04to01: scaleFONLL = 0.283/1.0; break; | |
274 | case k01to01: scaleFONLL = 0.192/1.0; break; | |
275 | case k01to04: scaleFONLL = 0.290/1.0; break; | |
276 | case k04to07: scaleFONLL = 0.291/1.0; break; | |
277 | case k04to08: scaleFONLL = (0.291+0.097)/1.0; break; | |
278 | } | |
279 | hFeedDownMCpt->Scale(scaleFONLL); | |
280 | hFeedDownMCptMax->Scale(scaleFONLL); | |
281 | hFeedDownMCptMin->Scale(scaleFONLL); | |
282 | } | |
283 | ||
284 | // | |
285 | // | |
286 | TFile * efffile = new TFile(efffilename,"read"); | |
287 | hDirectEffpt = (TH1D*)efffile->Get("hEffD"); | |
288 | hDirectEffpt->SetNameTitle("hDirectEffpt","direct acc x eff"); | |
289 | hFeedDownEffpt = (TH1D*)efffile->Get("hEffB"); | |
290 | hFeedDownEffpt->SetNameTitle("hFeedDownEffpt","feed-down acc x eff"); | |
291 | // | |
292 | // | |
293 | TFile * recofile = new TFile(recofilename,"read"); | |
294 | hRECpt = (TH1D*)recofile->Get(recohistoname); | |
295 | hRECpt->SetNameTitle("hRECpt","Reconstructed spectra"); | |
296 | ||
297 | // | |
298 | // Read the file of the EP resolution correction | |
299 | TFile *EPf=0; | |
300 | TH1D *hEPresolCorr=0; | |
301 | if(isRaavsEP>0.){ | |
302 | EPf = new TFile(epResolfile,"read"); | |
303 | if(isRaavsEP==kInPlane) hEPresolCorr = (TH1D*)EPf->Get("hCorrEPresol_InPlane"); | |
304 | else if(isRaavsEP==kOutOfPlane) hEPresolCorr = (TH1D*)EPf->Get("hCorrEPresol_OutOfPlane"); | |
305 | for(Int_t i=1; i<=hRECpt->GetNbinsX(); i++) { | |
306 | Double_t value = hRECpt->GetBinContent(i); | |
307 | Double_t error = hRECpt->GetBinError(i); | |
308 | Double_t pt = hRECpt->GetBinCenter(i); | |
309 | Int_t epbin = hEPresolCorr->FindBin( pt ); | |
310 | Double_t epcorr = hEPresolCorr->GetBinContent( epbin ); | |
311 | value = value*epcorr; | |
312 | error = error*epcorr; | |
313 | hRECpt->SetBinContent(i,value); | |
314 | hRECpt->SetBinError(i,error); | |
315 | } | |
316 | } | |
317 | ||
318 | // | |
319 | // Define the output histograms | |
320 | // | |
321 | TFile *out = new TFile(outfilename,"recreate"); | |
322 | // | |
323 | TH1D *histofc=0; | |
324 | TH1D *histofcMax=0; | |
325 | TH1D *histofcMin=0; | |
326 | TH1D *histoYieldCorr=0; | |
327 | TH1D *histoYieldCorrMax=0; | |
328 | TH1D *histoYieldCorrMin=0; | |
329 | TH1D *histoSigmaCorr=0; | |
330 | TH1D *histoSigmaCorrMax=0; | |
331 | TH1D *histoSigmaCorrMin=0; | |
332 | // | |
333 | TH2D *histofcRcb=0; | |
334 | TH1D *histofcRcb_px=0; | |
335 | TH2D *histoYieldCorrRcb=0; | |
336 | TH2D *histoSigmaCorrRcb=0; | |
337 | // | |
338 | TGraphAsymmErrors * gYieldCorr = 0; | |
339 | TGraphAsymmErrors * gSigmaCorr = 0; | |
340 | TGraphAsymmErrors * gFcExtreme = 0; | |
341 | TGraphAsymmErrors * gFcConservative = 0; | |
342 | TGraphAsymmErrors * gYieldCorrExtreme = 0; | |
343 | TGraphAsymmErrors * gSigmaCorrExtreme = 0; | |
344 | TGraphAsymmErrors * gYieldCorrConservative = 0; | |
345 | TGraphAsymmErrors * gSigmaCorrConservative = 0; | |
346 | // | |
347 | TNtuple * nSigma = 0; | |
348 | ||
349 | ||
350 | // | |
351 | // Main functionalities for the calculation | |
352 | // | |
353 | ||
354 | // Define and set the basic option flags | |
355 | AliHFPtSpectrum * spectra = new AliHFPtSpectrum("AliHFPtSpectrum","AliHFPtSpectrum",option); | |
356 | spectra->SetFeedDownCalculationOption(option); | |
357 | spectra->SetComputeAsymmetricUncertainties(asym); | |
358 | // Set flag on whether to additional PbPb Eloss hypothesis have to be computed | |
359 | spectra->SetComputeElossHypothesis(PbPbEloss); | |
360 | ||
361 | // Feed the input histograms | |
362 | // reconstructed spectra | |
363 | cout << " Setting the reconstructed spectrum,"; | |
364 | spectra->SetReconstructedSpectrum(hRECpt); | |
365 | // acceptance and efficiency corrections | |
366 | cout << " the efficiency,"; | |
367 | spectra->SetAccEffCorrection(hDirectEffpt,hFeedDownEffpt); | |
368 | // spectra->SetfIsStatUncEff(false); | |
369 | // option specific histos (theory) | |
370 | cout << " the theoretical spectra"; | |
371 | if(option==1){ | |
372 | spectra->SetMCptSpectra(hDirectMCpt,hFeedDownMCpt); | |
373 | if(asym) spectra->SetMCptDistributionsBounds(hDirectMCptMax,hDirectMCptMin,hFeedDownMCptMax,hFeedDownMCptMin); | |
374 | } | |
375 | else if(option==2){ | |
376 | spectra->SetFeedDownMCptSpectra(hFeedDownMCpt); | |
377 | if(asym) spectra->SetFeedDownMCptDistributionsBounds(hFeedDownMCptMax,hFeedDownMCptMin); | |
378 | } | |
379 | ||
380 | cout << " and the normalization" <<endl; | |
381 | // Set normalization factors (uncertainties set to 0. as example) | |
382 | spectra->SetNormalization(nevents,sigma); | |
383 | Double_t lumi = nevents / sigma ; | |
384 | Double_t lumiUnc = 0.04*lumi; // 10% uncertainty on the luminosity | |
385 | spectra->SetLuminosity(lumi,lumiUnc); | |
386 | Double_t effTrig = 1.0; | |
387 | spectra->SetTriggerEfficiency(effTrig,0.); | |
388 | if(isRaavsEP>0.) spectra->SetIsEventPlaneAnalysis(kTRUE); | |
389 | ||
390 | // Set the global uncertainties on the efficiencies (in percent) | |
391 | Double_t globalEffUnc = 0.15; | |
392 | Double_t globalBCEffRatioUnc = 0.15; | |
393 | spectra->SetAccEffPercentageUncertainty(globalEffUnc,globalBCEffRatioUnc); | |
394 | ||
395 | // Set if the yield is for particle+anti-particle or only one type | |
396 | spectra->SetIsParticlePlusAntiParticleYield(isParticlePlusAntiParticleYield); | |
397 | ||
398 | // Set the Tab parameter and uncertainties | |
399 | if ( (cc != kpp7) && (cc != kpp276) ) { | |
400 | spectra->SetTabParameter(tab,tabUnc); | |
401 | } | |
402 | ||
403 | // Do the calculations | |
404 | cout << " Doing the calculation... "<< endl; | |
405 | Double_t deltaY = 1.0; | |
406 | Double_t branchingRatioC = 1.0; | |
407 | Double_t branchingRatioBintoFinalDecay = 1.0; // this is relative to the input theoretical prediction | |
408 | spectra->ComputeHFPtSpectrum(deltaY,branchingRatioC,branchingRatioBintoFinalDecay); | |
409 | cout << " ended the calculation, getting the histograms back " << endl; | |
410 | ||
411 | // Set the systematics externally | |
412 | ||
413 | Bool_t combineFeedDown = true; | |
414 | AliHFSystErr *systematics = new AliHFSystErr(); | |
415 | if( cc==kpp276 ) { | |
416 | systematics->SetIsLowEnergy(true); | |
417 | } | |
418 | else if ( cc == kpPb0100 || cc == kpPb020 || cc == kpPb2040 || cc == kpPb4060 || cc == kpPb60100 ) { | |
419 | systematics->SetCollisionType(2); | |
420 | if(ccestimator==kV0A) { | |
421 | if(cc == kpPb020) systematics->SetCentrality("020V0A"); | |
422 | else if(cc == kpPb2040) systematics->SetCentrality("2040V0A"); | |
423 | else if(cc == kpPb4060) systematics->SetCentrality("4060V0A"); | |
424 | else if(cc == kpPb60100) systematics->SetCentrality("60100V0A"); | |
425 | } else if (ccestimator==kZNA) { | |
426 | if(cc == kpPb020) systematics->SetCentrality("020ZNA"); | |
427 | else if(cc == kpPb2040) systematics->SetCentrality("2040ZNA"); | |
428 | else if(cc == kpPb4060) systematics->SetCentrality("4060ZNA"); | |
429 | else if(cc == kpPb60100) systematics->SetCentrality("60100ZNA"); | |
430 | } else { | |
431 | if(!(cc == kpPb0100)) { | |
432 | cout <<" Error on the pPb options"<<endl; | |
433 | return; | |
434 | } | |
435 | } | |
436 | } | |
437 | // | |
438 | else if( cc!=kpp7 ) { | |
439 | systematics->SetCollisionType(1); | |
440 | if ( cc == k07half ) systematics->SetCentrality("07half"); | |
441 | else if ( cc == k010 ) systematics->SetCentrality("010"); | |
442 | else if ( cc == k1020 ) systematics->SetCentrality("1020"); | |
443 | else if ( cc == k020 ) systematics->SetCentrality("020"); | |
444 | else if ( cc == k2040 || cc == k2030 || cc == k3040 ) { | |
445 | systematics->SetCentrality("2040"); | |
446 | systematics->SetIsPbPb2010EnergyScan(true); | |
447 | } | |
448 | else if ( cc == k3050 ) { | |
449 | if (isRaavsEP == kPhiIntegrated) systematics->SetCentrality("4080"); | |
450 | else if (isRaavsEP == kInPlane) systematics->SetCentrality("3050InPlane"); | |
451 | else if (isRaavsEP == kOutOfPlane) systematics->SetCentrality("3050OutOfPlane"); | |
452 | } | |
453 | else if ( cc == k4060 || cc == k4050 || cc == k5060 ) systematics->SetCentrality("4060"); | |
454 | else if ( cc == k6080 ) systematics->SetCentrality("6080"); | |
455 | else if ( cc == k4080 ) systematics->SetCentrality("4080"); | |
456 | else { | |
457 | cout << " Systematics not yet implemented " << endl; | |
458 | return; | |
459 | } | |
460 | } else { systematics->SetCollisionType(0); } | |
461 | // | |
462 | systematics->Init(decay); | |
463 | spectra->ComputeSystUncertainties(systematics,combineFeedDown); | |
464 | ||
465 | // | |
466 | // Get the output histograms | |
467 | // | |
468 | // the corrected yield and cross-section | |
469 | histoYieldCorr = (TH1D*)spectra->GetHistoFeedDownCorrectedSpectrum(); | |
470 | histoSigmaCorr = (TH1D*)spectra->GetHistoCrossSectionFromYieldSpectrum(); | |
471 | histoYieldCorrMax = (TH1D*)spectra->GetHistoUpperLimitFeedDownCorrectedSpectrum(); | |
472 | histoYieldCorrMin = (TH1D*)spectra->GetHistoLowerLimitFeedDownCorrectedSpectrum(); | |
473 | histoSigmaCorrMax = (TH1D*)spectra->GetHistoUpperLimitCrossSectionFromYieldSpectrum(); | |
474 | histoSigmaCorrMin = (TH1D*)spectra->GetHistoLowerLimitCrossSectionFromYieldSpectrum(); | |
475 | histoYieldCorr->SetNameTitle("histoYieldCorr","corrected yield"); | |
476 | histoYieldCorrMax->SetNameTitle("histoYieldCorrMax","max corrected yield"); | |
477 | histoYieldCorrMin->SetNameTitle("histoYieldCorrMin","min corrected yield"); | |
478 | histoSigmaCorr->SetNameTitle("histoSigmaCorr","corrected invariant cross-section"); | |
479 | histoSigmaCorrMax->SetNameTitle("histoSigmaCorrMax","max corrected invariant cross-section"); | |
480 | histoSigmaCorrMin->SetNameTitle("histoSigmaCorrMin","min corrected invariant cross-section"); | |
481 | // the efficiencies | |
482 | if(!hDirectEffpt) hDirectEffpt = (TH1D*)spectra->GetDirectAccEffCorrection(); | |
483 | if(!hFeedDownEffpt) hFeedDownEffpt = (TH1D*)spectra->GetFeedDownAccEffCorrection(); | |
484 | // Get the PbPb Eloss hypothesis histograms | |
485 | if(PbPbEloss){ | |
486 | histofcRcb = spectra->GetHistoFeedDownCorrectionFcVsEloss(); | |
487 | histoYieldCorrRcb = spectra->GetHistoFeedDownCorrectedSpectrumVsEloss(); | |
488 | histoSigmaCorrRcb = spectra->GetHistoCrossSectionFromYieldSpectrumVsEloss(); | |
489 | histofcRcb->SetName("histofcRcb"); | |
490 | histoYieldCorrRcb->SetName("histoYieldCorrRcb"); | |
491 | histoSigmaCorrRcb->SetName("histoSigmaCorrRcb"); | |
492 | } | |
493 | ||
494 | // Get & Rename the TGraphs | |
495 | gSigmaCorr = spectra->GetCrossSectionFromYieldSpectrum(); | |
496 | gYieldCorr = spectra->GetFeedDownCorrectedSpectrum(); | |
497 | if (asym) { | |
498 | gSigmaCorrExtreme = spectra->GetCrossSectionFromYieldSpectrumExtreme(); | |
499 | gYieldCorrExtreme = spectra->GetFeedDownCorrectedSpectrumExtreme(); | |
500 | gSigmaCorrConservative = spectra->GetCrossSectionFromYieldSpectrumConservative(); | |
501 | gYieldCorrConservative = spectra->GetFeedDownCorrectedSpectrumConservative(); | |
502 | } | |
503 | ||
504 | // Get & Rename the TGraphs | |
505 | if (option==0){ | |
506 | gYieldCorr->SetNameTitle("gYieldCorr","gYieldCorr (uncorr)"); | |
507 | gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (uncorr)"); | |
508 | } | |
509 | if (option==1){ | |
510 | // fc histos | |
511 | histofc = (TH1D*)spectra->GetHistoFeedDownCorrectionFc(); | |
512 | histofcMax = (TH1D*)spectra->GetHistoUpperLimitFeedDownCorrectionFc(); | |
513 | histofcMin = (TH1D*)spectra->GetHistoLowerLimitFeedDownCorrectionFc(); | |
514 | histofc->SetNameTitle("histofc","fc correction factor"); | |
515 | histofcMax->SetNameTitle("histofcMax","max fc correction factor"); | |
516 | histofcMin->SetNameTitle("histofcMin","min fc correction factor"); | |
517 | if (asym) { | |
518 | gYieldCorr->SetNameTitle("gYieldCorr","gYieldCorr (by fc)"); | |
519 | gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (by fc)"); | |
520 | gFcExtreme = spectra->GetFeedDownCorrectionFcExtreme(); | |
521 | gFcExtreme->SetNameTitle("gFcExtreme","gFcExtreme"); | |
522 | gYieldCorrExtreme->SetNameTitle("gYieldCorrExtreme","Extreme gYieldCorr (by fc)"); | |
523 | gSigmaCorrExtreme->SetNameTitle("gSigmaCorrExtreme","Extreme gSigmaCorr (by fc)"); | |
524 | gFcConservative = spectra->GetFeedDownCorrectionFcConservative(); | |
525 | gFcConservative->SetNameTitle("gFcConservative","gFcConservative"); | |
526 | gYieldCorrConservative->SetNameTitle("gYieldCorrConservative","Conservative gYieldCorr (by fc)"); | |
527 | gSigmaCorrConservative->SetNameTitle("gSigmaCorrConservative","Conservative gSigmaCorr (by fc)"); | |
528 | } | |
529 | } | |
530 | if (option==2 && asym) { | |
531 | gYieldCorr->SetNameTitle("gYieldCorr","gYieldCorr (by Nb)"); | |
532 | gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (by Nb)"); | |
533 | gYieldCorrExtreme->SetNameTitle("gYieldCorrExtreme","Extreme gYieldCorr (by Nb)"); | |
534 | gSigmaCorrExtreme->SetNameTitle("gSigmaCorrExtreme","Extreme gSigmaCorr (by Nb)"); | |
535 | gYieldCorrConservative->SetNameTitle("gYieldCorrConservative","Conservative gYieldCorr (by Nb)"); | |
536 | gSigmaCorrConservative->SetNameTitle("gSigmaCorrConservative","Conservative gSigmaCorr (by Nb)"); | |
537 | gFcConservative = spectra->GetFeedDownCorrectionFcConservative(); | |
538 | gFcConservative->SetNameTitle("gFcConservative","gFcConservative"); | |
539 | } | |
540 | ||
541 | if(PbPbEloss){ | |
542 | nSigma = spectra->GetNtupleCrossSectionVsEloss(); | |
543 | } | |
544 | ||
545 | // | |
546 | // Now, plot the results ! :) | |
547 | // | |
548 | ||
549 | gROOT->SetStyle("Plain"); | |
550 | ||
551 | cout << " Drawing the results ! " << endl; | |
552 | ||
553 | // control plots | |
554 | if (option==1) { | |
555 | ||
556 | TCanvas *ceff = new TCanvas("ceff","efficiency drawing"); | |
557 | ceff->Divide(1,2); | |
558 | ceff->cd(1); | |
559 | hDirectEffpt->Draw(); | |
560 | ceff->cd(2); | |
561 | hFeedDownEffpt->Draw(); | |
562 | ceff->Update(); | |
563 | ||
564 | TCanvas *cTheoryRebin = new TCanvas("cTheoryRebin","control the theoretical spectra rebin"); | |
565 | cTheoryRebin->Divide(1,2); | |
566 | cTheoryRebin->cd(1); | |
567 | hDirectMCpt->Draw(""); | |
568 | TH1D *hDirectMCptRebin = (TH1D*)spectra->GetDirectTheoreticalSpectrum(); | |
569 | hDirectMCptRebin->SetLineColor(2); | |
570 | hDirectMCptRebin->Draw("same"); | |
571 | cTheoryRebin->cd(2); | |
572 | hFeedDownMCpt->Draw(""); | |
573 | TH1D *hFeedDownRebin = (TH1D*)spectra->GetFeedDownTheoreticalSpectrum(); | |
574 | hFeedDownRebin->SetLineColor(2); | |
575 | hFeedDownRebin->Draw("same"); | |
576 | cTheoryRebin->Update(); | |
577 | ||
578 | TCanvas *cTheoryRebinLimits = new TCanvas("cTheoryRebinLimits","control the theoretical spectra limits rebin"); | |
579 | cTheoryRebinLimits->Divide(1,2); | |
580 | cTheoryRebinLimits->cd(1); | |
581 | hDirectMCptMax->Draw(""); | |
582 | TH1D *hDirectMCptMaxRebin = (TH1D*)spectra->GetDirectTheoreticalUpperLimitSpectrum(); | |
583 | hDirectMCptMaxRebin->SetLineColor(2); | |
584 | hDirectMCptMaxRebin->Draw("same"); | |
585 | hDirectMCptMin->Draw("same"); | |
586 | TH1D *hDirectMCptMinRebin = (TH1D*)spectra->GetDirectTheoreticalLowerLimitSpectrum(); | |
587 | hDirectMCptMinRebin->SetLineColor(2); | |
588 | hDirectMCptMinRebin->Draw("same"); | |
589 | cTheoryRebinLimits->cd(2); | |
590 | hFeedDownMCptMax->Draw(""); | |
591 | TH1D *hFeedDownMaxRebin = (TH1D*)spectra->GetFeedDownTheoreticalUpperLimitSpectrum(); | |
592 | hFeedDownMaxRebin->SetLineColor(2); | |
593 | hFeedDownMaxRebin->Draw("same"); | |
594 | hFeedDownMCptMin->Draw("same"); | |
595 | TH1D *hFeedDownMinRebin = (TH1D*)spectra->GetFeedDownTheoreticalLowerLimitSpectrum(); | |
596 | hFeedDownMinRebin->SetLineColor(2); | |
597 | hFeedDownMinRebin->Draw("same"); | |
598 | cTheoryRebinLimits->Update(); | |
599 | } | |
600 | ||
601 | if (option==1) { | |
602 | ||
603 | TCanvas * cfc = new TCanvas("cfc","Fc"); | |
604 | histofcMax->Draw("c"); | |
605 | histofc->Draw("csame"); | |
606 | histofcMin->Draw("csame"); | |
607 | cfc->Update(); | |
608 | ||
609 | if (asym) { | |
610 | TH2F *histofcDraw= new TH2F("histofcDraw","histofc (for drawing)",100,0,33.25,100,0.01,1.25); | |
611 | histofcDraw->SetStats(0); | |
612 | histofcDraw->GetXaxis()->SetTitle("p_{T} [GeV]"); | |
613 | histofcDraw ->GetXaxis()->SetTitleSize(0.05); | |
614 | histofcDraw->GetXaxis()->SetTitleOffset(0.95); | |
615 | histofcDraw->GetYaxis()->SetTitle(" fc "); | |
616 | histofcDraw->GetYaxis()->SetTitleSize(0.05); | |
617 | ||
618 | if (gFcExtreme){ | |
619 | ||
620 | // for(Int_t item=0; item<gSigmaCorr->GetN(); item++){ | |
621 | // Double_t center=0., value=0.; | |
622 | // gFcExtreme->GetPoint(item,center,value); | |
623 | // Double_t highunc = gFcExtreme->GetErrorYhigh(item) / value ; | |
624 | // Double_t lowunc = gFcExtreme->GetErrorYlow(item) / value ; | |
625 | // cout << "Fc extreme: i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl; | |
626 | // } | |
627 | // for(Int_t item=0; item<gSigmaCorr->GetN(); item++){ | |
628 | // Double_t center=0., value=0.; | |
629 | // gFcConservative->GetPoint(item,center,value); | |
630 | // Double_t highunc = gFcConservative->GetErrorYhigh(item) / value ; | |
631 | // Double_t lowunc = gFcConservative->GetErrorYlow(item) / value ; | |
632 | // cout << "Fc conservative: i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl; | |
633 | // } | |
634 | TCanvas *cfcExtreme = new TCanvas("cfcExtreme","Extreme Asymmetric fc (TGraphAsymmErr)"); | |
635 | gFcExtreme->SetFillStyle(3006); | |
636 | gFcExtreme->SetLineWidth(3); | |
637 | gFcExtreme->SetMarkerStyle(20); | |
638 | gFcExtreme->SetFillColor(2); | |
639 | histofcDraw->Draw(); | |
640 | gFcExtreme->Draw("3same"); | |
641 | ||
642 | if(gFcConservative){ | |
643 | gFcConservative->SetFillStyle(3007); | |
644 | gFcConservative->SetFillColor(4); | |
645 | gFcConservative->Draw("3same"); | |
646 | } | |
647 | ||
648 | cfcExtreme->Update(); | |
649 | } | |
650 | } | |
651 | ||
652 | } | |
653 | ||
654 | // | |
655 | // Drawing the results (the raw-reconstructed, the expected, and the corrected spectra) | |
656 | // | |
657 | TCanvas * cresult = new TCanvas("cresult","corrected yields & sigma"); | |
658 | hDirectMCpt->SetMarkerStyle(20); | |
659 | hDirectMCpt->SetMarkerColor(4); | |
660 | hDirectMCpt->Draw("p"); | |
661 | histoSigmaCorr->SetMarkerStyle(21); | |
662 | histoSigmaCorr->SetMarkerColor(2); | |
663 | histoSigmaCorr->Draw("psame"); | |
664 | histoYieldCorr->SetMarkerStyle(22); | |
665 | histoYieldCorr->SetMarkerColor(6); | |
666 | histoYieldCorr->Draw("psame"); | |
667 | hRECpt->SetMarkerStyle(23); | |
668 | hRECpt->SetMarkerColor(3); | |
669 | hRECpt->Draw("psame"); | |
670 | cresult->SetLogy(); | |
671 | cresult->Update(); | |
672 | ||
673 | TCanvas * cresult2 = new TCanvas("cresult2","corrected yield & sigma"); | |
674 | histoSigmaCorr->SetMarkerStyle(21); | |
675 | histoSigmaCorr->SetMarkerColor(2); | |
676 | histoSigmaCorr->Draw("p"); | |
677 | histoYieldCorr->SetMarkerStyle(22); | |
678 | histoYieldCorr->SetMarkerColor(6); | |
679 | histoYieldCorr->Draw("psame"); | |
680 | hRECpt->SetMarkerStyle(23); | |
681 | hRECpt->SetMarkerColor(3); | |
682 | hRECpt->Draw("psame"); | |
683 | cresult2->SetLogy(); | |
684 | cresult2->Update(); | |
685 | ||
686 | ||
687 | if (asym) { | |
688 | ||
689 | TH2F *histoDraw = new TH2F("histoDraw","histo (for drawing)",100,0,33.25,100,50.,1e7); | |
690 | float max = 1.1*gYieldCorr->GetMaximum(); | |
691 | histoDraw->SetAxisRange(0.1,max,"Y"); | |
692 | histoDraw->SetStats(0); | |
693 | histoDraw->GetXaxis()->SetTitle("p_{T} [GeV]"); | |
694 | histoDraw->GetXaxis()->SetTitleSize(0.05); | |
695 | histoDraw->GetXaxis()->SetTitleOffset(0.95); | |
696 | histoDraw->GetYaxis()->SetTitle("#frac{d#N}{dp_{T}} |_{|y|<1} [L & trigger uncorr]"); | |
697 | histoDraw->GetYaxis()->SetTitleSize(0.05); | |
698 | TCanvas * cyieldAsym = new TCanvas("cyieldAsym","Asymmetric corrected yield (TGraphAsymmErr)"); | |
699 | gYieldCorr->SetFillStyle(3001); | |
700 | gYieldCorr->SetLineWidth(3); | |
701 | gYieldCorr->SetMarkerStyle(20); | |
702 | gYieldCorr->SetFillColor(3); | |
703 | histoDraw->Draw(); | |
704 | gYieldCorr->Draw("3LPsame"); | |
705 | gYieldCorr->Draw("Xsame"); | |
706 | cyieldAsym->SetLogy(); | |
707 | cyieldAsym->Update(); | |
708 | ||
709 | TCanvas * cyieldExtreme = new TCanvas("cyieldExtreme","Extreme Asymmetric corrected yield (TGraphAsymmErr)"); | |
710 | histoYieldCorr->Draw(); | |
711 | gYieldCorrExtreme->SetFillStyle(3002); | |
712 | gYieldCorrExtreme->SetLineWidth(3); | |
713 | gYieldCorrExtreme->SetMarkerStyle(20); | |
714 | gYieldCorrExtreme->SetFillColor(2); | |
715 | histoYieldCorr->Draw(); | |
716 | gYieldCorr->Draw("3same"); | |
717 | gYieldCorrExtreme->Draw("3same"); | |
718 | cyieldExtreme->SetLogy(); | |
719 | cyieldExtreme->Update(); | |
720 | ||
721 | TH2F *histo2Draw = new TH2F("histo2Draw","histo2 (for drawing)",100,0,33.25,100,50.,1e9); | |
722 | max = 1.1*gSigmaCorr->GetMaximum(); | |
723 | histo2Draw->SetAxisRange(0.1,max,"Y"); | |
724 | histo2Draw->SetStats(0); | |
725 | histo2Draw->GetXaxis()->SetTitle("p_{T} [GeV]"); | |
726 | histo2Draw->GetXaxis()->SetTitleSize(0.05); | |
727 | histo2Draw->GetXaxis()->SetTitleOffset(0.95); | |
728 | histo2Draw->GetYaxis()->SetTitle("#frac{1}{BR} #times #frac{d#sigma}{dp_{T}} |_{|y|<1}"); | |
729 | histo2Draw->GetYaxis()->SetTitleSize(0.05); | |
730 | TCanvas * csigmaAsym = new TCanvas("csigmaAsym","Asymmetric corrected sigma (TGraphAsymmErr)"); | |
731 | gSigmaCorr->SetFillStyle(3001); | |
732 | gSigmaCorr->SetLineWidth(3); | |
733 | gSigmaCorr->SetMarkerStyle(21); | |
734 | gSigmaCorr->SetFillColor(3); | |
735 | histo2Draw->Draw(); | |
736 | gSigmaCorr->Draw("3LPsame"); | |
737 | gSigmaCorr->Draw("Xsame"); | |
738 | csigmaAsym->SetLogy(); | |
739 | csigmaAsym->Update(); | |
740 | ||
741 | // cout << endl <<" Sytematics (stat approach) " <<endl; | |
742 | // for(Int_t item=0; item<gSigmaCorr->GetN(); item++){ | |
743 | // Double_t center=0., value=0.; | |
744 | // gSigmaCorr->GetPoint(item,center,value); | |
745 | // Double_t highunc = gSigmaCorr->GetErrorYhigh(item) / value ; | |
746 | // Double_t lowunc = gSigmaCorr->GetErrorYlow(item) / value ; | |
747 | // cout << "Sigma syst (stat), i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl; | |
748 | // } | |
749 | ||
750 | TCanvas * csigmaExtreme = new TCanvas("csigmaExtreme","Asymmetric extreme corrected sigma (TGraphAsymmErr)"); | |
751 | histoSigmaCorr->Draw(); | |
752 | gSigmaCorr->Draw("3Psame"); | |
753 | gSigmaCorrExtreme->SetFillStyle(3002); | |
754 | gSigmaCorrExtreme->SetLineWidth(3); | |
755 | gSigmaCorrExtreme->SetMarkerStyle(21); | |
756 | gSigmaCorrExtreme->SetFillColor(2); | |
757 | gSigmaCorrExtreme->Draw("3Psame"); | |
758 | csigmaExtreme->SetLogy(); | |
759 | csigmaExtreme->Update(); | |
760 | ||
761 | // cout << endl << " Sytematics (Extreme approach)" <<endl; | |
762 | // for(Int_t item=0; item<gSigmaCorrExtreme->GetN(); item++){ | |
763 | // Double_t center=0., value=0.; | |
764 | // gSigmaCorrExtreme->GetPoint(item,center,value); | |
765 | // Double_t highunc = gSigmaCorrExtreme->GetErrorYhigh(item) / value ; | |
766 | // Double_t lowunc = gSigmaCorrExtreme->GetErrorYlow(item) / value ; | |
767 | // cout << "Sigma syst (extreme) i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl; | |
768 | // } | |
769 | ||
770 | // cout << endl << " Sytematics (Conservative approach)" <<endl; | |
771 | // for(Int_t item=0; item<gSigmaCorrConservative->GetN(); item++){ | |
772 | // Double_t center=0., value=0.; | |
773 | // gSigmaCorrConservative->GetPoint(item,center,value); | |
774 | // Double_t highunc = gSigmaCorrConservative->GetErrorYhigh(item) / value ; | |
775 | // Double_t lowunc = gSigmaCorrConservative->GetErrorYlow(item) / value ; | |
776 | // cout << "Sigma syst (conservative) i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl; | |
777 | // } | |
778 | ||
779 | } | |
780 | ||
781 | // Draw the PbPb Eloss hypothesis histograms | |
782 | if(PbPbEloss){ | |
783 | AliHFPtSpectrum *CalcBins=NULL; | |
784 | gStyle->SetPalette(1); | |
785 | TCanvas *canvasfcRcb = new TCanvas("canvasfcRcb","fc vs pt vs Rcb"); | |
786 | // histofcRcb->Draw("cont4z"); | |
787 | histofcRcb->Draw("colz"); | |
788 | canvasfcRcb->Update(); | |
789 | canvasfcRcb->cd(2); | |
790 | TCanvas *canvasfcRcb1 = new TCanvas("canvasfcRcb1","fc vs pt vs Rcb=1"); | |
791 | histofcRcb_px = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px",40,40); | |
792 | histofcRcb_px->SetLineColor(2); | |
793 | if (option==1) { | |
794 | histofc->Draw(); | |
795 | histofcRcb_px->Draw("same"); | |
796 | } else histofcRcb_px->Draw(""); | |
797 | canvasfcRcb1->Update(); | |
798 | TCanvas *canvasfcRcb2 = new TCanvas("canvasfcRcb2","fc vs pt vs Rcb fixed Rcb"); | |
799 | Int_t bin0 = CalcBins->FindTH2YBin(histofcRcb,0.25); | |
800 | Int_t bin1 = CalcBins->FindTH2YBin(histofcRcb,0.5); | |
801 | Int_t bin2 = CalcBins->FindTH2YBin(histofcRcb,1.0); | |
802 | Int_t bin3 = CalcBins->FindTH2YBin(histofcRcb,1.5); | |
803 | Int_t bin4 = CalcBins->FindTH2YBin(histofcRcb,2.0); | |
804 | Int_t bin5 = CalcBins->FindTH2YBin(histofcRcb,3.0); | |
805 | Int_t bin6 = CalcBins->FindTH2YBin(histofcRcb,4.0); | |
806 | TH1D * histofcRcb_px0a = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px0a",bin0,bin0); | |
807 | TH1D * histofcRcb_px0 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px0",bin1,bin1); | |
808 | TH1D * histofcRcb_px1 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px1",bin2,bin2); | |
809 | TH1D * histofcRcb_px2 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px2",bin3,bin3); | |
810 | TH1D * histofcRcb_px3 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px3",bin4,bin4); | |
811 | TH1D * histofcRcb_px4 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px4",bin5,bin5); | |
812 | TH1D * histofcRcb_px5 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px5",bin6,bin6); | |
813 | if (option==1) { | |
814 | histofc->Draw(); | |
815 | // histofcRcb_px->Draw("same"); | |
816 | } else { | |
817 | // histofcRcb_px->Draw(""); | |
818 | histofcRcb_px0a->SetLineColor(2); | |
819 | histofcRcb_px0a->Draw(""); | |
820 | } | |
821 | histofcRcb_px0a->SetLineColor(2); | |
822 | histofcRcb_px0a->Draw("same"); | |
823 | histofcRcb_px0->SetLineColor(4); | |
824 | histofcRcb_px0->Draw("same"); | |
825 | histofcRcb_px1->SetLineColor(3); | |
826 | histofcRcb_px1->Draw("same"); | |
827 | histofcRcb_px2->SetLineColor(kCyan); | |
828 | histofcRcb_px2->Draw("same"); | |
829 | histofcRcb_px3->SetLineColor(kMagenta+1); | |
830 | histofcRcb_px3->Draw("same"); | |
831 | histofcRcb_px4->SetLineColor(kOrange+7); | |
832 | histofcRcb_px4->Draw("same"); | |
833 | histofcRcb_px5->SetLineColor(kGreen+3); | |
834 | histofcRcb_px5->Draw("same"); | |
835 | TLegend *legrcc = new TLegend(0.8,0.8,0.95,0.9); | |
836 | legrcc->SetFillColor(0); | |
837 | if (option==1) { | |
838 | legrcc->AddEntry(histofcRcb_px0a,"Rc/b=0.25","l"); | |
839 | legrcc->AddEntry(histofcRcb_px0,"Rc/b=0.5","l"); | |
840 | legrcc->AddEntry(histofcRcb_px1,"Rc/b=1.0","l"); | |
841 | legrcc->AddEntry(histofcRcb_px2,"Rc/b=1.5","l"); | |
842 | legrcc->AddEntry(histofcRcb_px3,"Rc/b=2.0","l"); | |
843 | legrcc->AddEntry(histofcRcb_px4,"Rc/b=3.0","l"); | |
844 | legrcc->AddEntry(histofcRcb_px5,"Rc/b=4.0","l"); | |
845 | }else{ | |
846 | legrcc->AddEntry(histofcRcb_px0a,"Rb=0.25","l"); | |
847 | legrcc->AddEntry(histofcRcb_px0,"Rb=0.5","l"); | |
848 | legrcc->AddEntry(histofcRcb_px1,"Rb=1.0","l"); | |
849 | legrcc->AddEntry(histofcRcb_px2,"Rb=1.5","l"); | |
850 | legrcc->AddEntry(histofcRcb_px3,"Rb=2.0","l"); | |
851 | legrcc->AddEntry(histofcRcb_px4,"Rb=3.0","l"); | |
852 | legrcc->AddEntry(histofcRcb_px5,"Rb=4.0","l"); | |
853 | } | |
854 | legrcc->Draw(); | |
855 | canvasfcRcb2->Update(); | |
856 | TCanvas *canvasYRcb = new TCanvas("canvasYRcb","corrected yield vs pt vs Rcb"); | |
857 | histoYieldCorrRcb->Draw("cont4z"); | |
858 | canvasYRcb->Update(); | |
859 | TCanvas *canvasSRcb = new TCanvas("canvasSRcb","sigma vs pt vs Rcb"); | |
860 | histoSigmaCorrRcb->Draw("cont4z"); | |
861 | canvasSRcb->Update(); | |
862 | TCanvas *canvasSRcb1 = new TCanvas("canvasSRcb1","sigma vs pt vs Rcb fixed Rcb"); | |
863 | TH1D * histoSigmaCorrRcb_px0a = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px0a",bin0,bin0); | |
864 | TH1D * histoSigmaCorrRcb_px0 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px0",bin1,bin1); | |
865 | TH1D * histoSigmaCorrRcb_px1 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px1",bin2,bin2); | |
866 | TH1D * histoSigmaCorrRcb_px2 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px2",bin3,bin3); | |
867 | TH1D * histoSigmaCorrRcb_px3 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px3",bin4,bin4); | |
868 | TH1D * histoSigmaCorrRcb_px4 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px4",bin5,bin5); | |
869 | TH1D * histoSigmaCorrRcb_px5 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px5",bin6,bin6); | |
870 | histoSigmaCorr->Draw(); | |
871 | histoSigmaCorrRcb_px0a->SetLineColor(2); | |
872 | histoSigmaCorrRcb_px0a->Draw("hsame"); | |
873 | histoSigmaCorrRcb_px0->SetLineColor(4); | |
874 | histoSigmaCorrRcb_px0->Draw("hsame"); | |
875 | histoSigmaCorrRcb_px1->SetLineColor(3); | |
876 | histoSigmaCorrRcb_px1->Draw("hsame"); | |
877 | histoSigmaCorrRcb_px2->SetLineColor(kCyan); | |
878 | histoSigmaCorrRcb_px2->Draw("hsame"); | |
879 | histoSigmaCorrRcb_px3->SetLineColor(kMagenta+1); | |
880 | histoSigmaCorrRcb_px3->Draw("hsame"); | |
881 | histoSigmaCorrRcb_px4->SetLineColor(kOrange+7); | |
882 | histoSigmaCorrRcb_px4->Draw("same"); | |
883 | histoSigmaCorrRcb_px5->SetLineColor(kGreen+3); | |
884 | histoSigmaCorrRcb_px5->Draw("same"); | |
885 | TLegend *legrcb = new TLegend(0.8,0.8,0.95,0.9); | |
886 | legrcb->SetFillColor(0); | |
887 | if (option==1) { | |
888 | legrcb->AddEntry(histoSigmaCorrRcb_px0a,"Rc/b=0.25","l"); | |
889 | legrcb->AddEntry(histoSigmaCorrRcb_px0,"Rc/b=0.5","l"); | |
890 | legrcb->AddEntry(histoSigmaCorrRcb_px1,"Rc/b=1.0","l"); | |
891 | legrcb->AddEntry(histoSigmaCorrRcb_px2,"Rc/b=1.5","l"); | |
892 | legrcb->AddEntry(histoSigmaCorrRcb_px3,"Rc/b=2.0","l"); | |
893 | legrcb->AddEntry(histoSigmaCorrRcb_px4,"Rc/b=3.0","l"); | |
894 | legrcb->AddEntry(histoSigmaCorrRcb_px5,"Rc/b=4.0","l"); | |
895 | }else{ | |
896 | legrcb->AddEntry(histoSigmaCorrRcb_px0a,"Rb=0.25","l"); | |
897 | legrcb->AddEntry(histoSigmaCorrRcb_px0,"Rb=0.5","l"); | |
898 | legrcb->AddEntry(histoSigmaCorrRcb_px1,"Rb=1.0","l"); | |
899 | legrcb->AddEntry(histoSigmaCorrRcb_px2,"Rb=1.5","l"); | |
900 | legrcb->AddEntry(histoSigmaCorrRcb_px3,"Rb=2.0","l"); | |
901 | legrcb->AddEntry(histoSigmaCorrRcb_px4,"Rb=3.0","l"); | |
902 | legrcb->AddEntry(histoSigmaCorrRcb_px5,"Rb=4.0","l"); | |
903 | } | |
904 | legrcb->Draw(); | |
905 | canvasSRcb1->Update(); | |
906 | } | |
907 | ||
908 | ||
909 | // | |
910 | // Write the histograms to the output file | |
911 | // | |
912 | cout << " Saving the results ! " << endl<< endl; | |
913 | ||
914 | out->cd(); | |
915 | // | |
916 | hDirectMCpt->Write(); hFeedDownMCpt->Write(); | |
917 | hDirectMCptMax->Write(); hDirectMCptMin->Write(); | |
918 | hFeedDownMCptMax->Write(); hFeedDownMCptMin->Write(); | |
919 | if(hDirectEffpt) hDirectEffpt->Write(); if(hFeedDownEffpt) hFeedDownEffpt->Write(); | |
920 | hRECpt->Write(); | |
921 | // | |
922 | histoYieldCorr->Write(); | |
923 | histoYieldCorrMax->Write(); histoYieldCorrMin->Write(); | |
924 | histoSigmaCorr->Write(); | |
925 | histoSigmaCorrMax->Write(); histoSigmaCorrMin->Write(); | |
926 | ||
927 | if(PbPbEloss){ | |
928 | histofcRcb->Write(); histofcRcb_px->Write(); | |
929 | histoYieldCorrRcb->Write(); | |
930 | histoSigmaCorrRcb->Write(); | |
931 | nSigma->Write(); | |
932 | } | |
933 | ||
934 | gYieldCorr->Write(); | |
935 | gSigmaCorr->Write(); | |
936 | if(asym){ | |
937 | if(gYieldCorrExtreme) gYieldCorrExtreme->Write(); | |
938 | if(gSigmaCorrExtreme) gSigmaCorrExtreme->Write(); | |
939 | if(gYieldCorrConservative) gYieldCorrConservative->Write(); | |
940 | if(gSigmaCorrConservative) gSigmaCorrConservative->Write(); | |
941 | if(asym && gFcConservative) gFcConservative->Write(); | |
942 | } | |
943 | ||
944 | if(option==1){ | |
945 | histofc->Write(); | |
946 | histofcMax->Write(); histofcMin->Write(); | |
947 | if(asym && gFcExtreme) gFcExtreme->Write(); | |
948 | } | |
949 | ||
950 | ||
951 | TH1D * hStatUncEffcSigma = spectra->GetDirectStatEffUncOnSigma(); | |
952 | TH1D * hStatUncEffbSigma = spectra->GetFeedDownStatEffUncOnSigma(); | |
953 | hStatUncEffcSigma->Write(); | |
954 | hStatUncEffbSigma->Write(); | |
955 | if(option!=0){ | |
956 | TH1D * hStatUncEffcFD = spectra->GetDirectStatEffUncOnFc(); | |
957 | TH1D * hStatUncEffbFD = spectra->GetFeedDownStatEffUncOnFc(); | |
958 | hStatUncEffcFD->Write(); | |
959 | hStatUncEffbFD->Write(); | |
960 | } | |
961 | ||
962 | // Draw the cross-section | |
963 | // spectra->DrawSpectrum(gPrediction); | |
964 | ||
965 | // out->Close(); | |
966 | ||
967 | } |