]> 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
Coverity fixes
[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
006faa18 41enum centrality{ kpp7, kpp276, k07half, kpPb0100, k010, k1020, k020, k2040, k2030, k3040, k4050, k3050, k5060, k4060, k6080, k4080, k5080, k80100, kpPb020, kpPb2040, kpPb4060, kpPb60100 };
42enum centestimator{ kV0M, kV0A, kZNA };
b188dc47 43enum BFDSubtrMethod { knone, kfc, kNb };
7a385c9e 44enum RaavsEP {kPhiIntegrated, kInPlane, kOutOfPlane};
89d2007e 45enum rapidity{ kdefault, k08to04, k07to04, k04to01, k01to01, k01to04, k04to07, k04to08 };
b188dc47 46
47void 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",
033a51ce 51 Int_t fdMethod=kNb, Double_t nevents=1.0, Double_t sigma=1.0, // sigma[pb]
006faa18 52 Bool_t isParticlePlusAntiParticleYield=true, Int_t cc=kpp7, Bool_t PbPbEloss=false,
53 Int_t ccestimator = kV0M,
bdf9f077 54 Int_t isRaavsEP=kPhiIntegrated,const char *epResolfile="",
55 Int_t rapiditySlice=kdefault) {
b188dc47 56
57
9f8b878f 58 gROOT->Macro("$ALICE_ROOT/PWGHF/vertexingHF/macros/LoadLibraries.C");
b188dc47 59
60 // Set if calculation considers asymmetric uncertainties or not
61 Bool_t asym = true;
62
53c91ed7 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)
f4024f34 65 Bool_t isD0Kpi = true;
b188dc47 66 Bool_t isDplusKpipi = false;
67 Bool_t isDstarD0pi = false;
f4024f34 68 Bool_t isDsKKpi = false;
4a277814 69 Bool_t isLctopKpi = false;
53c91ed7 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) {
b188dc47 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;
c7d86f5e 88 else if (fdMethod==knone) { option=0; asym=false; }
b188dc47 89 else option=3;
90
91 if (option>2) {
92 cout<< "Bad calculation option, should be <=2"<<endl;
93 return;
94 }
b188dc47 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.;
006faa18 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 }
b188dc47 134 }
5008ad9f 135
136 // pPb Glauber (A. Toia)
137 // https://twiki.cern.ch/twiki/bin/viewauth/ALICE/PACentStudies#Glauber_Calculations_with_sigma
006faa18 138 if( cc == kpPb0100 ){
5008ad9f 139 tab = 0.098334; tabUnc = 0.0070679;
006faa18 140 // A=207.2; B=1.;
5008ad9f 141 }
006faa18 142 else if( ccestimator == kV0A ){
143 if ( cc == kpPb020 ) {
4cb94a61 144 tab = 0.183; tabUnc = 0.006245;
006faa18 145 } else if ( cc == kpPb2040 ) {
4cb94a61 146 tab = 0.134; tabUnc = 0.004899;
006faa18 147 } else if ( cc == kpPb4060 ) {
4cb94a61 148 tab = 0.092; tabUnc = 0.004796;
006faa18 149 } else if ( cc == kpPb60100 ) {
4cb94a61 150 tab = 0.041; tabUnc = 0.008832;
006faa18 151 }
152 }
153 else if( ccestimator == kZNA ){
154 if ( cc == kpPb020 ) {
4cb94a61 155 tab = 0.164; tabUnc = 0.010724;
006faa18 156 } else if ( cc == kpPb2040 ) {
4cb94a61 157 tab = 0.137; tabUnc = 0.005099;
006faa18 158 } else if ( cc == kpPb4060 ) {
4cb94a61 159 tab = 0.1011; tabUnc = 0.006;
006faa18 160 } else if ( cc == kpPb60100 ) {
4cb94a61 161 tab = 0.0459; tabUnc = 0.003162;
006faa18 162 }
163 }
164
b188dc47 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
b188dc47 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 }
9e97ebdc 219 else if (isDsKKpi){
220 decay = 4;
1c996279 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");
9e97ebdc 227 }
4a277814 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 }
53c91ed7 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 }
b188dc47 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");
bdf9f077 253 //
254 // Scale FONLL inputs if we do the analysis in y-slices
255 //
256 if(rapiditySlice!=kdefault){
257 Double_t scaleFONLL = 1.0;
89d2007e 258 switch(rapiditySlice) {
934916bd 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;
89d2007e 266 }
bdf9f077 267 hDirectMCpt->Scale(scaleFONLL);
bdf9f077 268 hDirectMCptMax->Scale(scaleFONLL);
269 hDirectMCptMin->Scale(scaleFONLL);
89d2007e 270 switch(rapiditySlice) {
934916bd 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;
89d2007e 278 }
279 hFeedDownMCpt->Scale(scaleFONLL);
bdf9f077 280 hFeedDownMCptMax->Scale(scaleFONLL);
281 hFeedDownMCptMin->Scale(scaleFONLL);
282 }
283
b188dc47 284 //
285 //
286 TFile * efffile = new TFile(efffilename,"read");
d87868ee 287 hDirectEffpt = (TH1D*)efffile->Get("hEffD");
b188dc47 288 hDirectEffpt->SetNameTitle("hDirectEffpt","direct acc x eff");
d87868ee 289 hFeedDownEffpt = (TH1D*)efffile->Get("hEffB");
b188dc47 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
7a385c9e 297 //
298 // Read the file of the EP resolution correction
f4024f34 299 TFile *EPf=0;
300 TH1D *hEPresolCorr=0;
7a385c9e 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
b188dc47 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 //
b188dc47 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.);
7a385c9e 388 if(isRaavsEP>0.) spectra->SetIsEventPlaneAnalysis(kTRUE);
b188dc47 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
53c91ed7 412
b188dc47 413 Bool_t combineFeedDown = true;
414 AliHFSystErr *systematics = new AliHFSystErr();
415 if( cc==kpp276 ) {
416 systematics->SetIsLowEnergy(true);
d87868ee 417 }
e6eb9fad 418 else if ( cc == kpPb0100 || cc == kpPb020 || cc == kpPb2040 || cc == kpPb4060 || cc == kpPb60100 ) {
7feda2d4 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 }
d87868ee 436 }
437 //
438 else if( cc!=kpp7 ) {
b188dc47 439 systematics->SetCollisionType(1);
ddd86f95 440 if ( cc == k07half ) systematics->SetCentrality("07half");
a52e34c0 441 else if ( cc == k010 ) systematics->SetCentrality("010");
b188dc47 442 else if ( cc == k1020 ) systematics->SetCentrality("1020");
9f8b878f 443 else if ( cc == k020 ) systematics->SetCentrality("020");
6ae63a20 444 else if ( cc == k2040 || cc == k2030 || cc == k3040 ) {
b188dc47 445 systematics->SetCentrality("2040");
446 systematics->SetIsPbPb2010EnergyScan(true);
447 }
a52e34c0 448 else if ( cc == k3050 ) {
7a385c9e 449 if (isRaavsEP == kPhiIntegrated) systematics->SetCentrality("4080");
450 else if (isRaavsEP == kInPlane) systematics->SetCentrality("3050InPlane");
451 else if (isRaavsEP == kOutOfPlane) systematics->SetCentrality("3050OutOfPlane");
a52e34c0 452 }
6ae63a20 453 else if ( cc == k4060 || cc == k4050 || cc == k5060 ) systematics->SetCentrality("4060");
b188dc47 454 else if ( cc == k6080 ) systematics->SetCentrality("6080");
455 else if ( cc == k4080 ) systematics->SetCentrality("4080");
a52e34c0 456 else {
b188dc47 457 cout << " Systematics not yet implemented " << endl;
458 return;
459 }
460 } else { systematics->SetCollisionType(0); }
5008ad9f 461 //
66d514e1 462 systematics->Init(decay);
b188dc47 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
c7d86f5e 495 gSigmaCorr = spectra->GetCrossSectionFromYieldSpectrum();
496 gYieldCorr = spectra->GetFeedDownCorrectedSpectrum();
b188dc47 497 if (asym) {
b188dc47 498 gSigmaCorrExtreme = spectra->GetCrossSectionFromYieldSpectrumExtreme();
c7d86f5e 499 gYieldCorrExtreme = spectra->GetFeedDownCorrectedSpectrumExtreme();
b188dc47 500 gSigmaCorrConservative = spectra->GetCrossSectionFromYieldSpectrumConservative();
c7d86f5e 501 gYieldCorrConservative = spectra->GetFeedDownCorrectedSpectrumConservative();
b188dc47 502 }
503
504 // Get & Rename the TGraphs
c7d86f5e 505 if (option==0){
b188dc47 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){
d87868ee 783 AliHFPtSpectrum *CalcBins=NULL;
b188dc47 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
c7d86f5e 934 gYieldCorr->Write();
935 gSigmaCorr->Write();
b188dc47 936 if(asym){
b188dc47 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();
c7d86f5e 946 histofcMax->Write(); histofcMin->Write();
b188dc47 947 if(asym && gFcExtreme) gFcExtreme->Write();
948 }
949
950
951 TH1D * hStatUncEffcSigma = spectra->GetDirectStatEffUncOnSigma();
952 TH1D * hStatUncEffbSigma = spectra->GetFeedDownStatEffUncOnSigma();
c7d86f5e 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 }
b188dc47 961
962 // Draw the cross-section
963 // spectra->DrawSpectrum(gPrediction);
964
965 // out->Close();
966
967}