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