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