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