AliPhysics  7dff9e1 (7dff9e1)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
AliHFPtSpectrum.cxx
Go to the documentation of this file.
1 /**************************************************************************
2  * Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *
3  * *
4  * Author: The ALICE Off-line Project. *
5  * Contributors are mentioned in the code where appropriate. *
6  * *
7  * Permission to use, copy, modify and distribute this software and its *
8  * documentation strictly for non-commercial purposes is hereby granted *
9  * without fee, provided that the above copyright notice appears in all *
10  * copies and that both the copyright notice and this permission notice *
11  * appear in the supporting documentation. The authors make no claims *
12  * about the suitability of this software for any purpose. It is *
13  * provided "as is" without express or implied warranty. *
14  **************************************************************************/
15 
16 /* $Id$ */
17 
18 //***********************************************************************
19 // Class AliHFPtSpectrum
20 // Base class for feed-down corrections on heavy-flavour decays
21 // computes the cross-section via one of the three implemented methods:
22 // 0) Consider no feed-down prediction
23 // 1) Subtract the feed-down with the "fc" method
24 // Yield = Reco * fc; where fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) ) ;
25 // 2) Subtract the feed-down with the "Nb" method
26 // Yield = Reco - Feed-down (exact formula on the function implementation)
27 //
28 // (the corrected yields per bin are divided by the bin-width)
29 //
30 //
31 // In HIC you can also evaluate how the feed-down correction is influenced by an energy loss hypothesis:
32 // Raa(c-->D) / Raa(b-->D) defined here as Rcb for the "fc" method
33 // Raa(b-->D) defined here as Rb for the "Nb" method
34 //
35 // Author: Z.Conesa, zconesa@in2p3.fr
36 //***********************************************************************
37 
38 #include <Riostream.h>
39 
40 #include "TMath.h"
41 #include "TH1.h"
42 #include "TH1D.h"
43 #include "TH2.h"
44 #include "TH2D.h"
45 #include "TNtuple.h"
46 #include "TGraphAsymmErrors.h"
47 #include "TNamed.h"
48 #include "TCanvas.h"
49 #include "TLegend.h"
50 
51 //#include "AliLog.h"
52 #include "AliHFSystErr.h"
53 #include "AliHFPtSpectrum.h"
54 
56 ClassImp(AliHFPtSpectrum);
58 
59 
60 //_________________________________________________________________________________________________________
61 AliHFPtSpectrum::AliHFPtSpectrum(const char* name, const char* title, Int_t option):
62  TNamed(name,title),
63  fhDirectMCpt(NULL),
64  fhFeedDownMCpt(NULL),
65  fhDirectMCptMax(NULL),
66  fhDirectMCptMin(NULL),
67  fhFeedDownMCptMax(NULL),
68  fhFeedDownMCptMin(NULL),
69  fhDirectEffpt(NULL),
70  fhFeedDownEffpt(NULL),
71  fhRECpt(NULL),
72  fgRECSystematics(NULL),
73  fNevts(1),
74  fLuminosity(),
75  fTrigEfficiency(),
76  fGlobalEfficiencyUncertainties(),
77  fTab(),
78  fhFc(NULL),
79  fhFcMax(NULL),
80  fhFcMin(NULL),
81  fhFcRcb(NULL),
82  fgFcExtreme(NULL),
83  fgFcConservative(NULL),
84  fhYieldCorr(NULL),
85  fhYieldCorrMax(NULL),
86  fhYieldCorrMin(NULL),
87  fhYieldCorrRcb(NULL),
88  fgYieldCorr(NULL),
89  fgYieldCorrExtreme(NULL),
90  fgYieldCorrConservative(NULL),
91  fhSigmaCorr(NULL),
92  fhSigmaCorrMax(NULL),
93  fhSigmaCorrMin(NULL),
94  fhSigmaCorrDataSyst(NULL),
95  fhSigmaCorrRcb(NULL),
96  fgSigmaCorr(NULL),
97  fgSigmaCorrExtreme(NULL),
98  fgSigmaCorrConservative(NULL),
99  fnSigma(NULL),
100  fnHypothesis(NULL),
101  fCollisionType(0),
102  fFeedDownOption(option),
103  fAsymUncertainties(kTRUE),
104  fPbPbElossHypothesis(kFALSE),
105  fIsStatUncEff(kTRUE),
106  fParticleAntiParticle(2),
107  fIsEventPlane(kFALSE),
108  fnPtBins(0),
109  fPtBinLimits(NULL),
110  fPtBinWidths(NULL),
111  fhStatUncEffcSigma(NULL),
112  fhStatUncEffbSigma(NULL),
113  fhStatUncEffcFD(NULL),
114  fhStatUncEffbFD(NULL)
115 {
116  //
118  //
119 
120  fLuminosity[0]=1.; fLuminosity[1]=0.;
121  fTrigEfficiency[0]=1.; fTrigEfficiency[1]=0.;
122  fGlobalEfficiencyUncertainties[0]=0.; fGlobalEfficiencyUncertainties[1]=0.;
123  fTab[0]=1.; fTab[1]=0.;
124 
125 }
126 
127 //_________________________________________________________________________________________________________
128 AliHFPtSpectrum::AliHFPtSpectrum(const AliHFPtSpectrum &rhs):
129  TNamed(rhs),
130  fhDirectMCpt(rhs.fhDirectMCpt),
131  fhFeedDownMCpt(rhs.fhFeedDownMCpt),
132  fhDirectMCptMax(rhs.fhDirectMCptMax),
133  fhDirectMCptMin(rhs.fhDirectMCptMin),
134  fhFeedDownMCptMax(rhs.fhFeedDownMCptMax),
135  fhFeedDownMCptMin(rhs.fhFeedDownMCptMin),
136  fhDirectEffpt(rhs.fhDirectEffpt),
137  fhFeedDownEffpt(rhs.fhFeedDownEffpt),
138  fhRECpt(rhs.fhRECpt),
139  fgRECSystematics(rhs.fgRECSystematics),
140  fNevts(rhs.fNevts),
141  fLuminosity(),
142  fTrigEfficiency(),
143  fGlobalEfficiencyUncertainties(),
144  fTab(),
145  fhFc(rhs.fhFc),
146  fhFcMax(rhs.fhFcMax),
147  fhFcMin(rhs.fhFcMin),
148  fhFcRcb(rhs.fhFcRcb),
149  fgFcExtreme(rhs.fgFcExtreme),
150  fgFcConservative(rhs.fgFcConservative),
151  fhYieldCorr(rhs.fhYieldCorr),
152  fhYieldCorrMax(rhs.fhYieldCorrMax),
153  fhYieldCorrMin(rhs.fhYieldCorrMin),
154  fhYieldCorrRcb(rhs.fhYieldCorrRcb),
155  fgYieldCorr(rhs.fgYieldCorr),
156  fgYieldCorrExtreme(rhs.fgYieldCorrExtreme),
157  fgYieldCorrConservative(rhs.fgYieldCorrConservative),
158  fhSigmaCorr(rhs.fhSigmaCorr),
159  fhSigmaCorrMax(rhs.fhSigmaCorrMax),
160  fhSigmaCorrMin(rhs.fhSigmaCorrMin),
161  fhSigmaCorrDataSyst(rhs.fhSigmaCorrDataSyst),
162  fhSigmaCorrRcb(rhs.fhSigmaCorrRcb),
163  fgSigmaCorr(rhs.fgSigmaCorr),
164  fgSigmaCorrExtreme(rhs.fgSigmaCorrExtreme),
165  fgSigmaCorrConservative(rhs.fgSigmaCorrConservative),
166  fnSigma(rhs.fnSigma),
167  fnHypothesis(rhs.fnHypothesis),
168  fCollisionType(rhs.fCollisionType),
169  fFeedDownOption(rhs.fFeedDownOption),
170  fAsymUncertainties(rhs.fAsymUncertainties),
171  fPbPbElossHypothesis(rhs.fPbPbElossHypothesis),
172  fIsStatUncEff(rhs.fIsStatUncEff),
173  fParticleAntiParticle(rhs.fParticleAntiParticle),
174  fIsEventPlane(rhs.fIsEventPlane),
175  fnPtBins(rhs.fnPtBins),
176  fPtBinLimits(NULL),
177  fPtBinWidths(NULL),
178  fhStatUncEffcSigma(NULL),
179  fhStatUncEffbSigma(NULL),
180  fhStatUncEffcFD(NULL),
181  fhStatUncEffbFD(NULL)
182 {
183  //
185  //
186 
187  for(Int_t i=0; i<2; i++){
188  fLuminosity[i] = rhs.fLuminosity[i];
189  fTrigEfficiency[i] = rhs.fTrigEfficiency[i];
190  fGlobalEfficiencyUncertainties[i] = rhs.fGlobalEfficiencyUncertainties[i];
191  fTab[i] = rhs.fTab[i];
192  }
193 
194  fPtBinLimits = new Double_t[fnPtBins+1];
195  fPtBinWidths = new Double_t[fnPtBins];
196  for(Int_t i=0; i<fnPtBins; i++){
197  fPtBinLimits[i] = rhs.fPtBinLimits[i];
198  fPtBinWidths[i] = rhs.fPtBinWidths[i];
199  }
200  fPtBinLimits[fnPtBins] = rhs.fPtBinLimits[fnPtBins];
201 
202 }
203 
204 //_________________________________________________________________________________________________________
205 AliHFPtSpectrum &AliHFPtSpectrum::operator=(const AliHFPtSpectrum &source){
206  //
208  //
209 
210  if (&source == this) return *this;
211 
212  fhDirectMCpt = source.fhDirectMCpt;
213  fhFeedDownMCpt = source.fhFeedDownMCpt;
214  fhDirectMCptMax = source.fhDirectMCptMax;
215  fhDirectMCptMin = source.fhDirectMCptMin;
216  fhFeedDownMCptMax = source.fhFeedDownMCptMax;
217  fhFeedDownMCptMin = source.fhFeedDownMCptMin;
218  fhDirectEffpt = source.fhDirectEffpt;
219  fhFeedDownEffpt = source.fhFeedDownEffpt;
220  fhRECpt = source.fhRECpt;
221  fgRECSystematics = source.fgRECSystematics;
222  fNevts = source.fNevts;
223  fhFc = source.fhFc;
224  fhFcMax = source.fhFcMax;
225  fhFcMin = source.fhFcMin;
226  fhFcRcb = source.fhFcRcb;
227  fgFcExtreme = source.fgFcExtreme;
228  fgFcConservative = source.fgFcConservative;
229  fhYieldCorr = source.fhYieldCorr;
230  fhYieldCorrMax = source.fhYieldCorrMax;
231  fhYieldCorrMin = source.fhYieldCorrMin;
232  fhYieldCorrRcb = source.fhYieldCorrRcb;
233  fgYieldCorr = source.fgYieldCorr;
234  fgYieldCorrExtreme = source.fgYieldCorrExtreme;
235  fgYieldCorrConservative = source.fgYieldCorrConservative;
236  fhSigmaCorr = source.fhSigmaCorr;
237  fhSigmaCorrMax = source.fhSigmaCorrMax;
238  fhSigmaCorrMin = source.fhSigmaCorrMin;
239  fhSigmaCorrDataSyst = source.fhSigmaCorrDataSyst;
240  fhSigmaCorrRcb = source.fhSigmaCorrRcb;
241  fgSigmaCorr = source.fgSigmaCorr;
242  fgSigmaCorrExtreme = source.fgSigmaCorrExtreme;
243  fgSigmaCorrConservative = source.fgSigmaCorrConservative;
244  fnSigma = source.fnSigma;
245  fnHypothesis = source.fnHypothesis;
246  fCollisionType = source.fCollisionType;
247  fFeedDownOption = source.fFeedDownOption;
248  fAsymUncertainties = source.fAsymUncertainties;
249  fPbPbElossHypothesis = source.fPbPbElossHypothesis;
250  fIsStatUncEff = source.fIsStatUncEff;
251  fParticleAntiParticle = source.fParticleAntiParticle;
252  fIsEventPlane = source.fIsEventPlane;
253 
254  for(Int_t i=0; i<2; i++){
255  fLuminosity[i] = source.fLuminosity[i];
256  fTrigEfficiency[i] = source.fTrigEfficiency[i];
257  fGlobalEfficiencyUncertainties[i] = source.fGlobalEfficiencyUncertainties[i];
258  fTab[i] = source.fTab[i];
259  }
260 
261  fnPtBins = source.fnPtBins;
262  if(fPtBinLimits) delete fPtBinLimits;
263  if(fPtBinWidths) delete fPtBinWidths;
264  fPtBinLimits = new Double_t[fnPtBins+1];
265  fPtBinWidths = new Double_t[fnPtBins];
266  for(Int_t i=0; i<fnPtBins; i++){
267  fPtBinLimits[i] = source.fPtBinLimits[i];
268  fPtBinWidths[i] = source.fPtBinWidths[i];
269  }
270  fPtBinLimits[fnPtBins] = source.fPtBinLimits[fnPtBins];
271 
272  return *this;
273 }
274 
275 //_________________________________________________________________________________________________________
276 AliHFPtSpectrum::~AliHFPtSpectrum(){
277  //
279  //
280  if (fhDirectMCpt) delete fhDirectMCpt;
281  if (fhFeedDownMCpt) delete fhFeedDownMCpt;
282  if (fhDirectMCptMax) delete fhDirectMCptMax;
283  if (fhDirectMCptMin) delete fhDirectMCptMin;
284  if (fhFeedDownMCptMax) delete fhFeedDownMCptMax;
285  if (fhFeedDownMCptMin) delete fhFeedDownMCptMin;
286  if (fhDirectEffpt) delete fhDirectEffpt;
287  if (fhFeedDownEffpt) delete fhFeedDownEffpt;
288  if (fhRECpt) delete fhRECpt;
289  if (fgRECSystematics) delete fgRECSystematics;
290  if (fhFc) delete fhFc;
291  if (fhFcMax) delete fhFcMax;
292  if (fhFcMin) delete fhFcMin;
293  if (fhFcRcb) delete fhFcRcb;
294  if (fgFcExtreme) delete fgFcExtreme;
295  if (fgFcConservative) delete fgFcConservative;
296  if (fhYieldCorr) delete fhYieldCorr;
297  if (fhYieldCorrMax) delete fhYieldCorrMax;
298  if (fhYieldCorrMin) delete fhYieldCorrMin;
299  if (fhYieldCorrRcb) delete fhYieldCorrRcb;
300  if (fgYieldCorr) delete fgYieldCorr;
301  if (fgYieldCorrExtreme) delete fgYieldCorrExtreme;
302  if (fgYieldCorrConservative) delete fgYieldCorrConservative;
303  if (fhSigmaCorr) delete fhSigmaCorr;
304  if (fhSigmaCorrMax) delete fhSigmaCorrMax;
305  if (fhSigmaCorrMin) delete fhSigmaCorrMin;
306  if (fhSigmaCorrDataSyst) delete fhSigmaCorrDataSyst;
307  if (fgSigmaCorr) delete fgSigmaCorr;
308  if (fgSigmaCorrExtreme) delete fgSigmaCorrExtreme;
309  if (fgSigmaCorrConservative) delete fgSigmaCorrConservative;
310  if (fnSigma) delete fnSigma;
311  if (fnHypothesis) delete fnHypothesis;
312  if (fPtBinLimits) delete [] fPtBinLimits;
313  if (fPtBinWidths) delete [] fPtBinWidths;
314 }
315 
316 
317 //_________________________________________________________________________________________________________
318 TH1D * AliHFPtSpectrum::RebinTheoreticalSpectra(TH1D *hTheory, const char *name) {
323 
324  if (!hTheory || !fhRECpt) {
325  AliError("Feed-down or reconstructed spectra don't exist");
326  return NULL;
327  }
328 
329  //
330  // Get the reconstructed spectra bins & limits
331  Int_t nbinsMC = hTheory->GetNbinsX();
332 
333  // Check that the reconstructed spectra binning
334  // is larger than the theoretical one
335  Double_t thbinwidth = hTheory->GetBinWidth(1);
336  for (Int_t i=1; i<=fnPtBins; i++) {
337  if ( thbinwidth > fPtBinWidths[i-1] ) {
338  AliInfo(" Beware it seems that the reconstructed spectra has a smaller binning than the theoretical predictions !! ");
339  }
340  }
341 
342  //
343  // Define a new histogram with the real-data reconstructed spectrum binning
344  TH1D * hTheoryRebin = new TH1D(name," theoretical rebinned prediction",fnPtBins,fPtBinLimits);
345 
346  Double_t sum[fnPtBins], items[fnPtBins];
347  for (Int_t ibin=0; ibin<fnPtBins; ibin++) {
348  sum[ibin]=0.; items[ibin]=0.;
349  }
350  for (Int_t ibin=0; ibin<=nbinsMC; ibin++){
351 
352  for (Int_t ibinrec=0; ibinrec<fnPtBins; ibinrec++){
353  if (hTheory->GetBinCenter(ibin)>fPtBinLimits[ibinrec] &&
354  hTheory->GetBinCenter(ibin)<fPtBinLimits[ibinrec+1]){
355  sum[ibinrec]+=hTheory->GetBinContent(ibin);
356  items[ibinrec]+=1.;
357  }
358  }
359 
360  }
361 
362  // set the theoretical rebinned spectra to ( sum-bins / n-bins ) per new bin
363  for (Int_t ibinrec=0; ibinrec<fnPtBins; ibinrec++) {
364  hTheoryRebin->SetBinContent(ibinrec+1,sum[ibinrec]/items[ibinrec]);
365  }
366 
367  return (TH1D*)hTheoryRebin;
368 }
369 
370 //_________________________________________________________________________________________________________
371 void AliHFPtSpectrum::SetMCptSpectra(TH1D *hDirect, TH1D *hFeedDown){
376 
377  if (!hDirect || !hFeedDown || !fhRECpt) {
378  AliError("One or both (direct, feed-down) spectra or the reconstructed spectra don't exist");
379  return;
380  }
381 
382  Bool_t areconsistent = kTRUE;
383  areconsistent = CheckHistosConsistency(hDirect,hFeedDown);
384  if (!areconsistent) {
385  AliInfo("Histograms are not consistent (bin width, bounds)");
386  return;
387  }
388 
389  //
390  // Rebin the theoretical predictions to the reconstructed spectra binning
391  //
392  fhDirectMCpt = RebinTheoreticalSpectra(hDirect,"fhDirectMCpt");
393  fhDirectMCpt->SetNameTitle("fhDirectMCpt"," direct theoretical prediction");
394  fhFeedDownMCpt = RebinTheoreticalSpectra(hFeedDown,"fhFeedDownMCpt");
395  fhFeedDownMCpt->SetNameTitle("fhFeedDownMCpt"," feed-down theoretical prediction");
396 
397 }
398 
399 //_________________________________________________________________________________________________________
400 void AliHFPtSpectrum::SetFeedDownMCptSpectra(TH1D *hFeedDown){
405 
406  if (!hFeedDown || !fhRECpt) {
407  AliError("Feed-down or reconstructed spectra don't exist");
408  return;
409  }
410 
411  //
412  // Rebin the theoretical predictions to the reconstructed spectra binning
413  //
414  fhFeedDownMCpt = RebinTheoreticalSpectra(hFeedDown,"fhFeedDownMCpt");
415  fhFeedDownMCpt->SetNameTitle("fhFeedDownMCpt"," feed-down theoretical prediction");
416 
417 }
418 
419 //_________________________________________________________________________________________________________
420 void AliHFPtSpectrum::SetMCptDistributionsBounds(TH1D *hDirectMax, TH1D *hDirectMin, TH1D *hFeedDownMax, TH1D *hFeedDownMin){
426 
427  if (!hDirectMax || !hDirectMin || !hFeedDownMax|| !hFeedDownMin || !fhRECpt) {
428  AliError("One or all of the max/min direct/feed-down or the reconstructed spectra don't exist");
429  return;
430  }
431 
432  Bool_t areconsistent = kTRUE;
433  areconsistent &= CheckHistosConsistency(hDirectMax,hDirectMin);
434  areconsistent &= CheckHistosConsistency(hFeedDownMax,hFeedDownMin);
435  areconsistent &= CheckHistosConsistency(hDirectMax,hFeedDownMax);
436  if (!areconsistent) {
437  AliInfo("Histograms are not consistent (bin width, bounds)");
438  return;
439  }
440 
441 
442  //
443  // Rebin the theoretical predictions to the reconstructed spectra binning
444  //
445  fhDirectMCptMax = RebinTheoreticalSpectra(hDirectMax,"fhDirectMCptMax");
446  fhDirectMCptMax->SetNameTitle("fhDirectMCptMax"," maximum direct theoretical prediction");
447  fhDirectMCptMin = RebinTheoreticalSpectra(hDirectMin,"fhDirectMCptMin");
448  fhDirectMCptMin->SetNameTitle("fhDirectMCptMin"," minimum direct theoretical prediction");
449  fhFeedDownMCptMax = RebinTheoreticalSpectra(hFeedDownMax,"fhFeedDownMCptMax");
450  fhFeedDownMCptMax->SetNameTitle("fhFeedDownMCptMax"," maximum feed-down theoretical prediction");
451  fhFeedDownMCptMin = RebinTheoreticalSpectra(hFeedDownMin,"fhFeedDownMCptMin");
452  fhFeedDownMCptMin->SetNameTitle("fhFeedDownMCptMin"," minimum feed-down theoretical prediction");
453 
454 }
455 
456 //_________________________________________________________________________________________________________
457 void AliHFPtSpectrum::SetFeedDownMCptDistributionsBounds(TH1D *hFeedDownMax, TH1D *hFeedDownMin){
463 
464  if (!hFeedDownMax || !hFeedDownMin || !fhRECpt) {
465  AliError("One or all of the max/min direct/feed-down spectra don't exist");
466  return;
467  }
468 
469  Bool_t areconsistent = kTRUE;
470  areconsistent &= CheckHistosConsistency(hFeedDownMax,hFeedDownMin);
471  if (!areconsistent) {
472  AliInfo("Histograms are not consistent (bin width, bounds)");
473  return;
474  }
475 
476 
477  //
478  // Rebin the theoretical predictions to the reconstructed spectra binning
479  //
480  fhFeedDownMCptMax = RebinTheoreticalSpectra(hFeedDownMax,"fhFeedDownMCptMax");
481  fhFeedDownMCptMax->SetNameTitle("fhFeedDownMCptMax"," maximum feed-down theoretical prediction");
482  fhFeedDownMCptMin = RebinTheoreticalSpectra(hFeedDownMin,"fhFeedDownMCptMin");
483  fhFeedDownMCptMin->SetNameTitle("fhFeedDownMCptMin"," minimum feed-down theoretical prediction");
484 
485 }
486 
487 //_________________________________________________________________________________________________________
488 void AliHFPtSpectrum::SetDirectAccEffCorrection(TH1D *hDirectEff){
493 
494  if (!hDirectEff) {
495  AliError("The direct acceptance and efficiency corrections doesn't exist");
496  return;
497  }
498 
499  fhDirectEffpt = (TH1D*)hDirectEff->Clone();
500  fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance x efficiency correction");
501 }
502 
503 //_________________________________________________________________________________________________________
504 void AliHFPtSpectrum::SetAccEffCorrection(TH1D *hDirectEff, TH1D *hFeedDownEff){
509 
510  if (!hDirectEff || !hFeedDownEff) {
511  AliError("One or both (direct, feed-down) acceptance and efficiency corrections don't exist");
512  return;
513  }
514 
515  Bool_t areconsistent=kTRUE;
516  areconsistent = CheckHistosConsistency(hDirectEff,hFeedDownEff);
517  if (!areconsistent) {
518  AliInfo("Histograms are not consistent (bin width, bounds)");
519  return;
520  }
521 
522  fhDirectEffpt = (TH1D*)hDirectEff->Clone();
523  fhFeedDownEffpt = (TH1D*)hFeedDownEff->Clone();
524  fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance x efficiency correction");
525  fhFeedDownEffpt->SetNameTitle("fhFeedDownEffpt"," feed-down acceptance x efficiency correction");
526 }
527 
528 //_________________________________________________________________________________________________________
529 void AliHFPtSpectrum::SetReconstructedSpectrum(TH1D *hRec) {
533 
534  if (!hRec) {
535  AliError("The reconstructed spectrum doesn't exist");
536  return;
537  }
538 
539  fhRECpt = (TH1D*)hRec->Clone();
540  fhRECpt->SetNameTitle("fhRECpt"," reconstructed spectrum");
541 
542  //
543  // Evaluate the number of intervals and limits of the spectrum
544  //
545  fnPtBins = fhRECpt->GetNbinsX();
546  fPtBinLimits = new Double_t[fnPtBins+1];
547  fPtBinWidths = new Double_t[fnPtBins];
548  Double_t xlow=0., binwidth=0.;
549  for (Int_t i=1; i<=fnPtBins; i++) {
550  binwidth = fhRECpt->GetBinWidth(i);
551  xlow = fhRECpt->GetBinLowEdge(i);
552  fPtBinLimits[i-1] = xlow;
553  fPtBinWidths[i-1] = binwidth;
554  }
555  fPtBinLimits[fnPtBins] = xlow + binwidth;
556 }
557 
558 //_________________________________________________________________________________________________________
559 void AliHFPtSpectrum::SetReconstructedSpectrumSystematics(TGraphAsymmErrors *gRec) {
563 
564  // Check the compatibility with the reconstructed spectrum
565  Double_t gbinwidth = gRec->GetErrorXlow(1) + gRec->GetErrorXhigh(1) ;
566  Double_t hbinwidth = fhRECpt->GetBinWidth(1);
567  Double_t gxbincenter=0., gybincenter=0.;
568  gRec->GetPoint(1,gxbincenter,gybincenter);
569  Double_t hbincenter = fhRECpt->GetBinCenter(1);
570  if ( (gbinwidth != hbinwidth) || (gxbincenter!=hbincenter) ) {
571  AliError(" The reconstructed spectrum and its systematics don't seem compatible");
572  return;
573  }
574 
575  fgRECSystematics = gRec;
576 }
577 
578 //_________________________________________________________________________________________________________
579 void AliHFPtSpectrum::ComputeHFPtSpectrum(Double_t deltaY, Double_t branchingRatioC, Double_t branchingRatioBintoFinalDecay) {
593  //
594 
595  //
596  // First: Initialization
597  //
598  Bool_t areHistosOk = Initialize();
599  if (!areHistosOk) {
600  AliInfo(" Histos not properly initialized. Check : inconsistent binning ? missing histos ?");
601  return;
602  }
603  // Reset the statistical uncertainties on the efficiencies if needed
604  if(!fIsStatUncEff) ResetStatUncEff();
605 
606  //
607  // Second: Correct for feed-down
608  //
609  if (fFeedDownOption==1) {
610  // Compute the feed-down correction via fc-method
611  CalculateFeedDownCorrectionFc();
612  // Correct the yield for feed-down correction via fc-method
613  CalculateFeedDownCorrectedSpectrumFc();
614  }
615  else if (fFeedDownOption==2) {
616  // Correct the yield for feed-down correction via Nb-method
617  CalculateFeedDownCorrectedSpectrumNb(deltaY,branchingRatioBintoFinalDecay);
618  }
619  else if (fFeedDownOption==0) {
620  // If there is no need for feed-down correction,
621  // the "corrected" yield is equal to the raw yield
622  CalculateCorrectedSpectrumNoFeedDown();
623  }
624  else {
625  AliInfo(" Are you sure the feed-down correction option is right ?");
626  }
627 
628 
629  // Print out information
630  printf("\n\n Correcting the spectra with : \n luminosity = %2.2e +- %2.2e, trigger efficiency = %2.2e +- %2.2e, \n delta_y = %2.2f, BR_c = %2.2e, BR_b_decay = %2.2e \n %2.2f percent uncertainty on the efficiencies, and %2.2f percent uncertainty on the b/c efficiencies ratio \n\n",fLuminosity[0],fLuminosity[1],fTrigEfficiency[0],fTrigEfficiency[1],deltaY,branchingRatioC,branchingRatioBintoFinalDecay,fGlobalEfficiencyUncertainties[0],fGlobalEfficiencyUncertainties[1]);
631  if (fPbPbElossHypothesis) printf("\n\n The considered Tab is %4.2e +- %2.2e \n\n",fTab[0],fTab[1]);
632 
633  //
634  // Finally: Correct from yields to cross-section
635  //
636 
637  // declare the output histograms
638  fhSigmaCorr = new TH1D("fhSigmaCorr","corrected sigma",fnPtBins,fPtBinLimits);
639  fgSigmaCorr = new TGraphAsymmErrors(fnPtBins+1);
640 
641  fhSigmaCorrMax = new TH1D("fhSigmaCorrMax","max corrected sigma",fnPtBins,fPtBinLimits);
642  fhSigmaCorrMin = new TH1D("fhSigmaCorrMin","min corrected sigma",fnPtBins,fPtBinLimits);
643  fhSigmaCorrDataSyst = new TH1D("fhSigmaCorrDataSyst","data syst uncertainties on the corrected sigma",fnPtBins,fPtBinLimits);
644 
645  if (fPbPbElossHypothesis && fFeedDownOption==1) {
646  fhSigmaCorrRcb = new TH2D("fhSigmaCorrRcb","corrected sigma vs Rcb Eloss hypothesis; p_{T} [GeV/c] ; Rcb Eloss hypothesis ; #sigma",fnPtBins,fPtBinLimits,800,0.,4.);
647  fnSigma = new TNtuple("fnSigma"," Sigma ntuple calculation","pt:Signal:Rcb:fc:Yield:Sigma:SigmaStatUnc:SigmaMax:SigmaMin");
648  }
649  if (fPbPbElossHypothesis && fFeedDownOption==2) {
650  fhSigmaCorrRcb = new TH2D("fhSigmaCorrRcb","corrected sigma vs Rb Eloss hypothesis; p_{T} [GeV/c] ; Rb Eloss hypothesis ; #sigma",fnPtBins,fPtBinLimits,800,0.,4.);
651  fnSigma = new TNtuple("fnSigma"," Sigma ntuple calculation","pt:Signal:Rb:fc:Yield:Sigma:SigmaStatUnc:SigmaMax:SigmaMin");
652  }
653  // and the output TGraphAsymmErrors
654  if (fAsymUncertainties){
655  fgSigmaCorrExtreme = new TGraphAsymmErrors(fnPtBins+1);
656  fgSigmaCorrConservative = new TGraphAsymmErrors(fnPtBins+1);
657  }
658  fhStatUncEffcSigma = new TH1D("fhStatUncEffcSigma","direct charm stat unc on the cross section",fnPtBins,fPtBinLimits);
659  fhStatUncEffbSigma = new TH1D("fhStatUncEffbSigma","secondary charm stat unc on the cross section",fnPtBins,fPtBinLimits);
660 
661 
662  // protect against null denominator
663  if (deltaY==0. || fLuminosity[0]==0. || fTrigEfficiency[0]==0. || branchingRatioC==0.) {
664  AliError(" Hey you ! Why luminosity or trigger-efficiency or the c-BR or delta_y are set to zero ?! ");
665  return ;
666  }
667 
668  Double_t value=0, errValue=0, errvalueMax=0., errvalueMin=0.;
669  Double_t errvalueExtremeMax=0., errvalueExtremeMin=0.;
670  Double_t errvalueConservativeMax=0., errvalueConservativeMin=0.;
671  Double_t errvalueStatUncEffc=0.;
672  for(Int_t ibin=1; ibin<=fnPtBins; ibin++){
673 
674  // Variables initialization
675  value=0.; errValue=0.; errvalueMax=0.; errvalueMin=0.;
676  errvalueExtremeMax=0.; errvalueExtremeMin=0.;
677  errvalueConservativeMax=0.; errvalueConservativeMin=0.;
678  errvalueStatUncEffc=0.;
679 
680  // Sigma calculation
681  // Sigma = ( 1. / (lumi * delta_y * BR_c * ParticleAntiPartFactor * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
682  value = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0. && fhRECpt->GetBinContent(ibin)>0.) ?
683  ( fhYieldCorr->GetBinContent(ibin) / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
684  : 0. ;
685 
686  // Sigma statistical uncertainty:
687  // delta_sigma = sigma * sqrt ( (delta_spectra/spectra)^2 )
688  errValue = (value!=0.) ? value * (fhYieldCorr->GetBinError(ibin)/fhYieldCorr->GetBinContent(ibin)) : 0. ;
689 
690  // cout<< " x "<< fhRECpt->GetBinCenter(ibin) << " sigma " << value << " +- "<< errValue << " (stat)"<<endl;
691 
692  //
693  // Sigma systematic uncertainties
694  //
695  if (fAsymUncertainties && value>0.) {
696 
697  // (syst but feed-down) delta_sigma = sigma * sqrt ( (delta_spectra_syst/spectra)^2 +
698  // (delta_lumi/lumi)^2 + (delta_eff_trig/eff_trig)^2 + (delta_eff/eff)^2 + (global_eff)^2 )
699  errvalueMax = value * TMath::Sqrt( (fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin))*(fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin)) +
700  (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
701  (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) +
702  (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
703  fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] );
704  errvalueMin = value * TMath::Sqrt( (fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin))*(fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin)) +
705  (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
706  (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) +
707  (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
708  fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] );
709 
710  // Uncertainties from feed-down
711  // (feed-down syst) delta_sigma = sigma * sqrt ( (delta_spectra_fd/spectra_fd)^2 )
712  // extreme case
713  errvalueExtremeMax = value * (fgYieldCorrExtreme->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin));
714  errvalueExtremeMin = value * (fgYieldCorrExtreme->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin));
715  //
716  // conservative case
717  errvalueConservativeMax = value * (fgYieldCorrConservative->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin));
718  errvalueConservativeMin = value * (fgYieldCorrConservative->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin));
719 
720 
721  // stat unc of the efficiencies, separately
722  errvalueStatUncEffc = value * (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) ;
723 
724  }
725  else {
726  // protect against null denominator
727  errvalueMax = (value!=0.) ?
728  value * TMath::Sqrt( (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
729  (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) +
730  (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
731  fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] )
732  : 0. ;
733  errvalueMin = errvalueMax;
734  }
735 
736  //
737  // Fill the histograms
738  //
739  fhSigmaCorr->SetBinContent(ibin,value);
740  fhSigmaCorr->SetBinError(ibin,errValue);
741 
742  Double_t x = fhYieldCorr->GetBinCenter(ibin);
743  fgSigmaCorr->SetPoint(ibin,x,value); // i,x,y
744  fgSigmaCorr->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
745  fhSigmaCorrMax->SetBinContent(ibin,value+errvalueMax);
746  fhSigmaCorrMin->SetBinContent(ibin,value-errvalueMin);
747 
748  //
749  // Fill the histos and ntuple vs the Eloss hypothesis
750  //
751  if (fPbPbElossHypothesis) {
752 
753  // Loop over the Eloss hypothesis
754  if(!fnHypothesis) {
755  AliError("Error reading the fc hypothesis no ntuple, please check !!");
756  return ;
757  }
758  Int_t entriesHypo = fnHypothesis->GetEntries();
759  Float_t pt=0., Rhy=0., fc=0., fcMin=0., fcMax=0.;
760  fnHypothesis->SetBranchAddress("pt",&pt);
761  if (fFeedDownOption==2) fnHypothesis->SetBranchAddress("Rb",&Rhy);
762  else if (fFeedDownOption==1) fnHypothesis->SetBranchAddress("Rcb",&Rhy);
763  fnHypothesis->SetBranchAddress("fc",&fc);
764  fnHypothesis->SetBranchAddress("fcMax",&fcMax);
765  fnHypothesis->SetBranchAddress("fcMin",&fcMin);
766 
767  for (Int_t item=0; item<entriesHypo; item++) {
768 
769  fnHypothesis->GetEntry(item);
770  if ( TMath::Abs( pt - fhDirectEffpt->GetBinCenter(ibin) ) > 0.15 ) continue;
771 
772  Double_t yieldRcbvalue = (fhRECpt->GetBinContent(ibin) ) ? fhRECpt->GetBinContent(ibin) * fc : 0. ;
773  yieldRcbvalue /= fhRECpt->GetBinWidth(ibin) ;
774  Double_t yieldRcbvalueMax = (fhRECpt->GetBinContent(ibin) ) ? fhRECpt->GetBinContent(ibin) * fcMax : 0. ;
775  yieldRcbvalueMax /= fhRECpt->GetBinWidth(ibin) ;
776  Double_t yieldRcbvalueMin = (fhRECpt->GetBinContent(ibin) ) ? fhRECpt->GetBinContent(ibin) * fcMin : 0. ;
777  yieldRcbvalueMin /= fhRECpt->GetBinWidth(ibin) ;
778 
779  // Sigma calculation
780  // Sigma = ( 1. / (lumi * delta_y * BR_c * ParticleAntiPartFactor * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
781  Double_t sigmaRcbvalue = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)>0.) ?
782  ( yieldRcbvalue / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
783  : 0. ;
784  Double_t sigmaRcbvalueMax = (sigmaRcbvalue!=0.) ?
785  ( yieldRcbvalueMax / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
786  : 0. ;
787  Double_t sigmaRcbvalueMin = (sigmaRcbvalue!=0.) ?
788  ( yieldRcbvalueMin / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
789  : 0. ;
790  // Sigma statistical uncertainty:
791  // delta_sigma = sigma * sqrt ( (delta_spectra/spectra)^2 ) = sigma * ( delta_spectra / (spectra-corr * binwidth) )
792  Double_t sigmaRcbvalueStatUnc = (sigmaRcbvalue!=0.) ?
793  sigmaRcbvalue * ( fhRECpt->GetBinError(ibin) / ( yieldRcbvalue * fhRECpt->GetBinWidth(ibin) ) ) : 0. ;
794 
795  fhSigmaCorrRcb->Fill( fhSigmaCorr->GetBinCenter(ibin) , Rhy, sigmaRcbvalue );
796  // if(ibin==3)
797  // cout << " pt "<< fhRECpt->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<" fc-value="<< fhFcRcb->GetBinContent(ibin,rbin) <<", yield-fcRbvalue="<<yieldRcbvalue<<", sigma-fcRbvalue="<<sigmaRcbvalue<<endl;
798  fnSigma->Fill(fhRECpt->GetBinCenter(ibin), fhRECpt->GetBinContent(ibin),
799  Rhy, fc, yieldRcbvalue, sigmaRcbvalue, sigmaRcbvalueStatUnc, sigmaRcbvalueMax, sigmaRcbvalueMin );
800  }
801  }
802  //
803  // Fill the TGraphAsymmErrors
804  if (fAsymUncertainties) {
805  fgSigmaCorrExtreme->SetPoint(ibin,x,value); // i,x,y
806  fgSigmaCorrExtreme->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
807  fgSigmaCorrConservative->SetPoint(ibin,x,value); // i,x,y
808  fgSigmaCorrConservative->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),errvalueConservativeMin,errvalueConservativeMax); // i,xl,xh,yl,yh
809 
810  fhStatUncEffcSigma->SetBinContent(ibin,0.);
811  if(value>0.) fhStatUncEffcSigma->SetBinError(ibin,((errvalueStatUncEffc/value)*100.));
812  fhStatUncEffbSigma->SetBinContent(ibin,0.); fhStatUncEffbSigma->SetBinError(ibin,0.);
813  // cout << " pt "<< fhRECpt->GetBinCenter(ibin) <<" bin "<< ibin<<" stat-unc-c-sigma "<< errvalueStatUncEffc/value << endl;
814  }
815 
816  }
817 
818 }
819 
820 //_________________________________________________________________________________________________________
821 TH1D * AliHFPtSpectrum::EstimateEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco, const char *name) {
829 
830  if(!fhRECpt){
831  AliInfo("Hey, the reconstructed histogram was not set yet !");
832  return NULL;
833  }
834 
835  TH1D * hEfficiency = new TH1D(name," acceptance #times efficiency",fnPtBins,fPtBinLimits);
836 
837  Double_t *sumSimu=new Double_t[fnPtBins];
838  Double_t *sumReco=new Double_t[fnPtBins];
839  for (Int_t ibin=0; ibin<fnPtBins; ibin++){
840  sumSimu[ibin]=0.; sumReco[ibin]=0.;
841  }
842  for (Int_t ibin=0; ibin<=hSimu->GetNbinsX(); ibin++){
843 
844  for (Int_t ibinrec=0; ibinrec<fnPtBins; ibinrec++){
845  if ( hSimu->GetBinCenter(ibin)>fPtBinLimits[ibinrec] &&
846  hSimu->GetBinCenter(ibin)<fPtBinLimits[ibinrec+1] ) {
847  sumSimu[ibinrec]+=hSimu->GetBinContent(ibin);
848  }
849  if ( hReco->GetBinCenter(ibin)>fPtBinLimits[ibinrec] &&
850  hReco->GetBinCenter(ibin)<fPtBinLimits[ibinrec+1] ) {
851  sumReco[ibinrec]+=hReco->GetBinContent(ibin);
852  }
853  }
854 
855  }
856 
857 
858  // the efficiency is computed as reco/sim (in each bin)
859  // its uncertainty is err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
860  Double_t eff=0., erreff=0.;
861  for (Int_t ibinrec=0; ibinrec<fnPtBins; ibinrec++) {
862  if (sumSimu[ibinrec]!= 0. && sumReco[ibinrec]!=0.) {
863  eff = sumReco[ibinrec] / sumSimu[ibinrec] ;
864  // protection in case eff > 1.0
865  // test calculation (make the argument of the sqrt positive)
866  erreff = TMath::Sqrt( eff * TMath::Abs(1.0 - eff) ) / TMath::Sqrt( sumSimu[ibinrec] );
867  }
868  else { eff=0.0; erreff=0.; }
869  hEfficiency->SetBinContent(ibinrec+1,eff);
870  hEfficiency->SetBinError(ibinrec+1,erreff);
871  }
872 
873  delete [] sumSimu;
874  delete [] sumReco;
875 
876  return (TH1D*)hEfficiency;
877 }
878 
879 //_________________________________________________________________________________________________________
880 void AliHFPtSpectrum::EstimateAndSetDirectEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco) {
888 
889  if(!fhRECpt || !hSimu || !hReco){
890  AliError("Hey, the reconstructed histogram was not set yet !");
891  return;
892  }
893 
894  fhDirectEffpt = EstimateEfficiencyRecoBin(hSimu,hReco,"fhDirectEffpt");
895  fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance #times efficiency");
896 
897 }
898 
899 //_________________________________________________________________________________________________________
900 void AliHFPtSpectrum::EstimateAndSetFeedDownEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco) {
908 
909  if(!fhRECpt || !hSimu || !hReco){
910  AliError("Hey, the reconstructed histogram was not set yet !");
911  return;
912  }
913 
914  fhFeedDownEffpt = EstimateEfficiencyRecoBin(hSimu,hReco,"fhFeedDownEffpt");
915  fhFeedDownEffpt->SetNameTitle("fhFeedDownEffpt"," feed-down acceptance #times efficiency");
916 
917 }
918 
919 //_________________________________________________________________________________________________________
920 Bool_t AliHFPtSpectrum::Initialize(){
924 
925  if (fFeedDownOption==0) {
926  AliInfo("Getting ready for the corrections without feed-down consideration");
927  } else if (fFeedDownOption==1) {
928  AliInfo("Getting ready for the fc feed-down correction calculation");
929  } else if (fFeedDownOption==2) {
930  AliInfo("Getting ready for the Nb feed-down correction calculation");
931  } else { AliError("The calculation option must be <=2"); return kFALSE; }
932 
933  // Start checking the input histograms consistency
934  Bool_t areconsistent=kTRUE;
935 
936  // General checks
937  if (!fhDirectEffpt || !fhRECpt) {
938  AliError(" Reconstructed spectra and/or the Nc efficiency distributions are not defined");
939  return kFALSE;
940  }
941  areconsistent &= CheckHistosConsistency(fhRECpt,fhDirectEffpt);
942  if (!areconsistent) {
943  AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
944  return kFALSE;
945  }
946  if (fFeedDownOption==0) return kTRUE;
947 
948  //
949  // Common checks for options 1 (fc) & 2(Nb)
950  if (!fhFeedDownMCpt || !fhFeedDownEffpt) {
951  AliError(" Theoretical Nb and/or the Nb efficiency distributions are not defined");
952  return kFALSE;
953  }
954  areconsistent &= CheckHistosConsistency(fhRECpt,fhFeedDownMCpt);
955  areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownEffpt);
956  if (fAsymUncertainties) {
957  if (!fhFeedDownMCptMax || !fhFeedDownMCptMin) {
958  AliError(" Max/Min theoretical Nb distributions are not defined");
959  return kFALSE;
960  }
961  areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownMCptMax);
962  }
963  if (!areconsistent) {
964  AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
965  return kFALSE;
966  }
967  if (fFeedDownOption>1) return kTRUE;
968 
969  //
970  // Now checks for option 1 (fc correction)
971  if (!fhDirectMCpt) {
972  AliError("Theoretical Nc distributions is not defined");
973  return kFALSE;
974  }
975  areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhFeedDownMCpt);
976  areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectEffpt);
977  if (fAsymUncertainties) {
978  if (!fhDirectMCptMax || !fhDirectMCptMin) {
979  AliError(" Max/Min theoretical Nc distributions are not defined");
980  return kFALSE;
981  }
982  areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectMCptMax);
983  }
984  if (!areconsistent) {
985  AliInfo("Histograms required for fc correction are not consistent (bin width, bounds)");
986  return kFALSE;
987  }
988 
989  return kTRUE;
990 }
991 
992 //_________________________________________________________________________________________________________
993 Bool_t AliHFPtSpectrum::CheckHistosConsistency(TH1D *h1, TH1D *h2){
997 
998  if (!h1 || !h2) {
999  AliError("One or both histograms don't exist");
1000  return kFALSE;
1001  }
1002 
1003  Double_t binwidth1 = h1->GetBinWidth(1);
1004  Double_t binwidth2 = h2->GetBinWidth(1);
1005  Double_t min1 = h1->GetBinCenter(1) - (binwidth1/2.) ;
1006 // Double_t max1 = h1->GetBinCenter(nbins1) + (binwidth1/2.) ;
1007  Double_t min2 = h2->GetBinCenter(1) - (binwidth2/2.) ;
1008 // Double_t max2 = h2->GetBinCenter(nbins2) + (binwidth2/2.) ;
1009 
1010  if (binwidth1!=binwidth2) {
1011  AliInfo(" histograms with different bin width");
1012  return kFALSE;
1013  }
1014  if (min1!=min2) {
1015  AliInfo(" histograms with different minimum");
1016  return kFALSE;
1017  }
1018 // if (max1!=max2) {
1019 // AliInfo(" histograms with different maximum");
1020 // return kFALSE;
1021 // }
1022 
1023  return kTRUE;
1024 }
1025 
1026 //_________________________________________________________________________________________________________
1027 void AliHFPtSpectrum::CalculateCorrectedSpectrumNoFeedDown(){
1031 
1032 
1033  // declare the output histograms
1034  fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (no feed-down corr)",fnPtBins,fPtBinLimits);
1035  fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (no feed-down corr)",fnPtBins,fPtBinLimits);
1036  fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (no feed-down corr)",fnPtBins,fPtBinLimits);
1037  fgYieldCorr = new TGraphAsymmErrors(fnPtBins+1);
1038 
1039  //
1040  // Do the calculation
1041  //
1042  for (Int_t ibin=1; ibin<=fnPtBins; ibin++) {
1043  Double_t value = fhRECpt->GetBinContent(ibin) /fhRECpt->GetBinWidth(ibin);
1044  Double_t errvalue = fhRECpt->GetBinError(ibin) /fhRECpt->GetBinWidth(ibin);
1045  fhYieldCorr->SetBinContent(ibin,value);
1046  fhYieldCorr->SetBinError(ibin,errvalue);
1047  fhYieldCorrMax->SetBinContent(ibin,value);
1048  fhYieldCorrMax->SetBinError(ibin,errvalue);
1049  fhYieldCorrMin->SetBinContent(ibin,value);
1050  fhYieldCorrMin->SetBinError(ibin,errvalue);
1051  }
1052 
1053  fAsymUncertainties=kFALSE;
1054 }
1055 
1056 //_________________________________________________________________________________________________________
1057 void AliHFPtSpectrum::CalculateFeedDownCorrectionFc(){
1069  AliInfo("Calculating the feed-down correction factor (fc method)");
1070 
1071  // define the variables
1072  Double_t correction=1.;
1073  Double_t theoryRatio=1.;
1074  Double_t effRatio=1.;
1075  Double_t correctionExtremeA=1., correctionExtremeB=1.;
1076  Double_t theoryRatioExtremeA=1., theoryRatioExtremeB=1.;
1077  Double_t correctionConservativeA=1., correctionConservativeB=1.;
1078  Double_t theoryRatioConservativeA=1., theoryRatioConservativeB=1.;
1079  // Double_t correctionUnc=0.;
1080  Double_t correctionExtremeAUnc=0., correctionExtremeBUnc=0.;
1081  Double_t correctionConservativeAUnc=0., correctionConservativeBUnc=0.;
1082 
1083  // declare the output histograms
1084  fhFc = new TH1D("fhFc","fc correction factor",fnPtBins,fPtBinLimits);
1085  fhFcMax = new TH1D("fhFcMax","max fc correction factor",fnPtBins,fPtBinLimits);
1086  fhFcMin = new TH1D("fhFcMin","min fc correction factor",fnPtBins,fPtBinLimits);
1087  if(fPbPbElossHypothesis) {
1088  fhFcRcb = new TH2D("fhFcRcb","fc correction factor vs Rcb Eloss hypothesis; p_{T} [GeV/c] ; Rcb Eloss hypothesis ; fc correction",fnPtBins,fPtBinLimits,800,0.,4.);
1089  fnHypothesis = new TNtuple("fnHypothesis"," Feed-down correction vs hypothesis (fc)","pt:Rcb:fc:fcMax:fcMin");
1090  }
1091  // two local control histograms
1092  TH1D *hTheoryRatio = new TH1D("hTheoryRatio","Theoretical B-->D over c-->D (feed-down/direct) ratio",fnPtBins,fPtBinLimits);
1093  TH1D *hEffRatio = new TH1D("hEffRatio","Efficiency B-->D over c-->D (feed-down/direct) ratio",fnPtBins,fPtBinLimits);
1094  // and the output TGraphAsymmErrors
1095  if (fAsymUncertainties) {
1096  fgFcExtreme = new TGraphAsymmErrors(fnPtBins+1);
1097  fgFcExtreme->SetNameTitle("fgFcExtreme","fgFcExtreme");
1098  fgFcConservative = new TGraphAsymmErrors(fnPtBins+1);
1099  fgFcConservative->SetNameTitle("fgFcConservative","fgFcConservative");
1100  }
1101 
1102  fhStatUncEffcFD = new TH1D("fhStatUncEffcFD","direct charm stat unc on the feed-down correction",fnPtBins,fPtBinLimits);
1103  fhStatUncEffbFD = new TH1D("fhStatUncEffbFD","secondary charm stat unc on the feed-down correction",fnPtBins,fPtBinLimits);
1104  Double_t correctionConservativeAUncStatEffc=0., correctionConservativeBUncStatEffc=0.;
1105  Double_t correctionConservativeAUncStatEffb=0., correctionConservativeBUncStatEffb=0.;
1106 
1107  //
1108  // Compute fc
1109  //
1110  for (Int_t ibin=1; ibin<=fnPtBins; ibin++) {
1111 
1112  // Variables initialization
1113  correction=1.; theoryRatio=1.; effRatio=1.;
1114  correctionExtremeA=1.; correctionExtremeB=1.;
1115  theoryRatioExtremeA=1.; theoryRatioExtremeB=1.;
1116  correctionConservativeA=1.; correctionConservativeB=1.;
1117  theoryRatioConservativeA=1.; theoryRatioConservativeB=1.;
1118  // correctionUnc=0.;
1119  correctionExtremeAUnc=0.; correctionExtremeBUnc=0.;
1120  correctionConservativeAUnc=0.; correctionConservativeBUnc=0.;
1121  correctionConservativeAUncStatEffc=0.; correctionConservativeBUncStatEffc=0.;
1122  correctionConservativeAUncStatEffb=0.; correctionConservativeBUncStatEffb=0.;
1123 
1124  // theory_ratio = (N_b/N_c)
1125  theoryRatio = (fhDirectMCpt->GetBinContent(ibin)>0. && fhFeedDownMCpt->GetBinContent(ibin)>0.) ?
1126  fhFeedDownMCpt->GetBinContent(ibin) / fhDirectMCpt->GetBinContent(ibin) : 1.0 ;
1127 
1128  //
1129  // Calculate the uncertainty [ considering only the theoretical uncertainties on Nb & Nc for now !!! ]
1130  //
1131  // extreme A = direct-max, feed-down-min
1132  theoryRatioExtremeA = (fhDirectMCptMax->GetBinContent(ibin)>0. && fhFeedDownMCptMin->GetBinContent(ibin)>0.) ?
1133  fhFeedDownMCptMin->GetBinContent(ibin) / fhDirectMCptMax->GetBinContent(ibin) : 1.0 ;
1134  // extreme B = direct-min, feed-down-max
1135  theoryRatioExtremeB = (fhDirectMCptMin->GetBinContent(ibin)>0. && fhDirectMCptMax->GetBinContent(ibin)>0.) ?
1136  fhFeedDownMCptMax->GetBinContent(ibin) / fhDirectMCptMin->GetBinContent(ibin) : 1.0 ;
1137  // conservative A = direct-max, feed-down-max
1138  theoryRatioConservativeA = (fhDirectMCptMax->GetBinContent(ibin)>0. && fhFeedDownMCptMin->GetBinContent(ibin)>0.) ?
1139  fhFeedDownMCptMax->GetBinContent(ibin) / fhDirectMCptMax->GetBinContent(ibin) : 1.0 ;
1140  // conservative B = direct-min, feed-down-min
1141  theoryRatioConservativeB = (fhDirectMCptMin->GetBinContent(ibin)>0. && fhDirectMCptMax->GetBinContent(ibin)>0.) ?
1142  fhFeedDownMCptMin->GetBinContent(ibin) / fhDirectMCptMin->GetBinContent(ibin) : 1.0 ;
1143 
1144  // eff_ratio = (eff_b/eff_c)
1145  effRatio = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ?
1146  fhFeedDownEffpt->GetBinContent(ibin) / fhDirectEffpt->GetBinContent(ibin) : 1.0 ;
1147 
1148  // fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
1149  if( TMath::Abs(effRatio - 1.0)<0.0001 || TMath::Abs(theoryRatio - 1.0)<0.0001 ) {
1150  correction = 1.0;
1151  correctionExtremeA = 1.0;
1152  correctionExtremeB = 1.0;
1153  correctionConservativeA = 1.0;
1154  correctionConservativeB = 1.0;
1155  }
1156  else {
1157  correction = ( 1. / ( 1 + ( effRatio * theoryRatio ) ) );
1158  correctionExtremeA = ( 1. / ( 1 + ( effRatio * theoryRatioExtremeA ) ) );
1159  correctionExtremeB = ( 1. / ( 1 + ( effRatio * theoryRatioExtremeB ) ) );
1160  correctionConservativeA = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeA ) ) );
1161  correctionConservativeB = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeB ) ) );
1162  }
1163 
1164 
1165  // fc uncertainty from (eff_b/eff_c) = fc^2 * (N_b/N_c) * delta(eff_b/eff_c)
1166  // delta(eff_b/eff_c) is a percentage = effRatio * sqrt( fGlobalEfficiencyUncertainties[1]^2 + unc_eff_c ^2 + unc_eff_b ^2 )
1167  Double_t relEffUnc = TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] +
1168  (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
1169  (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))
1170  );
1171 
1172  // correctionUnc = correction*correction * theoryRatio * effRatio * relEffUnc;
1173  correctionExtremeAUnc = correctionExtremeA*correctionExtremeA * theoryRatioExtremeA * effRatio * relEffUnc;
1174  correctionExtremeBUnc = correctionExtremeB*correctionExtremeB * theoryRatioExtremeB * effRatio * relEffUnc;
1175 
1176  correctionConservativeAUnc = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA *effRatio * relEffUnc;
1177  //
1178  correctionConservativeAUncStatEffc = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA *effRatio *
1179  (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin));
1180  correctionConservativeAUncStatEffb = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA *effRatio *
1181  (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin));
1182 
1183  correctionConservativeBUnc = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB *effRatio * relEffUnc;
1184 
1185  correctionConservativeBUncStatEffb = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB *effRatio *
1186  (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin));
1187  correctionConservativeBUncStatEffc = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB *effRatio *
1188  (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin));
1189 
1190 
1191  // Fill in the histograms
1192  hTheoryRatio->SetBinContent(ibin,theoryRatio);
1193  hEffRatio->SetBinContent(ibin,effRatio);
1194  fhFc->SetBinContent(ibin,correction);
1195  //
1196  // Estimate how the result varies vs charm/beauty Eloss hypothesis
1197  //
1198  if ( correction>1.0e-16 && fPbPbElossHypothesis){
1199  // Loop over the Eloss hypothesis
1200  // Int_t rbin=0;
1201  for (Float_t rval=0.0025; rval<4.0; rval+=0.005){
1202  // Central fc with Eloss hypothesis
1203  Double_t correctionRcb = ( 1. / ( 1 + ( effRatio * theoryRatio * (1/rval) ) ) );
1204  fhFcRcb->Fill( fhFc->GetBinCenter(ibin) , rval, correctionRcb );
1205  // if(ibin==3){
1206  // cout << " pt "<< fhFc->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<", fc-Rcb-value="<<correctionRcb<<endl;
1207  // rbin++;
1208  // }
1209  // Upper / lower fc with up / low FONLL bands and Eloss hypothesis
1210  Double_t correctionConservativeARcb = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeA * (1/rval) ) ) );
1211  Double_t correctionConservativeBRcb = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeB * (1/rval) ) ) );
1212  Double_t correctionConservativeARcbUnc = correctionConservativeARcb*correctionConservativeARcb * theoryRatioConservativeA * (1/rval) *effRatio * relEffUnc;
1213  Double_t correctionConservativeBRcbUnc = correctionConservativeBRcb*correctionConservativeBRcb * theoryRatioConservativeB * (1/rval) *effRatio * relEffUnc;
1214  Double_t consvalRcb[4] = { correctionConservativeARcb - correctionConservativeARcbUnc, correctionConservativeARcb + correctionConservativeARcbUnc,
1215  correctionConservativeBRcb - correctionConservativeBRcbUnc, correctionConservativeBRcb + correctionConservativeBRcbUnc};
1216  Double_t uncConservativeRcbMin = correctionRcb - TMath::MinElement(4,consvalRcb);
1217  Double_t uncConservativeRcbMax = TMath::MaxElement(4,consvalRcb) - correctionRcb;
1218 // if(ibin==3)
1219 // cout << " pt="<<fhDirectEffpt->GetBinCenter(ibin)<<", hypo="<<rval<<", fc="<<correctionRcb<<" +"<<uncConservativeRcbMax<<" -"<<uncConservativeRcbMin<<endl;
1220  fnHypothesis->Fill( fhDirectEffpt->GetBinCenter(ibin), rval, correctionRcb, correctionRcb+uncConservativeRcbMax, correctionRcb-uncConservativeRcbMin);
1221  }
1222  }
1223  //
1224  // Fill the rest of (asymmetric) histograms
1225  //
1226  if (fAsymUncertainties) {
1227  Double_t x = fhDirectMCpt->GetBinCenter(ibin);
1228  Double_t val[4] = { correctionExtremeA + correctionExtremeAUnc, correctionExtremeA - correctionExtremeAUnc,
1229  correctionExtremeB + correctionExtremeBUnc, correctionExtremeB - correctionExtremeBUnc };
1230  Double_t uncExtremeMin = correction - TMath::MinElement(4,val);
1231  Double_t uncExtremeMax = TMath::MaxElement(4,val) - correction;
1232  fgFcExtreme->SetPoint(ibin,x,correction); // i,x,y
1233  fgFcExtreme->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),uncExtremeMin,uncExtremeMax); // i,xl,xh,yl,yh
1234  fhFcMax->SetBinContent(ibin,correction+uncExtremeMax);
1235  fhFcMin->SetBinContent(ibin,correction-uncExtremeMin);
1236  Double_t consval[4] = { correctionConservativeA - correctionConservativeAUnc, correctionConservativeA + correctionConservativeAUnc,
1237  correctionConservativeB - correctionConservativeBUnc, correctionConservativeB + correctionConservativeBUnc};
1238  Double_t uncConservativeMin = correction - TMath::MinElement(4,consval);
1239  Double_t uncConservativeMax = TMath::MaxElement(4,consval) - correction;
1240  fgFcConservative->SetPoint(ibin,x,correction); // i,x,y
1241  fgFcConservative->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),uncConservativeMin,uncConservativeMax); // i,xl,xh,yl,yh
1242  if( !(correction>0.) ){
1243  fgFcExtreme->SetPoint(ibin,x,0.); // i,x,y
1244  fgFcExtreme->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),0.,0.); // i,xl,xh,yl,yh
1245  fgFcConservative->SetPoint(ibin,x,0.); // i,x,y
1246  fgFcConservative->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),0.,0.); // i,xl,xh,yl,yh
1247  }
1248 
1249  Double_t valStatEffc[2] = { correctionConservativeAUncStatEffc/correctionConservativeA,
1250  correctionConservativeBUncStatEffc/correctionConservativeB };
1251  Double_t valStatEffb[2] = { correctionConservativeAUncStatEffb/correctionConservativeA,
1252  correctionConservativeBUncStatEffb/correctionConservativeB };
1253  Double_t uncConservativeStatEffc = TMath::MaxElement(2,valStatEffc);
1254  Double_t uncConservativeStatEffb = TMath::MaxElement(2,valStatEffb);
1255  fhStatUncEffcFD->SetBinContent(ibin,0.); fhStatUncEffcFD->SetBinError(ibin,uncConservativeStatEffc*100.);
1256  fhStatUncEffbFD->SetBinContent(ibin,0.); fhStatUncEffbFD->SetBinError(ibin,uncConservativeStatEffb*100.);
1257  // cout << " pt "<< fhStatUncEffcFD->GetBinCenter(ibin) <<" bin "<< ibin<<" fc-stat-c ="<<uncConservativeStatEffc<<" fc-stat-b ="<<uncConservativeStatEffb<<endl;
1258  }
1259 
1260  }
1261 
1262 }
1263 
1264 //_________________________________________________________________________________________________________
1265 void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumFc(){
1277 
1278  AliInfo(" Calculating the feed-down corrected spectrum (fc method)");
1279 
1280  if (!fhFc || !fhRECpt) {
1281  AliError(" Reconstructed or fc distributions are not defined");
1282  return;
1283  }
1284 
1285  Double_t value = 0., errvalue = 0., errvalueMax= 0., errvalueMin= 0.;
1286  Double_t valueExtremeMax= 0., valueExtremeMin= 0.;
1287  Double_t valueConservativeMax= 0., valueConservativeMin= 0.;
1288 
1289  // declare the output histograms
1290  fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by fc)",fnPtBins,fPtBinLimits);
1291  fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by fc)",fnPtBins,fPtBinLimits);
1292  fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by fc)",fnPtBins,fPtBinLimits);
1293  if(fPbPbElossHypothesis) fhYieldCorrRcb = new TH2D("fhYieldCorrRcb","corrected yield (by fc) vs Rcb Eloss hypothesis; p_{T} [GeV/c] ; Rcb Eloss hypothesis ; corrected yield",fnPtBins,fPtBinLimits,800,0.,4.);
1294  // and the output TGraphAsymmErrors
1295  fgYieldCorr = new TGraphAsymmErrors(fnPtBins+1);
1296  if (fAsymUncertainties){
1297  fgYieldCorrExtreme = new TGraphAsymmErrors(fnPtBins+1);
1298  fgYieldCorrConservative = new TGraphAsymmErrors(fnPtBins+1);
1299  }
1300 
1301  //
1302  // Do the calculation
1303  //
1304  for (Int_t ibin=1; ibin<=fnPtBins; ibin++) {
1305 
1306  // Variables initialization
1307  value = 0.; errvalue = 0.; errvalueMax= 0.; errvalueMin= 0.;
1308  valueExtremeMax= 0.; valueExtremeMin= 0.;
1309  valueConservativeMax= 0.; valueConservativeMin= 0.;
1310 
1311 
1312  // calculate the value
1313  // physics = reco * fc / bin-width
1314  value = (fhRECpt->GetBinContent(ibin) && fhFc->GetBinContent(ibin)) ?
1315  fhRECpt->GetBinContent(ibin) * fhFc->GetBinContent(ibin) : 0. ;
1316  value /= fhRECpt->GetBinWidth(ibin) ;
1317 
1318  // Statistical uncertainty
1319  // (stat) delta_physics = physics * sqrt ( (delta_reco/reco)^2 )
1320  errvalue = (value!=0. && fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) ?
1321  value * (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) : 0. ;
1322 
1323  // Calculate the systematic uncertainties
1324  // (syst but feed-down) delta_physics = physics * sqrt ( (delta_reco_syst/reco)^2 )
1325  // (feed-down syst) delta_physics = physics * sqrt ( (delta_fc/fc)^2 )
1326  //
1327  // Protect against null denominator. If so, define uncertainty as null
1328  if (fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) {
1329 
1330  if (fAsymUncertainties) {
1331 
1332  // Systematics but feed-down
1333  if (fgRECSystematics) {
1334  errvalueMax = value * ( fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinContent(ibin) );
1335  errvalueMin = value * ( fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinContent(ibin) );
1336  }
1337  else { errvalueMax = 0.; errvalueMin = 0.; }
1338 
1339  // Extreme feed-down systematics
1340  valueExtremeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcExtreme->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1341  valueExtremeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcExtreme->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1342 
1343  // Conservative feed-down systematics
1344  valueConservativeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcConservative->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1345  valueConservativeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcConservative->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1346 
1347  }
1348 
1349  }
1350  else { errvalueMax = 0.; errvalueMin = 0.; }
1351 
1352  //
1353  // Fill in the histograms
1354  //
1355  fhYieldCorr->SetBinContent(ibin,value);
1356  fhYieldCorr->SetBinError(ibin,errvalue);
1357  //
1358  // Fill the histos and ntuple vs the Eloss hypothesis
1359  //
1360  if (fPbPbElossHypothesis) {
1361  // Loop over the Eloss hypothesis
1362  for (Float_t rval=0.0025; rval<4.0; rval+=0.005){
1363  Int_t rbin = FindTH2YBin(fhYieldCorrRcb,rval);
1364  Double_t fcRcbvalue = fhFcRcb->GetBinContent(ibin,rbin);
1365  // physics = reco * fcRcb / bin-width
1366  Double_t Rcbvalue = (fhRECpt->GetBinContent(ibin) && fcRcbvalue) ?
1367  fhRECpt->GetBinContent(ibin) * fcRcbvalue : 0. ;
1368  Rcbvalue /= fhRECpt->GetBinWidth(ibin) ;
1369  fhYieldCorrRcb->Fill( fhYieldCorr->GetBinCenter(ibin) , rval, Rcbvalue );
1370  // cout << " pt "<< fhRECpt->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<" fc-fcRbvalue="<<fcRcbvalue<<", yield="<<Rcbvalue<<endl;
1371  }
1372  }
1373  if (fAsymUncertainties) {
1374  Double_t center = fhYieldCorr->GetBinCenter(ibin);
1375  fgYieldCorr->SetPoint(ibin,center,value); // i,x,y
1376  fgYieldCorr->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
1377  fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax);
1378  fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
1379  fgYieldCorrExtreme->SetPoint(ibin,center,value);
1380  fgYieldCorrExtreme->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),value-valueExtremeMin,valueExtremeMax-value);
1381  fgYieldCorrConservative->SetPoint(ibin,center,value);
1382  fgYieldCorrConservative->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),value-valueConservativeMin,valueConservativeMax-value);
1383  }
1384 
1385  }
1386 
1387 }
1388 
1389 //_________________________________________________________________________________________________________
1390 void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumNb(Double_t deltaY, Double_t branchingRatioBintoFinalDecay) {
1404  AliInfo("Calculating the feed-down correction factor and spectrum (Nb method)");
1405 
1406  Double_t value = 0., errvalue = 0., errvalueMax = 0., errvalueMin = 0., kfactor = 0.;
1407  Double_t errvalueExtremeMax = 0., errvalueExtremeMin = 0.;
1408 
1409  // declare the output histograms
1410  fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by Nb)",fnPtBins,fPtBinLimits);
1411  fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by Nb)",fnPtBins,fPtBinLimits);
1412  fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by Nb)",fnPtBins,fPtBinLimits);
1413  if(fPbPbElossHypothesis) {
1414  fhFcRcb = new TH2D("fhFcRcb","fc correction factor (Nb method) vs Rb Eloss hypothesis; p_{T} [GeV/c] ; Rb Eloss hypothesis ; fc correction",fnPtBins,fPtBinLimits,800,0.,4.);
1415  fhYieldCorrRcb = new TH2D("fhYieldCorrRcb","corrected yield (by Nb) vs Rb Eloss hypothesis; p_{T} [GeV/c] ; Rb Eloss hypothesis ; corrected yield",fnPtBins,fPtBinLimits,800,0.,4.);
1416  fnHypothesis = new TNtuple("fnHypothesis"," Feed-down correction vs hypothesis (Nb)","pt:Rb:fc:fcMax:fcMin");
1417  }
1418  // and the output TGraphAsymmErrors
1419  fgYieldCorr = new TGraphAsymmErrors(fnPtBins+1);
1420  if (fAsymUncertainties){
1421  fgYieldCorrExtreme = new TGraphAsymmErrors(fnPtBins+1);
1422  fgYieldCorrConservative = new TGraphAsymmErrors(fnPtBins+1);
1423  // Define fc-conservative
1424  fgFcConservative = new TGraphAsymmErrors(fnPtBins+1);
1425  AliInfo(" Beware the conservative & extreme uncertainties are equal by definition !");
1426  }
1427 
1428  // variables to define fc-conservative
1429  double correction=0, correctionMax=0., correctionMin=0.;
1430 
1431  fhStatUncEffcFD = new TH1D("fhStatUncEffcFD","direct charm stat unc on the feed-down correction",fnPtBins,fPtBinLimits);
1432  fhStatUncEffbFD = new TH1D("fhStatUncEffbFD","secondary charm stat unc on the feed-down correction",fnPtBins,fPtBinLimits);
1433  Double_t correctionUncStatEffc=0.;
1434  Double_t correctionUncStatEffb=0.;
1435 
1436 
1437  //
1438  // Do the calculation
1439  //
1440  for (Int_t ibin=1; ibin<=fnPtBins; ibin++) {
1441 
1442  // Calculate the value
1443  // physics = [ reco - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) ] / bin-width
1444  // In HIC : physics = [ reco - ( Tab * Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th * Rb ) ] / bin-width
1445  //
1446  //
1447  Double_t frac = 1.0, errfrac =0.;
1448 
1449  // Variables initialization
1450  value = 0.; errvalue = 0.; errvalueMax = 0.; errvalueMin = 0.; kfactor = 0.;
1451  errvalueExtremeMax = 0.; errvalueExtremeMin = 0.;
1452  correction=0; correctionMax=0.; correctionMin=0.;
1453  correctionUncStatEffc=0.; correctionUncStatEffb=0.;
1454 
1455  if(fPbPbElossHypothesis || fCollisionType>0) {
1456  frac = fTab[0]*fNevts;
1457  if(fIsEventPlane) frac = frac/2.0;
1458  errfrac = frac * TMath::Sqrt( (fTab[1]/fTab[0])*(fTab[1]/fTab[0]) + (1/fNevts) );
1459  } else {
1460  frac = fLuminosity[0];
1461  errfrac = fLuminosity[1];
1462  }
1463  printf("Tab=%e events=%d frac=%f\n",fTab[0],(Int_t)fNevts,frac);
1464  value = ( fhRECpt->GetBinContent(ibin)>0. && fhRECpt->GetBinContent(ibin)!=0. &&
1465  fhFeedDownMCpt->GetBinContent(ibin)>0. && fhFeedDownEffpt->GetBinContent(ibin)>0. ) ?
1466  fhRECpt->GetBinContent(ibin) - frac*(deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) * fhRECpt->GetBinWidth(ibin) )
1467  : 0. ;
1468  printf("%d raw=%f after=%f \n",ibin,fhRECpt->GetBinContent(ibin),value);
1469  value /= fhRECpt->GetBinWidth(ibin);
1470  if (value<0.) value =0.;
1471 
1472  // Statistical uncertainty: delta_physics = sqrt ( (delta_reco)^2 ) / bin-width
1473  errvalue = (value!=0. && fhRECpt->GetBinError(ibin) && fhRECpt->GetBinError(ibin)!=0.) ?
1474  fhRECpt->GetBinError(ibin) : 0.;
1475  errvalue /= fhRECpt->GetBinWidth(ibin);
1476 
1477  // Correction (fc) : Estimate of the relative amount feed-down subtracted
1478  // correction = [ 1 - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th)/reco ]
1479  // in HIC: correction = [ 1 - ( Tab * Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th)/reco ]
1480  correction = (value>0.) ?
1481  1 - (frac*deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) * fhRECpt->GetBinWidth(ibin) ) / fhRECpt->GetBinContent(ibin) : 0. ;
1482  if (correction<0.) correction = 0.;
1483 
1484  // cout << " pt="<< fhRECpt->GetBinCenter(ibin) << " rec="<< fhRECpt->GetBinContent(ibin) << ", corr="<<correction<<" = 1 - "<< (frac*deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) * fhRECpt->GetBinWidth(ibin) ) / fhRECpt->GetBinContent(ibin) << endl;
1485 
1486  // Systematic uncertainties
1487  // (syst but feed-down) delta_physics = sqrt ( (delta_reco_syst)^2 ) / bin-width
1488  // (feed-down syst) delta_physics = sqrt ( (k*delta_lumi/lumi)^2 + (k*delta_eff_trig/eff_trig)^2
1489  // + (k*delta_Nb/Nb)^2 + (k*delta_eff/eff)^2 + (k*global_eff_ratio)^2 ) / bin-width
1490  // where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th * bin-width
1491  kfactor = frac*deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) * fhRECpt->GetBinWidth(ibin) ;
1492  //
1493  if (fAsymUncertainties && value>0.) {
1494  Double_t nb = fhFeedDownMCpt->GetBinContent(ibin);
1495  Double_t nbDmax = fhFeedDownMCptMax->GetBinContent(ibin) - fhFeedDownMCpt->GetBinContent(ibin);
1496  Double_t nbDmin = fhFeedDownMCpt->GetBinContent(ibin) - fhFeedDownMCptMin->GetBinContent(ibin);
1497 
1498  // Systematics but feed-down
1499  if (fgRECSystematics){
1500  errvalueMax = fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinWidth(ibin) ;
1501  errvalueMin = fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinWidth(ibin);
1502  }
1503  else { errvalueMax = 0.; errvalueMin = 0.; }
1504 
1505  // Feed-down systematics
1506  // min value with the maximum Nb
1507  Double_t errCom = ( (kfactor*errfrac/frac)*(kfactor*errfrac/frac) ) +
1508  ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1509  ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
1510  ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) ) ;
1511  errvalueExtremeMin = TMath::Sqrt( errCom + ( (kfactor*nbDmax/nb)*(kfactor*nbDmax/nb) ) ) / fhRECpt->GetBinWidth(ibin);
1512  // max value with the minimum Nb
1513  errvalueExtremeMax = TMath::Sqrt( errCom + ( (kfactor*nbDmin/nb)*(kfactor*nbDmin/nb) ) ) / fhRECpt->GetBinWidth(ibin);
1514 
1515  // Correction systematics (fc)
1516  // min value with the maximum Nb
1517  correctionMin = TMath::Sqrt( errCom + ( (kfactor*nbDmax/nb)*(kfactor*nbDmax/nb) ) ) / fhRECpt->GetBinContent(ibin) ;
1518  // max value with the minimum Nb
1519  correctionMax = TMath::Sqrt( errCom + ( (kfactor*nbDmin/nb)*(kfactor*nbDmin/nb) ) ) / fhRECpt->GetBinContent(ibin) ;
1520  //
1521  correctionUncStatEffb = TMath::Sqrt( ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) )
1522  ) / fhRECpt->GetBinContent(ibin) ;
1523  correctionUncStatEffc = 0.;
1524  }
1525  else{ // Don't consider Nb uncertainty in this case [ to be tested!!! ]
1526  errvalueExtremeMax = TMath::Sqrt( ( (kfactor*errfrac/frac)*(kfactor*errfrac/frac) ) +
1527  ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1528  ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
1529  ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1530  ) / fhRECpt->GetBinWidth(ibin);
1531  errvalueExtremeMin = errvalueExtremeMax ;
1532  }
1533 
1534 
1535  // fill in histograms
1536  fhYieldCorr->SetBinContent(ibin,value);
1537  fhYieldCorr->SetBinError(ibin,errvalue);
1538  //
1539  // Estimate how the result varies vs charm/beauty Eloss hypothesis
1540  //
1541  if ( correction>1.0e-16 && fPbPbElossHypothesis){
1542  // Loop over the Eloss hypothesis
1543  // Int_t rbin=0;
1544  for (Float_t rval=0.0025; rval<4.0; rval+=0.005){
1545  // correction = [ 1 - (Tab *Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th *binwidth )* (rval) /reco ]
1546  Double_t fcRcbvalue = 1 - (frac*deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin)*fhRECpt->GetBinWidth(ibin) * rval ) / fhRECpt->GetBinContent(ibin) ;
1547  // cout << " rval="<<rval <<" , fc="<<fcRcbvalue<<" ";
1548  if(fcRcbvalue<1.e-16) fcRcbvalue=0.;
1549  fhFcRcb->Fill( fhRECpt->GetBinCenter(ibin) , rval, fcRcbvalue );
1550  // physics = reco * fcRcb / bin-width
1551  Double_t Rcbvalue = (fhRECpt->GetBinContent(ibin) && fcRcbvalue) ?
1552  fhRECpt->GetBinContent(ibin) * fcRcbvalue : 0. ;
1553  Rcbvalue /= fhRECpt->GetBinWidth(ibin) ;
1554  fhYieldCorrRcb->Fill( fhYieldCorr->GetBinCenter(ibin) , rval, Rcbvalue );
1555  // if(ibin==3){
1556  // cout << " pt "<< fhFcRcb->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<", fc-Rb-value="<< fcRcbvalue << ", yield-Rb-value="<< Rcbvalue <<endl;
1557  // cout << " pt "<< fhFcRcb->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", fc-Rb-value="<< fcRcbvalue << ", yield-Rb-value="<< Rcbvalue <<endl;
1558  // rbin++;
1559  // }
1560  Double_t correctionMaxRcb = correctionMax*rval;
1561  Double_t correctionMinRcb = correctionMin*rval;
1562  fnHypothesis->Fill( fhYieldCorr->GetBinCenter(ibin), rval, fcRcbvalue, fcRcbvalue + correctionMaxRcb, fcRcbvalue - correctionMinRcb);
1563 // if(ibin==3){
1564 // cout << " pt="<< fhFcRcb->GetBinCenter(ibin) <<", rval="<<rval<<", fc="<<fcRcbvalue<<" +"<<correctionMaxRcb<<" -"<<correctionMinRcb<<endl;}
1565  }
1566  }
1567  //
1568  // Fill the rest of (asymmetric) histograms
1569  //
1570  if (fAsymUncertainties) {
1571  Double_t x = fhYieldCorr->GetBinCenter(ibin);
1572  fgYieldCorr->SetPoint(ibin,x,value); // i,x,y
1573  fgYieldCorr->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
1574  fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax);
1575  fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
1576  fgYieldCorrExtreme->SetPoint(ibin,x,value); // i,x,y
1577  fgYieldCorrExtreme->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
1578  fgYieldCorrConservative->SetPoint(ibin,x,value); // i,x,y
1579  fgYieldCorrConservative->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
1580  // cout << " bin " << ibin << ", correction " << correction << ", min correction unc " << correctionMin << ", max correction unc " << correctionMax << endl;
1581  if(correction>0.){
1582  fgFcConservative->SetPoint(ibin,x,correction);
1583  fgFcConservative->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),correctionMin,correctionMax);
1584 
1585  fhStatUncEffbFD->SetBinContent(ibin,0.); fhStatUncEffbFD->SetBinError(ibin,correctionUncStatEffb/correction*100.);
1586  fhStatUncEffcFD->SetBinContent(ibin,0.); fhStatUncEffcFD->SetBinError(ibin,correctionUncStatEffc/correction*100.);
1587  // cout << " pt "<< fhStatUncEffcFD->GetBinCenter(ibin) <<" bin "<< ibin<<" fc-stat-c ="<< correctionUncStatEffc/correction <<" fc-stat-b ="<< correctionUncStatEffb/correction <<endl;
1588  }
1589  else{
1590  fgFcConservative->SetPoint(ibin,x,0.);
1591  fgFcConservative->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),0.,0.);
1592  }
1593  }
1594 
1595  }
1596 
1597 }
1598 
1599 
1600 //_________________________________________________________________________________________________________
1601 void AliHFPtSpectrum::ComputeSystUncertainties(AliHFSystErr *systematics, Bool_t combineFeedDown) {
1607 
1608  // Estimate the feed-down uncertainty in percentage
1609  Int_t nentries = 0;
1610  TGraphAsymmErrors *grErrFeeddown = 0;
1611  Double_t x=0., y=0., errx=0., erryl=0., erryh=0;
1612  if(fFeedDownOption!=0) {
1613  nentries = fgSigmaCorrConservative->GetN();
1614  grErrFeeddown = new TGraphAsymmErrors(nentries);
1615  for(Int_t i=0; i<nentries; i++) {
1616  x=0.; y=0.; errx=0.; erryl=0.; erryh=0.;
1617  fgSigmaCorrConservative->GetPoint(i,x,y);
1618  if(y>0.){
1619  errx = fgSigmaCorrConservative->GetErrorXlow(i) ;
1620  erryl = fgSigmaCorrConservative->GetErrorYlow(i) / y ;
1621  erryh = fgSigmaCorrConservative->GetErrorYhigh(i) / y ;
1622  }
1623  // cout << " x "<< x << " +- "<<errx<<" , y "<<y<<" + "<<erryh<<" - "<<erryl<<endl;
1624  grErrFeeddown->SetPoint(i,x,0.);
1625  grErrFeeddown->SetPointError(i,errx,errx,erryl,erryh); //i, xl, xh, yl, yh
1626  }
1627  }
1628 
1629  // Draw all the systematics independently
1630  systematics->DrawErrors(grErrFeeddown);
1631 
1632  // Set the sigma systematic uncertainties
1633  // possibly combine with the feed-down uncertainties
1634  Double_t errylcomb=0., erryhcomb=0;
1635  for(Int_t i=1; i<nentries; i++) {
1636  fgSigmaCorr->GetPoint(i,x,y);
1637  errx = grErrFeeddown->GetErrorXlow(i) ;
1638  erryl = grErrFeeddown->GetErrorYlow(i);
1639  erryh = grErrFeeddown->GetErrorYhigh(i);
1640  if (combineFeedDown) {
1641  errylcomb = systematics->GetTotalSystErr(x,erryl) * y ;
1642  erryhcomb = systematics->GetTotalSystErr(x,erryh) * y ;
1643  } else {
1644  errylcomb = systematics->GetTotalSystErr(x) * y ;
1645  erryhcomb = systematics->GetTotalSystErr(x) * y ;
1646  }
1647  fgSigmaCorr->SetPointError(i,errx,errx,errylcomb,erryhcomb);
1648  //
1649  fhSigmaCorrDataSyst->SetBinContent(i,y);
1650  erryl = systematics->GetTotalSystErr(x) * y ;
1651  fhSigmaCorrDataSyst->SetBinError(i,erryl);
1652  }
1653 
1654 }
1655 
1656 
1657 //_________________________________________________________________________________________________________
1658 void AliHFPtSpectrum::DrawSpectrum(TGraphAsymmErrors *gPrediction) {
1662 
1663  TCanvas *csigma = new TCanvas("csigma","Draw the corrected cross-section & the prediction");
1664  csigma->SetFillColor(0);
1665  gPrediction->GetXaxis()->SetTitleSize(0.05);
1666  gPrediction->GetXaxis()->SetTitleOffset(0.95);
1667  gPrediction->GetYaxis()->SetTitleSize(0.05);
1668  gPrediction->GetYaxis()->SetTitleOffset(0.95);
1669  gPrediction->GetXaxis()->SetTitle("p_{T} [GeV]");
1670  gPrediction->GetYaxis()->SetTitle("BR #times #frac{d#sigma}{dp_{T}} |_{|y|<0.5} [pb/GeV]");
1671  gPrediction->SetLineColor(kGreen+2);
1672  gPrediction->SetLineWidth(3);
1673  gPrediction->SetFillColor(kGreen+1);
1674  gPrediction->Draw("3CA");
1675  fgSigmaCorr->SetLineColor(kRed);
1676  fgSigmaCorr->SetLineWidth(1);
1677  fgSigmaCorr->SetFillColor(kRed);
1678  fgSigmaCorr->SetFillStyle(0);
1679  fgSigmaCorr->Draw("2");
1680  fhSigmaCorr->SetMarkerColor(kRed);
1681  fhSigmaCorr->Draw("esame");
1682  csigma->SetLogy();
1683  TLegend * leg = new TLegend(0.7,0.75,0.87,0.5);
1684  leg->SetBorderSize(0);
1685  leg->SetLineColor(0);
1686  leg->SetFillColor(0);
1687  leg->SetTextFont(42);
1688  leg->AddEntry(gPrediction,"FONLL ","fl");
1689  leg->AddEntry(fhSigmaCorr,"data stat. unc.","pl");
1690  leg->AddEntry(fgSigmaCorr,"data syst. unc.","f");
1691  leg->Draw();
1692  csigma->Draw();
1693 
1694 }
1695 
1696 //_________________________________________________________________________________________________________
1697 TH1D * AliHFPtSpectrum::ReweightHisto(TH1D *hToReweight, TH1D *hReference){
1703 
1704  // check histograms consistency
1705  Bool_t areconsistent=kTRUE;
1706  areconsistent &= CheckHistosConsistency(hToReweight,hReference);
1707  if (!areconsistent) {
1708  AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
1709  return NULL;
1710  }
1711 
1712  // define a new empty histogram
1713  TH1D *hReweighted = (TH1D*)hToReweight->Clone("hReweighted");
1714  hReweighted->Reset();
1715  Double_t weight=1.0;
1716  Double_t yvalue=1.0;
1717  Double_t integralRef = hReference->Integral();
1718  Double_t integralH = hToReweight->Integral();
1719 
1720  // now reweight the spectra
1721  //
1722  // the weight is the relative probability of the given pt bin in the reference histo
1723  // divided by its relative probability (to normalize it) on the histo to re-weight
1724  for (Int_t i=0; i<=hToReweight->GetNbinsX(); i++) {
1725  weight = (hReference->GetBinContent(i)/integralRef) / (hToReweight->GetBinContent(i)/integralH) ;
1726  yvalue = hToReweight->GetBinContent(i);
1727  hReweighted->SetBinContent(i,yvalue*weight);
1728  }
1729 
1730  return (TH1D*)hReweighted;
1731 }
1732 
1733 //_________________________________________________________________________________________________________
1734 TH1D * AliHFPtSpectrum::ReweightRecHisto(TH1D *hRecToReweight, TH1D *hMCToReweight, TH1D *hMCReference){
1740 
1741  // check histograms consistency
1742  Bool_t areconsistent=kTRUE;
1743  areconsistent &= CheckHistosConsistency(hMCToReweight,hMCReference);
1744  areconsistent &= CheckHistosConsistency(hRecToReweight,hMCReference);
1745  if (!areconsistent) {
1746  AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
1747  return NULL;
1748  }
1749 
1750  // define a new empty histogram
1751  TH1D *hReweighted = (TH1D*)hMCToReweight->Clone("hReweighted");
1752  hReweighted->Reset();
1753  TH1D *hRecReweighted = (TH1D*)hRecToReweight->Clone("hRecReweighted");
1754  hRecReweighted->Reset();
1755  Double_t weight=1.0;
1756  Double_t yvalue=1.0, yrecvalue=1.0;
1757  Double_t integralRef = hMCReference->Integral();
1758  Double_t integralH = hMCToReweight->Integral();
1759 
1760  // now reweight the spectra
1761  //
1762  // the weight is the relative probability of the given pt bin
1763  // that should be applied in the MC histo to get the reference histo shape
1764  // Probabilities are properly normalized.
1765  for (Int_t i=0; i<=hMCToReweight->GetNbinsX(); i++) {
1766  weight = (hMCReference->GetBinContent(i)/integralRef) / (hMCToReweight->GetBinContent(i)/integralH) ;
1767  yvalue = hMCToReweight->GetBinContent(i);
1768  hReweighted->SetBinContent(i,yvalue*weight);
1769  yrecvalue = hRecToReweight->GetBinContent(i);
1770  hRecReweighted->SetBinContent(i,yrecvalue*weight);
1771  }
1772 
1773  return (TH1D*)hRecReweighted;
1774 }
1775 
1776 
1777 
1778 //_________________________________________________________________________________________________________
1779 Int_t AliHFPtSpectrum::FindTH2YBin(TH2D *histo, Float_t yvalue){
1783 
1784  Int_t nbins = histo->GetNbinsY();
1785  Int_t ybin=0;
1786  for (int j=0; j<=nbins; j++) {
1787  Float_t value = histo->GetYaxis()->GetBinCenter(j);
1788  Float_t width = histo->GetYaxis()->GetBinWidth(j);
1789  // if( TMath::Abs(yvalue-value)<= width/2. ) {
1790  if( TMath::Abs(yvalue-value)<= width ) {
1791  ybin =j;
1792  // cout <<" value "<<value << ", yval "<< yvalue<<", bin width "<<width/2.<< " y ="<<ybin<<endl;
1793  break;
1794  }
1795  }
1796 
1797  return ybin;
1798 }
1799 
1800 //_________________________________________________________________________________________________________
1801 void AliHFPtSpectrum::ResetStatUncEff(){
1802 
1803  Int_t entries = fhDirectEffpt->GetNbinsX();
1804  for(Int_t i=0; i<=entries; i++){
1805  fhDirectEffpt->SetBinError(i,0.);
1806  fhFeedDownEffpt->SetBinError(i,0.);
1807  }
1808 
1809 }
AliHFPtSpectrum::fhSigmaCorrRcb
TH2D * fhSigmaCorrRcb
Corrected cross-section (syst. unc. from data only)
Definition: AliHFPtSpectrum.h:287
AliHFPtSpectrum::fgYieldCorr
TGraphAsymmErrors * fgYieldCorr
Corrected yield (stat unc. only) vs the Ratio(c/b eloss)
Definition: AliHFPtSpectrum.h:280
AliHFPtSpectrum::SetMCptDistributionsBounds
void SetMCptDistributionsBounds(TH1D *hDirectMax, TH1D *hDirectMin, TH1D *hFeedDownMax, TH1D *hFeedDownMin)
Set the theoretical direct & feeddown pt spectrum upper and lower bounds.
Definition: AliHFPtSpectrum.cxx:420
AliHFPtSpectrum::fhStatUncEffbFD
TH1D * fhStatUncEffbFD
Uncertainty on the feed-down correction due to the prompt efficiency statistical uncertainty.
Definition: AliHFPtSpectrum.h:311
AliHFPtSpectrum::fhStatUncEffcFD
TH1D * fhStatUncEffcFD
Uncertainty on the cross-section due to the feed-down efficiency statistical uncertainty.
Definition: AliHFPtSpectrum.h:310
ClassImp
ClassImp(AliAnalysisTaskTriggerRates) AliAnalysisTaskTriggerRates
Definition: AliAnalysisTaskTriggerRates.cxx:34
AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumFc
void CalculateFeedDownCorrectedSpectrumFc()
Correct the yield for feed-down correction via fc-method.
Definition: AliHFPtSpectrum.cxx:1265
AliHFPtSpectrum::fgYieldCorrExtreme
TGraphAsymmErrors * fgYieldCorrExtreme
Corrected yield as TGraphAsymmErrors (syst but feed-down)
Definition: AliHFPtSpectrum.h:281
title
const char * title
Definition: MakeQAPdf.C:26
AliHFPtSpectrum::DrawSpectrum
void DrawSpectrum(TGraphAsymmErrors *gPrediction)
Drawing the corrected spectrum comparing to theoretical prediction.
Definition: AliHFPtSpectrum.cxx:1658
AliHFPtSpectrum::fPtBinWidths
Double_t * fPtBinWidths
Definition: AliHFPtSpectrum.h:305
AliHFPtSpectrum::fNevts
Int_t fNevts
all reconstructed D Systematic uncertainties
Definition: AliHFPtSpectrum.h:261
AliHFPtSpectrum::fgYieldCorrConservative
TGraphAsymmErrors * fgYieldCorrConservative
Extreme corrected yield as TGraphAsymmErrors (syst from feed-down)
Definition: AliHFPtSpectrum.h:282
AliHFPtSpectrum::CalculateFeedDownCorrectionFc
void CalculateFeedDownCorrectionFc()
Compute the feed-down correction via fc-method.
Definition: AliHFPtSpectrum.cxx:1057
AliHFPtSpectrum::fnHypothesis
TNtuple * fnHypothesis
Ntuple of the calculation vs the Ratio(c/b eloss)
Definition: AliHFPtSpectrum.h:293
AliHFPtSpectrum::SetDirectAccEffCorrection
void SetDirectAccEffCorrection(TH1D *hDirectEff)
Set the acceptance and efficiency corrections for direct.
Definition: AliHFPtSpectrum.cxx:488
AliHFPtSpectrum::EstimateAndSetDirectEfficiencyRecoBin
void EstimateAndSetDirectEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco)
Definition: AliHFPtSpectrum.cxx:880
AliHFPtSpectrum::fhDirectMCptMin
TH1D * fhDirectMCptMin
Input MC maximum c–>D spectra.
Definition: AliHFPtSpectrum.h:251
AliHFPtSpectrum::fhYieldCorrMax
TH1D * fhYieldCorrMax
Corrected yield (stat unc. only)
Definition: AliHFPtSpectrum.h:277
AliHFPtSpectrum::EstimateAndSetFeedDownEfficiencyRecoBin
void EstimateAndSetFeedDownEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco)
Definition: AliHFPtSpectrum.cxx:900
AliHFPtSpectrum::fhSigmaCorrMin
TH1D * fhSigmaCorrMin
Maximum corrected cross-section.
Definition: AliHFPtSpectrum.h:285
AliHFPtSpectrum::fhYieldCorrRcb
TH2D * fhYieldCorrRcb
Minimum corrected yield.
Definition: AliHFPtSpectrum.h:279
AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumNb
void CalculateFeedDownCorrectedSpectrumNb(Double_t deltaY, Double_t branchingRatioBintoFinalDecay)
Correct the yield for feed-down correction via Nb-method.
Definition: AliHFPtSpectrum.cxx:1390
AliHFPtSpectrum::SetFeedDownMCptDistributionsBounds
void SetFeedDownMCptDistributionsBounds(TH1D *hFeedDownMax, TH1D *hFeedDownMin)
Set the theoretical feeddown pt spectrum upper and lower bounds.
Definition: AliHFPtSpectrum.cxx:457
AliHFPtSpectrum::fPbPbElossHypothesis
Bool_t fPbPbElossHypothesis
flag: asymmetric uncertainties are (1) or not (0) considered
Definition: AliHFPtSpectrum.h:298
AliHFPtSpectrum::fAsymUncertainties
Bool_t fAsymUncertainties
feed-down correction flag: 0=none, 1=fc, 2=Nb
Definition: AliHFPtSpectrum.h:297
AliHFPtSpectrum::ReweightRecHisto
TH1D * ReweightRecHisto(TH1D *hRecToReweight, TH1D *hMCToReweight, TH1D *hMCReference)
to reweight the reco-histos: hRecToReweight is reweighted as hReference/hMCToReweight ...
Definition: AliHFPtSpectrum.cxx:1734
AliHFPtSpectrum::fhDirectMCptMax
TH1D * fhDirectMCptMax
Input MC b–>D spectra.
Definition: AliHFPtSpectrum.h:250
AliHFPtSpectrum::fnPtBins
Int_t fnPtBins
flag : when the analysis is done for In/Out of plane, divide the B-prediction by two ...
Definition: AliHFPtSpectrum.h:303
AliHFPtSpectrum::fhFcRcb
TH2D * fhFcRcb
Minimum fc histo.
Definition: AliHFPtSpectrum.h:273
AliHFPtSpectrum::fhFeedDownEffpt
TH1D * fhFeedDownEffpt
c–>D Acceptance and efficiency correction
Definition: AliHFPtSpectrum.h:255
AliHFPtSpectrum::fhDirectEffpt
TH1D * fhDirectEffpt
Input MC minimum b–>D spectra.
Definition: AliHFPtSpectrum.h:254
AliHFPtSpectrum::ResetStatUncEff
void ResetStatUncEff()
Reset stat unc on the efficiencies.
Definition: AliHFPtSpectrum.cxx:1801
AliHFPtSpectrum::Initialize
Bool_t Initialize()
Initialization.
Definition: AliHFPtSpectrum.cxx:920
AliHFPtSpectrum::fhFeedDownMCptMin
TH1D * fhFeedDownMCptMin
Input MC maximum b–>D spectra.
Definition: AliHFPtSpectrum.h:253
AliHFPtSpectrum::ComputeSystUncertainties
void ComputeSystUncertainties(AliHFSystErr *systematics, Bool_t combineFeedDown)
Definition: AliHFPtSpectrum.cxx:1601
AliHFPtSpectrum::fgFcExtreme
TGraphAsymmErrors * fgFcExtreme
Correction histo fc vs the Ratio(c/b eloss)
Definition: AliHFPtSpectrum.h:274
AliHFPtSpectrum::fhSigmaCorrDataSyst
TH1D * fhSigmaCorrDataSyst
Minimum corrected cross-section.
Definition: AliHFPtSpectrum.h:286
AliHFPtSpectrum::fnSigma
TNtuple * fnSigma
Conservative corrected cross-section as TGraphAsymmErrors (syst from feed-down)
Definition: AliHFPtSpectrum.h:292
AliHFPtSpectrum::fFeedDownOption
Int_t fFeedDownOption
0=pp, 1=Pb-Pb, 2=p-Pb
Definition: AliHFPtSpectrum.h:296
AliHFPtSpectrum::fhFeedDownMCpt
TH1D * fhFeedDownMCpt
Input MC c–>D spectra.
Definition: AliHFPtSpectrum.h:249
AliHFPtSpectrum::SetReconstructedSpectrum
void SetReconstructedSpectrum(TH1D *hRec)
Set the reconstructed spectrum.
Definition: AliHFPtSpectrum.cxx:529
AliHFPtSpectrum::fgSigmaCorrConservative
TGraphAsymmErrors * fgSigmaCorrConservative
Extreme corrected cross-section as TGraphAsymmErrors (syst from feed-down)
Definition: AliHFPtSpectrum.h:290
AliHFSystErr::GetTotalSystErr
Double_t GetTotalSystErr(Double_t pt, Double_t feeddownErr=0) const
Definition: AliHFSystErr.cxx:4605
AliHFPtSpectrum::fhFcMin
TH1D * fhFcMin
Maximum fc histo.
Definition: AliHFPtSpectrum.h:272
AliHFPtSpectrum::FindTH2YBin
Int_t FindTH2YBin(TH2D *histo, Float_t yvalue)
Functionality to find the y-axis bin of a TH2 for a given y-value.
Definition: AliHFPtSpectrum.cxx:1779
AliHFPtSpectrum::fhFeedDownMCptMax
TH1D * fhFeedDownMCptMax
Input MC minimum c–>D spectra.
Definition: AliHFPtSpectrum.h:252
AliHFPtSpectrum::SetReconstructedSpectrumSystematics
void SetReconstructedSpectrumSystematics(TGraphAsymmErrors *gRec)
Definition: AliHFPtSpectrum.cxx:559
AliHFPtSpectrum::fhSigmaCorr
TH1D * fhSigmaCorr
Conservative corrected yield as TGraphAsymmErrors (syst from feed-down)
Definition: AliHFPtSpectrum.h:283
AliHFPtSpectrum::fTab
Double_t fTab[2]
uncertainties on the efficiency [0]=c, b, [1]=b/c
Definition: AliHFPtSpectrum.h:265
AliHFPtSpectrum::fgSigmaCorr
TGraphAsymmErrors * fgSigmaCorr
Corrected cross-section (stat unc. only) vs the Ratio(c/b eloss)
Definition: AliHFPtSpectrum.h:288
AliHFPtSpectrum::operator=
AliHFPtSpectrum & operator=(const AliHFPtSpectrum &source)
Assignment operator.
Definition: AliHFPtSpectrum.cxx:205
AliHFPtSpectrum::fhYieldCorrMin
TH1D * fhYieldCorrMin
Maximum corrected yield.
Definition: AliHFPtSpectrum.h:278
AliHFPtSpectrum::fCollisionType
Int_t fCollisionType
Ntuple of the calculation vs the Ratio(c/b eloss)
Definition: AliHFPtSpectrum.h:294
AliHFPtSpectrum::CalculateCorrectedSpectrumNoFeedDown
void CalculateCorrectedSpectrumNoFeedDown()
Definition: AliHFPtSpectrum.cxx:1027
AliHFPtSpectrum::fhYieldCorr
TH1D * fhYieldCorr
Extreme correction as TGraphAsymmErrors.
Definition: AliHFPtSpectrum.h:276
AliHFPtSpectrum::SetFeedDownMCptSpectra
void SetFeedDownMCptSpectra(TH1D *hFeedDown)
Set the theoretical feeddown pt spectrum.
Definition: AliHFPtSpectrum.cxx:400
AliHFPtSpectrum::~AliHFPtSpectrum
virtual ~AliHFPtSpectrum()
Destructor.
Definition: AliHFPtSpectrum.cxx:276
AliHFPtSpectrum::RebinTheoreticalSpectra
TH1D * RebinTheoreticalSpectra(TH1D *hTheory, const char *name)
Function to rebin the theoretical spectra in the data-reconstructed spectra binning.
Definition: AliHFPtSpectrum.cxx:318
AliHFPtSpectrum::fgRECSystematics
TGraphAsymmErrors * fgRECSystematics
all reconstructed D
Definition: AliHFPtSpectrum.h:258
AliHFPtSpectrum::CheckHistosConsistency
Bool_t CheckHistosConsistency(TH1D *h1, TH1D *h2)
Check histograms consistency function.
Definition: AliHFPtSpectrum.cxx:993
AliHFPtSpectrum::fhStatUncEffbSigma
TH1D * fhStatUncEffbSigma
Uncertainty on the cross-section due to the prompt efficiency statistical uncertainty.
Definition: AliHFPtSpectrum.h:309
AliHFPtSpectrum::ComputeHFPtSpectrum
void ComputeHFPtSpectrum(Double_t deltaY=1.0, Double_t branchingRatioC=1.0, Double_t branchingRatioBintoFinalDecay=1.0)
Definition: AliHFPtSpectrum.cxx:579
AliHFPtSpectrum::fgSigmaCorrExtreme
TGraphAsymmErrors * fgSigmaCorrExtreme
Corrected cross-section as TGraphAsymmErrors (syst but feed-down)
Definition: AliHFPtSpectrum.h:289
AliHFPtSpectrum::fTrigEfficiency
Double_t fTrigEfficiency[2]
analyzed luminosity & uncertainty
Definition: AliHFPtSpectrum.h:263
AliHFPtSpectrum::SetAccEffCorrection
void SetAccEffCorrection(TH1D *hDirectEff, TH1D *hFeedDownEff)
Set the acceptance and efficiency corrections for direct & feeddown.
Definition: AliHFPtSpectrum.cxx:504
AliHFPtSpectrum::fhRECpt
TH1D * fhRECpt
b–>D Acceptance and efficiency correction
Definition: AliHFPtSpectrum.h:256
AliHFPtSpectrum::ReweightHisto
TH1D * ReweightHisto(TH1D *hToReweight, TH1D *hReference)
Definition: AliHFPtSpectrum.cxx:1697
AliHFSystErr.h
AliHFPtSpectrum::fhFcMax
TH1D * fhFcMax
Correction histo fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
Definition: AliHFPtSpectrum.h:271
AliHFPtSpectrum::fGlobalEfficiencyUncertainties
Double_t fGlobalEfficiencyUncertainties[2]
trigger efficiency & uncertainty
Definition: AliHFPtSpectrum.h:264
nbins
const Int_t nbins
Definition: AverageDmesonRaa.C:63
AliHFSystErr::DrawErrors
void DrawErrors(TGraphAsymmErrors *grErrFeeddown=0) const
Definition: AliHFSystErr.cxx:4625
AliHFPtSpectrum::fgFcConservative
TGraphAsymmErrors * fgFcConservative
Extreme correction as TGraphAsymmErrors.
Definition: AliHFPtSpectrum.h:275
AliHFPtSpectrum::fParticleAntiParticle
Int_t fParticleAntiParticle
flag : consider (1) or not (0) the stat unc on the efficiencies
Definition: AliHFPtSpectrum.h:300
AliHFPtSpectrum::AliHFPtSpectrum
AliHFPtSpectrum(const char *name="AliHFPtSpectrum", const char *title="HF feed down correction class", Int_t option=1)
Constructor.
Definition: AliHFPtSpectrum.cxx:61
AliHFPtSpectrum::EstimateEfficiencyRecoBin
TH1D * EstimateEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco, const char *name)
Function to estimate the efficiency in the data-reconstructed spectra binning.
Definition: AliHFPtSpectrum.cxx:821
AliHFPtSpectrum::fLuminosity
Double_t fLuminosity[2]
nb of analyzed events
Definition: AliHFPtSpectrum.h:262
AliHFPtSpectrum::SetMCptSpectra
void SetMCptSpectra(TH1D *hDirect, TH1D *hFeedDown)
Definition: AliHFPtSpectrum.cxx:371
AliHFPtSpectrum::fhFc
TH1D * fhFc
Tab parameter and its uncertainty.
Definition: AliHFPtSpectrum.h:270
AliHFPtSpectrum.h
AliHFPtSpectrum::fPtBinLimits
Double_t * fPtBinLimits
number of pt bins
Definition: AliHFPtSpectrum.h:304
AliHFPtSpectrum::fIsStatUncEff
Bool_t fIsStatUncEff
flag: whether to do estimates vs Ratio(c/b eloss) hypothesis
Definition: AliHFPtSpectrum.h:299
AliHFPtSpectrum::fhSigmaCorrMax
TH1D * fhSigmaCorrMax
Corrected cross-section (stat unc. only)
Definition: AliHFPtSpectrum.h:284
AliHFPtSpectrum::fIsEventPlane
Bool_t fIsEventPlane
1: only one sign, 2: yield is for particle+anti-particle
Definition: AliHFPtSpectrum.h:301