AliPhysics  vAN-20151012 (2287573)
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Macros
HFPtSpectrum.C
Go to the documentation of this file.
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 [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
39 //
40 
47 
48 void HFPtSpectrum ( const char *mcfilename="FeedDownCorrectionMC.root",
49  const char *efffilename="Efficiencies.root",
50  const char *recofilename="Reconstructed.root", const char *recohistoname="hRawSpectrumD0",
51  const char *outfilename="HFPtSpectrum.root",
52  Int_t fdMethod=kNb, Double_t nevents=1.0, Double_t sigma=1.0, // sigma[pb]
53  Bool_t isParticlePlusAntiParticleYield=true, Int_t cc=kpp7, Bool_t PbPbEloss=false,
54  Int_t ccestimator = kV0M,
55  Int_t isRaavsEP=kPhiIntegrated,const char *epResolfile="",
56  Int_t rapiditySlice=kdefault,
57  Int_t analysisSpeciality=kTopological) {
58 
59 
60  // gROOT->Macro("$ALICE_PHYSICS/PWGHF/vertexingHF/macros/LoadLibraries.C");
61 
62  // Set if calculation considers asymmetric uncertainties or not
63  Bool_t asym = true;
64 
65  // Set the meson/baryon and decay
66  // (only D0 -> K pi, D+--> K pi pi, D* --> D0 pi, D+s -->KKpi, Lc+ --> pKpi & Lc+ --> pK0S implemented here)
67  Bool_t isD0Kpi = true;
68  Bool_t isDplusKpipi = false;
69  Bool_t isDstarD0pi = false;
70  Bool_t isDsKKpi = false;
71  Bool_t isLctopKpi = false;
72  Bool_t isLcK0Sp = false;
73 
74  Int_t shouldBeOne=0;
75  if(isD0Kpi) shouldBeOne++;
76  if(isDplusKpipi) shouldBeOne++;
77  if(isDstarD0pi) shouldBeOne++;
78  if(isDsKKpi) shouldBeOne++;
79  if(isLctopKpi) shouldBeOne++;
80  if(isLcK0Sp) shouldBeOne++;
81 
82  if (shouldBeOne!=1) {
83  cout << "Sorry, can not deal with more than one correction at the same time"<<endl;
84  return;
85  }
86 
87  Int_t option=3;
88  if (fdMethod==kfc) option=1;
89  else if (fdMethod==kNb) option=2;
90  else if (fdMethod==knone) { option=0; asym=false; }
91  else option=3;
92 
93  if (option>2) {
94  cout<< "Bad calculation option, should be <=2"<<endl;
95  return;
96  }
97 
98 
99  //
100  // Defining the Tab values for the given centrality class
101  // https://twiki.cern.ch/twiki/bin/viewauth/ALICE/CentStudies
102  //
103  Double_t tab = 1., tabUnc = 0.;
104  if( ccestimator == kV0M ) {
105  if ( cc == k07half ) {
106  tab = 24.81; tabUnc = 0.8037;
107  } else if ( cc == k010 ) {
108  tab = 23.48; tabUnc = 0.97;
109  } else if ( cc == k1020 ) {
110  tab = 14.4318; tabUnc = 0.5733;
111  } else if ( cc == k020 ) {
112  tab = 18.93; tabUnc = 0.74;
113  } else if ( cc == k2040 ) {
114  tab = 6.86; tabUnc = 0.28;
115  } else if ( cc == k2030 ) {
116  tab = 8.73769; tabUnc = 0.370219;
117  } else if ( cc == k3040 ) {
118  tab = 5.02755; tabUnc = 0.22099;
119  } else if ( cc == k4050 ) {
120  tab = 2.68327; tabUnc = 0.137073;
121  } else if ( cc == k3050 ) {
122  tab = 3.87011; tabUnc = 0.183847;
123  } else if ( cc == k4060 ) {
124  tab = 2.00; tabUnc= 0.11;
125  } else if ( cc == k4080 ) {
126  tab = 1.20451; tabUnc = 0.071843;
127  } else if ( cc == k5060 ) {
128  tab = 1.32884; tabUnc = 0.0929536;
129  } else if ( cc == k6080 ) {
130  tab = 0.419; tabUnc = 0.033;
131  } else if ( cc == k5080 ) {
132  tab = 0.719; tabUnc = 0.054;
133  } else if ( cc == k80100 ){
134  tab = 0.0690; tabUnc = 0.0062;
135  }
136  }
137 
138  // pPb Glauber (A. Toia)
139  // https://twiki.cern.ch/twiki/bin/viewauth/ALICE/PACentStudies#Glauber_Calculations_with_sigma
140  if( cc == kpPb0100 ){
141  tab = 0.098334; tabUnc = 0.0070679;
142  // A=207.2; B=1.;
143  }
144  else if( ccestimator == kV0A ){
145  if ( cc == kpPb020 ) {
146  tab = 0.183; tabUnc = 0.006245;
147  } else if ( cc == kpPb2040 ) {
148  tab = 0.134; tabUnc = 0.004899;
149  } else if ( cc == kpPb4060 ) {
150  tab = 0.092; tabUnc = 0.004796;
151  } else if ( cc == kpPb60100 ) {
152  tab = 0.041; tabUnc = 0.008832;
153  }
154  }
155  else if( ccestimator == kZNA ){
156  if ( cc == kpPb020 ) {
157  tab = 0.164; tabUnc = 0.010724;
158  } else if ( cc == kpPb2040 ) {
159  tab = 0.137; tabUnc = 0.005099;
160  } else if ( cc == kpPb4060 ) {
161  tab = 0.1011; tabUnc = 0.006;
162  } else if ( cc == kpPb60100 ) {
163  tab = 0.0459; tabUnc = 0.003162;
164  }
165  }
166  else if( ccestimator == kCL1 ){
167  if ( cc == kpPb020 ) {
168  tab = 0.19; tabUnc = 0.007;
169  } else if ( cc == kpPb2040 ) {
170  tab = 0.136; tabUnc = 0.005;
171  } else if ( cc == kpPb4060 ) {
172  tab = 0.088; tabUnc = 0.005;
173  } else if ( cc == kpPb60100 ) {
174  tab = 0.0369; tabUnc = 0.0085;
175  }
176  }
177 
178  tab *= 1e-9; // to pass from mb^{-1} to pb^{-1}
179  tabUnc *= 1e-9;
180 
181 
182 
183  //
184  // Get the histograms from the files
185  //
186  TH1D *hDirectMCpt=0; // Input MC c-->D spectra
187  TH1D *hFeedDownMCpt=0; // Input MC b-->D spectra
188  TH1D *hDirectMCptMax=0; // Input MC maximum c-->D spectra
189  TH1D *hDirectMCptMin=0; // Input MC minimum c-->D spectra
190  TH1D *hFeedDownMCptMax=0; // Input MC maximum b-->D spectra
191  TH1D *hFeedDownMCptMin=0; // Input MC minimum b-->D spectra
192  // TGraphAsymmErrors *gPrediction=0; // Input MC c-->D spectra
193  TH1D *hDirectEffpt=0; // c-->D Acceptance and efficiency correction
194  TH1D *hFeedDownEffpt=0; // b-->D Acceptance and efficiency correction
195  TH1D *hRECpt=0; // all reconstructed D
196 
197  //
198  // Define/Get the input histograms
199  //
200  Int_t decay=0;
201  TFile * mcfile = new TFile(mcfilename,"read");
202  if (isD0Kpi){
203  decay = 1;
204  hDirectMCpt = (TH1D*)mcfile->Get("hD0Kpipred_central");
205  hFeedDownMCpt = (TH1D*)mcfile->Get("hD0KpifromBpred_central_corr");
206  hDirectMCptMax = (TH1D*)mcfile->Get("hD0Kpipred_max");
207  hDirectMCptMin = (TH1D*)mcfile->Get("hD0Kpipred_min");
208  hFeedDownMCptMax = (TH1D*)mcfile->Get("hD0KpifromBpred_max_corr");
209  hFeedDownMCptMin = (TH1D*)mcfile->Get("hD0KpifromBpred_min_corr");
210  // gPrediction = (TGraphAsymmErrors*)mcfile->Get("D0Kpiprediction");
211  }
212  else if (isDplusKpipi){
213  decay = 2;
214  hDirectMCpt = (TH1D*)mcfile->Get("hDpluskpipipred_central");
215  hFeedDownMCpt = (TH1D*)mcfile->Get("hDpluskpipifromBpred_central_corr");
216  hDirectMCptMax = (TH1D*)mcfile->Get("hDpluskpipipred_max");
217  hDirectMCptMin = (TH1D*)mcfile->Get("hDpluskpipipred_min");
218  hFeedDownMCptMax = (TH1D*)mcfile->Get("hDpluskpipifromBpred_max_corr");
219  hFeedDownMCptMin = (TH1D*)mcfile->Get("hDpluskpipifromBpred_min_corr");
220  // gPrediction = (TGraphAsymmErrors*)mcfile->Get("Dpluskpipiprediction");
221  }
222  else if(isDstarD0pi){
223  decay = 3;
224  hDirectMCpt = (TH1D*)mcfile->Get("hDstarD0pipred_central");
225  hFeedDownMCpt = (TH1D*)mcfile->Get("hDstarD0pifromBpred_central_corr");
226  hDirectMCptMax = (TH1D*)mcfile->Get("hDstarD0pipred_max");
227  hDirectMCptMin = (TH1D*)mcfile->Get("hDstarD0pipred_min");
228  hFeedDownMCptMax = (TH1D*)mcfile->Get("hDstarD0pifromBpred_max_corr");
229  hFeedDownMCptMin = (TH1D*)mcfile->Get("hDstarD0pifromBpred_min_corr");
230  // gPrediction = (TGraphAsymmErrors*)mcfile->Get("DstarD0piprediction");
231  }
232  else if (isDsKKpi){
233  decay = 4;
234  hDirectMCpt = (TH1D*)mcfile->Get("hDsPhipitoKkpipred_central");
235  hFeedDownMCpt = (TH1D*)mcfile->Get("hDsPhipitoKkpifromBpred_central_corr");
236  hDirectMCptMax = (TH1D*)mcfile->Get("hDsPhipitoKkpipred_max");
237  hDirectMCptMin = (TH1D*)mcfile->Get("hDsPhipitoKkpipred_min");
238  hFeedDownMCptMax = (TH1D*)mcfile->Get("hDsPhipitoKkpifromBpred_max_corr");
239  hFeedDownMCptMin = (TH1D*)mcfile->Get("hDsPhipitoKkpifromBpred_min_corr");
240  }
241  else if (isLctopKpi){
242  decay = 5;
243  hDirectMCpt = (TH1D*)mcfile->Get("hLcpkpipred_central");
244  hFeedDownMCpt = (TH1D*)mcfile->Get("hLcpkpifromBpred_central_corr");
245  hDirectMCptMax = (TH1D*)mcfile->Get("hLcpkpipred_max");
246  hDirectMCptMin = (TH1D*)mcfile->Get("hLcpkpipred_min");
247  hFeedDownMCptMax = (TH1D*)mcfile->Get("hLcpkpifromBpred_max_corr");
248  hFeedDownMCptMin = (TH1D*)mcfile->Get("hLcpkpifromBpred_min_corr");
249  }
250  else if (isLcK0Sp){
251  decay = 6;
252  hDirectMCpt = (TH1D*)mcfile->Get("hLcK0sppred_central");
253  hFeedDownMCpt = (TH1D*)mcfile->Get("hLcK0spfromBpred_central_corr");
254  hDirectMCptMax = (TH1D*)mcfile->Get("hLcK0sppred_max");
255  hDirectMCptMin = (TH1D*)mcfile->Get("hLcK0sppred_min");
256  hFeedDownMCptMax = (TH1D*)mcfile->Get("hLcK0spfromBpred_max_corr");
257  hFeedDownMCptMin = (TH1D*)mcfile->Get("hLcK0spfromBpred_min_corr");
258  }
259  //
260  hDirectMCpt->SetNameTitle("hDirectMCpt","direct MC spectra");
261  hFeedDownMCpt->SetNameTitle("hFeedDownMCpt","feed-down MC spectra");
262  hDirectMCptMax->SetNameTitle("hDirectMCptMax","max direct MC spectra");
263  hDirectMCptMin->SetNameTitle("hDirectMCptMin","min direct MC spectra");
264  hFeedDownMCptMax->SetNameTitle("hFeedDownMCptMax","max feed-down MC spectra");
265  hFeedDownMCptMin->SetNameTitle("hFeedDownMCptMin","min feed-down MC spectra");
266  //
267  // Scale FONLL inputs if we do the analysis in y-slices
268  //
269  if(rapiditySlice!=kdefault){
270  Double_t scaleFONLL = 1.0;
271  switch(rapiditySlice) {
272  case k08to04: scaleFONLL = (0.093+0.280)/1.0; break;
273  case k07to04: scaleFONLL = 0.280/1.0; break;
274  case k04to01: scaleFONLL = 0.284/1.0; break;
275  case k01to01: scaleFONLL = 0.191/1.0; break;
276  case k01to04: scaleFONLL = 0.288/1.0; break;
277  case k04to07: scaleFONLL = 0.288/1.0; break;
278  case k04to08: scaleFONLL = (0.288+0.096)/1.0; break;
279  case k01to05: scaleFONLL = (0.4)/1.0; break;
280  }
281  hDirectMCpt->Scale(scaleFONLL);
282  hDirectMCptMax->Scale(scaleFONLL);
283  hDirectMCptMin->Scale(scaleFONLL);
284  switch(rapiditySlice) {
285  case k08to04: scaleFONLL = (0.089+0.274)/1.0; break;
286  case k07to04: scaleFONLL = 0.274/1.0; break;
287  case k04to01: scaleFONLL = 0.283/1.0; break;
288  case k01to01: scaleFONLL = 0.192/1.0; break;
289  case k01to04: scaleFONLL = 0.290/1.0; break;
290  case k04to07: scaleFONLL = 0.291/1.0; break;
291  case k04to08: scaleFONLL = (0.291+0.097)/1.0; break;
292  case k01to05: scaleFONLL = (0.4)/1.0; break;
293  }
294  hFeedDownMCpt->Scale(scaleFONLL);
295  hFeedDownMCptMax->Scale(scaleFONLL);
296  hFeedDownMCptMin->Scale(scaleFONLL);
297  }
298 
299  //
300  //
301  TFile * efffile = new TFile(efffilename,"read");
302  hDirectEffpt = (TH1D*)efffile->Get("hEffD");
303  hDirectEffpt->SetNameTitle("hDirectEffpt","direct acc x eff");
304  hFeedDownEffpt = (TH1D*)efffile->Get("hEffB");
305  hFeedDownEffpt->SetNameTitle("hFeedDownEffpt","feed-down acc x eff");
306  //
307  //
308  TFile * recofile = new TFile(recofilename,"read");
309  hRECpt = (TH1D*)recofile->Get(recohistoname);
310  hRECpt->SetNameTitle("hRECpt","Reconstructed spectra");
311  //
312  // Read the file of the EP resolution correction
313  TFile *EPf=0;
314  TH1D *hEPresolCorr=0;
315  if(isRaavsEP>0.){
316  EPf = new TFile(epResolfile,"read");
317  if(isRaavsEP==kInPlane) hEPresolCorr = (TH1D*)EPf->Get("hCorrEPresol_InPlane");
318  else if(isRaavsEP==kOutOfPlane) hEPresolCorr = (TH1D*)EPf->Get("hCorrEPresol_OutOfPlane");
319  for(Int_t i=1; i<=hRECpt->GetNbinsX(); i++) {
320  Double_t value = hRECpt->GetBinContent(i);
321  Double_t error = hRECpt->GetBinError(i);
322  Double_t pt = hRECpt->GetBinCenter(i);
323  Int_t epbin = hEPresolCorr->FindBin( pt );
324  Double_t epcorr = hEPresolCorr->GetBinContent( epbin );
325  value = value*epcorr;
326  error = error*epcorr;
327  hRECpt->SetBinContent(i,value);
328  hRECpt->SetBinError(i,error);
329  }
330  }
331 
332  //
333  // Define the output histograms
334  //
335  TFile *out = new TFile(outfilename,"recreate");
336  //
337  TH1D *histofc=0;
338  TH1D *histofcMax=0;
339  TH1D *histofcMin=0;
340  TH1D *histoYieldCorr=0;
341  TH1D *histoYieldCorrMax=0;
342  TH1D *histoYieldCorrMin=0;
343  TH1D *histoSigmaCorr=0;
344  TH1D *histoSigmaCorrMax=0;
345  TH1D *histoSigmaCorrMin=0;
346  //
347  TH2D *histofcRcb=0;
348  TH1D *histofcRcb_px=0;
349  TH2D *histoYieldCorrRcb=0;
350  TH2D *histoSigmaCorrRcb=0;
351  //
352  TGraphAsymmErrors * gYieldCorr = 0;
353  TGraphAsymmErrors * gSigmaCorr = 0;
354  TGraphAsymmErrors * gFcExtreme = 0;
355  TGraphAsymmErrors * gFcConservative = 0;
356  TGraphAsymmErrors * gYieldCorrExtreme = 0;
357  TGraphAsymmErrors * gSigmaCorrExtreme = 0;
358  TGraphAsymmErrors * gYieldCorrConservative = 0;
359  TGraphAsymmErrors * gSigmaCorrConservative = 0;
360  //
361  TNtuple * nSigma = 0;
362 
363 
364  //
365  // Main functionalities for the calculation
366  //
367 
368  // Define and set the basic option flags
369  AliHFPtSpectrum * spectra = new AliHFPtSpectrum("AliHFPtSpectrum","AliHFPtSpectrum",option);
370  spectra->SetFeedDownCalculationOption(option);
371  spectra->SetComputeAsymmetricUncertainties(asym);
372  // Set flag on whether to additional PbPb Eloss hypothesis have to be computed
373  spectra->SetComputeElossHypothesis(PbPbEloss);
374 
375  // Feed the input histograms
376  // reconstructed spectra
377  cout << " Setting the reconstructed spectrum,";
378  spectra->SetReconstructedSpectrum(hRECpt);
379  // acceptance and efficiency corrections
380  cout << " the efficiency,";
381  spectra->SetAccEffCorrection(hDirectEffpt,hFeedDownEffpt);
382  // spectra->SetfIsStatUncEff(false);
383  // option specific histos (theory)
384  cout << " the theoretical spectra";
385  if(option==1){
386  spectra->SetMCptSpectra(hDirectMCpt,hFeedDownMCpt);
387  if(asym) spectra->SetMCptDistributionsBounds(hDirectMCptMax,hDirectMCptMin,hFeedDownMCptMax,hFeedDownMCptMin);
388  }
389  else if(option==2){
390  spectra->SetFeedDownMCptSpectra(hFeedDownMCpt);
391  if(asym) spectra->SetFeedDownMCptDistributionsBounds(hFeedDownMCptMax,hFeedDownMCptMin);
392  }
393 
394  cout << " and the normalization" <<endl;
395  // Set normalization factors (uncertainties set to 0. as example)
396  spectra->SetNormalization(nevents,sigma);
397  Double_t lumi = nevents / sigma ;
398  Double_t lumiUnc = 0.04*lumi; // 10% uncertainty on the luminosity
399  spectra->SetLuminosity(lumi,lumiUnc);
400  Double_t effTrig = 1.0;
401  spectra->SetTriggerEfficiency(effTrig,0.);
402  if(isRaavsEP>0.) spectra->SetIsEventPlaneAnalysis(kTRUE);
403 
404  // Set the global uncertainties on the efficiencies (in percent)
405  Double_t globalEffUnc = 0.15;
406  Double_t globalBCEffRatioUnc = 0.15;
407  spectra->SetAccEffPercentageUncertainty(globalEffUnc,globalBCEffRatioUnc);
408 
409  // Set if the yield is for particle+anti-particle or only one type
410  spectra->SetIsParticlePlusAntiParticleYield(isParticlePlusAntiParticleYield);
411 
412  // Set the Tab parameter and uncertainties
413  if ( (cc != kpp7) && (cc != kpp276) ) {
414  spectra->SetTabParameter(tab,tabUnc);
415  }
416  if ( cc == kpPb0100 || cc == kpPb020 || cc == kpPb2040 || cc == kpPb4060 || cc == kpPb60100 ) {
417  spectra->SetCollisionType(2);
418  } else if ( !( cc==kpp7 || cc==kpp276 ) ) {
419  spectra->SetCollisionType(1);
420  }
421 
422  // Do the calculations
423  cout << " Doing the calculation... "<< endl;
424  Double_t deltaY = 1.0;
425  Double_t branchingRatioC = 1.0;
426  Double_t branchingRatioBintoFinalDecay = 1.0; // this is relative to the input theoretical prediction
427  spectra->ComputeHFPtSpectrum(deltaY,branchingRatioC,branchingRatioBintoFinalDecay);
428  cout << " ended the calculation, getting the histograms back " << endl;
429 
430  // Set the systematics externally
431 
432  Bool_t combineFeedDown = true;
433  AliHFSystErr *systematics = new AliHFSystErr();
434  if( cc==kpp276 ) {
435  systematics->SetIsLowEnergy(true);
436  }
437  else if ( cc == kpPb0100 || cc == kpPb020 || cc == kpPb2040 || cc == kpPb4060 || cc == kpPb60100 ) {
438  systematics->SetCollisionType(2);
439  // Rapidity slices
440  if(rapiditySlice!=kdefault){
441  systematics->SetIspPb2011RapidityScan(true);
442  TString rapidity="";
443  switch(rapiditySlice) {
444  case k08to04: rapidity="0804"; break;
445  case k07to04: rapidity="0804"; break;
446  case k04to01: rapidity="0401"; break;
447  case k01to01: rapidity="0101"; break;
448  case k01to04: rapidity="0104"; break;
449  case k04to07: rapidity="0408"; break;
450  case k04to08: rapidity="0408"; break;
451  case k01to05: rapidity="0401"; break;
452  }
453  systematics->SetRapidity(rapidity);
454  }
455  // Centrality slices
456  if(ccestimator==kV0A) {
457  if(cc == kpPb020) systematics->SetCentrality("020V0A");
458  else if(cc == kpPb2040) systematics->SetCentrality("2040V0A");
459  else if(cc == kpPb4060) systematics->SetCentrality("4060V0A");
460  else if(cc == kpPb60100) systematics->SetCentrality("60100V0A");
461  } else if (ccestimator==kZNA) {
462  if(cc == kpPb020) systematics->SetCentrality("020ZNA");
463  else if(cc == kpPb2040) systematics->SetCentrality("2040ZNA");
464  else if(cc == kpPb4060) systematics->SetCentrality("4060ZNA");
465  else if(cc == kpPb60100) systematics->SetCentrality("60100ZNA");
466  } else if (ccestimator==kCL1) {
467  if(cc == kpPb020) systematics->SetCentrality("020CL1");
468  else if(cc == kpPb2040) systematics->SetCentrality("2040CL1");
469  else if(cc == kpPb4060) systematics->SetCentrality("4060CL1");
470  else if(cc == kpPb60100) systematics->SetCentrality("60100CL1");
471  } else {
472  if(!(cc == kpPb0100)) {
473  cout <<" Error on the pPb options"<<endl;
474  return;
475  }
476  }
477  }
478  //
479  else if( cc!=kpp7 ) {
480  systematics->SetCollisionType(1);
481  if ( cc == k07half ) systematics->SetCentrality("07half");
482  else if ( cc == k010 ) systematics->SetCentrality("010");
483  else if ( cc == k1020 ) systematics->SetCentrality("1020");
484  else if ( cc == k020 ) systematics->SetCentrality("020");
485  else if ( cc == k2040 || cc == k2030 || cc == k3040 ) {
486  systematics->SetCentrality("2040");
487  systematics->SetIsPbPb2010EnergyScan(true);
488  }
489  else if ( cc == k3050 ) {
490  if (isRaavsEP == kPhiIntegrated) systematics->SetCentrality("4080");
491  else if (isRaavsEP == kInPlane) systematics->SetCentrality("3050InPlane");
492  else if (isRaavsEP == kOutOfPlane) systematics->SetCentrality("3050OutOfPlane");
493  }
494  else if ( cc == k4060 || cc == k4050 || cc == k5060 ) systematics->SetCentrality("4060");
495  else if ( cc == k6080 ) systematics->SetCentrality("6080");
496  else if ( cc == k4080 ) systematics->SetCentrality("4080");
497  else {
498  cout << " Systematics not yet implemented " << endl;
499  return;
500  }
501  } else { systematics->SetCollisionType(0); }
502  if(analysisSpeciality==kLowPt){
503  systematics->SetIsLowPtAnalysis(true);
504  }
505  //
506  systematics->Init(decay);
507 
508  spectra->ComputeSystUncertainties(systematics,combineFeedDown);
509 
510  //
511  // Get the output histograms
512  //
513  // the corrected yield and cross-section
514  histoYieldCorr = (TH1D*)spectra->GetHistoFeedDownCorrectedSpectrum();
515  histoSigmaCorr = (TH1D*)spectra->GetHistoCrossSectionFromYieldSpectrum();
516  histoYieldCorrMax = (TH1D*)spectra->GetHistoUpperLimitFeedDownCorrectedSpectrum();
517  histoYieldCorrMin = (TH1D*)spectra->GetHistoLowerLimitFeedDownCorrectedSpectrum();
518  histoSigmaCorrMax = (TH1D*)spectra->GetHistoUpperLimitCrossSectionFromYieldSpectrum();
519  histoSigmaCorrMin = (TH1D*)spectra->GetHistoLowerLimitCrossSectionFromYieldSpectrum();
520  histoYieldCorr->SetNameTitle("histoYieldCorr","corrected yield");
521  histoYieldCorrMax->SetNameTitle("histoYieldCorrMax","max corrected yield");
522  histoYieldCorrMin->SetNameTitle("histoYieldCorrMin","min corrected yield");
523  histoSigmaCorr->SetNameTitle("histoSigmaCorr","corrected invariant cross-section");
524  histoSigmaCorrMax->SetNameTitle("histoSigmaCorrMax","max corrected invariant cross-section");
525  histoSigmaCorrMin->SetNameTitle("histoSigmaCorrMin","min corrected invariant cross-section");
526  // the efficiencies
527  if(!hDirectEffpt) hDirectEffpt = (TH1D*)spectra->GetDirectAccEffCorrection();
528  if(!hFeedDownEffpt) hFeedDownEffpt = (TH1D*)spectra->GetFeedDownAccEffCorrection();
529  // Get the PbPb Eloss hypothesis histograms
530  if(PbPbEloss){
531  histofcRcb = spectra->GetHistoFeedDownCorrectionFcVsEloss();
532  histoYieldCorrRcb = spectra->GetHistoFeedDownCorrectedSpectrumVsEloss();
533  histoSigmaCorrRcb = spectra->GetHistoCrossSectionFromYieldSpectrumVsEloss();
534  histofcRcb->SetName("histofcRcb");
535  histoYieldCorrRcb->SetName("histoYieldCorrRcb");
536  histoSigmaCorrRcb->SetName("histoSigmaCorrRcb");
537  }
538 
539  // Get & Rename the TGraphs
540  gSigmaCorr = spectra->GetCrossSectionFromYieldSpectrum();
541  gYieldCorr = spectra->GetFeedDownCorrectedSpectrum();
542  if (asym) {
543  gSigmaCorrExtreme = spectra->GetCrossSectionFromYieldSpectrumExtreme();
544  gYieldCorrExtreme = spectra->GetFeedDownCorrectedSpectrumExtreme();
545  gSigmaCorrConservative = spectra->GetCrossSectionFromYieldSpectrumConservative();
546  gYieldCorrConservative = spectra->GetFeedDownCorrectedSpectrumConservative();
547  }
548 
549  // Get & Rename the TGraphs
550  if (option==0){
551  gYieldCorr->SetNameTitle("gYieldCorr","gYieldCorr (uncorr)");
552  gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (uncorr)");
553  }
554  if (option==1){
555  // fc histos
556  histofc = (TH1D*)spectra->GetHistoFeedDownCorrectionFc();
557  histofcMax = (TH1D*)spectra->GetHistoUpperLimitFeedDownCorrectionFc();
558  histofcMin = (TH1D*)spectra->GetHistoLowerLimitFeedDownCorrectionFc();
559  histofc->SetNameTitle("histofc","fc correction factor");
560  histofcMax->SetNameTitle("histofcMax","max fc correction factor");
561  histofcMin->SetNameTitle("histofcMin","min fc correction factor");
562  if (asym) {
563  gYieldCorr->SetNameTitle("gYieldCorr","gYieldCorr (by fc)");
564  gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (by fc)");
565  gFcExtreme = spectra->GetFeedDownCorrectionFcExtreme();
566  gFcExtreme->SetNameTitle("gFcExtreme","gFcExtreme");
567  gYieldCorrExtreme->SetNameTitle("gYieldCorrExtreme","Extreme gYieldCorr (by fc)");
568  gSigmaCorrExtreme->SetNameTitle("gSigmaCorrExtreme","Extreme gSigmaCorr (by fc)");
569  gFcConservative = spectra->GetFeedDownCorrectionFcConservative();
570  gFcConservative->SetNameTitle("gFcConservative","gFcConservative");
571  gYieldCorrConservative->SetNameTitle("gYieldCorrConservative","Conservative gYieldCorr (by fc)");
572  gSigmaCorrConservative->SetNameTitle("gSigmaCorrConservative","Conservative gSigmaCorr (by fc)");
573  }
574  }
575  if (option==2 && asym) {
576  gYieldCorr->SetNameTitle("gYieldCorr","gYieldCorr (by Nb)");
577  gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (by Nb)");
578  gYieldCorrExtreme->SetNameTitle("gYieldCorrExtreme","Extreme gYieldCorr (by Nb)");
579  gSigmaCorrExtreme->SetNameTitle("gSigmaCorrExtreme","Extreme gSigmaCorr (by Nb)");
580  gYieldCorrConservative->SetNameTitle("gYieldCorrConservative","Conservative gYieldCorr (by Nb)");
581  gSigmaCorrConservative->SetNameTitle("gSigmaCorrConservative","Conservative gSigmaCorr (by Nb)");
582  gFcConservative = spectra->GetFeedDownCorrectionFcConservative();
583  gFcConservative->SetNameTitle("gFcConservative","gFcConservative");
584  }
585 
586  if(PbPbEloss){
587  nSigma = spectra->GetNtupleCrossSectionVsEloss();
588  }
589 
590  //
591  // Now, plot the results ! :)
592  //
593 
594  gROOT->SetStyle("Plain");
595 
596  cout << " Drawing the results ! " << endl;
597 
598  // control plots
599  if (option==1) {
600 
601  TCanvas *ceff = new TCanvas("ceff","efficiency drawing");
602  ceff->Divide(1,2);
603  ceff->cd(1);
604  hDirectEffpt->Draw();
605  ceff->cd(2);
606  hFeedDownEffpt->Draw();
607  ceff->Update();
608 
609  TCanvas *cTheoryRebin = new TCanvas("cTheoryRebin","control the theoretical spectra rebin");
610  cTheoryRebin->Divide(1,2);
611  cTheoryRebin->cd(1);
612  hDirectMCpt->Draw("");
613  TH1D *hDirectMCptRebin = (TH1D*)spectra->GetDirectTheoreticalSpectrum();
614  hDirectMCptRebin->SetLineColor(2);
615  hDirectMCptRebin->Draw("same");
616  cTheoryRebin->cd(2);
617  hFeedDownMCpt->Draw("");
618  TH1D *hFeedDownRebin = (TH1D*)spectra->GetFeedDownTheoreticalSpectrum();
619  hFeedDownRebin->SetLineColor(2);
620  hFeedDownRebin->Draw("same");
621  cTheoryRebin->Update();
622 
623  TCanvas *cTheoryRebinLimits = new TCanvas("cTheoryRebinLimits","control the theoretical spectra limits rebin");
624  cTheoryRebinLimits->Divide(1,2);
625  cTheoryRebinLimits->cd(1);
626  hDirectMCptMax->Draw("");
627  TH1D *hDirectMCptMaxRebin = (TH1D*)spectra->GetDirectTheoreticalUpperLimitSpectrum();
628  hDirectMCptMaxRebin->SetLineColor(2);
629  hDirectMCptMaxRebin->Draw("same");
630  hDirectMCptMin->Draw("same");
631  TH1D *hDirectMCptMinRebin = (TH1D*)spectra->GetDirectTheoreticalLowerLimitSpectrum();
632  hDirectMCptMinRebin->SetLineColor(2);
633  hDirectMCptMinRebin->Draw("same");
634  cTheoryRebinLimits->cd(2);
635  hFeedDownMCptMax->Draw("");
636  TH1D *hFeedDownMaxRebin = (TH1D*)spectra->GetFeedDownTheoreticalUpperLimitSpectrum();
637  hFeedDownMaxRebin->SetLineColor(2);
638  hFeedDownMaxRebin->Draw("same");
639  hFeedDownMCptMin->Draw("same");
640  TH1D *hFeedDownMinRebin = (TH1D*)spectra->GetFeedDownTheoreticalLowerLimitSpectrum();
641  hFeedDownMinRebin->SetLineColor(2);
642  hFeedDownMinRebin->Draw("same");
643  cTheoryRebinLimits->Update();
644  }
645 
646  if (option==1) {
647 
648  TCanvas * cfc = new TCanvas("cfc","Fc");
649  histofcMax->Draw("c");
650  histofc->Draw("csame");
651  histofcMin->Draw("csame");
652  cfc->Update();
653 
654  if (asym) {
655  TH2F *histofcDraw= new TH2F("histofcDraw","histofc (for drawing)",100,0,33.25,100,0.01,1.25);
656  histofcDraw->SetStats(0);
657  histofcDraw->GetXaxis()->SetTitle("p_{T} [GeV]");
658  histofcDraw ->GetXaxis()->SetTitleSize(0.05);
659  histofcDraw->GetXaxis()->SetTitleOffset(0.95);
660  histofcDraw->GetYaxis()->SetTitle(" fc ");
661  histofcDraw->GetYaxis()->SetTitleSize(0.05);
662 
663  if (gFcExtreme){
664 
665 // for(Int_t item=0; item<gSigmaCorr->GetN(); item++){
666 // Double_t center=0., value=0.;
667 // gFcExtreme->GetPoint(item,center,value);
668 // Double_t highunc = gFcExtreme->GetErrorYhigh(item) / value ;
669 // Double_t lowunc = gFcExtreme->GetErrorYlow(item) / value ;
670 // cout << "Fc extreme: i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
671 // }
672 // for(Int_t item=0; item<gSigmaCorr->GetN(); item++){
673 // Double_t center=0., value=0.;
674 // gFcConservative->GetPoint(item,center,value);
675 // Double_t highunc = gFcConservative->GetErrorYhigh(item) / value ;
676 // Double_t lowunc = gFcConservative->GetErrorYlow(item) / value ;
677 // cout << "Fc conservative: i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
678 // }
679  TCanvas *cfcExtreme = new TCanvas("cfcExtreme","Extreme Asymmetric fc (TGraphAsymmErr)");
680  gFcExtreme->SetFillStyle(3006);
681  gFcExtreme->SetLineWidth(3);
682  gFcExtreme->SetMarkerStyle(20);
683  gFcExtreme->SetFillColor(2);
684  histofcDraw->Draw();
685  gFcExtreme->Draw("3same");
686 
687  if(gFcConservative){
688  gFcConservative->SetFillStyle(3007);
689  gFcConservative->SetFillColor(4);
690  gFcConservative->Draw("3same");
691  }
692 
693  cfcExtreme->Update();
694  }
695  }
696 
697  }
698 
699  //
700  // Drawing the results (the raw-reconstructed, the expected, and the corrected spectra)
701  //
702  TCanvas * cresult = new TCanvas("cresult","corrected yields & sigma");
703  hDirectMCpt->SetMarkerStyle(20);
704  hDirectMCpt->SetMarkerColor(4);
705  hDirectMCpt->Draw("p");
706  histoSigmaCorr->SetMarkerStyle(21);
707  histoSigmaCorr->SetMarkerColor(2);
708  histoSigmaCorr->Draw("psame");
709  histoYieldCorr->SetMarkerStyle(22);
710  histoYieldCorr->SetMarkerColor(6);
711  histoYieldCorr->Draw("psame");
712  hRECpt->SetMarkerStyle(23);
713  hRECpt->SetMarkerColor(3);
714  hRECpt->Draw("psame");
715  cresult->SetLogy();
716  cresult->Update();
717 
718  TCanvas * cresult2 = new TCanvas("cresult2","corrected yield & sigma");
719  histoSigmaCorr->SetMarkerStyle(21);
720  histoSigmaCorr->SetMarkerColor(2);
721  histoSigmaCorr->Draw("p");
722  histoYieldCorr->SetMarkerStyle(22);
723  histoYieldCorr->SetMarkerColor(6);
724  histoYieldCorr->Draw("psame");
725  hRECpt->SetMarkerStyle(23);
726  hRECpt->SetMarkerColor(3);
727  hRECpt->Draw("psame");
728  cresult2->SetLogy();
729  cresult2->Update();
730 
731 
732  if (asym) {
733 
734  TH2F *histoDraw = new TH2F("histoDraw","histo (for drawing)",100,0,33.25,100,50.,1e7);
735  float max = 1.1*gYieldCorr->GetMaximum();
736  histoDraw->SetAxisRange(0.1,max,"Y");
737  histoDraw->SetStats(0);
738  histoDraw->GetXaxis()->SetTitle("p_{T} [GeV]");
739  histoDraw->GetXaxis()->SetTitleSize(0.05);
740  histoDraw->GetXaxis()->SetTitleOffset(0.95);
741  histoDraw->GetYaxis()->SetTitle("#frac{d#N}{dp_{T}} |_{|y|<1} [L & trigger uncorr]");
742  histoDraw->GetYaxis()->SetTitleSize(0.05);
743  TCanvas * cyieldAsym = new TCanvas("cyieldAsym","Asymmetric corrected yield (TGraphAsymmErr)");
744  gYieldCorr->SetFillStyle(3001);
745  gYieldCorr->SetLineWidth(3);
746  gYieldCorr->SetMarkerStyle(20);
747  gYieldCorr->SetFillColor(3);
748  histoDraw->Draw();
749  gYieldCorr->Draw("3LPsame");
750  gYieldCorr->Draw("Xsame");
751  cyieldAsym->SetLogy();
752  cyieldAsym->Update();
753 
754  TCanvas * cyieldExtreme = new TCanvas("cyieldExtreme","Extreme Asymmetric corrected yield (TGraphAsymmErr)");
755  histoYieldCorr->Draw();
756  gYieldCorrExtreme->SetFillStyle(3002);
757  gYieldCorrExtreme->SetLineWidth(3);
758  gYieldCorrExtreme->SetMarkerStyle(20);
759  gYieldCorrExtreme->SetFillColor(2);
760  histoYieldCorr->Draw();
761  gYieldCorr->Draw("3same");
762  gYieldCorrExtreme->Draw("3same");
763  cyieldExtreme->SetLogy();
764  cyieldExtreme->Update();
765 
766  TH2F *histo2Draw = new TH2F("histo2Draw","histo2 (for drawing)",100,0,33.25,100,50.,1e9);
767  max = 1.1*gSigmaCorr->GetMaximum();
768  histo2Draw->SetAxisRange(0.1,max,"Y");
769  histo2Draw->SetStats(0);
770  histo2Draw->GetXaxis()->SetTitle("p_{T} [GeV]");
771  histo2Draw->GetXaxis()->SetTitleSize(0.05);
772  histo2Draw->GetXaxis()->SetTitleOffset(0.95);
773  histo2Draw->GetYaxis()->SetTitle("#frac{1}{BR} #times #frac{d#sigma}{dp_{T}} |_{|y|<1}");
774  histo2Draw->GetYaxis()->SetTitleSize(0.05);
775  TCanvas * csigmaAsym = new TCanvas("csigmaAsym","Asymmetric corrected sigma (TGraphAsymmErr)");
776  gSigmaCorr->SetFillStyle(3001);
777  gSigmaCorr->SetLineWidth(3);
778  gSigmaCorr->SetMarkerStyle(21);
779  gSigmaCorr->SetFillColor(3);
780  histo2Draw->Draw();
781  gSigmaCorr->Draw("3LPsame");
782  gSigmaCorr->Draw("Xsame");
783  csigmaAsym->SetLogy();
784  csigmaAsym->Update();
785 
786 // cout << endl <<" Sytematics (stat approach) " <<endl;
787 // for(Int_t item=0; item<gSigmaCorr->GetN(); item++){
788 // Double_t center=0., value=0.;
789 // gSigmaCorr->GetPoint(item,center,value);
790 // Double_t highunc = gSigmaCorr->GetErrorYhigh(item) / value ;
791 // Double_t lowunc = gSigmaCorr->GetErrorYlow(item) / value ;
792 // cout << "Sigma syst (stat), i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
793 // }
794 
795  TCanvas * csigmaExtreme = new TCanvas("csigmaExtreme","Asymmetric extreme corrected sigma (TGraphAsymmErr)");
796  histoSigmaCorr->Draw();
797  gSigmaCorr->Draw("3Psame");
798  gSigmaCorrExtreme->SetFillStyle(3002);
799  gSigmaCorrExtreme->SetLineWidth(3);
800  gSigmaCorrExtreme->SetMarkerStyle(21);
801  gSigmaCorrExtreme->SetFillColor(2);
802  gSigmaCorrExtreme->Draw("3Psame");
803  csigmaExtreme->SetLogy();
804  csigmaExtreme->Update();
805 
806 // cout << endl << " Sytematics (Extreme approach)" <<endl;
807 // for(Int_t item=0; item<gSigmaCorrExtreme->GetN(); item++){
808 // Double_t center=0., value=0.;
809 // gSigmaCorrExtreme->GetPoint(item,center,value);
810 // Double_t highunc = gSigmaCorrExtreme->GetErrorYhigh(item) / value ;
811 // Double_t lowunc = gSigmaCorrExtreme->GetErrorYlow(item) / value ;
812 // cout << "Sigma syst (extreme) i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
813 // }
814 
815 // cout << endl << " Sytematics (Conservative approach)" <<endl;
816 // for(Int_t item=0; item<gSigmaCorrConservative->GetN(); item++){
817 // Double_t center=0., value=0.;
818 // gSigmaCorrConservative->GetPoint(item,center,value);
819 // Double_t highunc = gSigmaCorrConservative->GetErrorYhigh(item) / value ;
820 // Double_t lowunc = gSigmaCorrConservative->GetErrorYlow(item) / value ;
821 // cout << "Sigma syst (conservative) i=" << item << ", center=" << center <<", value=" << value << " high unc=" << highunc*100 << "%, low unc=" << lowunc*100 << "%"<<endl;
822 // }
823 
824  }
825 
826  // Draw the PbPb Eloss hypothesis histograms
827  if(PbPbEloss){
828  AliHFPtSpectrum *CalcBins=NULL;
829  gStyle->SetPalette(1);
830  TCanvas *canvasfcRcb = new TCanvas("canvasfcRcb","fc vs pt vs Rcb");
831  // histofcRcb->Draw("cont4z");
832  histofcRcb->Draw("colz");
833  canvasfcRcb->Update();
834  canvasfcRcb->cd(2);
835  TCanvas *canvasfcRcb1 = new TCanvas("canvasfcRcb1","fc vs pt vs Rcb=1");
836  histofcRcb_px = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px",40,40);
837  histofcRcb_px->SetLineColor(2);
838  if (option==1) {
839  histofc->Draw();
840  histofcRcb_px->Draw("same");
841  } else histofcRcb_px->Draw("");
842  canvasfcRcb1->Update();
843  TCanvas *canvasfcRcb2 = new TCanvas("canvasfcRcb2","fc vs pt vs Rcb fixed Rcb");
844  Int_t bin0 = CalcBins->FindTH2YBin(histofcRcb,0.25);
845  Int_t bin1 = CalcBins->FindTH2YBin(histofcRcb,0.5);
846  Int_t bin2 = CalcBins->FindTH2YBin(histofcRcb,1.0);
847  Int_t bin3 = CalcBins->FindTH2YBin(histofcRcb,1.5);
848  Int_t bin4 = CalcBins->FindTH2YBin(histofcRcb,2.0);
849  Int_t bin5 = CalcBins->FindTH2YBin(histofcRcb,3.0);
850  Int_t bin6 = CalcBins->FindTH2YBin(histofcRcb,4.0);
851  TH1D * histofcRcb_px0a = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px0a",bin0,bin0);
852  TH1D * histofcRcb_px0 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px0",bin1,bin1);
853  TH1D * histofcRcb_px1 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px1",bin2,bin2);
854  TH1D * histofcRcb_px2 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px2",bin3,bin3);
855  TH1D * histofcRcb_px3 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px3",bin4,bin4);
856  TH1D * histofcRcb_px4 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px4",bin5,bin5);
857  TH1D * histofcRcb_px5 = (TH1D*)histofcRcb->ProjectionX("histofcRcb_px5",bin6,bin6);
858  if (option==1) {
859  histofc->Draw();
860  // histofcRcb_px->Draw("same");
861  } else {
862  // histofcRcb_px->Draw("");
863  histofcRcb_px0a->SetLineColor(2);
864  histofcRcb_px0a->Draw("");
865  }
866  histofcRcb_px0a->SetLineColor(2);
867  histofcRcb_px0a->Draw("same");
868  histofcRcb_px0->SetLineColor(4);
869  histofcRcb_px0->Draw("same");
870  histofcRcb_px1->SetLineColor(3);
871  histofcRcb_px1->Draw("same");
872  histofcRcb_px2->SetLineColor(kCyan);
873  histofcRcb_px2->Draw("same");
874  histofcRcb_px3->SetLineColor(kMagenta+1);
875  histofcRcb_px3->Draw("same");
876  histofcRcb_px4->SetLineColor(kOrange+7);
877  histofcRcb_px4->Draw("same");
878  histofcRcb_px5->SetLineColor(kGreen+3);
879  histofcRcb_px5->Draw("same");
880  TLegend *legrcc = new TLegend(0.8,0.8,0.95,0.9);
881  legrcc->SetFillColor(0);
882  if (option==1) {
883  legrcc->AddEntry(histofcRcb_px0a,"Rc/b=0.25","l");
884  legrcc->AddEntry(histofcRcb_px0,"Rc/b=0.5","l");
885  legrcc->AddEntry(histofcRcb_px1,"Rc/b=1.0","l");
886  legrcc->AddEntry(histofcRcb_px2,"Rc/b=1.5","l");
887  legrcc->AddEntry(histofcRcb_px3,"Rc/b=2.0","l");
888  legrcc->AddEntry(histofcRcb_px4,"Rc/b=3.0","l");
889  legrcc->AddEntry(histofcRcb_px5,"Rc/b=4.0","l");
890  }else{
891  legrcc->AddEntry(histofcRcb_px0a,"Rb=0.25","l");
892  legrcc->AddEntry(histofcRcb_px0,"Rb=0.5","l");
893  legrcc->AddEntry(histofcRcb_px1,"Rb=1.0","l");
894  legrcc->AddEntry(histofcRcb_px2,"Rb=1.5","l");
895  legrcc->AddEntry(histofcRcb_px3,"Rb=2.0","l");
896  legrcc->AddEntry(histofcRcb_px4,"Rb=3.0","l");
897  legrcc->AddEntry(histofcRcb_px5,"Rb=4.0","l");
898  }
899  legrcc->Draw();
900  canvasfcRcb2->Update();
901  TCanvas *canvasYRcb = new TCanvas("canvasYRcb","corrected yield vs pt vs Rcb");
902  histoYieldCorrRcb->Draw("cont4z");
903  canvasYRcb->Update();
904  TCanvas *canvasSRcb = new TCanvas("canvasSRcb","sigma vs pt vs Rcb");
905  histoSigmaCorrRcb->Draw("cont4z");
906  canvasSRcb->Update();
907  TCanvas *canvasSRcb1 = new TCanvas("canvasSRcb1","sigma vs pt vs Rcb fixed Rcb");
908  TH1D * histoSigmaCorrRcb_px0a = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px0a",bin0,bin0);
909  TH1D * histoSigmaCorrRcb_px0 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px0",bin1,bin1);
910  TH1D * histoSigmaCorrRcb_px1 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px1",bin2,bin2);
911  TH1D * histoSigmaCorrRcb_px2 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px2",bin3,bin3);
912  TH1D * histoSigmaCorrRcb_px3 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px3",bin4,bin4);
913  TH1D * histoSigmaCorrRcb_px4 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px4",bin5,bin5);
914  TH1D * histoSigmaCorrRcb_px5 = (TH1D*)histoSigmaCorrRcb->ProjectionX("histoSigmaCorrRcb_px5",bin6,bin6);
915  histoSigmaCorr->Draw();
916  histoSigmaCorrRcb_px0a->SetLineColor(2);
917  histoSigmaCorrRcb_px0a->Draw("hsame");
918  histoSigmaCorrRcb_px0->SetLineColor(4);
919  histoSigmaCorrRcb_px0->Draw("hsame");
920  histoSigmaCorrRcb_px1->SetLineColor(3);
921  histoSigmaCorrRcb_px1->Draw("hsame");
922  histoSigmaCorrRcb_px2->SetLineColor(kCyan);
923  histoSigmaCorrRcb_px2->Draw("hsame");
924  histoSigmaCorrRcb_px3->SetLineColor(kMagenta+1);
925  histoSigmaCorrRcb_px3->Draw("hsame");
926  histoSigmaCorrRcb_px4->SetLineColor(kOrange+7);
927  histoSigmaCorrRcb_px4->Draw("same");
928  histoSigmaCorrRcb_px5->SetLineColor(kGreen+3);
929  histoSigmaCorrRcb_px5->Draw("same");
930  TLegend *legrcb = new TLegend(0.8,0.8,0.95,0.9);
931  legrcb->SetFillColor(0);
932  if (option==1) {
933  legrcb->AddEntry(histoSigmaCorrRcb_px0a,"Rc/b=0.25","l");
934  legrcb->AddEntry(histoSigmaCorrRcb_px0,"Rc/b=0.5","l");
935  legrcb->AddEntry(histoSigmaCorrRcb_px1,"Rc/b=1.0","l");
936  legrcb->AddEntry(histoSigmaCorrRcb_px2,"Rc/b=1.5","l");
937  legrcb->AddEntry(histoSigmaCorrRcb_px3,"Rc/b=2.0","l");
938  legrcb->AddEntry(histoSigmaCorrRcb_px4,"Rc/b=3.0","l");
939  legrcb->AddEntry(histoSigmaCorrRcb_px5,"Rc/b=4.0","l");
940  }else{
941  legrcb->AddEntry(histoSigmaCorrRcb_px0a,"Rb=0.25","l");
942  legrcb->AddEntry(histoSigmaCorrRcb_px0,"Rb=0.5","l");
943  legrcb->AddEntry(histoSigmaCorrRcb_px1,"Rb=1.0","l");
944  legrcb->AddEntry(histoSigmaCorrRcb_px2,"Rb=1.5","l");
945  legrcb->AddEntry(histoSigmaCorrRcb_px3,"Rb=2.0","l");
946  legrcb->AddEntry(histoSigmaCorrRcb_px4,"Rb=3.0","l");
947  legrcb->AddEntry(histoSigmaCorrRcb_px5,"Rb=4.0","l");
948  }
949  legrcb->Draw();
950  canvasSRcb1->Update();
951  }
952 
953 
954  //
955  // Write the histograms to the output file
956  //
957  cout << " Saving the results ! " << endl<< endl;
958 
959  out->cd();
960  //
961  hDirectMCpt->Write(); hFeedDownMCpt->Write();
962  hDirectMCptMax->Write(); hDirectMCptMin->Write();
963  hFeedDownMCptMax->Write(); hFeedDownMCptMin->Write();
964  if(hDirectEffpt) hDirectEffpt->Write(); if(hFeedDownEffpt) hFeedDownEffpt->Write();
965  hRECpt->Write();
966  //
967  histoYieldCorr->Write();
968  histoYieldCorrMax->Write(); histoYieldCorrMin->Write();
969  histoSigmaCorr->Write();
970  histoSigmaCorrMax->Write(); histoSigmaCorrMin->Write();
971 
972  if(PbPbEloss){
973  histofcRcb->Write(); histofcRcb_px->Write();
974  histoYieldCorrRcb->Write();
975  histoSigmaCorrRcb->Write();
976  nSigma->Write();
977  }
978 
979  gYieldCorr->Write();
980  gSigmaCorr->Write();
981  if(asym){
982  if(gYieldCorrExtreme) gYieldCorrExtreme->Write();
983  if(gSigmaCorrExtreme) gSigmaCorrExtreme->Write();
984  if(gYieldCorrConservative) gYieldCorrConservative->Write();
985  if(gSigmaCorrConservative) gSigmaCorrConservative->Write();
986  if(asym && gFcConservative) gFcConservative->Write();
987  }
988 
989  if(option==1){
990  histofc->Write();
991  histofcMax->Write(); histofcMin->Write();
992  if(asym && gFcExtreme) gFcExtreme->Write();
993  }
994 
995 
996  TH1D * hStatUncEffcSigma = spectra->GetDirectStatEffUncOnSigma();
997  TH1D * hStatUncEffbSigma = spectra->GetFeedDownStatEffUncOnSigma();
998  hStatUncEffcSigma->Write();
999  hStatUncEffbSigma->Write();
1000  if(option!=0){
1001  TH1D * hStatUncEffcFD = spectra->GetDirectStatEffUncOnFc();
1002  TH1D * hStatUncEffbFD = spectra->GetFeedDownStatEffUncOnFc();
1003  hStatUncEffcFD->Write();
1004  hStatUncEffbFD->Write();
1005  }
1006 
1007  // Draw the cross-section
1008  // spectra->DrawSpectrum(gPrediction);
1009 
1010  // out->Close();
1011 
1012 }
particularity
Definition: HFPtSpectrum.C:46
void SetMCptDistributionsBounds(TH1D *hDirectMax, TH1D *hDirectMin, TH1D *hFeedDownMax, TH1D *hFeedDownMin)
Set the theoretical direct & feeddown pt spectrum upper and lower bounds.
void SetIsLowEnergy(Bool_t flag)
Definition: AliHFSystErr.h:64
void SetFeedDownCalculationOption(Int_t option)
Set the calculation option flag for feed-down correction: 0=none, 1=fc , 2=Nb.
TH1D * GetDirectTheoreticalLowerLimitSpectrum() const
TGraphAsymmErrors * GetFeedDownCorrectionFcExtreme() const
Return the TGraphAsymmErrors of the feed-down correction (extreme systematics)
TGraphAsymmErrors * GetFeedDownCorrectedSpectrumConservative() const
Return the TGraphAsymmErrors of the yield after feed-down correction (feed-down conservative systemat...
BFDSubtrMethod
Definition: HFPtSpectrum.C:43
TH1D * GetHistoFeedDownCorrectionFc() const
Return the histogram of the feed-down correction.
TH1D * GetHistoLowerLimitCrossSectionFromYieldSpectrum() const
centrality
void SetCentrality(TString centrality)
Definition: AliHFSystErr.h:60
TGraphAsymmErrors * GetFeedDownCorrectionFcConservative() const
Return the TGraphAsymmErrors of the feed-down correction (conservative systematics) ...
TH1D * GetHistoUpperLimitFeedDownCorrectedSpectrum() const
Return the histogram of the yield after feed-down correction bounds.
TH1D * GetFeedDownTheoreticalLowerLimitSpectrum() const
void SetIsPbPb2010EnergyScan(Bool_t flag)
Definition: AliHFSystErr.h:73
TH1D * GetHistoUpperLimitCrossSectionFromYieldSpectrum() const
Return the equivalent invariant cross-section histogram bounds.
TGraphAsymmErrors * GetCrossSectionFromYieldSpectrumConservative() const
Return the equivalent invariant cross-section TGraphAsymmErrors (feed-down conservative systematics) ...
TH1D * GetHistoLowerLimitFeedDownCorrectedSpectrum() const
void SetFeedDownMCptDistributionsBounds(TH1D *hFeedDownMax, TH1D *hFeedDownMin)
Set the theoretical feeddown pt spectrum upper and lower bounds.
TH1D * GetFeedDownStatEffUncOnSigma() const
void SetTriggerEfficiency(Double_t efficiency, Double_t unc)
Set the trigger efficiency and its uncertainty.
TH1D * GetDirectAccEffCorrection() const
Return the acceptance and efficiency corrections (rebinned if needed)
Double_t * sigma
TH1D * GetHistoUpperLimitFeedDownCorrectionFc() const
Return the histograms of the feed-down correction bounds.
void SetComputeAsymmetricUncertainties(Bool_t flag)
Set if the calculation has to consider asymmetric uncertaInt_ties or not.
RaavsEP
Definition: HFPtSpectrum.C:44
centestimator
TH2D * GetHistoFeedDownCorrectedSpectrumVsEloss() const
Return the histogram of the yield after feed-down correction vs the Ratio(c/b eloss) ...
TH1D * GetDirectTheoreticalSpectrum() const
TGraphAsymmErrors * GetFeedDownCorrectedSpectrumExtreme() const
Return the TGraphAsymmErrors of the yield after feed-down correction (feed-down extreme systematics) ...
TNtuple * GetNtupleCrossSectionVsEloss()
Return the ntuple of the calculation vs the Ratio(c/b eloss)
void SetAccEffPercentageUncertainty(Double_t globalEffUnc, Double_t globalBCEffRatioUnc)
Set global acceptance x efficiency correction uncertainty (in percentages)
void ComputeSystUncertainties(AliHFSystErr *systematics, Bool_t combineFeedDown)
TH1D * GetFeedDownAccEffCorrection() const
TH1D * GetHistoFeedDownCorrectedSpectrum() const
Return the histogram of the yield after feed-down correction.
TGraphAsymmErrors * GetCrossSectionFromYieldSpectrum() const
Return the equivalent invariant cross-section TGraphAsymmErrors (systematics but feed-down) ...
TH1D * GetDirectTheoreticalUpperLimitSpectrum() const
rapidity
Definition: HFPtSpectrum.C:45
void SetNormalization(Double_t normalization)
Set the normalization factors.
void HFPtSpectrum(const char *mcfilename="FeedDownCorrectionMC.root", const char *efffilename="Efficiencies.root", const char *recofilename="Reconstructed.root", const char *recohistoname="hRawSpectrumD0", const char *outfilename="HFPtSpectrum.root", Int_t fdMethod=kNb, Double_t nevents=1.0, Double_t sigma=1.0, Bool_t isParticlePlusAntiParticleYield=true, Int_t cc=kpp7, Bool_t PbPbEloss=false, Int_t ccestimator=kV0M, Int_t isRaavsEP=kPhiIntegrated, const char *epResolfile="", Int_t rapiditySlice=kdefault, Int_t analysisSpeciality=kTopological)
Definition: HFPtSpectrum.C:48
TGraphAsymmErrors * GetCrossSectionFromYieldSpectrumExtreme() const
Return the equivalent invariant cross-section TGraphAsymmErrors (feed-down extreme systematics) ...
void SetReconstructedSpectrum(TH1D *hRec)
Set the reconstructed spectrum.
TH1D * GetHistoLowerLimitFeedDownCorrectionFc() const
Int_t FindTH2YBin(TH2D *histo, Float_t yvalue)
Functionality to find the y-axis bin of a TH2 for a given y-value.
void SetIsLowPtAnalysis(Bool_t flag)
Definition: AliHFSystErr.h:68
void SetFeedDownMCptSpectra(TH1D *hFeedDown)
Set the theoretical feeddown pt spectrum.
void Init(Int_t decay)
Function to initialize the variables/histograms.
TH1D * GetDirectStatEffUncOnSigma() const
void SetIspPb2011RapidityScan(Bool_t flag)
Definition: AliHFSystErr.h:83
TH1D * GetHistoCrossSectionFromYieldSpectrum() const
Return the equivalent invariant cross-section histogram.
void SetComputeElossHypothesis(Bool_t flag)
Set if the calculation has to consider Ratio(c/b eloss) hypothesis.
void ComputeHFPtSpectrum(Double_t deltaY=1.0, Double_t branchingRatioC=1.0, Double_t branchingRatioBintoFinalDecay=1.0)
void SetCollisionType(Int_t ct)
void SetAccEffCorrection(TH1D *hDirectEff, TH1D *hFeedDownEff)
Set the acceptance and efficiency corrections for direct & feeddown.
TH1D * GetFeedDownTheoreticalSpectrum() const
Int_t nevents[nsamples]
void SetCollisionType(Int_t type)
Definition: AliHFSystErr.h:51
void SetIsEventPlaneAnalysis(Bool_t flag)
TH1D * GetFeedDownTheoreticalUpperLimitSpectrum() const
void SetIsParticlePlusAntiParticleYield(Bool_t flag)
Set if the yield is for particle plus anti-particle or not.
TGraphAsymmErrors * GetFeedDownCorrectedSpectrum() const
Return the TGraphAsymmErrors of the yield after feed-down correction (systematics but feed-down) ...
void SetRapidity(TString rapidity)
Settings of rapidity ranges for pPb 0-100% CC.
Definition: AliHFSystErr.h:79
void SetMCptSpectra(TH1D *hDirect, TH1D *hFeedDown)
TH2D * GetHistoCrossSectionFromYieldSpectrumVsEloss() const
TH2D * GetHistoFeedDownCorrectionFcVsEloss() const
Return the histogram of the feed-down correction times the Ratio(c/b eloss)
void SetLuminosity(Double_t luminosity, Double_t unc)
Set the luminosity and its uncertainty.
TH1D * GetFeedDownStatEffUncOnFc() const
TH1D * GetDirectStatEffUncOnFc() const
Histograms to keep track of the influence of the efficiencies statistical uncertainty on the feed-dow...
void SetTabParameter(Double_t tabvalue, Double_t uncertainty)
Set the Tab parameter and its uncertainty.