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