1 #if !defined(__CINT__) || defined(__MAKECINT__)
8 #include "TGraphAsymmErrors.h"
13 #include "AliHFSystErr.h"
14 #include "AliHFPtSpectrum.h"
20 // Macro to use the AliHFPtSpectrum class
21 // to compute the feed-down corrections for heavy-flavor
23 // Z.Conesa, September 2010 (zconesa@in2p3.fr)
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
41 enum centrality{ kpp7, kpp276, k010, k1020, k020, k2040, k4060, k6080, k4080, k80100 };
42 enum BFDSubtrMethod { knone, kfc, kNb };
44 void HFPtSpectrum ( const char *mcfilename="FeedDownCorrectionMC.root",
45 const char *efffilename="Efficiencies.root",
46 const char *recofilename="Reconstructed.root", const char *recohistoname="hRawSpectrumD0",
47 const char *outfilename="HFPtSpectrum.root",
48 Int_t fdMethod=kNb, Double_t nevents=1.0, Double_t sigma=1.0, // sigma[nb]
49 Bool_t isParticlePlusAntiParticleYield=true, Int_t cc=kpp7, Bool_t PbPbEloss=false ) {
52 gROOT->Macro("$ALICE_ROOT/PWGHF/vertexingHF/macros/LoadLibraries.C");
54 // Set if calculation considers asymmetric uncertainties or not
57 // Set the meson and decay
58 // (only D0 -> K pi, D+--> K pi pi & D* --> D0 pi & D+s -->KKpi implemented here)
59 Bool_t isD0Kpi = true;
60 Bool_t isDplusKpipi = false;
61 Bool_t isDstarD0pi = false;
62 Bool_t isDsKKpi = false;
63 if (isD0Kpi && isDplusKpipi && isDstarD0pi && isDsKKpi) {
64 cout << "Sorry, can not deal with more than one correction at the same time"<<endl;
69 if (fdMethod==kfc) option=1;
70 else if (fdMethod==kNb) option=2;
71 else if (fdMethod==knone) option=0;
75 cout<< "Bad calculation option, should be <=2"<<endl;
78 if (option==0) asym = false;
82 // Defining the Tab values for the given centrality class
83 // https://twiki.cern.ch/twiki/bin/viewauth/ALICE/CentStudies
85 Double_t tab = 1., tabUnc = 0.;
87 tab = 23.48; tabUnc = 0.97;
88 } else if ( cc == k1020 ) {
89 tab = 14.4318; tabUnc = 0.5733;
90 } else if ( cc == k020 ) {
91 tab = 18.93; tabUnc = 0.74;
92 } else if ( cc == k2040 ) {
93 tab = 6.86; tabUnc = 0.28;
94 } else if ( cc == k4060 ) {
95 tab = 2.00; tabUnc= 0.11;
96 } else if ( cc == k4080 ) {
97 tab = 1.20451; tabUnc = 0.071843;
98 } else if ( cc == k6080 ) {
99 tab = 0.419; tabUnc = 0.033;
100 } else if ( cc == k80100 ){
101 tab = 0.0690; tabUnc = 0.0062;
103 tab *= 1e-9; // to pass from mb^{-1} to pb^{-1}
109 // Get the histograms from the files
111 TH1D *hDirectMCpt=0; // Input MC c-->D spectra
112 TH1D *hFeedDownMCpt=0; // Input MC b-->D spectra
113 TH1D *hDirectMCptMax=0; // Input MC maximum c-->D spectra
114 TH1D *hDirectMCptMin=0; // Input MC minimum c-->D spectra
115 TH1D *hFeedDownMCptMax=0; // Input MC maximum b-->D spectra
116 TH1D *hFeedDownMCptMin=0; // Input MC minimum b-->D spectra
117 // TGraphAsymmErrors *gPrediction=0; // Input MC c-->D spectra
118 TH1D *hDirectEffpt=0; // c-->D Acceptance and efficiency correction
119 TH1D *hFeedDownEffpt=0; // b-->D Acceptance and efficiency correction
120 TH1D *hRECpt=0; // all reconstructed D
123 // Define/Get the input histograms
126 TFile * mcfile = new TFile(mcfilename,"read");
129 hDirectMCpt = (TH1D*)mcfile->Get("hD0Kpipred_central");
130 hFeedDownMCpt = (TH1D*)mcfile->Get("hD0KpifromBpred_central_corr");
131 hDirectMCptMax = (TH1D*)mcfile->Get("hD0Kpipred_max");
132 hDirectMCptMin = (TH1D*)mcfile->Get("hD0Kpipred_min");
133 hFeedDownMCptMax = (TH1D*)mcfile->Get("hD0KpifromBpred_max_corr");
134 hFeedDownMCptMin = (TH1D*)mcfile->Get("hD0KpifromBpred_min_corr");
135 // gPrediction = (TGraphAsymmErrors*)mcfile->Get("D0Kpiprediction");
137 else if (isDplusKpipi){
139 hDirectMCpt = (TH1D*)mcfile->Get("hDpluskpipipred_central");
140 hFeedDownMCpt = (TH1D*)mcfile->Get("hDpluskpipifromBpred_central_corr");
141 hDirectMCptMax = (TH1D*)mcfile->Get("hDpluskpipipred_max");
142 hDirectMCptMin = (TH1D*)mcfile->Get("hDpluskpipipred_min");
143 hFeedDownMCptMax = (TH1D*)mcfile->Get("hDpluskpipifromBpred_max_corr");
144 hFeedDownMCptMin = (TH1D*)mcfile->Get("hDpluskpipifromBpred_min_corr");
145 // gPrediction = (TGraphAsymmErrors*)mcfile->Get("Dpluskpipiprediction");
147 else if(isDstarD0pi){
149 hDirectMCpt = (TH1D*)mcfile->Get("hDstarD0pipred_central");
150 hFeedDownMCpt = (TH1D*)mcfile->Get("hDstarD0pifromBpred_central_corr");
151 hDirectMCptMax = (TH1D*)mcfile->Get("hDstarD0pipred_max");
152 hDirectMCptMin = (TH1D*)mcfile->Get("hDstarD0pipred_min");
153 hFeedDownMCptMax = (TH1D*)mcfile->Get("hDstarD0pifromBpred_max_corr");
154 hFeedDownMCptMin = (TH1D*)mcfile->Get("hDstarD0pifromBpred_min_corr");
155 // gPrediction = (TGraphAsymmErrors*)mcfile->Get("DstarD0piprediction");
159 hDirectMCpt = (TH1D*)mcfile->Get("hDsKkpipred_central");
160 hFeedDownMCpt = (TH1D*)mcfile->Get("hDsKkpifromBpred_central");
161 hDirectMCptMax = (TH1D*)mcfile->Get("hDsKkpipred_max");
162 hDirectMCptMin = (TH1D*)mcfile->Get("hDsKkpipred_min");
163 hFeedDownMCptMax = (TH1D*)mcfile->Get("hDsKkpifromBpred_max");
164 hFeedDownMCptMin = (TH1D*)mcfile->Get("hDsKkpifromBpred_min");
167 hDirectMCpt->SetNameTitle("hDirectMCpt","direct MC spectra");
168 hFeedDownMCpt->SetNameTitle("hFeedDownMCpt","feed-down MC spectra");
169 hDirectMCptMax->SetNameTitle("hDirectMCptMax","max direct MC spectra");
170 hDirectMCptMin->SetNameTitle("hDirectMCptMin","min direct MC spectra");
171 hFeedDownMCptMax->SetNameTitle("hFeedDownMCptMax","max feed-down MC spectra");
172 hFeedDownMCptMin->SetNameTitle("hFeedDownMCptMin","min feed-down MC spectra");
175 TFile * efffile = new TFile(efffilename,"read");
176 hDirectEffpt = (TH1D*)efffile->Get("hDirectEffpt");
177 hDirectEffpt->SetNameTitle("hDirectEffpt","direct acc x eff");
178 hFeedDownEffpt = (TH1D*)efffile->Get("hFeedDownEffpt");
179 hFeedDownEffpt->SetNameTitle("hFeedDownEffpt","feed-down acc x eff");
182 TFile * recofile = new TFile(recofilename,"read");
183 hRECpt = (TH1D*)recofile->Get(recohistoname);
184 hRECpt->SetNameTitle("hRECpt","Reconstructed spectra");
187 // Define the output histograms
189 TFile *out = new TFile(outfilename,"recreate");
194 TH1D *histoYieldCorr=0;
195 TH1D *histoYieldCorrMax=0;
196 TH1D *histoYieldCorrMin=0;
197 TH1D *histoSigmaCorr=0;
198 TH1D *histoSigmaCorrMax=0;
199 TH1D *histoSigmaCorrMin=0;
202 TH1D *histofcRcb_px=0;
203 TH2D *histoYieldCorrRcb=0;
204 TH2D *histoSigmaCorrRcb=0;
206 TGraphAsymmErrors * gYieldCorr = 0;
207 TGraphAsymmErrors * gSigmaCorr = 0;
208 TGraphAsymmErrors * gFcExtreme = 0;
209 TGraphAsymmErrors * gFcConservative = 0;
210 TGraphAsymmErrors * gYieldCorrExtreme = 0;
211 TGraphAsymmErrors * gSigmaCorrExtreme = 0;
212 TGraphAsymmErrors * gYieldCorrConservative = 0;
213 TGraphAsymmErrors * gSigmaCorrConservative = 0;
215 TNtuple * nSigma = 0;
219 // Main functionalities for the calculation
222 // Define and set the basic option flags
223 AliHFPtSpectrum * spectra = new AliHFPtSpectrum("AliHFPtSpectrum","AliHFPtSpectrum",option);
224 spectra->SetFeedDownCalculationOption(option);
225 spectra->SetComputeAsymmetricUncertainties(asym);
226 // Set flag on whether to additional PbPb Eloss hypothesis have to be computed
227 spectra->SetComputeElossHypothesis(PbPbEloss);
229 // Feed the input histograms
230 // reconstructed spectra
231 cout << " Setting the reconstructed spectrum,";
232 spectra->SetReconstructedSpectrum(hRECpt);
233 // acceptance and efficiency corrections
234 cout << " the efficiency,";
235 spectra->SetAccEffCorrection(hDirectEffpt,hFeedDownEffpt);
236 // spectra->SetfIsStatUncEff(false);
237 // option specific histos (theory)
238 cout << " the theoretical spectra";
240 spectra->SetMCptSpectra(hDirectMCpt,hFeedDownMCpt);
241 if(asym) spectra->SetMCptDistributionsBounds(hDirectMCptMax,hDirectMCptMin,hFeedDownMCptMax,hFeedDownMCptMin);
244 spectra->SetFeedDownMCptSpectra(hFeedDownMCpt);
245 if(asym) spectra->SetFeedDownMCptDistributionsBounds(hFeedDownMCptMax,hFeedDownMCptMin);
248 cout << " and the normalization" <<endl;
249 // Set normalization factors (uncertainties set to 0. as example)
250 spectra->SetNormalization(nevents,sigma);
251 Double_t lumi = nevents / sigma ;
252 Double_t lumiUnc = 0.04*lumi; // 10% uncertainty on the luminosity
253 spectra->SetLuminosity(lumi,lumiUnc);
254 Double_t effTrig = 1.0;
255 spectra->SetTriggerEfficiency(effTrig,0.);
257 // Set the global uncertainties on the efficiencies (in percent)
258 Double_t globalEffUnc = 0.15;
259 Double_t globalBCEffRatioUnc = 0.15;
260 spectra->SetAccEffPercentageUncertainty(globalEffUnc,globalBCEffRatioUnc);
262 // Set if the yield is for particle+anti-particle or only one type
263 spectra->SetIsParticlePlusAntiParticleYield(isParticlePlusAntiParticleYield);
265 // Set the Tab parameter and uncertainties
266 if ( (cc != kpp7) && (cc != kpp276) ) {
267 spectra->SetTabParameter(tab,tabUnc);
270 // Do the calculations
271 cout << " Doing the calculation... "<< endl;
272 Double_t deltaY = 1.0;
273 Double_t branchingRatioC = 1.0;
274 Double_t branchingRatioBintoFinalDecay = 1.0; // this is relative to the input theoretical prediction
275 spectra->ComputeHFPtSpectrum(deltaY,branchingRatioC,branchingRatioBintoFinalDecay);
276 cout << " ended the calculation, getting the histograms back " << endl;
278 // Set the systematics externally
279 Bool_t combineFeedDown = true;
280 AliHFSystErr *systematics = new AliHFSystErr();
282 systematics->SetIsLowEnergy(true);
283 } else if( cc!=kpp7 ) {
284 systematics->SetCollisionType(1);
285 if ( cc == k010 ) systematics->SetCentrality("010");
286 else if ( cc == k1020 ) systematics->SetCentrality("1020");
287 else if ( cc == k020 ) systematics->SetCentrality("020");
288 else if ( cc == k2040 ) {
289 systematics->SetCentrality("2040");
290 systematics->SetIsPbPb2010EnergyScan(true);
292 else if ( cc == k4060 ) systematics->SetCentrality("4060");
293 else if ( cc == k6080 ) systematics->SetCentrality("6080");
294 else if ( cc == k4080 ) systematics->SetCentrality("4080");
296 cout << " Systematics not yet implemented " << endl;
299 } else { systematics->SetCollisionType(0); }
300 systematics->Init(decay);
301 spectra->ComputeSystUncertainties(systematics,combineFeedDown);
304 // Get the output histograms
306 // the corrected yield and cross-section
307 histoYieldCorr = (TH1D*)spectra->GetHistoFeedDownCorrectedSpectrum();
308 histoSigmaCorr = (TH1D*)spectra->GetHistoCrossSectionFromYieldSpectrum();
309 histoYieldCorrMax = (TH1D*)spectra->GetHistoUpperLimitFeedDownCorrectedSpectrum();
310 histoYieldCorrMin = (TH1D*)spectra->GetHistoLowerLimitFeedDownCorrectedSpectrum();
311 histoSigmaCorrMax = (TH1D*)spectra->GetHistoUpperLimitCrossSectionFromYieldSpectrum();
312 histoSigmaCorrMin = (TH1D*)spectra->GetHistoLowerLimitCrossSectionFromYieldSpectrum();
313 histoYieldCorr->SetNameTitle("histoYieldCorr","corrected yield");
314 histoYieldCorrMax->SetNameTitle("histoYieldCorrMax","max corrected yield");
315 histoYieldCorrMin->SetNameTitle("histoYieldCorrMin","min corrected yield");
316 histoSigmaCorr->SetNameTitle("histoSigmaCorr","corrected invariant cross-section");
317 histoSigmaCorrMax->SetNameTitle("histoSigmaCorrMax","max corrected invariant cross-section");
318 histoSigmaCorrMin->SetNameTitle("histoSigmaCorrMin","min corrected invariant cross-section");
320 if(!hDirectEffpt) hDirectEffpt = (TH1D*)spectra->GetDirectAccEffCorrection();
321 if(!hFeedDownEffpt) hFeedDownEffpt = (TH1D*)spectra->GetFeedDownAccEffCorrection();
322 // Get the PbPb Eloss hypothesis histograms
324 histofcRcb = spectra->GetHistoFeedDownCorrectionFcVsEloss();
325 histoYieldCorrRcb = spectra->GetHistoFeedDownCorrectedSpectrumVsEloss();
326 histoSigmaCorrRcb = spectra->GetHistoCrossSectionFromYieldSpectrumVsEloss();
327 histofcRcb->SetName("histofcRcb");
328 histoYieldCorrRcb->SetName("histoYieldCorrRcb");
329 histoSigmaCorrRcb->SetName("histoSigmaCorrRcb");
332 // Get & Rename the TGraphs
334 gSigmaCorr = spectra->GetCrossSectionFromYieldSpectrum();
335 gYieldCorr = spectra->GetFeedDownCorrectedSpectrum();
336 gSigmaCorrExtreme = spectra->GetCrossSectionFromYieldSpectrumExtreme();
337 gYieldCorrExtreme = spectra->GetFeedDownCorrectedSpectrumExtreme();
338 gSigmaCorrConservative = spectra->GetCrossSectionFromYieldSpectrumConservative();
339 gYieldCorrConservative = spectra->GetFeedDownCorrectedSpectrumConservative();
342 // Get & Rename the TGraphs
343 if (option==0 && asym){
344 gYieldCorr->SetNameTitle("gYieldCorr","gYieldCorr (uncorr)");
345 gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (uncorr)");
349 histofc = (TH1D*)spectra->GetHistoFeedDownCorrectionFc();
350 histofcMax = (TH1D*)spectra->GetHistoUpperLimitFeedDownCorrectionFc();
351 histofcMin = (TH1D*)spectra->GetHistoLowerLimitFeedDownCorrectionFc();
352 histofc->SetNameTitle("histofc","fc correction factor");
353 histofcMax->SetNameTitle("histofcMax","max fc correction factor");
354 histofcMin->SetNameTitle("histofcMin","min fc correction factor");
356 gYieldCorr->SetNameTitle("gYieldCorr","gYieldCorr (by fc)");
357 gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (by fc)");
358 gFcExtreme = spectra->GetFeedDownCorrectionFcExtreme();
359 gFcExtreme->SetNameTitle("gFcExtreme","gFcExtreme");
360 gYieldCorrExtreme->SetNameTitle("gYieldCorrExtreme","Extreme gYieldCorr (by fc)");
361 gSigmaCorrExtreme->SetNameTitle("gSigmaCorrExtreme","Extreme gSigmaCorr (by fc)");
362 gFcConservative = spectra->GetFeedDownCorrectionFcConservative();
363 gFcConservative->SetNameTitle("gFcConservative","gFcConservative");
364 gYieldCorrConservative->SetNameTitle("gYieldCorrConservative","Conservative gYieldCorr (by fc)");
365 gSigmaCorrConservative->SetNameTitle("gSigmaCorrConservative","Conservative gSigmaCorr (by fc)");
368 if (option==2 && asym) {
369 gYieldCorr->SetNameTitle("gYieldCorr","gYieldCorr (by Nb)");
370 gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (by Nb)");
371 gYieldCorrExtreme->SetNameTitle("gYieldCorrExtreme","Extreme gYieldCorr (by Nb)");
372 gSigmaCorrExtreme->SetNameTitle("gSigmaCorrExtreme","Extreme gSigmaCorr (by Nb)");
373 gYieldCorrConservative->SetNameTitle("gYieldCorrConservative","Conservative gYieldCorr (by Nb)");
374 gSigmaCorrConservative->SetNameTitle("gSigmaCorrConservative","Conservative gSigmaCorr (by Nb)");
375 gFcConservative = spectra->GetFeedDownCorrectionFcConservative();
376 gFcConservative->SetNameTitle("gFcConservative","gFcConservative");
380 nSigma = spectra->GetNtupleCrossSectionVsEloss();
384 // Now, plot the results ! :)
387 gROOT->SetStyle("Plain");
389 cout << " Drawing the results ! " << endl;
394 TCanvas *ceff = new TCanvas("ceff","efficiency drawing");
397 hDirectEffpt->Draw();
399 hFeedDownEffpt->Draw();
402 TCanvas *cTheoryRebin = new TCanvas("cTheoryRebin","control the theoretical spectra rebin");
403 cTheoryRebin->Divide(1,2);
405 hDirectMCpt->Draw("");
406 TH1D *hDirectMCptRebin = (TH1D*)spectra->GetDirectTheoreticalSpectrum();
407 hDirectMCptRebin->SetLineColor(2);
408 hDirectMCptRebin->Draw("same");
410 hFeedDownMCpt->Draw("");
411 TH1D *hFeedDownRebin = (TH1D*)spectra->GetFeedDownTheoreticalSpectrum();
412 hFeedDownRebin->SetLineColor(2);
413 hFeedDownRebin->Draw("same");
414 cTheoryRebin->Update();
416 TCanvas *cTheoryRebinLimits = new TCanvas("cTheoryRebinLimits","control the theoretical spectra limits rebin");
417 cTheoryRebinLimits->Divide(1,2);
418 cTheoryRebinLimits->cd(1);
419 hDirectMCptMax->Draw("");
420 TH1D *hDirectMCptMaxRebin = (TH1D*)spectra->GetDirectTheoreticalUpperLimitSpectrum();
421 hDirectMCptMaxRebin->SetLineColor(2);
422 hDirectMCptMaxRebin->Draw("same");
423 hDirectMCptMin->Draw("same");
424 TH1D *hDirectMCptMinRebin = (TH1D*)spectra->GetDirectTheoreticalLowerLimitSpectrum();
425 hDirectMCptMinRebin->SetLineColor(2);
426 hDirectMCptMinRebin->Draw("same");
427 cTheoryRebinLimits->cd(2);
428 hFeedDownMCptMax->Draw("");
429 TH1D *hFeedDownMaxRebin = (TH1D*)spectra->GetFeedDownTheoreticalUpperLimitSpectrum();
430 hFeedDownMaxRebin->SetLineColor(2);
431 hFeedDownMaxRebin->Draw("same");
432 hFeedDownMCptMin->Draw("same");
433 TH1D *hFeedDownMinRebin = (TH1D*)spectra->GetFeedDownTheoreticalLowerLimitSpectrum();
434 hFeedDownMinRebin->SetLineColor(2);
435 hFeedDownMinRebin->Draw("same");
436 cTheoryRebinLimits->Update();
441 TCanvas * cfc = new TCanvas("cfc","Fc");
442 histofcMax->Draw("c");
443 histofc->Draw("csame");
444 histofcMin->Draw("csame");
448 TH2F *histofcDraw= new TH2F("histofcDraw","histofc (for drawing)",100,0,33.25,100,0.01,1.25);
449 histofcDraw->SetStats(0);
450 histofcDraw->GetXaxis()->SetTitle("p_{T} [GeV]");
451 histofcDraw ->GetXaxis()->SetTitleSize(0.05);
452 histofcDraw->GetXaxis()->SetTitleOffset(0.95);
453 histofcDraw->GetYaxis()->SetTitle(" fc ");
454 histofcDraw->GetYaxis()->SetTitleSize(0.05);
458 // for(Int_t item=0; item<gSigmaCorr->GetN(); item++){
459 // Double_t center=0., value=0.;
460 // gFcExtreme->GetPoint(item,center,value);
461 // Double_t highunc = gFcExtreme->GetErrorYhigh(item) / value ;
462 // Double_t lowunc = gFcExtreme->GetErrorYlow(item) / value ;
463 // cout << "Fc extreme: i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
465 // for(Int_t item=0; item<gSigmaCorr->GetN(); item++){
466 // Double_t center=0., value=0.;
467 // gFcConservative->GetPoint(item,center,value);
468 // Double_t highunc = gFcConservative->GetErrorYhigh(item) / value ;
469 // Double_t lowunc = gFcConservative->GetErrorYlow(item) / value ;
470 // cout << "Fc conservative: i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
472 TCanvas *cfcExtreme = new TCanvas("cfcExtreme","Extreme Asymmetric fc (TGraphAsymmErr)");
473 gFcExtreme->SetFillStyle(3006);
474 gFcExtreme->SetLineWidth(3);
475 gFcExtreme->SetMarkerStyle(20);
476 gFcExtreme->SetFillColor(2);
478 gFcExtreme->Draw("3same");
481 gFcConservative->SetFillStyle(3007);
482 gFcConservative->SetFillColor(4);
483 gFcConservative->Draw("3same");
486 cfcExtreme->Update();
493 // Drawing the results (the raw-reconstructed, the expected, and the corrected spectra)
495 TCanvas * cresult = new TCanvas("cresult","corrected yields & sigma");
496 hDirectMCpt->SetMarkerStyle(20);
497 hDirectMCpt->SetMarkerColor(4);
498 hDirectMCpt->Draw("p");
499 histoSigmaCorr->SetMarkerStyle(21);
500 histoSigmaCorr->SetMarkerColor(2);
501 histoSigmaCorr->Draw("psame");
502 histoYieldCorr->SetMarkerStyle(22);
503 histoYieldCorr->SetMarkerColor(6);
504 histoYieldCorr->Draw("psame");
505 hRECpt->SetMarkerStyle(23);
506 hRECpt->SetMarkerColor(3);
507 hRECpt->Draw("psame");
511 TCanvas * cresult2 = new TCanvas("cresult2","corrected yield & sigma");
512 histoSigmaCorr->SetMarkerStyle(21);
513 histoSigmaCorr->SetMarkerColor(2);
514 histoSigmaCorr->Draw("p");
515 histoYieldCorr->SetMarkerStyle(22);
516 histoYieldCorr->SetMarkerColor(6);
517 histoYieldCorr->Draw("psame");
518 hRECpt->SetMarkerStyle(23);
519 hRECpt->SetMarkerColor(3);
520 hRECpt->Draw("psame");
527 TH2F *histoDraw = new TH2F("histoDraw","histo (for drawing)",100,0,33.25,100,50.,1e7);
528 float max = 1.1*gYieldCorr->GetMaximum();
529 histoDraw->SetAxisRange(0.1,max,"Y");
530 histoDraw->SetStats(0);
531 histoDraw->GetXaxis()->SetTitle("p_{T} [GeV]");
532 histoDraw->GetXaxis()->SetTitleSize(0.05);
533 histoDraw->GetXaxis()->SetTitleOffset(0.95);
534 histoDraw->GetYaxis()->SetTitle("#frac{d#N}{dp_{T}} |_{|y|<1} [L & trigger uncorr]");
535 histoDraw->GetYaxis()->SetTitleSize(0.05);
536 TCanvas * cyieldAsym = new TCanvas("cyieldAsym","Asymmetric corrected yield (TGraphAsymmErr)");
537 gYieldCorr->SetFillStyle(3001);
538 gYieldCorr->SetLineWidth(3);
539 gYieldCorr->SetMarkerStyle(20);
540 gYieldCorr->SetFillColor(3);
542 gYieldCorr->Draw("3LPsame");
543 gYieldCorr->Draw("Xsame");
544 cyieldAsym->SetLogy();
545 cyieldAsym->Update();
547 TCanvas * cyieldExtreme = new TCanvas("cyieldExtreme","Extreme Asymmetric corrected yield (TGraphAsymmErr)");
548 histoYieldCorr->Draw();
549 gYieldCorrExtreme->SetFillStyle(3002);
550 gYieldCorrExtreme->SetLineWidth(3);
551 gYieldCorrExtreme->SetMarkerStyle(20);
552 gYieldCorrExtreme->SetFillColor(2);
553 histoYieldCorr->Draw();
554 gYieldCorr->Draw("3same");
555 gYieldCorrExtreme->Draw("3same");
556 cyieldExtreme->SetLogy();
557 cyieldExtreme->Update();
559 TH2F *histo2Draw = new TH2F("histo2Draw","histo2 (for drawing)",100,0,33.25,100,50.,1e9);
560 max = 1.1*gSigmaCorr->GetMaximum();
561 histo2Draw->SetAxisRange(0.1,max,"Y");
562 histo2Draw->SetStats(0);
563 histo2Draw->GetXaxis()->SetTitle("p_{T} [GeV]");
564 histo2Draw->GetXaxis()->SetTitleSize(0.05);
565 histo2Draw->GetXaxis()->SetTitleOffset(0.95);
566 histo2Draw->GetYaxis()->SetTitle("#frac{1}{BR} #times #frac{d#sigma}{dp_{T}} |_{|y|<1}");
567 histo2Draw->GetYaxis()->SetTitleSize(0.05);
568 TCanvas * csigmaAsym = new TCanvas("csigmaAsym","Asymmetric corrected sigma (TGraphAsymmErr)");
569 gSigmaCorr->SetFillStyle(3001);
570 gSigmaCorr->SetLineWidth(3);
571 gSigmaCorr->SetMarkerStyle(21);
572 gSigmaCorr->SetFillColor(3);
574 gSigmaCorr->Draw("3LPsame");
575 gSigmaCorr->Draw("Xsame");
576 csigmaAsym->SetLogy();
577 csigmaAsym->Update();
579 // cout << endl <<" Sytematics (stat approach) " <<endl;
580 // for(Int_t item=0; item<gSigmaCorr->GetN(); item++){
581 // Double_t center=0., value=0.;
582 // gSigmaCorr->GetPoint(item,center,value);
583 // Double_t highunc = gSigmaCorr->GetErrorYhigh(item) / value ;
584 // Double_t lowunc = gSigmaCorr->GetErrorYlow(item) / value ;
585 // cout << "Sigma syst (stat), i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
588 TCanvas * csigmaExtreme = new TCanvas("csigmaExtreme","Asymmetric extreme corrected sigma (TGraphAsymmErr)");
589 histoSigmaCorr->Draw();
590 gSigmaCorr->Draw("3Psame");
591 gSigmaCorrExtreme->SetFillStyle(3002);
592 gSigmaCorrExtreme->SetLineWidth(3);
593 gSigmaCorrExtreme->SetMarkerStyle(21);
594 gSigmaCorrExtreme->SetFillColor(2);
595 gSigmaCorrExtreme->Draw("3Psame");
596 csigmaExtreme->SetLogy();
597 csigmaExtreme->Update();
599 // cout << endl << " Sytematics (Extreme approach)" <<endl;
600 // for(Int_t item=0; item<gSigmaCorrExtreme->GetN(); item++){
601 // Double_t center=0., value=0.;
602 // gSigmaCorrExtreme->GetPoint(item,center,value);
603 // Double_t highunc = gSigmaCorrExtreme->GetErrorYhigh(item) / value ;
604 // Double_t lowunc = gSigmaCorrExtreme->GetErrorYlow(item) / value ;
605 // cout << "Sigma syst (extreme) i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
608 // cout << endl << " Sytematics (Conservative approach)" <<endl;
609 // for(Int_t item=0; item<gSigmaCorrConservative->GetN(); item++){
610 // Double_t center=0., value=0.;
611 // gSigmaCorrConservative->GetPoint(item,center,value);
612 // Double_t highunc = gSigmaCorrConservative->GetErrorYhigh(item) / value ;
613 // Double_t lowunc = gSigmaCorrConservative->GetErrorYlow(item) / value ;
614 // cout << "Sigma syst (conservative) i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
619 // Draw the PbPb Eloss hypothesis histograms
621 AliHFPtSpectrum *CalcBins;
622 gStyle->SetPalette(1);
623 TCanvas *canvasfcRcb = new TCanvas("canvasfcRcb","fc vs pt vs Rcb");
624 // histofcRcb->Draw("cont4z");
625 histofcRcb->Draw("colz");
626 canvasfcRcb->Update();
628 TCanvas *canvasfcRcb1 = new TCanvas("canvasfcRcb1","fc vs pt vs Rcb=1");
629 histofcRcb_px = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px",40,40);
630 histofcRcb_px->SetLineColor(2);
633 histofcRcb_px->Draw("same");
634 } else histofcRcb_px->Draw("");
635 canvasfcRcb1->Update();
636 TCanvas *canvasfcRcb2 = new TCanvas("canvasfcRcb2","fc vs pt vs Rcb fixed Rcb");
637 Int_t bin0 = CalcBins->FindTH2YBin(histofcRcb,0.25);
638 Int_t bin1 = CalcBins->FindTH2YBin(histofcRcb,0.5);
639 Int_t bin2 = CalcBins->FindTH2YBin(histofcRcb,1.0);
640 Int_t bin3 = CalcBins->FindTH2YBin(histofcRcb,1.5);
641 Int_t bin4 = CalcBins->FindTH2YBin(histofcRcb,2.0);
642 Int_t bin5 = CalcBins->FindTH2YBin(histofcRcb,3.0);
643 Int_t bin6 = CalcBins->FindTH2YBin(histofcRcb,4.0);
644 TH1D * histofcRcb_px0a = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px0a",bin0,bin0);
645 TH1D * histofcRcb_px0 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px0",bin1,bin1);
646 TH1D * histofcRcb_px1 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px1",bin2,bin2);
647 TH1D * histofcRcb_px2 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px2",bin3,bin3);
648 TH1D * histofcRcb_px3 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px3",bin4,bin4);
649 TH1D * histofcRcb_px4 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px4",bin5,bin5);
650 TH1D * histofcRcb_px5 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px5",bin6,bin6);
653 // histofcRcb_px->Draw("same");
655 // histofcRcb_px->Draw("");
656 histofcRcb_px0a->SetLineColor(2);
657 histofcRcb_px0a->Draw("");
659 histofcRcb_px0a->SetLineColor(2);
660 histofcRcb_px0a->Draw("same");
661 histofcRcb_px0->SetLineColor(4);
662 histofcRcb_px0->Draw("same");
663 histofcRcb_px1->SetLineColor(3);
664 histofcRcb_px1->Draw("same");
665 histofcRcb_px2->SetLineColor(kCyan);
666 histofcRcb_px2->Draw("same");
667 histofcRcb_px3->SetLineColor(kMagenta+1);
668 histofcRcb_px3->Draw("same");
669 histofcRcb_px4->SetLineColor(kOrange+7);
670 histofcRcb_px4->Draw("same");
671 histofcRcb_px5->SetLineColor(kGreen+3);
672 histofcRcb_px5->Draw("same");
673 TLegend *legrcc = new TLegend(0.8,0.8,0.95,0.9);
674 legrcc->SetFillColor(0);
676 legrcc->AddEntry(histofcRcb_px0a,"Rc/b=0.25","l");
677 legrcc->AddEntry(histofcRcb_px0,"Rc/b=0.5","l");
678 legrcc->AddEntry(histofcRcb_px1,"Rc/b=1.0","l");
679 legrcc->AddEntry(histofcRcb_px2,"Rc/b=1.5","l");
680 legrcc->AddEntry(histofcRcb_px3,"Rc/b=2.0","l");
681 legrcc->AddEntry(histofcRcb_px4,"Rc/b=3.0","l");
682 legrcc->AddEntry(histofcRcb_px5,"Rc/b=4.0","l");
684 legrcc->AddEntry(histofcRcb_px0a,"Rb=0.25","l");
685 legrcc->AddEntry(histofcRcb_px0,"Rb=0.5","l");
686 legrcc->AddEntry(histofcRcb_px1,"Rb=1.0","l");
687 legrcc->AddEntry(histofcRcb_px2,"Rb=1.5","l");
688 legrcc->AddEntry(histofcRcb_px3,"Rb=2.0","l");
689 legrcc->AddEntry(histofcRcb_px4,"Rb=3.0","l");
690 legrcc->AddEntry(histofcRcb_px5,"Rb=4.0","l");
693 canvasfcRcb2->Update();
694 TCanvas *canvasYRcb = new TCanvas("canvasYRcb","corrected yield vs pt vs Rcb");
695 histoYieldCorrRcb->Draw("cont4z");
696 canvasYRcb->Update();
697 TCanvas *canvasSRcb = new TCanvas("canvasSRcb","sigma vs pt vs Rcb");
698 histoSigmaCorrRcb->Draw("cont4z");
699 canvasSRcb->Update();
700 TCanvas *canvasSRcb1 = new TCanvas("canvasSRcb1","sigma vs pt vs Rcb fixed Rcb");
701 TH1D * histoSigmaCorrRcb_px0a = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px0a",bin0,bin0);
702 TH1D * histoSigmaCorrRcb_px0 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px0",bin1,bin1);
703 TH1D * histoSigmaCorrRcb_px1 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px1",bin2,bin2);
704 TH1D * histoSigmaCorrRcb_px2 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px2",bin3,bin3);
705 TH1D * histoSigmaCorrRcb_px3 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px3",bin4,bin4);
706 TH1D * histoSigmaCorrRcb_px4 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px4",bin5,bin5);
707 TH1D * histoSigmaCorrRcb_px5 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px5",bin6,bin6);
708 histoSigmaCorr->Draw();
709 histoSigmaCorrRcb_px0a->SetLineColor(2);
710 histoSigmaCorrRcb_px0a->Draw("hsame");
711 histoSigmaCorrRcb_px0->SetLineColor(4);
712 histoSigmaCorrRcb_px0->Draw("hsame");
713 histoSigmaCorrRcb_px1->SetLineColor(3);
714 histoSigmaCorrRcb_px1->Draw("hsame");
715 histoSigmaCorrRcb_px2->SetLineColor(kCyan);
716 histoSigmaCorrRcb_px2->Draw("hsame");
717 histoSigmaCorrRcb_px3->SetLineColor(kMagenta+1);
718 histoSigmaCorrRcb_px3->Draw("hsame");
719 histoSigmaCorrRcb_px4->SetLineColor(kOrange+7);
720 histoSigmaCorrRcb_px4->Draw("same");
721 histoSigmaCorrRcb_px5->SetLineColor(kGreen+3);
722 histoSigmaCorrRcb_px5->Draw("same");
723 TLegend *legrcb = new TLegend(0.8,0.8,0.95,0.9);
724 legrcb->SetFillColor(0);
726 legrcb->AddEntry(histoSigmaCorrRcb_px0a,"Rc/b=0.25","l");
727 legrcb->AddEntry(histoSigmaCorrRcb_px0,"Rc/b=0.5","l");
728 legrcb->AddEntry(histoSigmaCorrRcb_px1,"Rc/b=1.0","l");
729 legrcb->AddEntry(histoSigmaCorrRcb_px2,"Rc/b=1.5","l");
730 legrcb->AddEntry(histoSigmaCorrRcb_px3,"Rc/b=2.0","l");
731 legrcb->AddEntry(histoSigmaCorrRcb_px4,"Rc/b=3.0","l");
732 legrcb->AddEntry(histoSigmaCorrRcb_px5,"Rc/b=4.0","l");
734 legrcb->AddEntry(histoSigmaCorrRcb_px0a,"Rb=0.25","l");
735 legrcb->AddEntry(histoSigmaCorrRcb_px0,"Rb=0.5","l");
736 legrcb->AddEntry(histoSigmaCorrRcb_px1,"Rb=1.0","l");
737 legrcb->AddEntry(histoSigmaCorrRcb_px2,"Rb=1.5","l");
738 legrcb->AddEntry(histoSigmaCorrRcb_px3,"Rb=2.0","l");
739 legrcb->AddEntry(histoSigmaCorrRcb_px4,"Rb=3.0","l");
740 legrcb->AddEntry(histoSigmaCorrRcb_px5,"Rb=4.0","l");
743 canvasSRcb1->Update();
748 // Write the histograms to the output file
750 cout << " Saving the results ! " << endl<< endl;
754 hDirectMCpt->Write(); hFeedDownMCpt->Write();
755 hDirectMCptMax->Write(); hDirectMCptMin->Write();
756 hFeedDownMCptMax->Write(); hFeedDownMCptMin->Write();
757 if(hDirectEffpt) hDirectEffpt->Write(); if(hFeedDownEffpt) hFeedDownEffpt->Write();
760 histoYieldCorr->Write();
761 histoYieldCorrMax->Write(); histoYieldCorrMin->Write();
762 histoSigmaCorr->Write();
763 histoSigmaCorrMax->Write(); histoSigmaCorrMin->Write();
766 histofcRcb->Write(); histofcRcb_px->Write();
767 histoYieldCorrRcb->Write();
768 histoSigmaCorrRcb->Write();
775 if(gYieldCorrExtreme) gYieldCorrExtreme->Write();
776 if(gSigmaCorrExtreme) gSigmaCorrExtreme->Write();
777 if(gYieldCorrConservative) gYieldCorrConservative->Write();
778 if(gSigmaCorrConservative) gSigmaCorrConservative->Write();
779 if(asym && gFcConservative) gFcConservative->Write();
784 histofcMax->Write(); histofcMin->Write();
785 if(asym && gFcExtreme) gFcExtreme->Write();
789 TH1D * hStatUncEffcSigma = spectra->GetDirectStatEffUncOnSigma();
790 TH1D * hStatUncEffbSigma = spectra->GetFeedDownStatEffUncOnSigma();
791 TH1D * hStatUncEffcFD = spectra->GetDirectStatEffUncOnFc();
792 TH1D * hStatUncEffbFD = spectra->GetFeedDownStatEffUncOnFc();
793 hStatUncEffcSigma->Write();
794 hStatUncEffbSigma->Write();
795 hStatUncEffcFD->Write();
796 hStatUncEffbFD->Write();
798 // Draw the cross-section
799 // spectra->DrawSpectrum(gPrediction);