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 [pb]
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, k07half, kpPb0100, k010, k1020, k020, k2040, k2030, k3040, k4050, k3050, k5060, k4060, k6080, k4080, k80100 };
42 enum BFDSubtrMethod { knone, kfc, kNb };
43 enum RaavsEP {kPhiIntegrated, kInPlane, kOutOfPlane};
44 enum rapidity{ kdefault, km01to01, k01to04, k04to07, k04to08 };
46 void HFPtSpectrum ( const char *mcfilename="FeedDownCorrectionMC.root",
47 const char *efffilename="Efficiencies.root",
48 const char *recofilename="Reconstructed.root", const char *recohistoname="hRawSpectrumD0",
49 const char *outfilename="HFPtSpectrum.root",
50 Int_t fdMethod=kNb, Double_t nevents=1.0, Double_t sigma=1.0, // sigma[pb]
51 Bool_t isParticlePlusAntiParticleYield=true, Int_t cc=kpp7, Bool_t PbPbEloss=false,
52 Int_t isRaavsEP=kPhiIntegrated,const char *epResolfile="",
53 Int_t rapiditySlice=kdefault) {
56 gROOT->Macro("$ALICE_ROOT/PWGHF/vertexingHF/macros/LoadLibraries.C");
58 // Set if calculation considers asymmetric uncertainties or not
61 // Set the meson and decay
62 // (only D0 -> K pi, D+--> K pi pi & D* --> D0 pi & D+s -->KKpi implemented here)
63 Bool_t isD0Kpi = false;
64 Bool_t isDplusKpipi = false;
65 Bool_t isDstarD0pi = false;
66 Bool_t isDsKKpi = true;
67 Bool_t isLctopKpi = false;
68 if (isD0Kpi && isDplusKpipi && isDstarD0pi && isDsKKpi) {
69 cout << "Sorry, can not deal with more than one correction at the same time"<<endl;
74 if (fdMethod==kfc) option=1;
75 else if (fdMethod==kNb) option=2;
76 else if (fdMethod==knone) { option=0; asym=false; }
80 cout<< "Bad calculation option, should be <=2"<<endl;
86 // Defining the Tab values for the given centrality class
87 // https://twiki.cern.ch/twiki/bin/viewauth/ALICE/CentStudies
89 Double_t tab = 1., tabUnc = 0.;
90 if ( cc == k07half ) {
91 tab = 24.81; tabUnc = 0.8037;
92 } else if ( cc == k010 ) {
93 tab = 23.48; tabUnc = 0.97;
94 } else if ( cc == k1020 ) {
95 tab = 14.4318; tabUnc = 0.5733;
96 } else if ( cc == k020 ) {
97 tab = 18.93; tabUnc = 0.74;
98 } else if ( cc == k2040 ) {
99 tab = 6.86; tabUnc = 0.28;
100 } else if ( cc == k2030 ) {
101 tab = 8.73769; tabUnc = 0.370219;
102 } else if ( cc == k3040 ) {
103 tab = 5.02755; tabUnc = 0.22099;
104 } else if ( cc == k4050 ) {
105 tab = 2.68327; tabUnc = 0.137073;
106 } else if ( cc == k3050 ) {
107 tab = 3.87011; tabUnc = 0.183847;
108 } else if ( cc == k4060 ) {
109 tab = 2.00; tabUnc= 0.11;
110 } else if ( cc == k4080 ) {
111 tab = 1.20451; tabUnc = 0.071843;
112 } else if ( cc == k5060 ) {
113 tab = 1.32884; tabUnc = 0.0929536;
114 } else if ( cc == k6080 ) {
115 tab = 0.419; tabUnc = 0.033;
116 } else if ( cc == k80100 ){
117 tab = 0.0690; tabUnc = 0.0062;
120 // pPb Glauber (A. Toia)
121 // https://twiki.cern.ch/twiki/bin/viewauth/ALICE/PACentStudies#Glauber_Calculations_with_sigma
122 else if( cc == kpPb0100 ){
123 tab = 0.098334; tabUnc = 0.0070679;
125 tab *= 1e-9; // to pass from mb^{-1} to pb^{-1}
131 // Get the histograms from the files
133 TH1D *hDirectMCpt=0; // Input MC c-->D spectra
134 TH1D *hFeedDownMCpt=0; // Input MC b-->D spectra
135 TH1D *hDirectMCptMax=0; // Input MC maximum c-->D spectra
136 TH1D *hDirectMCptMin=0; // Input MC minimum c-->D spectra
137 TH1D *hFeedDownMCptMax=0; // Input MC maximum b-->D spectra
138 TH1D *hFeedDownMCptMin=0; // Input MC minimum b-->D spectra
139 // TGraphAsymmErrors *gPrediction=0; // Input MC c-->D spectra
140 TH1D *hDirectEffpt=0; // c-->D Acceptance and efficiency correction
141 TH1D *hFeedDownEffpt=0; // b-->D Acceptance and efficiency correction
142 TH1D *hRECpt=0; // all reconstructed D
145 // Define/Get the input histograms
148 TFile * mcfile = new TFile(mcfilename,"read");
151 hDirectMCpt = (TH1D*)mcfile->Get("hD0Kpipred_central");
152 hFeedDownMCpt = (TH1D*)mcfile->Get("hD0KpifromBpred_central_corr");
153 hDirectMCptMax = (TH1D*)mcfile->Get("hD0Kpipred_max");
154 hDirectMCptMin = (TH1D*)mcfile->Get("hD0Kpipred_min");
155 hFeedDownMCptMax = (TH1D*)mcfile->Get("hD0KpifromBpred_max_corr");
156 hFeedDownMCptMin = (TH1D*)mcfile->Get("hD0KpifromBpred_min_corr");
157 // gPrediction = (TGraphAsymmErrors*)mcfile->Get("D0Kpiprediction");
159 else if (isDplusKpipi){
161 hDirectMCpt = (TH1D*)mcfile->Get("hDpluskpipipred_central");
162 hFeedDownMCpt = (TH1D*)mcfile->Get("hDpluskpipifromBpred_central_corr");
163 hDirectMCptMax = (TH1D*)mcfile->Get("hDpluskpipipred_max");
164 hDirectMCptMin = (TH1D*)mcfile->Get("hDpluskpipipred_min");
165 hFeedDownMCptMax = (TH1D*)mcfile->Get("hDpluskpipifromBpred_max_corr");
166 hFeedDownMCptMin = (TH1D*)mcfile->Get("hDpluskpipifromBpred_min_corr");
167 // gPrediction = (TGraphAsymmErrors*)mcfile->Get("Dpluskpipiprediction");
169 else if(isDstarD0pi){
171 hDirectMCpt = (TH1D*)mcfile->Get("hDstarD0pipred_central");
172 hFeedDownMCpt = (TH1D*)mcfile->Get("hDstarD0pifromBpred_central_corr");
173 hDirectMCptMax = (TH1D*)mcfile->Get("hDstarD0pipred_max");
174 hDirectMCptMin = (TH1D*)mcfile->Get("hDstarD0pipred_min");
175 hFeedDownMCptMax = (TH1D*)mcfile->Get("hDstarD0pifromBpred_max_corr");
176 hFeedDownMCptMin = (TH1D*)mcfile->Get("hDstarD0pifromBpred_min_corr");
177 // gPrediction = (TGraphAsymmErrors*)mcfile->Get("DstarD0piprediction");
181 hDirectMCpt = (TH1D*)mcfile->Get("hDsPhipitoKkpipred_central");
182 hFeedDownMCpt = (TH1D*)mcfile->Get("hDsPhipitoKkpifromBpred_central_corr");
183 hDirectMCptMax = (TH1D*)mcfile->Get("hDsPhipitoKkpipred_max");
184 hDirectMCptMin = (TH1D*)mcfile->Get("hDsPhipitoKkpipred_min");
185 hFeedDownMCptMax = (TH1D*)mcfile->Get("hDsPhipitoKkpifromBpred_max_corr");
186 hFeedDownMCptMin = (TH1D*)mcfile->Get("hDsPhipitoKkpifromBpred_min_corr");
188 else if (isLctopKpi){
190 hDirectMCpt = (TH1D*)mcfile->Get("hLcpkpipred_central");
191 hFeedDownMCpt = (TH1D*)mcfile->Get("hLcpkpifromBpred_central_corr");
192 hDirectMCptMax = (TH1D*)mcfile->Get("hLcpkpipred_max");
193 hDirectMCptMin = (TH1D*)mcfile->Get("hLcpkpipred_min");
194 hFeedDownMCptMax = (TH1D*)mcfile->Get("hLcpkpifromBpred_max_corr");
195 hFeedDownMCptMin = (TH1D*)mcfile->Get("hLcpkpifromBpred_min_corr");
198 hDirectMCpt->SetNameTitle("hDirectMCpt","direct MC spectra");
199 hFeedDownMCpt->SetNameTitle("hFeedDownMCpt","feed-down MC spectra");
200 hDirectMCptMax->SetNameTitle("hDirectMCptMax","max direct MC spectra");
201 hDirectMCptMin->SetNameTitle("hDirectMCptMin","min direct MC spectra");
202 hFeedDownMCptMax->SetNameTitle("hFeedDownMCptMax","max feed-down MC spectra");
203 hFeedDownMCptMin->SetNameTitle("hFeedDownMCptMin","min feed-down MC spectra");
205 // Scale FONLL inputs if we do the analysis in y-slices
207 if(rapiditySlice!=kdefault){
208 Double_t scaleFONLL = 1.0;
209 if(rapiditySlice==km01to01) scaleFONLL = 0.2/1.0;
210 else if (rapiditySlice==k01to04) scaleFONLL = 0.3/1.0;
211 else if (rapiditySlice==k04to07) scaleFONLL = 0.1/1.0;
212 else if (rapiditySlice==k04to08) scaleFONLL = 0.1/1.0;
213 hDirectMCpt->Scale(scaleFONLL);
214 hFeedDownMCpt->Scale(scaleFONLL);
215 hDirectMCptMax->Scale(scaleFONLL);
216 hDirectMCptMin->Scale(scaleFONLL);
217 hFeedDownMCptMax->Scale(scaleFONLL);
218 hFeedDownMCptMin->Scale(scaleFONLL);
223 TFile * efffile = new TFile(efffilename,"read");
224 hDirectEffpt = (TH1D*)efffile->Get("hEffD");
225 hDirectEffpt->SetNameTitle("hDirectEffpt","direct acc x eff");
226 hFeedDownEffpt = (TH1D*)efffile->Get("hEffB");
227 hFeedDownEffpt->SetNameTitle("hFeedDownEffpt","feed-down acc x eff");
230 TFile * recofile = new TFile(recofilename,"read");
231 hRECpt = (TH1D*)recofile->Get(recohistoname);
232 hRECpt->SetNameTitle("hRECpt","Reconstructed spectra");
235 // Read the file of the EP resolution correction
239 EPf = new TFile(epResolfile,"read");
240 if(isRaavsEP==kInPlane) hEPresolCorr = (TH1D*)EPf->Get("hCorrEPresol_InPlane");
241 else if(isRaavsEP==kOutOfPlane) hEPresolCorr = (TH1D*)EPf->Get("hCorrEPresol_OutOfPlane");
242 for(Int_t i=1; i<=hRECpt->GetNbinsX(); i++) {
243 Double_t value = hRECpt->GetBinContent(i);
244 Double_t error = hRECpt->GetBinError(i);
245 Double_t pt = hRECpt->GetBinCenter(i);
246 Int_t epbin = hEPresolCorr->FindBin( pt );
247 Double_t epcorr = hEPresolCorr->GetBinContent( epbin );
248 value = value*epcorr;
249 error = error*epcorr;
250 hRECpt->SetBinContent(i,value);
251 hRECpt->SetBinError(i,error);
256 // Define the output histograms
258 TFile *out = new TFile(outfilename,"recreate");
263 TH1D *histoYieldCorr=0;
264 TH1D *histoYieldCorrMax=0;
265 TH1D *histoYieldCorrMin=0;
266 TH1D *histoSigmaCorr=0;
267 TH1D *histoSigmaCorrMax=0;
268 TH1D *histoSigmaCorrMin=0;
271 TH1D *histofcRcb_px=0;
272 TH2D *histoYieldCorrRcb=0;
273 TH2D *histoSigmaCorrRcb=0;
275 TGraphAsymmErrors * gYieldCorr = 0;
276 TGraphAsymmErrors * gSigmaCorr = 0;
277 TGraphAsymmErrors * gFcExtreme = 0;
278 TGraphAsymmErrors * gFcConservative = 0;
279 TGraphAsymmErrors * gYieldCorrExtreme = 0;
280 TGraphAsymmErrors * gSigmaCorrExtreme = 0;
281 TGraphAsymmErrors * gYieldCorrConservative = 0;
282 TGraphAsymmErrors * gSigmaCorrConservative = 0;
284 TNtuple * nSigma = 0;
288 // Main functionalities for the calculation
291 // Define and set the basic option flags
292 AliHFPtSpectrum * spectra = new AliHFPtSpectrum("AliHFPtSpectrum","AliHFPtSpectrum",option);
293 spectra->SetFeedDownCalculationOption(option);
294 spectra->SetComputeAsymmetricUncertainties(asym);
295 // Set flag on whether to additional PbPb Eloss hypothesis have to be computed
296 spectra->SetComputeElossHypothesis(PbPbEloss);
298 // Feed the input histograms
299 // reconstructed spectra
300 cout << " Setting the reconstructed spectrum,";
301 spectra->SetReconstructedSpectrum(hRECpt);
302 // acceptance and efficiency corrections
303 cout << " the efficiency,";
304 spectra->SetAccEffCorrection(hDirectEffpt,hFeedDownEffpt);
305 // spectra->SetfIsStatUncEff(false);
306 // option specific histos (theory)
307 cout << " the theoretical spectra";
309 spectra->SetMCptSpectra(hDirectMCpt,hFeedDownMCpt);
310 if(asym) spectra->SetMCptDistributionsBounds(hDirectMCptMax,hDirectMCptMin,hFeedDownMCptMax,hFeedDownMCptMin);
313 spectra->SetFeedDownMCptSpectra(hFeedDownMCpt);
314 if(asym) spectra->SetFeedDownMCptDistributionsBounds(hFeedDownMCptMax,hFeedDownMCptMin);
317 cout << " and the normalization" <<endl;
318 // Set normalization factors (uncertainties set to 0. as example)
319 spectra->SetNormalization(nevents,sigma);
320 Double_t lumi = nevents / sigma ;
321 Double_t lumiUnc = 0.04*lumi; // 10% uncertainty on the luminosity
322 spectra->SetLuminosity(lumi,lumiUnc);
323 Double_t effTrig = 1.0;
324 spectra->SetTriggerEfficiency(effTrig,0.);
325 if(isRaavsEP>0.) spectra->SetIsEventPlaneAnalysis(kTRUE);
327 // Set the global uncertainties on the efficiencies (in percent)
328 Double_t globalEffUnc = 0.15;
329 Double_t globalBCEffRatioUnc = 0.15;
330 spectra->SetAccEffPercentageUncertainty(globalEffUnc,globalBCEffRatioUnc);
332 // Set if the yield is for particle+anti-particle or only one type
333 spectra->SetIsParticlePlusAntiParticleYield(isParticlePlusAntiParticleYield);
335 // Set the Tab parameter and uncertainties
336 if ( (cc != kpp7) && (cc != kpp276) ) {
337 spectra->SetTabParameter(tab,tabUnc);
340 // Do the calculations
341 cout << " Doing the calculation... "<< endl;
342 Double_t deltaY = 1.0;
343 Double_t branchingRatioC = 1.0;
344 Double_t branchingRatioBintoFinalDecay = 1.0; // this is relative to the input theoretical prediction
345 spectra->ComputeHFPtSpectrum(deltaY,branchingRatioC,branchingRatioBintoFinalDecay);
346 cout << " ended the calculation, getting the histograms back " << endl;
348 // Set the systematics externally
349 Bool_t combineFeedDown = true;
350 AliHFSystErr *systematics = new AliHFSystErr();
352 systematics->SetIsLowEnergy(true);
354 else if ( cc == kpPb0100 ){
355 systematics->SetCollisionType(0);
356 cout <<endl<<" Beware pPb systematics not yet implemented, using pp at 7 TeV !!"<<endl<<endl;
359 else if( cc!=kpp7 ) {
360 systematics->SetCollisionType(1);
361 if ( cc == k07half ) systematics->SetCentrality("07half");
362 else if ( cc == k010 ) systematics->SetCentrality("010");
363 else if ( cc == k1020 ) systematics->SetCentrality("1020");
364 else if ( cc == k020 ) systematics->SetCentrality("020");
365 else if ( cc == k2040 || cc == k2030 || cc == k3040 ) {
366 systematics->SetCentrality("2040");
367 systematics->SetIsPbPb2010EnergyScan(true);
369 else if ( cc == k3050 ) {
370 if (isRaavsEP == kPhiIntegrated) systematics->SetCentrality("4080");
371 else if (isRaavsEP == kInPlane) systematics->SetCentrality("3050InPlane");
372 else if (isRaavsEP == kOutOfPlane) systematics->SetCentrality("3050OutOfPlane");
374 else if ( cc == k4060 || cc == k4050 || cc == k5060 ) systematics->SetCentrality("4060");
375 else if ( cc == k6080 ) systematics->SetCentrality("6080");
376 else if ( cc == k4080 ) systematics->SetCentrality("4080");
378 cout << " Systematics not yet implemented " << endl;
381 } else { systematics->SetCollisionType(0); }
383 systematics->Init(decay);
384 spectra->ComputeSystUncertainties(systematics,combineFeedDown);
387 // Get the output histograms
389 // the corrected yield and cross-section
390 histoYieldCorr = (TH1D*)spectra->GetHistoFeedDownCorrectedSpectrum();
391 histoSigmaCorr = (TH1D*)spectra->GetHistoCrossSectionFromYieldSpectrum();
392 histoYieldCorrMax = (TH1D*)spectra->GetHistoUpperLimitFeedDownCorrectedSpectrum();
393 histoYieldCorrMin = (TH1D*)spectra->GetHistoLowerLimitFeedDownCorrectedSpectrum();
394 histoSigmaCorrMax = (TH1D*)spectra->GetHistoUpperLimitCrossSectionFromYieldSpectrum();
395 histoSigmaCorrMin = (TH1D*)spectra->GetHistoLowerLimitCrossSectionFromYieldSpectrum();
396 histoYieldCorr->SetNameTitle("histoYieldCorr","corrected yield");
397 histoYieldCorrMax->SetNameTitle("histoYieldCorrMax","max corrected yield");
398 histoYieldCorrMin->SetNameTitle("histoYieldCorrMin","min corrected yield");
399 histoSigmaCorr->SetNameTitle("histoSigmaCorr","corrected invariant cross-section");
400 histoSigmaCorrMax->SetNameTitle("histoSigmaCorrMax","max corrected invariant cross-section");
401 histoSigmaCorrMin->SetNameTitle("histoSigmaCorrMin","min corrected invariant cross-section");
403 if(!hDirectEffpt) hDirectEffpt = (TH1D*)spectra->GetDirectAccEffCorrection();
404 if(!hFeedDownEffpt) hFeedDownEffpt = (TH1D*)spectra->GetFeedDownAccEffCorrection();
405 // Get the PbPb Eloss hypothesis histograms
407 histofcRcb = spectra->GetHistoFeedDownCorrectionFcVsEloss();
408 histoYieldCorrRcb = spectra->GetHistoFeedDownCorrectedSpectrumVsEloss();
409 histoSigmaCorrRcb = spectra->GetHistoCrossSectionFromYieldSpectrumVsEloss();
410 histofcRcb->SetName("histofcRcb");
411 histoYieldCorrRcb->SetName("histoYieldCorrRcb");
412 histoSigmaCorrRcb->SetName("histoSigmaCorrRcb");
415 // Get & Rename the TGraphs
416 gSigmaCorr = spectra->GetCrossSectionFromYieldSpectrum();
417 gYieldCorr = spectra->GetFeedDownCorrectedSpectrum();
419 gSigmaCorrExtreme = spectra->GetCrossSectionFromYieldSpectrumExtreme();
420 gYieldCorrExtreme = spectra->GetFeedDownCorrectedSpectrumExtreme();
421 gSigmaCorrConservative = spectra->GetCrossSectionFromYieldSpectrumConservative();
422 gYieldCorrConservative = spectra->GetFeedDownCorrectedSpectrumConservative();
425 // Get & Rename the TGraphs
427 gYieldCorr->SetNameTitle("gYieldCorr","gYieldCorr (uncorr)");
428 gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (uncorr)");
432 histofc = (TH1D*)spectra->GetHistoFeedDownCorrectionFc();
433 histofcMax = (TH1D*)spectra->GetHistoUpperLimitFeedDownCorrectionFc();
434 histofcMin = (TH1D*)spectra->GetHistoLowerLimitFeedDownCorrectionFc();
435 histofc->SetNameTitle("histofc","fc correction factor");
436 histofcMax->SetNameTitle("histofcMax","max fc correction factor");
437 histofcMin->SetNameTitle("histofcMin","min fc correction factor");
439 gYieldCorr->SetNameTitle("gYieldCorr","gYieldCorr (by fc)");
440 gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (by fc)");
441 gFcExtreme = spectra->GetFeedDownCorrectionFcExtreme();
442 gFcExtreme->SetNameTitle("gFcExtreme","gFcExtreme");
443 gYieldCorrExtreme->SetNameTitle("gYieldCorrExtreme","Extreme gYieldCorr (by fc)");
444 gSigmaCorrExtreme->SetNameTitle("gSigmaCorrExtreme","Extreme gSigmaCorr (by fc)");
445 gFcConservative = spectra->GetFeedDownCorrectionFcConservative();
446 gFcConservative->SetNameTitle("gFcConservative","gFcConservative");
447 gYieldCorrConservative->SetNameTitle("gYieldCorrConservative","Conservative gYieldCorr (by fc)");
448 gSigmaCorrConservative->SetNameTitle("gSigmaCorrConservative","Conservative gSigmaCorr (by fc)");
451 if (option==2 && asym) {
452 gYieldCorr->SetNameTitle("gYieldCorr","gYieldCorr (by Nb)");
453 gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (by Nb)");
454 gYieldCorrExtreme->SetNameTitle("gYieldCorrExtreme","Extreme gYieldCorr (by Nb)");
455 gSigmaCorrExtreme->SetNameTitle("gSigmaCorrExtreme","Extreme gSigmaCorr (by Nb)");
456 gYieldCorrConservative->SetNameTitle("gYieldCorrConservative","Conservative gYieldCorr (by Nb)");
457 gSigmaCorrConservative->SetNameTitle("gSigmaCorrConservative","Conservative gSigmaCorr (by Nb)");
458 gFcConservative = spectra->GetFeedDownCorrectionFcConservative();
459 gFcConservative->SetNameTitle("gFcConservative","gFcConservative");
463 nSigma = spectra->GetNtupleCrossSectionVsEloss();
467 // Now, plot the results ! :)
470 gROOT->SetStyle("Plain");
472 cout << " Drawing the results ! " << endl;
477 TCanvas *ceff = new TCanvas("ceff","efficiency drawing");
480 hDirectEffpt->Draw();
482 hFeedDownEffpt->Draw();
485 TCanvas *cTheoryRebin = new TCanvas("cTheoryRebin","control the theoretical spectra rebin");
486 cTheoryRebin->Divide(1,2);
488 hDirectMCpt->Draw("");
489 TH1D *hDirectMCptRebin = (TH1D*)spectra->GetDirectTheoreticalSpectrum();
490 hDirectMCptRebin->SetLineColor(2);
491 hDirectMCptRebin->Draw("same");
493 hFeedDownMCpt->Draw("");
494 TH1D *hFeedDownRebin = (TH1D*)spectra->GetFeedDownTheoreticalSpectrum();
495 hFeedDownRebin->SetLineColor(2);
496 hFeedDownRebin->Draw("same");
497 cTheoryRebin->Update();
499 TCanvas *cTheoryRebinLimits = new TCanvas("cTheoryRebinLimits","control the theoretical spectra limits rebin");
500 cTheoryRebinLimits->Divide(1,2);
501 cTheoryRebinLimits->cd(1);
502 hDirectMCptMax->Draw("");
503 TH1D *hDirectMCptMaxRebin = (TH1D*)spectra->GetDirectTheoreticalUpperLimitSpectrum();
504 hDirectMCptMaxRebin->SetLineColor(2);
505 hDirectMCptMaxRebin->Draw("same");
506 hDirectMCptMin->Draw("same");
507 TH1D *hDirectMCptMinRebin = (TH1D*)spectra->GetDirectTheoreticalLowerLimitSpectrum();
508 hDirectMCptMinRebin->SetLineColor(2);
509 hDirectMCptMinRebin->Draw("same");
510 cTheoryRebinLimits->cd(2);
511 hFeedDownMCptMax->Draw("");
512 TH1D *hFeedDownMaxRebin = (TH1D*)spectra->GetFeedDownTheoreticalUpperLimitSpectrum();
513 hFeedDownMaxRebin->SetLineColor(2);
514 hFeedDownMaxRebin->Draw("same");
515 hFeedDownMCptMin->Draw("same");
516 TH1D *hFeedDownMinRebin = (TH1D*)spectra->GetFeedDownTheoreticalLowerLimitSpectrum();
517 hFeedDownMinRebin->SetLineColor(2);
518 hFeedDownMinRebin->Draw("same");
519 cTheoryRebinLimits->Update();
524 TCanvas * cfc = new TCanvas("cfc","Fc");
525 histofcMax->Draw("c");
526 histofc->Draw("csame");
527 histofcMin->Draw("csame");
531 TH2F *histofcDraw= new TH2F("histofcDraw","histofc (for drawing)",100,0,33.25,100,0.01,1.25);
532 histofcDraw->SetStats(0);
533 histofcDraw->GetXaxis()->SetTitle("p_{T} [GeV]");
534 histofcDraw ->GetXaxis()->SetTitleSize(0.05);
535 histofcDraw->GetXaxis()->SetTitleOffset(0.95);
536 histofcDraw->GetYaxis()->SetTitle(" fc ");
537 histofcDraw->GetYaxis()->SetTitleSize(0.05);
541 // for(Int_t item=0; item<gSigmaCorr->GetN(); item++){
542 // Double_t center=0., value=0.;
543 // gFcExtreme->GetPoint(item,center,value);
544 // Double_t highunc = gFcExtreme->GetErrorYhigh(item) / value ;
545 // Double_t lowunc = gFcExtreme->GetErrorYlow(item) / value ;
546 // cout << "Fc extreme: i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
548 // for(Int_t item=0; item<gSigmaCorr->GetN(); item++){
549 // Double_t center=0., value=0.;
550 // gFcConservative->GetPoint(item,center,value);
551 // Double_t highunc = gFcConservative->GetErrorYhigh(item) / value ;
552 // Double_t lowunc = gFcConservative->GetErrorYlow(item) / value ;
553 // cout << "Fc conservative: i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
555 TCanvas *cfcExtreme = new TCanvas("cfcExtreme","Extreme Asymmetric fc (TGraphAsymmErr)");
556 gFcExtreme->SetFillStyle(3006);
557 gFcExtreme->SetLineWidth(3);
558 gFcExtreme->SetMarkerStyle(20);
559 gFcExtreme->SetFillColor(2);
561 gFcExtreme->Draw("3same");
564 gFcConservative->SetFillStyle(3007);
565 gFcConservative->SetFillColor(4);
566 gFcConservative->Draw("3same");
569 cfcExtreme->Update();
576 // Drawing the results (the raw-reconstructed, the expected, and the corrected spectra)
578 TCanvas * cresult = new TCanvas("cresult","corrected yields & sigma");
579 hDirectMCpt->SetMarkerStyle(20);
580 hDirectMCpt->SetMarkerColor(4);
581 hDirectMCpt->Draw("p");
582 histoSigmaCorr->SetMarkerStyle(21);
583 histoSigmaCorr->SetMarkerColor(2);
584 histoSigmaCorr->Draw("psame");
585 histoYieldCorr->SetMarkerStyle(22);
586 histoYieldCorr->SetMarkerColor(6);
587 histoYieldCorr->Draw("psame");
588 hRECpt->SetMarkerStyle(23);
589 hRECpt->SetMarkerColor(3);
590 hRECpt->Draw("psame");
594 TCanvas * cresult2 = new TCanvas("cresult2","corrected yield & sigma");
595 histoSigmaCorr->SetMarkerStyle(21);
596 histoSigmaCorr->SetMarkerColor(2);
597 histoSigmaCorr->Draw("p");
598 histoYieldCorr->SetMarkerStyle(22);
599 histoYieldCorr->SetMarkerColor(6);
600 histoYieldCorr->Draw("psame");
601 hRECpt->SetMarkerStyle(23);
602 hRECpt->SetMarkerColor(3);
603 hRECpt->Draw("psame");
610 TH2F *histoDraw = new TH2F("histoDraw","histo (for drawing)",100,0,33.25,100,50.,1e7);
611 float max = 1.1*gYieldCorr->GetMaximum();
612 histoDraw->SetAxisRange(0.1,max,"Y");
613 histoDraw->SetStats(0);
614 histoDraw->GetXaxis()->SetTitle("p_{T} [GeV]");
615 histoDraw->GetXaxis()->SetTitleSize(0.05);
616 histoDraw->GetXaxis()->SetTitleOffset(0.95);
617 histoDraw->GetYaxis()->SetTitle("#frac{d#N}{dp_{T}} |_{|y|<1} [L & trigger uncorr]");
618 histoDraw->GetYaxis()->SetTitleSize(0.05);
619 TCanvas * cyieldAsym = new TCanvas("cyieldAsym","Asymmetric corrected yield (TGraphAsymmErr)");
620 gYieldCorr->SetFillStyle(3001);
621 gYieldCorr->SetLineWidth(3);
622 gYieldCorr->SetMarkerStyle(20);
623 gYieldCorr->SetFillColor(3);
625 gYieldCorr->Draw("3LPsame");
626 gYieldCorr->Draw("Xsame");
627 cyieldAsym->SetLogy();
628 cyieldAsym->Update();
630 TCanvas * cyieldExtreme = new TCanvas("cyieldExtreme","Extreme Asymmetric corrected yield (TGraphAsymmErr)");
631 histoYieldCorr->Draw();
632 gYieldCorrExtreme->SetFillStyle(3002);
633 gYieldCorrExtreme->SetLineWidth(3);
634 gYieldCorrExtreme->SetMarkerStyle(20);
635 gYieldCorrExtreme->SetFillColor(2);
636 histoYieldCorr->Draw();
637 gYieldCorr->Draw("3same");
638 gYieldCorrExtreme->Draw("3same");
639 cyieldExtreme->SetLogy();
640 cyieldExtreme->Update();
642 TH2F *histo2Draw = new TH2F("histo2Draw","histo2 (for drawing)",100,0,33.25,100,50.,1e9);
643 max = 1.1*gSigmaCorr->GetMaximum();
644 histo2Draw->SetAxisRange(0.1,max,"Y");
645 histo2Draw->SetStats(0);
646 histo2Draw->GetXaxis()->SetTitle("p_{T} [GeV]");
647 histo2Draw->GetXaxis()->SetTitleSize(0.05);
648 histo2Draw->GetXaxis()->SetTitleOffset(0.95);
649 histo2Draw->GetYaxis()->SetTitle("#frac{1}{BR} #times #frac{d#sigma}{dp_{T}} |_{|y|<1}");
650 histo2Draw->GetYaxis()->SetTitleSize(0.05);
651 TCanvas * csigmaAsym = new TCanvas("csigmaAsym","Asymmetric corrected sigma (TGraphAsymmErr)");
652 gSigmaCorr->SetFillStyle(3001);
653 gSigmaCorr->SetLineWidth(3);
654 gSigmaCorr->SetMarkerStyle(21);
655 gSigmaCorr->SetFillColor(3);
657 gSigmaCorr->Draw("3LPsame");
658 gSigmaCorr->Draw("Xsame");
659 csigmaAsym->SetLogy();
660 csigmaAsym->Update();
662 // cout << endl <<" Sytematics (stat approach) " <<endl;
663 // for(Int_t item=0; item<gSigmaCorr->GetN(); item++){
664 // Double_t center=0., value=0.;
665 // gSigmaCorr->GetPoint(item,center,value);
666 // Double_t highunc = gSigmaCorr->GetErrorYhigh(item) / value ;
667 // Double_t lowunc = gSigmaCorr->GetErrorYlow(item) / value ;
668 // cout << "Sigma syst (stat), i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
671 TCanvas * csigmaExtreme = new TCanvas("csigmaExtreme","Asymmetric extreme corrected sigma (TGraphAsymmErr)");
672 histoSigmaCorr->Draw();
673 gSigmaCorr->Draw("3Psame");
674 gSigmaCorrExtreme->SetFillStyle(3002);
675 gSigmaCorrExtreme->SetLineWidth(3);
676 gSigmaCorrExtreme->SetMarkerStyle(21);
677 gSigmaCorrExtreme->SetFillColor(2);
678 gSigmaCorrExtreme->Draw("3Psame");
679 csigmaExtreme->SetLogy();
680 csigmaExtreme->Update();
682 // cout << endl << " Sytematics (Extreme approach)" <<endl;
683 // for(Int_t item=0; item<gSigmaCorrExtreme->GetN(); item++){
684 // Double_t center=0., value=0.;
685 // gSigmaCorrExtreme->GetPoint(item,center,value);
686 // Double_t highunc = gSigmaCorrExtreme->GetErrorYhigh(item) / value ;
687 // Double_t lowunc = gSigmaCorrExtreme->GetErrorYlow(item) / value ;
688 // cout << "Sigma syst (extreme) i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
691 // cout << endl << " Sytematics (Conservative approach)" <<endl;
692 // for(Int_t item=0; item<gSigmaCorrConservative->GetN(); item++){
693 // Double_t center=0., value=0.;
694 // gSigmaCorrConservative->GetPoint(item,center,value);
695 // Double_t highunc = gSigmaCorrConservative->GetErrorYhigh(item) / value ;
696 // Double_t lowunc = gSigmaCorrConservative->GetErrorYlow(item) / value ;
697 // cout << "Sigma syst (conservative) i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
702 // Draw the PbPb Eloss hypothesis histograms
704 AliHFPtSpectrum *CalcBins=NULL;
705 gStyle->SetPalette(1);
706 TCanvas *canvasfcRcb = new TCanvas("canvasfcRcb","fc vs pt vs Rcb");
707 // histofcRcb->Draw("cont4z");
708 histofcRcb->Draw("colz");
709 canvasfcRcb->Update();
711 TCanvas *canvasfcRcb1 = new TCanvas("canvasfcRcb1","fc vs pt vs Rcb=1");
712 histofcRcb_px = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px",40,40);
713 histofcRcb_px->SetLineColor(2);
716 histofcRcb_px->Draw("same");
717 } else histofcRcb_px->Draw("");
718 canvasfcRcb1->Update();
719 TCanvas *canvasfcRcb2 = new TCanvas("canvasfcRcb2","fc vs pt vs Rcb fixed Rcb");
720 Int_t bin0 = CalcBins->FindTH2YBin(histofcRcb,0.25);
721 Int_t bin1 = CalcBins->FindTH2YBin(histofcRcb,0.5);
722 Int_t bin2 = CalcBins->FindTH2YBin(histofcRcb,1.0);
723 Int_t bin3 = CalcBins->FindTH2YBin(histofcRcb,1.5);
724 Int_t bin4 = CalcBins->FindTH2YBin(histofcRcb,2.0);
725 Int_t bin5 = CalcBins->FindTH2YBin(histofcRcb,3.0);
726 Int_t bin6 = CalcBins->FindTH2YBin(histofcRcb,4.0);
727 TH1D * histofcRcb_px0a = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px0a",bin0,bin0);
728 TH1D * histofcRcb_px0 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px0",bin1,bin1);
729 TH1D * histofcRcb_px1 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px1",bin2,bin2);
730 TH1D * histofcRcb_px2 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px2",bin3,bin3);
731 TH1D * histofcRcb_px3 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px3",bin4,bin4);
732 TH1D * histofcRcb_px4 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px4",bin5,bin5);
733 TH1D * histofcRcb_px5 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px5",bin6,bin6);
736 // histofcRcb_px->Draw("same");
738 // histofcRcb_px->Draw("");
739 histofcRcb_px0a->SetLineColor(2);
740 histofcRcb_px0a->Draw("");
742 histofcRcb_px0a->SetLineColor(2);
743 histofcRcb_px0a->Draw("same");
744 histofcRcb_px0->SetLineColor(4);
745 histofcRcb_px0->Draw("same");
746 histofcRcb_px1->SetLineColor(3);
747 histofcRcb_px1->Draw("same");
748 histofcRcb_px2->SetLineColor(kCyan);
749 histofcRcb_px2->Draw("same");
750 histofcRcb_px3->SetLineColor(kMagenta+1);
751 histofcRcb_px3->Draw("same");
752 histofcRcb_px4->SetLineColor(kOrange+7);
753 histofcRcb_px4->Draw("same");
754 histofcRcb_px5->SetLineColor(kGreen+3);
755 histofcRcb_px5->Draw("same");
756 TLegend *legrcc = new TLegend(0.8,0.8,0.95,0.9);
757 legrcc->SetFillColor(0);
759 legrcc->AddEntry(histofcRcb_px0a,"Rc/b=0.25","l");
760 legrcc->AddEntry(histofcRcb_px0,"Rc/b=0.5","l");
761 legrcc->AddEntry(histofcRcb_px1,"Rc/b=1.0","l");
762 legrcc->AddEntry(histofcRcb_px2,"Rc/b=1.5","l");
763 legrcc->AddEntry(histofcRcb_px3,"Rc/b=2.0","l");
764 legrcc->AddEntry(histofcRcb_px4,"Rc/b=3.0","l");
765 legrcc->AddEntry(histofcRcb_px5,"Rc/b=4.0","l");
767 legrcc->AddEntry(histofcRcb_px0a,"Rb=0.25","l");
768 legrcc->AddEntry(histofcRcb_px0,"Rb=0.5","l");
769 legrcc->AddEntry(histofcRcb_px1,"Rb=1.0","l");
770 legrcc->AddEntry(histofcRcb_px2,"Rb=1.5","l");
771 legrcc->AddEntry(histofcRcb_px3,"Rb=2.0","l");
772 legrcc->AddEntry(histofcRcb_px4,"Rb=3.0","l");
773 legrcc->AddEntry(histofcRcb_px5,"Rb=4.0","l");
776 canvasfcRcb2->Update();
777 TCanvas *canvasYRcb = new TCanvas("canvasYRcb","corrected yield vs pt vs Rcb");
778 histoYieldCorrRcb->Draw("cont4z");
779 canvasYRcb->Update();
780 TCanvas *canvasSRcb = new TCanvas("canvasSRcb","sigma vs pt vs Rcb");
781 histoSigmaCorrRcb->Draw("cont4z");
782 canvasSRcb->Update();
783 TCanvas *canvasSRcb1 = new TCanvas("canvasSRcb1","sigma vs pt vs Rcb fixed Rcb");
784 TH1D * histoSigmaCorrRcb_px0a = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px0a",bin0,bin0);
785 TH1D * histoSigmaCorrRcb_px0 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px0",bin1,bin1);
786 TH1D * histoSigmaCorrRcb_px1 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px1",bin2,bin2);
787 TH1D * histoSigmaCorrRcb_px2 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px2",bin3,bin3);
788 TH1D * histoSigmaCorrRcb_px3 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px3",bin4,bin4);
789 TH1D * histoSigmaCorrRcb_px4 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px4",bin5,bin5);
790 TH1D * histoSigmaCorrRcb_px5 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px5",bin6,bin6);
791 histoSigmaCorr->Draw();
792 histoSigmaCorrRcb_px0a->SetLineColor(2);
793 histoSigmaCorrRcb_px0a->Draw("hsame");
794 histoSigmaCorrRcb_px0->SetLineColor(4);
795 histoSigmaCorrRcb_px0->Draw("hsame");
796 histoSigmaCorrRcb_px1->SetLineColor(3);
797 histoSigmaCorrRcb_px1->Draw("hsame");
798 histoSigmaCorrRcb_px2->SetLineColor(kCyan);
799 histoSigmaCorrRcb_px2->Draw("hsame");
800 histoSigmaCorrRcb_px3->SetLineColor(kMagenta+1);
801 histoSigmaCorrRcb_px3->Draw("hsame");
802 histoSigmaCorrRcb_px4->SetLineColor(kOrange+7);
803 histoSigmaCorrRcb_px4->Draw("same");
804 histoSigmaCorrRcb_px5->SetLineColor(kGreen+3);
805 histoSigmaCorrRcb_px5->Draw("same");
806 TLegend *legrcb = new TLegend(0.8,0.8,0.95,0.9);
807 legrcb->SetFillColor(0);
809 legrcb->AddEntry(histoSigmaCorrRcb_px0a,"Rc/b=0.25","l");
810 legrcb->AddEntry(histoSigmaCorrRcb_px0,"Rc/b=0.5","l");
811 legrcb->AddEntry(histoSigmaCorrRcb_px1,"Rc/b=1.0","l");
812 legrcb->AddEntry(histoSigmaCorrRcb_px2,"Rc/b=1.5","l");
813 legrcb->AddEntry(histoSigmaCorrRcb_px3,"Rc/b=2.0","l");
814 legrcb->AddEntry(histoSigmaCorrRcb_px4,"Rc/b=3.0","l");
815 legrcb->AddEntry(histoSigmaCorrRcb_px5,"Rc/b=4.0","l");
817 legrcb->AddEntry(histoSigmaCorrRcb_px0a,"Rb=0.25","l");
818 legrcb->AddEntry(histoSigmaCorrRcb_px0,"Rb=0.5","l");
819 legrcb->AddEntry(histoSigmaCorrRcb_px1,"Rb=1.0","l");
820 legrcb->AddEntry(histoSigmaCorrRcb_px2,"Rb=1.5","l");
821 legrcb->AddEntry(histoSigmaCorrRcb_px3,"Rb=2.0","l");
822 legrcb->AddEntry(histoSigmaCorrRcb_px4,"Rb=3.0","l");
823 legrcb->AddEntry(histoSigmaCorrRcb_px5,"Rb=4.0","l");
826 canvasSRcb1->Update();
831 // Write the histograms to the output file
833 cout << " Saving the results ! " << endl<< endl;
837 hDirectMCpt->Write(); hFeedDownMCpt->Write();
838 hDirectMCptMax->Write(); hDirectMCptMin->Write();
839 hFeedDownMCptMax->Write(); hFeedDownMCptMin->Write();
840 if(hDirectEffpt) hDirectEffpt->Write(); if(hFeedDownEffpt) hFeedDownEffpt->Write();
843 histoYieldCorr->Write();
844 histoYieldCorrMax->Write(); histoYieldCorrMin->Write();
845 histoSigmaCorr->Write();
846 histoSigmaCorrMax->Write(); histoSigmaCorrMin->Write();
849 histofcRcb->Write(); histofcRcb_px->Write();
850 histoYieldCorrRcb->Write();
851 histoSigmaCorrRcb->Write();
858 if(gYieldCorrExtreme) gYieldCorrExtreme->Write();
859 if(gSigmaCorrExtreme) gSigmaCorrExtreme->Write();
860 if(gYieldCorrConservative) gYieldCorrConservative->Write();
861 if(gSigmaCorrConservative) gSigmaCorrConservative->Write();
862 if(asym && gFcConservative) gFcConservative->Write();
867 histofcMax->Write(); histofcMin->Write();
868 if(asym && gFcExtreme) gFcExtreme->Write();
872 TH1D * hStatUncEffcSigma = spectra->GetDirectStatEffUncOnSigma();
873 TH1D * hStatUncEffbSigma = spectra->GetFeedDownStatEffUncOnSigma();
874 hStatUncEffcSigma->Write();
875 hStatUncEffbSigma->Write();
877 TH1D * hStatUncEffcFD = spectra->GetDirectStatEffUncOnFc();
878 TH1D * hStatUncEffbFD = spectra->GetFeedDownStatEffUncOnFc();
879 hStatUncEffcFD->Write();
880 hStatUncEffbFD->Write();
883 // Draw the cross-section
884 // spectra->DrawSpectrum(gPrediction);