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