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