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