AliPhysics  e6d2b2b (e6d2b2b)
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 
55 using std::cout;
56 using std::endl;
57 
59 ClassImp(AliHFPtSpectrum);
61 
62 
63 //_________________________________________________________________________________________________________
64 AliHFPtSpectrum::AliHFPtSpectrum(const char* name, const char* title, Int_t option):
65  TNamed(name,title),
66  fhDirectMCpt(NULL),
67  fhFeedDownMCpt(NULL),
68  fhDirectMCptMax(NULL),
69  fhDirectMCptMin(NULL),
70  fhFeedDownMCptMax(NULL),
71  fhFeedDownMCptMin(NULL),
72  fhDirectEffpt(NULL),
73  fhFeedDownEffpt(NULL),
74  fhRECpt(NULL),
75  fgRECSystematics(NULL),
76  fNevts(1),
77  fLuminosity(),
78  fTrigEfficiency(),
79  fGlobalEfficiencyUncertainties(),
80  fGlobalEfficiencyPtDependent(false),
81  fTab(),
82  fSystematics(NULL),
83  fhFc(NULL),
84  fhFcMax(NULL),
85  fhFcMin(NULL),
86  fhFcRcb(NULL),
87  fgFcExtreme(NULL),
88  fgFcConservative(NULL),
89  fhYieldCorr(NULL),
90  fhYieldCorrMax(NULL),
91  fhYieldCorrMin(NULL),
92  fhYieldCorrRcb(NULL),
93  fgYieldCorr(NULL),
94  fgYieldCorrExtreme(NULL),
95  fgYieldCorrConservative(NULL),
96  fhSigmaCorr(NULL),
97  fhSigmaCorrMax(NULL),
98  fhSigmaCorrMin(NULL),
99  fhSigmaCorrDataSyst(NULL),
100  fhSigmaCorrRcb(NULL),
101  fgSigmaCorr(NULL),
102  fgSigmaCorrExtreme(NULL),
103  fgSigmaCorrConservative(NULL),
104  fnSigma(NULL),
105  fnHypothesis(NULL),
106  fCollisionType(0),
107  fFeedDownOption(option),
108  fAsymUncertainties(kTRUE),
109  fPbPbElossHypothesis(kFALSE),
110  fIsStatUncEff(kTRUE),
111  fParticleAntiParticle(2),
112  fIsEventPlane(kFALSE),
113  fnPtBins(0),
114  fPtBinLimits(NULL),
115  fPtBinWidths(NULL),
116  fhStatUncEffcSigma(NULL),
117  fhStatUncEffbSigma(NULL),
118  fhStatUncEffcFD(NULL),
119  fhStatUncEffbFD(NULL)
120 {
121  //
123  //
124 
125  fLuminosity[0]=1.; fLuminosity[1]=0.;
126  fTrigEfficiency[0]=1.; fTrigEfficiency[1]=0.;
128  fTab[0]=1.; fTab[1]=0.;
129 
130 }
131 
132 //_________________________________________________________________________________________________________
134  TNamed(rhs),
143  fhRECpt(rhs.fhRECpt),
145  fNevts(rhs.fNevts),
146  fLuminosity(),
147  fTrigEfficiency(),
150  fTab(),
152  fhFc(rhs.fhFc),
153  fhFcMax(rhs.fhFcMax),
154  fhFcMin(rhs.fhFcMin),
155  fhFcRcb(rhs.fhFcRcb),
173  fnSigma(rhs.fnSigma),
182  fnPtBins(rhs.fnPtBins),
183  fPtBinLimits(NULL),
184  fPtBinWidths(NULL),
185  fhStatUncEffcSigma(NULL),
186  fhStatUncEffbSigma(NULL),
187  fhStatUncEffcFD(NULL),
188  fhStatUncEffbFD(NULL)
189 {
190  //
192  //
193 
194  for(Int_t i=0; i<2; i++){
195  fLuminosity[i] = rhs.fLuminosity[i];
196  fTrigEfficiency[i] = rhs.fTrigEfficiency[i];
198  fTab[i] = rhs.fTab[i];
199  }
200 
201  fPtBinLimits = new Double_t[fnPtBins+1];
203  for(Int_t i=0; i<fnPtBins; i++){
204  fPtBinLimits[i] = rhs.fPtBinLimits[i];
205  fPtBinWidths[i] = rhs.fPtBinWidths[i];
206  }
208 
209 }
210 
211 //_________________________________________________________________________________________________________
213  //
215  //
216 
217  if (&source == this) return *this;
218 
219  fhDirectMCpt = source.fhDirectMCpt;
225  fhDirectEffpt = source.fhDirectEffpt;
227  fhRECpt = source.fhRECpt;
229  fNevts = source.fNevts;
231  fSystematics = source.fSystematics;
232  fhFc = source.fhFc;
233  fhFcMax = source.fhFcMax;
234  fhFcMin = source.fhFcMin;
235  fhFcRcb = source.fhFcRcb;
236  fgFcExtreme = source.fgFcExtreme;
238  fhYieldCorr = source.fhYieldCorr;
242  fgYieldCorr = source.fgYieldCorr;
245  fhSigmaCorr = source.fhSigmaCorr;
250  fgSigmaCorr = source.fgSigmaCorr;
253  fnSigma = source.fnSigma;
254  fnHypothesis = source.fnHypothesis;
259  fIsStatUncEff = source.fIsStatUncEff;
261  fIsEventPlane = source.fIsEventPlane;
262 
263  for(Int_t i=0; i<2; i++){
264  fLuminosity[i] = source.fLuminosity[i];
265  fTrigEfficiency[i] = source.fTrigEfficiency[i];
267  fTab[i] = source.fTab[i];
268  }
269 
270  fnPtBins = source.fnPtBins;
271  if(fPtBinLimits) delete fPtBinLimits;
272  if(fPtBinWidths) delete fPtBinWidths;
273  fPtBinLimits = new Double_t[fnPtBins+1];
275  for(Int_t i=0; i<fnPtBins; i++){
276  fPtBinLimits[i] = source.fPtBinLimits[i];
277  fPtBinWidths[i] = source.fPtBinWidths[i];
278  }
280 
281  return *this;
282 }
283 
284 //_________________________________________________________________________________________________________
286  //
288  //
289  if (fhDirectMCpt) delete fhDirectMCpt;
290  if (fhFeedDownMCpt) delete fhFeedDownMCpt;
291  if (fhDirectMCptMax) delete fhDirectMCptMax;
292  if (fhDirectMCptMin) delete fhDirectMCptMin;
295  if (fhDirectEffpt) delete fhDirectEffpt;
296  if (fhFeedDownEffpt) delete fhFeedDownEffpt;
297  if (fhRECpt) delete fhRECpt;
299  if (fhFc) delete fhFc;
300  if (fhFcMax) delete fhFcMax;
301  if (fhFcMin) delete fhFcMin;
302  if (fhFcRcb) delete fhFcRcb;
303  if (fgFcExtreme) delete fgFcExtreme;
305  if (fhYieldCorr) delete fhYieldCorr;
306  if (fhYieldCorrMax) delete fhYieldCorrMax;
307  if (fhYieldCorrMin) delete fhYieldCorrMin;
308  if (fhYieldCorrRcb) delete fhYieldCorrRcb;
309  if (fgYieldCorr) delete fgYieldCorr;
312  if (fhSigmaCorr) delete fhSigmaCorr;
313  if (fhSigmaCorrMax) delete fhSigmaCorrMax;
314  if (fhSigmaCorrMin) delete fhSigmaCorrMin;
316  if (fgSigmaCorr) delete fgSigmaCorr;
319  if (fnSigma) delete fnSigma;
320  if (fnHypothesis) delete fnHypothesis;
321  if (fPtBinLimits) delete [] fPtBinLimits;
322  if (fPtBinWidths) delete [] fPtBinWidths;
323 }
324 
325 
326 //_________________________________________________________________________________________________________
327 TH1D * AliHFPtSpectrum::RebinTheoreticalSpectra(TH1D *hTheory, const char *name) {
332 
333  if (!hTheory || !fhRECpt) {
334  AliError("Feed-down or reconstructed spectra don't exist");
335  return NULL;
336  }
337 
338  //
339  // Get the reconstructed spectra bins & limits
340  Int_t nbinsMC = hTheory->GetNbinsX();
341 
342  // Check that the reconstructed spectra binning
343  // is larger than the theoretical one
344  Double_t thbinwidth = hTheory->GetBinWidth(1);
345  for (Int_t i=1; i<=fnPtBins; i++) {
346  if ( thbinwidth > fPtBinWidths[i-1] ) {
347  AliInfo(" Beware it seems that the reconstructed spectra has a smaller binning than the theoretical predictions !! ");
348  }
349  }
350 
351  //
352  // Define a new histogram with the real-data reconstructed spectrum binning
353  TH1D * hTheoryRebin = new TH1D(name," theoretical rebinned prediction",fnPtBins,fPtBinLimits);
354 
355  Double_t sum[fnPtBins], items[fnPtBins];
356  for (Int_t ibin=0; ibin<fnPtBins; ibin++) {
357  sum[ibin]=0.; items[ibin]=0.;
358  }
359  for (Int_t ibin=0; ibin<=nbinsMC; ibin++){
360 
361  for (Int_t ibinrec=0; ibinrec<fnPtBins; ibinrec++){
362  if (hTheory->GetBinCenter(ibin)>fPtBinLimits[ibinrec] &&
363  hTheory->GetBinCenter(ibin)<fPtBinLimits[ibinrec+1]){
364  sum[ibinrec]+=hTheory->GetBinContent(ibin);
365  items[ibinrec]+=1.;
366  }
367  }
368 
369  }
370 
371  // set the theoretical rebinned spectra to ( sum-bins / n-bins ) per new bin
372  for (Int_t ibinrec=0; ibinrec<fnPtBins; ibinrec++) {
373  hTheoryRebin->SetBinContent(ibinrec+1,sum[ibinrec]/items[ibinrec]);
374  }
375 
376  return (TH1D*)hTheoryRebin;
377 }
378 
379 //_________________________________________________________________________________________________________
380 void AliHFPtSpectrum::SetMCptSpectra(TH1D *hDirect, TH1D *hFeedDown){
385 
386  if (!hDirect || !hFeedDown || !fhRECpt) {
387  AliError("One or both (direct, feed-down) spectra or the reconstructed spectra don't exist");
388  return;
389  }
390 
391  Bool_t areconsistent = kTRUE;
392  areconsistent = CheckHistosConsistency(hDirect,hFeedDown);
393  if (!areconsistent) {
394  AliInfo("Histograms are not consistent (bin width, bounds)");
395  return;
396  }
397 
398  //
399  // Rebin the theoretical predictions to the reconstructed spectra binning
400  //
401  fhDirectMCpt = RebinTheoreticalSpectra(hDirect,"fhDirectMCpt");
402  fhDirectMCpt->SetNameTitle("fhDirectMCpt"," direct theoretical prediction");
403  fhFeedDownMCpt = RebinTheoreticalSpectra(hFeedDown,"fhFeedDownMCpt");
404  fhFeedDownMCpt->SetNameTitle("fhFeedDownMCpt"," feed-down theoretical prediction");
405 
406 }
407 
408 //_________________________________________________________________________________________________________
414 
415  if (!hFeedDown || !fhRECpt) {
416  AliError("Feed-down or reconstructed spectra don't exist");
417  return;
418  }
419 
420  //
421  // Rebin the theoretical predictions to the reconstructed spectra binning
422  //
423  fhFeedDownMCpt = RebinTheoreticalSpectra(hFeedDown,"fhFeedDownMCpt");
424  fhFeedDownMCpt->SetNameTitle("fhFeedDownMCpt"," feed-down theoretical prediction");
425 
426 }
427 
428 //_________________________________________________________________________________________________________
429 void AliHFPtSpectrum::SetMCptDistributionsBounds(TH1D *hDirectMax, TH1D *hDirectMin, TH1D *hFeedDownMax, TH1D *hFeedDownMin){
435 
436  if (!hDirectMax || !hDirectMin || !hFeedDownMax|| !hFeedDownMin || !fhRECpt) {
437  AliError("One or all of the max/min direct/feed-down or the reconstructed spectra don't exist");
438  return;
439  }
440 
441  Bool_t areconsistent = kTRUE;
442  areconsistent &= CheckHistosConsistency(hDirectMax,hDirectMin);
443  areconsistent &= CheckHistosConsistency(hFeedDownMax,hFeedDownMin);
444  areconsistent &= CheckHistosConsistency(hDirectMax,hFeedDownMax);
445  if (!areconsistent) {
446  AliInfo("Histograms are not consistent (bin width, bounds)");
447  return;
448  }
449 
450 
451  //
452  // Rebin the theoretical predictions to the reconstructed spectra binning
453  //
454  fhDirectMCptMax = RebinTheoreticalSpectra(hDirectMax,"fhDirectMCptMax");
455  fhDirectMCptMax->SetNameTitle("fhDirectMCptMax"," maximum direct theoretical prediction");
456  fhDirectMCptMin = RebinTheoreticalSpectra(hDirectMin,"fhDirectMCptMin");
457  fhDirectMCptMin->SetNameTitle("fhDirectMCptMin"," minimum direct theoretical prediction");
458  fhFeedDownMCptMax = RebinTheoreticalSpectra(hFeedDownMax,"fhFeedDownMCptMax");
459  fhFeedDownMCptMax->SetNameTitle("fhFeedDownMCptMax"," maximum feed-down theoretical prediction");
460  fhFeedDownMCptMin = RebinTheoreticalSpectra(hFeedDownMin,"fhFeedDownMCptMin");
461  fhFeedDownMCptMin->SetNameTitle("fhFeedDownMCptMin"," minimum feed-down theoretical prediction");
462 
463 }
464 
465 //_________________________________________________________________________________________________________
472 
473  if (!hFeedDownMax || !hFeedDownMin || !fhRECpt) {
474  AliError("One or all of the max/min direct/feed-down spectra don't exist");
475  return;
476  }
477 
478  Bool_t areconsistent = kTRUE;
479  areconsistent &= CheckHistosConsistency(hFeedDownMax,hFeedDownMin);
480  if (!areconsistent) {
481  AliInfo("Histograms are not consistent (bin width, bounds)");
482  return;
483  }
484 
485 
486  //
487  // Rebin the theoretical predictions to the reconstructed spectra binning
488  //
489  fhFeedDownMCptMax = RebinTheoreticalSpectra(hFeedDownMax,"fhFeedDownMCptMax");
490  fhFeedDownMCptMax->SetNameTitle("fhFeedDownMCptMax"," maximum feed-down theoretical prediction");
491  fhFeedDownMCptMin = RebinTheoreticalSpectra(hFeedDownMin,"fhFeedDownMCptMin");
492  fhFeedDownMCptMin->SetNameTitle("fhFeedDownMCptMin"," minimum feed-down theoretical prediction");
493 
494 }
495 
496 //_________________________________________________________________________________________________________
502 
503  if (!hDirectEff) {
504  AliError("The direct acceptance and efficiency corrections doesn't exist");
505  return;
506  }
507 
508  fhDirectEffpt = (TH1D*)hDirectEff->Clone();
509  fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance x efficiency correction");
510 }
511 
512 //_________________________________________________________________________________________________________
513 void AliHFPtSpectrum::SetAccEffCorrection(TH1D *hDirectEff, TH1D *hFeedDownEff){
518 
519  if (!hDirectEff || !hFeedDownEff) {
520  AliError("One or both (direct, feed-down) acceptance and efficiency corrections don't exist");
521  return;
522  }
523 
524  Bool_t areconsistent=kTRUE;
525  areconsistent = CheckHistosConsistency(hDirectEff,hFeedDownEff);
526  if (!areconsistent) {
527  AliInfo("Histograms are not consistent (bin width, bounds)");
528  return;
529  }
530 
531  fhDirectEffpt = (TH1D*)hDirectEff->Clone();
532  fhFeedDownEffpt = (TH1D*)hFeedDownEff->Clone();
533  fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance x efficiency correction");
534  fhFeedDownEffpt->SetNameTitle("fhFeedDownEffpt"," feed-down acceptance x efficiency correction");
535 }
536 
537 //_________________________________________________________________________________________________________
542 
543  if (!hRec) {
544  AliError("The reconstructed spectrum doesn't exist");
545  return;
546  }
547 
548  fhRECpt = (TH1D*)hRec->Clone();
549  fhRECpt->SetNameTitle("fhRECpt"," reconstructed spectrum");
550 
551  //
552  // Evaluate the number of intervals and limits of the spectrum
553  //
554  fnPtBins = fhRECpt->GetNbinsX();
555  fPtBinLimits = new Double_t[fnPtBins+1];
557  Double_t xlow=0., binwidth=0.;
558  for (Int_t i=1; i<=fnPtBins; i++) {
559  binwidth = fhRECpt->GetBinWidth(i);
560  xlow = fhRECpt->GetBinLowEdge(i);
561  fPtBinLimits[i-1] = xlow;
562  fPtBinWidths[i-1] = binwidth;
563  }
564  fPtBinLimits[fnPtBins] = xlow + binwidth;
565 }
566 
567 //_________________________________________________________________________________________________________
572 
573  // Check the compatibility with the reconstructed spectrum
574  Double_t gbinwidth = gRec->GetErrorXlow(1) + gRec->GetErrorXhigh(1) ;
575  Double_t hbinwidth = fhRECpt->GetBinWidth(1);
576  Double_t gxbincenter=0., gybincenter=0.;
577  gRec->GetPoint(1,gxbincenter,gybincenter);
578  Double_t hbincenter = fhRECpt->GetBinCenter(1);
579  if ( (gbinwidth != hbinwidth) || (gxbincenter!=hbincenter) ) {
580  AliError(" The reconstructed spectrum and its systematics don't seem compatible");
581  return;
582  }
583 
584  fgRECSystematics = gRec;
585 }
586 
587 //_________________________________________________________________________________________________________
588 void AliHFPtSpectrum::ComputeHFPtSpectrum(Double_t deltaY, Double_t branchingRatioC, Double_t branchingRatioBintoFinalDecay) {
602  //
603 
604  //
605  // First: Initialization
606  //
607  Bool_t areHistosOk = Initialize();
608  if (!areHistosOk) {
609  AliInfo(" Histos not properly initialized. Check : inconsistent binning ? missing histos ?");
610  return;
611  }
612  // Reset the statistical uncertainties on the efficiencies if needed
614 
615  //
616  // Second: Correct for feed-down
617  //
618  if (fFeedDownOption==1) {
619  // Compute the feed-down correction via fc-method
621  // Correct the yield for feed-down correction via fc-method
623  }
624  else if (fFeedDownOption==2) {
625  // Correct the yield for feed-down correction via Nb-method
626  CalculateFeedDownCorrectedSpectrumNb(deltaY,branchingRatioBintoFinalDecay);
627  }
628  else if (fFeedDownOption==0) {
629  // If there is no need for feed-down correction,
630  // the "corrected" yield is equal to the raw yield
632  }
633  else {
634  AliInfo(" Are you sure the feed-down correction option is right ?");
635  }
636 
637 
638  // Print out information
639  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 usage of pt-dependent efficiency uncertainty for Nb uncertainy calculation? %1.0d \n\n",fLuminosity[0],fLuminosity[1],fTrigEfficiency[0],fTrigEfficiency[1],deltaY,branchingRatioC,branchingRatioBintoFinalDecay,fGlobalEfficiencyUncertainties[0],fGlobalEfficiencyUncertainties[1],fGlobalEfficiencyPtDependent);
640  if (fPbPbElossHypothesis) printf("\n\n The considered Tab is %4.2e +- %2.2e \n\n",fTab[0],fTab[1]);
641 
642  //
643  // Finally: Correct from yields to cross-section
644  //
645 
646  // declare the output histograms
647  fhSigmaCorr = new TH1D("fhSigmaCorr","corrected sigma",fnPtBins,fPtBinLimits);
649 
650  fhSigmaCorrMax = new TH1D("fhSigmaCorrMax","max corrected sigma",fnPtBins,fPtBinLimits);
651  fhSigmaCorrMin = new TH1D("fhSigmaCorrMin","min corrected sigma",fnPtBins,fPtBinLimits);
652  fhSigmaCorrDataSyst = new TH1D("fhSigmaCorrDataSyst","data syst uncertainties on the corrected sigma",fnPtBins,fPtBinLimits);
653 
655  fhSigmaCorrRcb = new TH2D("fhSigmaCorrRcb","corrected sigma vs Rcb Eloss hypothesis; p_{T} [GeV/c] ; Rcb Eloss hypothesis ; #sigma",fnPtBins,fPtBinLimits,800,0.,4.);
656  fnSigma = new TNtuple("fnSigma"," Sigma ntuple calculation","pt:Signal:Rcb:fc:Yield:Sigma:SigmaStatUnc:SigmaMax:SigmaMin");
657  }
659  fhSigmaCorrRcb = new TH2D("fhSigmaCorrRcb","corrected sigma vs Rb Eloss hypothesis; p_{T} [GeV/c] ; Rb Eloss hypothesis ; #sigma",fnPtBins,fPtBinLimits,800,0.,4.);
660  fnSigma = new TNtuple("fnSigma"," Sigma ntuple calculation","pt:Signal:Rb:fc:Yield:Sigma:SigmaStatUnc:SigmaMax:SigmaMin");
661  }
662  // and the output TGraphAsymmErrors
663  if (fAsymUncertainties){
666  }
667  fhStatUncEffcSigma = new TH1D("fhStatUncEffcSigma","direct charm stat unc on the cross section",fnPtBins,fPtBinLimits);
668  fhStatUncEffbSigma = new TH1D("fhStatUncEffbSigma","secondary charm stat unc on the cross section",fnPtBins,fPtBinLimits);
669 
670 
671  // protect against null denominator
672  if (deltaY==0. || fLuminosity[0]==0. || fTrigEfficiency[0]==0. || branchingRatioC==0.) {
673  AliError(" Hey you ! Why luminosity or trigger-efficiency or the c-BR or delta_y are set to zero ?! ");
674  return ;
675  }
676 
677  Double_t value=0, errValue=0, errvalueMax=0., errvalueMin=0.;
678  Double_t errvalueExtremeMax=0., errvalueExtremeMin=0.;
679  Double_t errvalueConservativeMax=0., errvalueConservativeMin=0.;
680  Double_t errvalueStatUncEffc=0.;
681  for(Int_t ibin=1; ibin<=fnPtBins; ibin++){
682 
683  // Variables initialization
684  value=0.; errValue=0.; errvalueMax=0.; errvalueMin=0.;
685  errvalueExtremeMax=0.; errvalueExtremeMin=0.;
686  errvalueConservativeMax=0.; errvalueConservativeMin=0.;
687  errvalueStatUncEffc=0.;
688 
689  // Sigma calculation
690  // Sigma = ( 1. / (lumi * delta_y * BR_c * ParticleAntiPartFactor * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
691  value = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0. && fhRECpt->GetBinContent(ibin)>0.) ?
692  ( fhYieldCorr->GetBinContent(ibin) / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
693  : 0. ;
694 
695  // Sigma statistical uncertainty:
696  // delta_sigma = sigma * sqrt ( (delta_spectra/spectra)^2 )
697  errValue = (value!=0.) ? value * (fhYieldCorr->GetBinError(ibin)/fhYieldCorr->GetBinContent(ibin)) : 0. ;
698 
699  // cout<< " x "<< fhRECpt->GetBinCenter(ibin) << " sigma " << value << " +- "<< errValue << " (stat)"<<endl;
700 
701  //
702  // Sigma systematic uncertainties
703  //
704  if (fAsymUncertainties && value>0.) {
705 
706  // (syst but feed-down) delta_sigma = sigma * sqrt ( (delta_spectra_syst/spectra)^2 +
707  // (delta_lumi/lumi)^2 + (delta_eff_trig/eff_trig)^2 + (delta_eff/eff)^2 + (global_eff)^2 )
708  errvalueMax = value * TMath::Sqrt( (fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin))*(fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin)) +
711  (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
713  errvalueMin = value * TMath::Sqrt( (fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin))*(fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin)) +
716  (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
717  fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] );
718 
719  // Uncertainties from feed-down
720  // (feed-down syst) delta_sigma = sigma * sqrt ( (delta_spectra_fd/spectra_fd)^2 )
721  // extreme case
722  errvalueExtremeMax = value * (fgYieldCorrExtreme->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin));
723  errvalueExtremeMin = value * (fgYieldCorrExtreme->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin));
724  //
725  // conservative case
726  errvalueConservativeMax = value * (fgYieldCorrConservative->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin));
727  errvalueConservativeMin = value * (fgYieldCorrConservative->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin));
728 
729 
730  // stat unc of the efficiencies, separately
731  errvalueStatUncEffc = value * (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) ;
732 
733  }
734  else {
735  // protect against null denominator
736  errvalueMax = (value!=0.) ?
737  value * TMath::Sqrt( (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
739  (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
741  : 0. ;
742  errvalueMin = errvalueMax;
743  }
744  // if(fFeedDownOption==0 && fSystematics) {
745  // errvalueMax = value * fSystematics->GetTotalSystErr(fhRECpt->GetBinCenter(ibin));
746  // errvalueMin = errvalueMax;
747  // }
748  //
749  // Fill the histograms
750  //
751  fhSigmaCorr->SetBinContent(ibin,value);
752  fhSigmaCorr->SetBinError(ibin,errValue);
753 
754  Double_t x = fhYieldCorr->GetBinCenter(ibin);
755  fgSigmaCorr->SetPoint(ibin,x,value); // i,x,y
756  fgSigmaCorr->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
757  fhSigmaCorrMax->SetBinContent(ibin,value+errvalueMax);
758  fhSigmaCorrMin->SetBinContent(ibin,value-errvalueMin);
759 
760  //
761  // Fill the histos and ntuple vs the Eloss hypothesis
762  //
763  if (fPbPbElossHypothesis) {
764 
765  // Loop over the Eloss hypothesis
766  if(!fnHypothesis) {
767  AliError("Error reading the fc hypothesis no ntuple, please check !!");
768  return ;
769  }
770  Int_t entriesHypo = fnHypothesis->GetEntries();
771  Float_t pt=0., Rhy=0., fc=0., fcMin=0., fcMax=0.;
772  fnHypothesis->SetBranchAddress("pt",&pt);
773  if (fFeedDownOption==2) fnHypothesis->SetBranchAddress("Rb",&Rhy);
774  else if (fFeedDownOption==1) fnHypothesis->SetBranchAddress("Rcb",&Rhy);
775  fnHypothesis->SetBranchAddress("fc",&fc);
776  fnHypothesis->SetBranchAddress("fcMax",&fcMax);
777  fnHypothesis->SetBranchAddress("fcMin",&fcMin);
778 
779  for (Int_t item=0; item<entriesHypo; item++) {
780 
781  fnHypothesis->GetEntry(item);
782  if ( TMath::Abs( pt - fhDirectEffpt->GetBinCenter(ibin) ) > 0.15 ) continue;
783 
784  Double_t yieldRcbvalue = (fhRECpt->GetBinContent(ibin) ) ? fhRECpt->GetBinContent(ibin) * fc : 0. ;
785  yieldRcbvalue /= fhRECpt->GetBinWidth(ibin) ;
786  Double_t yieldRcbvalueMax = (fhRECpt->GetBinContent(ibin) ) ? fhRECpt->GetBinContent(ibin) * fcMax : 0. ;
787  yieldRcbvalueMax /= fhRECpt->GetBinWidth(ibin) ;
788  Double_t yieldRcbvalueMin = (fhRECpt->GetBinContent(ibin) ) ? fhRECpt->GetBinContent(ibin) * fcMin : 0. ;
789  yieldRcbvalueMin /= fhRECpt->GetBinWidth(ibin) ;
790 
791  // Sigma calculation
792  // Sigma = ( 1. / (lumi * delta_y * BR_c * ParticleAntiPartFactor * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
793  Double_t sigmaRcbvalue = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)>0.) ?
794  ( yieldRcbvalue / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
795  : 0. ;
796  Double_t sigmaRcbvalueMax = (sigmaRcbvalue!=0.) ?
797  ( yieldRcbvalueMax / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
798  : 0. ;
799  Double_t sigmaRcbvalueMin = (sigmaRcbvalue!=0.) ?
800  ( yieldRcbvalueMin / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
801  : 0. ;
802  // Sigma statistical uncertainty:
803  // delta_sigma = sigma * sqrt ( (delta_spectra/spectra)^2 ) = sigma * ( delta_spectra / (spectra-corr * binwidth) )
804  Double_t sigmaRcbvalueStatUnc = (sigmaRcbvalue!=0.) ?
805  sigmaRcbvalue * ( fhRECpt->GetBinError(ibin) / ( yieldRcbvalue * fhRECpt->GetBinWidth(ibin) ) ) : 0. ;
806 
807  fhSigmaCorrRcb->Fill( fhSigmaCorr->GetBinCenter(ibin) , Rhy, sigmaRcbvalue );
808  // if(ibin==3)
809  // cout << " pt "<< fhRECpt->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<" fc-value="<< fhFcRcb->GetBinContent(ibin,rbin) <<", yield-fcRbvalue="<<yieldRcbvalue<<", sigma-fcRbvalue="<<sigmaRcbvalue<<endl;
810  fnSigma->Fill(fhRECpt->GetBinCenter(ibin), fhRECpt->GetBinContent(ibin),
811  Rhy, fc, yieldRcbvalue, sigmaRcbvalue, sigmaRcbvalueStatUnc, sigmaRcbvalueMax, sigmaRcbvalueMin );
812  }
813  }
814  //
815  // Fill the TGraphAsymmErrors
816  if (fAsymUncertainties) {
817  fgSigmaCorrExtreme->SetPoint(ibin,x,value); // i,x,y
818  fgSigmaCorrExtreme->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
819  fgSigmaCorrConservative->SetPoint(ibin,x,value); // i,x,y
820  fgSigmaCorrConservative->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),errvalueConservativeMin,errvalueConservativeMax); // i,xl,xh,yl,yh
821 
822  fhStatUncEffcSigma->SetBinContent(ibin,0.);
823  if(value>0.) fhStatUncEffcSigma->SetBinError(ibin,((errvalueStatUncEffc/value)*100.));
824  fhStatUncEffbSigma->SetBinContent(ibin,0.); fhStatUncEffbSigma->SetBinError(ibin,0.);
825  // cout << " pt "<< fhRECpt->GetBinCenter(ibin) <<" bin "<< ibin<<" stat-unc-c-sigma "<< errvalueStatUncEffc/value << endl;
826  }
827 
828  }
829 
830 }
831 
832 //_________________________________________________________________________________________________________
833 TH1D * AliHFPtSpectrum::EstimateEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco, const char *name) {
841 
842  if(!fhRECpt){
843  AliInfo("Hey, the reconstructed histogram was not set yet !");
844  return NULL;
845  }
846 
847  TH1D * hEfficiency = new TH1D(name," acceptance #times efficiency",fnPtBins,fPtBinLimits);
848 
849  Double_t *sumSimu=new Double_t[fnPtBins];
850  Double_t *sumReco=new Double_t[fnPtBins];
851  for (Int_t ibin=0; ibin<fnPtBins; ibin++){
852  sumSimu[ibin]=0.; sumReco[ibin]=0.;
853  }
854  for (Int_t ibin=0; ibin<=hSimu->GetNbinsX(); ibin++){
855 
856  for (Int_t ibinrec=0; ibinrec<fnPtBins; ibinrec++){
857  if ( hSimu->GetBinCenter(ibin)>fPtBinLimits[ibinrec] &&
858  hSimu->GetBinCenter(ibin)<fPtBinLimits[ibinrec+1] ) {
859  sumSimu[ibinrec]+=hSimu->GetBinContent(ibin);
860  }
861  if ( hReco->GetBinCenter(ibin)>fPtBinLimits[ibinrec] &&
862  hReco->GetBinCenter(ibin)<fPtBinLimits[ibinrec+1] ) {
863  sumReco[ibinrec]+=hReco->GetBinContent(ibin);
864  }
865  }
866 
867  }
868 
869 
870  // the efficiency is computed as reco/sim (in each bin)
871  // its uncertainty is err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
872  Double_t eff=0., erreff=0.;
873  for (Int_t ibinrec=0; ibinrec<fnPtBins; ibinrec++) {
874  if (sumSimu[ibinrec]!= 0. && sumReco[ibinrec]!=0.) {
875  eff = sumReco[ibinrec] / sumSimu[ibinrec] ;
876  // protection in case eff > 1.0
877  // test calculation (make the argument of the sqrt positive)
878  erreff = TMath::Sqrt( eff * TMath::Abs(1.0 - eff) ) / TMath::Sqrt( sumSimu[ibinrec] );
879  }
880  else { eff=0.0; erreff=0.; }
881  hEfficiency->SetBinContent(ibinrec+1,eff);
882  hEfficiency->SetBinError(ibinrec+1,erreff);
883  }
884 
885  delete [] sumSimu;
886  delete [] sumReco;
887 
888  return (TH1D*)hEfficiency;
889 }
890 
891 //_________________________________________________________________________________________________________
900 
901  if(!fhRECpt || !hSimu || !hReco){
902  AliError("Hey, the reconstructed histogram was not set yet !");
903  return;
904  }
905 
906  fhDirectEffpt = EstimateEfficiencyRecoBin(hSimu,hReco,"fhDirectEffpt");
907  fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance #times efficiency");
908 
909 }
910 
911 //_________________________________________________________________________________________________________
920 
921  if(!fhRECpt || !hSimu || !hReco){
922  AliError("Hey, the reconstructed histogram was not set yet !");
923  return;
924  }
925 
926  fhFeedDownEffpt = EstimateEfficiencyRecoBin(hSimu,hReco,"fhFeedDownEffpt");
927  fhFeedDownEffpt->SetNameTitle("fhFeedDownEffpt"," feed-down acceptance #times efficiency");
928 
929 }
930 
931 //_________________________________________________________________________________________________________
936 
937  if (fFeedDownOption==0) {
938  AliInfo("Getting ready for the corrections without feed-down consideration");
939  } else if (fFeedDownOption==1) {
940  AliInfo("Getting ready for the fc feed-down correction calculation");
941  } else if (fFeedDownOption==2) {
942  AliInfo("Getting ready for the Nb feed-down correction calculation");
943  } else { AliError("The calculation option must be <=2"); return kFALSE; }
944 
945  // Start checking the input histograms consistency
946  Bool_t areconsistent=kTRUE;
947 
948  // General checks
949  if (!fhDirectEffpt || !fhRECpt) {
950  AliError(" Reconstructed spectra and/or the Nc efficiency distributions are not defined");
951  return kFALSE;
952  }
953  areconsistent &= CheckHistosConsistency(fhRECpt,fhDirectEffpt);
954  if (!areconsistent) {
955  AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
956  return kFALSE;
957  }
958  if (fFeedDownOption==0) return kTRUE;
959 
960  //
961  // Common checks for options 1 (fc) & 2(Nb)
962  if (!fhFeedDownMCpt || !fhFeedDownEffpt) {
963  AliError(" Theoretical Nb and/or the Nb efficiency distributions are not defined");
964  return kFALSE;
965  }
968  if (fAsymUncertainties) {
970  AliError(" Max/Min theoretical Nb distributions are not defined");
971  return kFALSE;
972  }
974  }
975  if (!areconsistent) {
976  AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
977  return kFALSE;
978  }
979  if (fFeedDownOption>1) return kTRUE;
980 
981  //
982  // Now checks for option 1 (fc correction)
983  if (!fhDirectMCpt) {
984  AliError("Theoretical Nc distributions is not defined");
985  return kFALSE;
986  }
989  if (fAsymUncertainties) {
990  if (!fhDirectMCptMax || !fhDirectMCptMin) {
991  AliError(" Max/Min theoretical Nc distributions are not defined");
992  return kFALSE;
993  }
995  }
996  if (!areconsistent) {
997  AliInfo("Histograms required for fc correction are not consistent (bin width, bounds)");
998  return kFALSE;
999  }
1000 
1001  return kTRUE;
1002 }
1003 
1004 //_________________________________________________________________________________________________________
1009 
1010  if (!h1 || !h2) {
1011  AliError("One or both histograms don't exist");
1012  return kFALSE;
1013  }
1014 
1015  Double_t binwidth1 = h1->GetBinWidth(1);
1016  Double_t binwidth2 = h2->GetBinWidth(1);
1017  Double_t min1 = h1->GetBinCenter(1) - (binwidth1/2.) ;
1018 // Double_t max1 = h1->GetBinCenter(nbins1) + (binwidth1/2.) ;
1019  Double_t min2 = h2->GetBinCenter(1) - (binwidth2/2.) ;
1020 // Double_t max2 = h2->GetBinCenter(nbins2) + (binwidth2/2.) ;
1021 
1022  if (binwidth1!=binwidth2) {
1023  AliInfo(" histograms with different bin width");
1024  return kFALSE;
1025  }
1026  if (min1!=min2) {
1027  AliInfo(" histograms with different minimum");
1028  return kFALSE;
1029  }
1030 // if (max1!=max2) {
1031 // AliInfo(" histograms with different maximum");
1032 // return kFALSE;
1033 // }
1034 
1035  return kTRUE;
1036 }
1037 
1038 //_________________________________________________________________________________________________________
1043 
1044 
1045  // declare the output histograms
1046  fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (no feed-down corr)",fnPtBins,fPtBinLimits);
1047  fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (no feed-down corr)",fnPtBins,fPtBinLimits);
1048  fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (no feed-down corr)",fnPtBins,fPtBinLimits);
1050 
1051  //
1052  // Do the calculation
1053  //
1054  for (Int_t ibin=1; ibin<=fnPtBins; ibin++) {
1055  Double_t value = fhRECpt->GetBinContent(ibin) /fhRECpt->GetBinWidth(ibin);
1056  Double_t errvalue = fhRECpt->GetBinError(ibin) /fhRECpt->GetBinWidth(ibin);
1057  fhYieldCorr->SetBinContent(ibin,value);
1058  fhYieldCorr->SetBinError(ibin,errvalue);
1059  fhYieldCorrMax->SetBinContent(ibin,value);
1060  fhYieldCorrMax->SetBinError(ibin,errvalue);
1061  fhYieldCorrMin->SetBinContent(ibin,value);
1062  fhYieldCorrMin->SetBinError(ibin,errvalue);
1063  }
1064 
1065  fAsymUncertainties=kFALSE;
1066 }
1067 
1068 //_________________________________________________________________________________________________________
1081  AliInfo("Calculating the feed-down correction factor (fc method)");
1082 
1083  // define the variables
1084  Double_t correction=1.;
1085  Double_t theoryRatio=1.;
1086  Double_t effRatio=1.;
1087  Double_t correctionExtremeA=1., correctionExtremeB=1.;
1088  Double_t theoryRatioExtremeA=1., theoryRatioExtremeB=1.;
1089  Double_t correctionConservativeA=1., correctionConservativeB=1.;
1090  Double_t theoryRatioConservativeA=1., theoryRatioConservativeB=1.;
1091  // Double_t correctionUnc=0.;
1092  Double_t correctionExtremeAUnc=0., correctionExtremeBUnc=0.;
1093  Double_t correctionConservativeAUnc=0., correctionConservativeBUnc=0.;
1094 
1095  // declare the output histograms
1096  fhFc = new TH1D("fhFc","fc correction factor",fnPtBins,fPtBinLimits);
1097  fhFcMax = new TH1D("fhFcMax","max fc correction factor",fnPtBins,fPtBinLimits);
1098  fhFcMin = new TH1D("fhFcMin","min fc correction factor",fnPtBins,fPtBinLimits);
1099  if(fPbPbElossHypothesis) {
1100  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.);
1101  fnHypothesis = new TNtuple("fnHypothesis"," Feed-down correction vs hypothesis (fc)","pt:Rcb:fc:fcMax:fcMin");
1102  }
1103  // two local control histograms
1104  TH1D *hTheoryRatio = new TH1D("hTheoryRatio","Theoretical B-->D over c-->D (feed-down/direct) ratio",fnPtBins,fPtBinLimits);
1105  TH1D *hEffRatio = new TH1D("hEffRatio","Efficiency B-->D over c-->D (feed-down/direct) ratio",fnPtBins,fPtBinLimits);
1106  // and the output TGraphAsymmErrors
1107  if (fAsymUncertainties) {
1109  fgFcExtreme->SetNameTitle("fgFcExtreme","fgFcExtreme");
1111  fgFcConservative->SetNameTitle("fgFcConservative","fgFcConservative");
1112  }
1113 
1114  fhStatUncEffcFD = new TH1D("fhStatUncEffcFD","direct charm stat unc on the feed-down correction",fnPtBins,fPtBinLimits);
1115  fhStatUncEffbFD = new TH1D("fhStatUncEffbFD","secondary charm stat unc on the feed-down correction",fnPtBins,fPtBinLimits);
1116  Double_t correctionConservativeAUncStatEffc=0., correctionConservativeBUncStatEffc=0.;
1117  Double_t correctionConservativeAUncStatEffb=0., correctionConservativeBUncStatEffb=0.;
1118 
1119  //
1120  // Compute fc
1121  //
1122  for (Int_t ibin=1; ibin<=fnPtBins; ibin++) {
1123 
1124  // Variables initialization
1125  correction=1.; theoryRatio=1.; effRatio=1.;
1126  correctionExtremeA=1.; correctionExtremeB=1.;
1127  theoryRatioExtremeA=1.; theoryRatioExtremeB=1.;
1128  correctionConservativeA=1.; correctionConservativeB=1.;
1129  theoryRatioConservativeA=1.; theoryRatioConservativeB=1.;
1130  // correctionUnc=0.;
1131  correctionExtremeAUnc=0.; correctionExtremeBUnc=0.;
1132  correctionConservativeAUnc=0.; correctionConservativeBUnc=0.;
1133  correctionConservativeAUncStatEffc=0.; correctionConservativeBUncStatEffc=0.;
1134  correctionConservativeAUncStatEffb=0.; correctionConservativeBUncStatEffb=0.;
1135 
1136  // theory_ratio = (N_b/N_c)
1137  theoryRatio = (fhDirectMCpt->GetBinContent(ibin)>0. && fhFeedDownMCpt->GetBinContent(ibin)>0.) ?
1138  fhFeedDownMCpt->GetBinContent(ibin) / fhDirectMCpt->GetBinContent(ibin) : 1.0 ;
1139 
1140  //
1141  // Calculate the uncertainty [ considering only the theoretical uncertainties on Nb & Nc for now !!! ]
1142  //
1143  // extreme A = direct-max, feed-down-min
1144  theoryRatioExtremeA = (fhDirectMCptMax->GetBinContent(ibin)>0. && fhFeedDownMCptMin->GetBinContent(ibin)>0.) ?
1145  fhFeedDownMCptMin->GetBinContent(ibin) / fhDirectMCptMax->GetBinContent(ibin) : 1.0 ;
1146  // extreme B = direct-min, feed-down-max
1147  theoryRatioExtremeB = (fhDirectMCptMin->GetBinContent(ibin)>0. && fhDirectMCptMax->GetBinContent(ibin)>0.) ?
1148  fhFeedDownMCptMax->GetBinContent(ibin) / fhDirectMCptMin->GetBinContent(ibin) : 1.0 ;
1149  // conservative A = direct-max, feed-down-max
1150  theoryRatioConservativeA = (fhDirectMCptMax->GetBinContent(ibin)>0. && fhFeedDownMCptMin->GetBinContent(ibin)>0.) ?
1151  fhFeedDownMCptMax->GetBinContent(ibin) / fhDirectMCptMax->GetBinContent(ibin) : 1.0 ;
1152  // conservative B = direct-min, feed-down-min
1153  theoryRatioConservativeB = (fhDirectMCptMin->GetBinContent(ibin)>0. && fhDirectMCptMax->GetBinContent(ibin)>0.) ?
1154  fhFeedDownMCptMin->GetBinContent(ibin) / fhDirectMCptMin->GetBinContent(ibin) : 1.0 ;
1155 
1156  // eff_ratio = (eff_b/eff_c)
1157  effRatio = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ?
1158  fhFeedDownEffpt->GetBinContent(ibin) / fhDirectEffpt->GetBinContent(ibin) : 1.0 ;
1159 
1160  // fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
1161  if( TMath::Abs(effRatio - 1.0)<0.0001 || TMath::Abs(theoryRatio - 1.0)<0.0001 ) {
1162  correction = 1.0;
1163  correctionExtremeA = 1.0;
1164  correctionExtremeB = 1.0;
1165  correctionConservativeA = 1.0;
1166  correctionConservativeB = 1.0;
1167  }
1168  else {
1169  correction = ( 1. / ( 1 + ( effRatio * theoryRatio ) ) );
1170  correctionExtremeA = ( 1. / ( 1 + ( effRatio * theoryRatioExtremeA ) ) );
1171  correctionExtremeB = ( 1. / ( 1 + ( effRatio * theoryRatioExtremeB ) ) );
1172  correctionConservativeA = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeA ) ) );
1173  correctionConservativeB = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeB ) ) );
1174  }
1175 
1176 
1177 
1178  // uncertainy on the feed-down efficiency, which impacts the feed-down estimate
1179  // pt-dependent estimate is 50% of the cut variation uncertainty from HFSystErr
1180  Double_t x = fhDirectEffpt->GetBinCenter(ibin);
1181  Double_t efficiencyUnc = fGlobalEfficiencyUncertainties[1];
1183 // cout <<endl<<" mine test: gl-unc="<<fGlobalEfficiencyUncertainties[1]<<" cutvar-unc="<<CalculateEfficiencyPtDepedentUncertainty(x,true)<<" ==> effUnc="<<efficiencyUnc<<endl;
1184 
1185  // fc uncertainty from (eff_b/eff_c) = fc^2 * (N_b/N_c) * delta(eff_b/eff_c)
1186  // delta(eff_b/eff_c) is a percentage = effRatio * sqrt( fGlobalEfficiencyUncertainties[1]^2 + unc_eff_c ^2 + unc_eff_b ^2 )
1187  Double_t relEffUnc = TMath::Sqrt( efficiencyUnc*efficiencyUnc +
1188  (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
1189  (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))
1190  );
1191 
1192  // correctionUnc = correction*correction * theoryRatio * effRatio * relEffUnc;
1193  correctionExtremeAUnc = correctionExtremeA*correctionExtremeA * theoryRatioExtremeA * effRatio * relEffUnc;
1194  correctionExtremeBUnc = correctionExtremeB*correctionExtremeB * theoryRatioExtremeB * effRatio * relEffUnc;
1195 
1196  correctionConservativeAUnc = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA *effRatio * relEffUnc;
1197  //
1198  correctionConservativeAUncStatEffc = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA *effRatio *
1199  (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin));
1200  correctionConservativeAUncStatEffb = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA *effRatio *
1201  (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin));
1202 
1203  correctionConservativeBUnc = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB *effRatio * relEffUnc;
1204 
1205  correctionConservativeBUncStatEffb = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB *effRatio *
1206  (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin));
1207  correctionConservativeBUncStatEffc = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB *effRatio *
1208  (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin));
1209 
1210 
1211  // Fill in the histograms
1212  hTheoryRatio->SetBinContent(ibin,theoryRatio);
1213  hEffRatio->SetBinContent(ibin,effRatio);
1214  fhFc->SetBinContent(ibin,correction);
1215  //
1216  // Estimate how the result varies vs charm/beauty Eloss hypothesis
1217  //
1218  if ( correction>1.0e-16 && fPbPbElossHypothesis){
1219  // Loop over the Eloss hypothesis
1220  // Int_t rbin=0;
1221  for (Float_t rval=0.0025; rval<4.0; rval+=0.005){
1222  // Central fc with Eloss hypothesis
1223  Double_t correctionRcb = ( 1. / ( 1 + ( effRatio * theoryRatio * (1/rval) ) ) );
1224  fhFcRcb->Fill( fhFc->GetBinCenter(ibin) , rval, correctionRcb );
1225  // if(ibin==3){
1226  // cout << " pt "<< fhFc->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<", fc-Rcb-value="<<correctionRcb<<endl;
1227  // rbin++;
1228  // }
1229  // Upper / lower fc with up / low FONLL bands and Eloss hypothesis
1230  Double_t correctionConservativeARcb = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeA * (1/rval) ) ) );
1231  Double_t correctionConservativeBRcb = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeB * (1/rval) ) ) );
1232  Double_t correctionConservativeARcbUnc = correctionConservativeARcb*correctionConservativeARcb * theoryRatioConservativeA * (1/rval) *effRatio * relEffUnc;
1233  Double_t correctionConservativeBRcbUnc = correctionConservativeBRcb*correctionConservativeBRcb * theoryRatioConservativeB * (1/rval) *effRatio * relEffUnc;
1234  Double_t consvalRcb[4] = { correctionConservativeARcb - correctionConservativeARcbUnc, correctionConservativeARcb + correctionConservativeARcbUnc,
1235  correctionConservativeBRcb - correctionConservativeBRcbUnc, correctionConservativeBRcb + correctionConservativeBRcbUnc};
1236  Double_t uncConservativeRcbMin = correctionRcb - TMath::MinElement(4,consvalRcb);
1237  Double_t uncConservativeRcbMax = TMath::MaxElement(4,consvalRcb) - correctionRcb;
1238 // if(ibin==3)
1239 // cout << " pt="<<fhDirectEffpt->GetBinCenter(ibin)<<", hypo="<<rval<<", fc="<<correctionRcb<<" +"<<uncConservativeRcbMax<<" -"<<uncConservativeRcbMin<<endl;
1240  fnHypothesis->Fill( fhDirectEffpt->GetBinCenter(ibin), rval, correctionRcb, correctionRcb+uncConservativeRcbMax, correctionRcb-uncConservativeRcbMin);
1241  }
1242  }
1243  //
1244  // Fill the rest of (asymmetric) histograms
1245  //
1246  if (fAsymUncertainties) {
1247  Double_t x = fhDirectMCpt->GetBinCenter(ibin);
1248  Double_t val[4] = { correctionExtremeA + correctionExtremeAUnc, correctionExtremeA - correctionExtremeAUnc,
1249  correctionExtremeB + correctionExtremeBUnc, correctionExtremeB - correctionExtremeBUnc };
1250  Double_t uncExtremeMin = correction - TMath::MinElement(4,val);
1251  Double_t uncExtremeMax = TMath::MaxElement(4,val) - correction;
1252  fgFcExtreme->SetPoint(ibin,x,correction); // i,x,y
1253  fgFcExtreme->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),uncExtremeMin,uncExtremeMax); // i,xl,xh,yl,yh
1254  fhFcMax->SetBinContent(ibin,correction+uncExtremeMax);
1255  fhFcMin->SetBinContent(ibin,correction-uncExtremeMin);
1256  Double_t consval[4] = { correctionConservativeA - correctionConservativeAUnc, correctionConservativeA + correctionConservativeAUnc,
1257  correctionConservativeB - correctionConservativeBUnc, correctionConservativeB + correctionConservativeBUnc};
1258  Double_t uncConservativeMin = correction - TMath::MinElement(4,consval);
1259  Double_t uncConservativeMax = TMath::MaxElement(4,consval) - correction;
1260  fgFcConservative->SetPoint(ibin,x,correction); // i,x,y
1261  fgFcConservative->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),uncConservativeMin,uncConservativeMax); // i,xl,xh,yl,yh
1262  if( !(correction>0.) ){
1263  fgFcExtreme->SetPoint(ibin,x,0.); // i,x,y
1264  fgFcExtreme->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),0.,0.); // i,xl,xh,yl,yh
1265  fgFcConservative->SetPoint(ibin,x,0.); // i,x,y
1266  fgFcConservative->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),0.,0.); // i,xl,xh,yl,yh
1267  }
1268 
1269  Double_t valStatEffc[2] = { correctionConservativeAUncStatEffc/correctionConservativeA,
1270  correctionConservativeBUncStatEffc/correctionConservativeB };
1271  Double_t valStatEffb[2] = { correctionConservativeAUncStatEffb/correctionConservativeA,
1272  correctionConservativeBUncStatEffb/correctionConservativeB };
1273  Double_t uncConservativeStatEffc = TMath::MaxElement(2,valStatEffc);
1274  Double_t uncConservativeStatEffb = TMath::MaxElement(2,valStatEffb);
1275  fhStatUncEffcFD->SetBinContent(ibin,0.); fhStatUncEffcFD->SetBinError(ibin,uncConservativeStatEffc*100.);
1276  fhStatUncEffbFD->SetBinContent(ibin,0.); fhStatUncEffbFD->SetBinError(ibin,uncConservativeStatEffb*100.);
1277  // cout << " pt "<< fhStatUncEffcFD->GetBinCenter(ibin) <<" bin "<< ibin<<" fc-stat-c ="<<uncConservativeStatEffc<<" fc-stat-b ="<<uncConservativeStatEffb<<endl;
1278  }
1279 
1280  }
1281 
1282 }
1283 
1284 //_________________________________________________________________________________________________________
1297 
1298  AliInfo(" Calculating the feed-down corrected spectrum (fc method)");
1299 
1300  if (!fhFc || !fhRECpt) {
1301  AliError(" Reconstructed or fc distributions are not defined");
1302  return;
1303  }
1304 
1305  Double_t value = 0., errvalue = 0., errvalueMax= 0., errvalueMin= 0.;
1306  Double_t valueExtremeMax= 0., valueExtremeMin= 0.;
1307  Double_t valueConservativeMax= 0., valueConservativeMin= 0.;
1308 
1309  // declare the output histograms
1310  fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by fc)",fnPtBins,fPtBinLimits);
1311  fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by fc)",fnPtBins,fPtBinLimits);
1312  fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by fc)",fnPtBins,fPtBinLimits);
1313  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.);
1314  // and the output TGraphAsymmErrors
1316  if (fAsymUncertainties){
1319  }
1320 
1321  //
1322  // Do the calculation
1323  //
1324  for (Int_t ibin=1; ibin<=fnPtBins; ibin++) {
1325 
1326  // Variables initialization
1327  value = 0.; errvalue = 0.; errvalueMax= 0.; errvalueMin= 0.;
1328  valueExtremeMax= 0.; valueExtremeMin= 0.;
1329  valueConservativeMax= 0.; valueConservativeMin= 0.;
1330 
1331 
1332  // calculate the value
1333  // physics = reco * fc / bin-width
1334  value = (fhRECpt->GetBinContent(ibin) && fhFc->GetBinContent(ibin)) ?
1335  fhRECpt->GetBinContent(ibin) * fhFc->GetBinContent(ibin) : 0. ;
1336  value /= fhRECpt->GetBinWidth(ibin) ;
1337 
1338  // Statistical uncertainty
1339  // (stat) delta_physics = physics * sqrt ( (delta_reco/reco)^2 )
1340  errvalue = (value!=0. && fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) ?
1341  value * (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) : 0. ;
1342 
1343  // Calculate the systematic uncertainties
1344  // (syst but feed-down) delta_physics = physics * sqrt ( (delta_reco_syst/reco)^2 )
1345  // (feed-down syst) delta_physics = physics * sqrt ( (delta_fc/fc)^2 )
1346  //
1347  // Protect against null denominator. If so, define uncertainty as null
1348  if (fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) {
1349 
1350  if (fAsymUncertainties) {
1351 
1352  // Systematics but feed-down
1353  if (fgRECSystematics) {
1354  errvalueMax = value * ( fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinContent(ibin) );
1355  errvalueMin = value * ( fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinContent(ibin) );
1356  }
1357  else { errvalueMax = 0.; errvalueMin = 0.; }
1358 
1359  // Extreme feed-down systematics
1360  valueExtremeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcExtreme->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1361  valueExtremeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcExtreme->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1362 
1363  // Conservative feed-down systematics
1364  valueConservativeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcConservative->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1365  valueConservativeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcConservative->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1366 
1367  }
1368 
1369  }
1370  else { errvalueMax = 0.; errvalueMin = 0.; }
1371 
1372  //
1373  // Fill in the histograms
1374  //
1375  fhYieldCorr->SetBinContent(ibin,value);
1376  fhYieldCorr->SetBinError(ibin,errvalue);
1377  //
1378  // Fill the histos and ntuple vs the Eloss hypothesis
1379  //
1380  if (fPbPbElossHypothesis) {
1381  // Loop over the Eloss hypothesis
1382  for (Float_t rval=0.0025; rval<4.0; rval+=0.005){
1383  Int_t rbin = FindTH2YBin(fhYieldCorrRcb,rval);
1384  Double_t fcRcbvalue = fhFcRcb->GetBinContent(ibin,rbin);
1385  // physics = reco * fcRcb / bin-width
1386  Double_t Rcbvalue = (fhRECpt->GetBinContent(ibin) && fcRcbvalue) ?
1387  fhRECpt->GetBinContent(ibin) * fcRcbvalue : 0. ;
1388  Rcbvalue /= fhRECpt->GetBinWidth(ibin) ;
1389  fhYieldCorrRcb->Fill( fhYieldCorr->GetBinCenter(ibin) , rval, Rcbvalue );
1390  // cout << " pt "<< fhRECpt->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<" fc-fcRbvalue="<<fcRcbvalue<<", yield="<<Rcbvalue<<endl;
1391  }
1392  }
1393  if (fAsymUncertainties) {
1394  Double_t center = fhYieldCorr->GetBinCenter(ibin);
1395  fgYieldCorr->SetPoint(ibin,center,value); // i,x,y
1396  fgYieldCorr->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
1397  fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax);
1398  fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
1399  fgYieldCorrExtreme->SetPoint(ibin,center,value);
1400  fgYieldCorrExtreme->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),value-valueExtremeMin,valueExtremeMax-value);
1401  fgYieldCorrConservative->SetPoint(ibin,center,value);
1402  fgYieldCorrConservative->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),value-valueConservativeMin,valueConservativeMax-value);
1403  }
1404 
1405  }
1406 
1407 }
1408 
1409 //_________________________________________________________________________________________________________
1410 void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumNb(Double_t deltaY, Double_t branchingRatioBintoFinalDecay) {
1424  AliInfo("Calculating the feed-down correction factor and spectrum (Nb method)");
1425 
1426  Double_t value = 0., errvalue = 0., errvalueMax = 0., errvalueMin = 0., kfactor = 0.;
1427  Double_t errvalueExtremeMax = 0., errvalueExtremeMin = 0.;
1428 
1429  // declare the output histograms
1430  fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by Nb)",fnPtBins,fPtBinLimits);
1431  fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by Nb)",fnPtBins,fPtBinLimits);
1432  fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by Nb)",fnPtBins,fPtBinLimits);
1433  if(fPbPbElossHypothesis) {
1434  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.);
1435  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.);
1436  fnHypothesis = new TNtuple("fnHypothesis"," Feed-down correction vs hypothesis (Nb)","pt:Rb:fc:fcMax:fcMin");
1437  }
1438  // and the output TGraphAsymmErrors
1440  if (fAsymUncertainties){
1443  // Define fc-conservative
1445  AliInfo(" Beware the conservative & extreme uncertainties are equal by definition !");
1446  }
1447 
1448  // variables to define fc-conservative
1449  double correction=0, correctionMax=0., correctionMin=0.;
1450 
1451  fhStatUncEffcFD = new TH1D("fhStatUncEffcFD","direct charm stat unc on the feed-down correction",fnPtBins,fPtBinLimits);
1452  fhStatUncEffbFD = new TH1D("fhStatUncEffbFD","secondary charm stat unc on the feed-down correction",fnPtBins,fPtBinLimits);
1453  Double_t correctionUncStatEffc=0.;
1454  Double_t correctionUncStatEffb=0.;
1455 
1456 
1457  //
1458  // Do the calculation
1459  //
1460  for (Int_t ibin=1; ibin<=fnPtBins; ibin++) {
1461 
1462  // Calculate the value
1463  // physics = [ reco - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) ] / bin-width
1464  // In HIC : physics = [ reco - ( Tab * Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th * Rb ) ] / bin-width
1465  //
1466  //
1467  Double_t frac = 1.0, errfrac =0.;
1468 
1469  // Variables initialization
1470  value = 0.; errvalue = 0.; errvalueMax = 0.; errvalueMin = 0.; kfactor = 0.;
1471  errvalueExtremeMax = 0.; errvalueExtremeMin = 0.;
1472  correction=0; correctionMax=0.; correctionMin=0.;
1473  correctionUncStatEffc=0.; correctionUncStatEffb=0.;
1474 
1475  // uncertainy on the feed-down efficiency, which impacts the feed-down estimate
1476  Double_t x = fhYieldCorr->GetBinCenter(ibin);
1477  Double_t efficiencyUnc = fGlobalEfficiencyUncertainties[0];
1479 
1481  frac = fTab[0]*fNevts;
1482  if(fIsEventPlane) frac = frac/2.0;
1483  errfrac = frac * TMath::Sqrt( (fTab[1]/fTab[0])*(fTab[1]/fTab[0]) + (1/fNevts) );
1484  } else {
1485  frac = fLuminosity[0];
1486  errfrac = fLuminosity[1];
1487  }
1488  // printf("Tab=%e events=%d frac=%f\n",fTab[0],(Int_t)fNevts,frac);
1489  value = ( fhRECpt->GetBinContent(ibin)>0. && fhRECpt->GetBinContent(ibin)!=0. &&
1490  fhFeedDownMCpt->GetBinContent(ibin)>0. && fhFeedDownEffpt->GetBinContent(ibin)>0. ) ?
1491  fhRECpt->GetBinContent(ibin) - frac*(deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) * fhRECpt->GetBinWidth(ibin) )
1492  : 0. ;
1493  // printf("%d raw=%f after=%f \n",ibin,fhRECpt->GetBinContent(ibin),value);
1494  value /= fhRECpt->GetBinWidth(ibin);
1495  if (value<0.) value =0.;
1496 
1497  // Statistical uncertainty: delta_physics = sqrt ( (delta_reco)^2 ) / bin-width
1498  errvalue = (value!=0. && fhRECpt->GetBinError(ibin) && fhRECpt->GetBinError(ibin)!=0.) ?
1499  fhRECpt->GetBinError(ibin) : 0.;
1500  errvalue /= fhRECpt->GetBinWidth(ibin);
1501 
1502  // Correction (fc) : Estimate of the relative amount feed-down subtracted
1503  // correction = [ 1 - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th)/reco ]
1504  // in HIC: correction = [ 1 - ( Tab * Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th)/reco ]
1505  correction = (value>0.) ?
1506  1 - (frac*deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) * fhRECpt->GetBinWidth(ibin) ) / fhRECpt->GetBinContent(ibin) : 0. ;
1507  if (correction<0.) correction = 0.;
1508 
1509  // 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;
1510 
1511  // Systematic uncertainties
1512  // (syst but feed-down) delta_physics = sqrt ( (delta_reco_syst)^2 ) / bin-width
1513  // (feed-down syst) delta_physics = sqrt ( (k*delta_lumi/lumi)^2 + (k*delta_eff_trig/eff_trig)^2
1514  // + (k*delta_Nb/Nb)^2 + (k*delta_eff/eff)^2 + (k*global_eff_ratio)^2 ) / bin-width
1515  // where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th * bin-width
1516  kfactor = frac*deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) * fhRECpt->GetBinWidth(ibin) ;
1517  //
1518  if (fAsymUncertainties && value>0.) {
1519  Double_t nb = fhFeedDownMCpt->GetBinContent(ibin);
1520  Double_t nbDmax = fhFeedDownMCptMax->GetBinContent(ibin) - fhFeedDownMCpt->GetBinContent(ibin);
1521  Double_t nbDmin = fhFeedDownMCpt->GetBinContent(ibin) - fhFeedDownMCptMin->GetBinContent(ibin);
1522 
1523  // Systematics but feed-down
1524  if (fgRECSystematics){
1525  errvalueMax = fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinWidth(ibin) ;
1526  errvalueMin = fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinWidth(ibin);
1527  }
1528  else { errvalueMax = 0.; errvalueMin = 0.; }
1529 
1530  // Feed-down systematics
1531  // min value with the maximum Nb
1532  Double_t errCom = ( (kfactor*errfrac/frac)*(kfactor*errfrac/frac) ) +
1533  ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1534  ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
1535  ( (kfactor*efficiencyUnc)*(kfactor*efficiencyUnc) ) ;
1536  errvalueExtremeMin = TMath::Sqrt( errCom + ( (kfactor*nbDmax/nb)*(kfactor*nbDmax/nb) ) ) / fhRECpt->GetBinWidth(ibin);
1537  // max value with the minimum Nb
1538  errvalueExtremeMax = TMath::Sqrt( errCom + ( (kfactor*nbDmin/nb)*(kfactor*nbDmin/nb) ) ) / fhRECpt->GetBinWidth(ibin);
1539 
1540  // Correction systematics (fc)
1541  // min value with the maximum Nb
1542  correctionMin = TMath::Sqrt( errCom + ( (kfactor*nbDmax/nb)*(kfactor*nbDmax/nb) ) ) / fhRECpt->GetBinContent(ibin) ;
1543  // max value with the minimum Nb
1544  correctionMax = TMath::Sqrt( errCom + ( (kfactor*nbDmin/nb)*(kfactor*nbDmin/nb) ) ) / fhRECpt->GetBinContent(ibin) ;
1545  //
1546  correctionUncStatEffb = TMath::Sqrt( ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) )
1547  ) / fhRECpt->GetBinContent(ibin) ;
1548  correctionUncStatEffc = 0.;
1549  }
1550  else{ // Don't consider Nb uncertainty in this case [ to be tested!!! ]
1551  errvalueExtremeMax = TMath::Sqrt( ( (kfactor*errfrac/frac)*(kfactor*errfrac/frac) ) +
1552  ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1553  ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
1554  ( (kfactor*efficiencyUnc)*(kfactor*efficiencyUnc) )
1555  ) / fhRECpt->GetBinWidth(ibin);
1556  errvalueExtremeMin = errvalueExtremeMax ;
1557  }
1558 
1559 
1560  // fill in histograms
1561  fhYieldCorr->SetBinContent(ibin,value);
1562  fhYieldCorr->SetBinError(ibin,errvalue);
1563  //
1564  // Estimate how the result varies vs charm/beauty Eloss hypothesis
1565  //
1566  if ( correction>1.0e-16 && fPbPbElossHypothesis){
1567  // Loop over the Eloss hypothesis
1568  // Int_t rbin=0;
1569  for (Float_t rval=0.0025; rval<4.0; rval+=0.005){
1570  // correction = [ 1 - (Tab *Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th *binwidth )* (rval) /reco ]
1571  Double_t fcRcbvalue = 1 - (frac*deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin)*fhRECpt->GetBinWidth(ibin) * rval ) / fhRECpt->GetBinContent(ibin) ;
1572  // cout << " rval="<<rval <<" , fc="<<fcRcbvalue<<" ";
1573  if(fcRcbvalue<1.e-16) fcRcbvalue=0.;
1574  fhFcRcb->Fill( fhRECpt->GetBinCenter(ibin) , rval, fcRcbvalue );
1575  // physics = reco * fcRcb / bin-width
1576  Double_t Rcbvalue = (fhRECpt->GetBinContent(ibin) && fcRcbvalue) ?
1577  fhRECpt->GetBinContent(ibin) * fcRcbvalue : 0. ;
1578  Rcbvalue /= fhRECpt->GetBinWidth(ibin) ;
1579  fhYieldCorrRcb->Fill( fhYieldCorr->GetBinCenter(ibin) , rval, Rcbvalue );
1580  // if(ibin==3){
1581  // cout << " pt "<< fhFcRcb->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<", fc-Rb-value="<< fcRcbvalue << ", yield-Rb-value="<< Rcbvalue <<endl;
1582  // cout << " pt "<< fhFcRcb->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", fc-Rb-value="<< fcRcbvalue << ", yield-Rb-value="<< Rcbvalue <<endl;
1583  // rbin++;
1584  // }
1585  Double_t correctionMaxRcb = correctionMax*rval;
1586  Double_t correctionMinRcb = correctionMin*rval;
1587  fnHypothesis->Fill( fhYieldCorr->GetBinCenter(ibin), rval, fcRcbvalue, fcRcbvalue + correctionMaxRcb, fcRcbvalue - correctionMinRcb);
1588 // if(ibin==3){
1589 // cout << " pt="<< fhFcRcb->GetBinCenter(ibin) <<", rval="<<rval<<", fc="<<fcRcbvalue<<" +"<<correctionMaxRcb<<" -"<<correctionMinRcb<<endl;}
1590  }
1591  }
1592  //
1593  // Fill the rest of (asymmetric) histograms
1594  //
1595  if (fAsymUncertainties) {
1596 // Double_t x = fhYieldCorr->GetBinCenter(ibin);
1597  fgYieldCorr->SetPoint(ibin,x,value); // i,x,y
1598  fgYieldCorr->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
1599  fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax);
1600  fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
1601  fgYieldCorrExtreme->SetPoint(ibin,x,value); // i,x,y
1602  fgYieldCorrExtreme->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
1603  fgYieldCorrConservative->SetPoint(ibin,x,value); // i,x,y
1604  fgYieldCorrConservative->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
1605  // cout << " bin " << ibin << ", correction " << correction << ", min correction unc " << correctionMin << ", max correction unc " << correctionMax << endl;
1606  if(correction>0.){
1607  fgFcConservative->SetPoint(ibin,x,correction);
1608  fgFcConservative->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),correctionMin,correctionMax);
1609 
1610  fhStatUncEffbFD->SetBinContent(ibin,0.); fhStatUncEffbFD->SetBinError(ibin,correctionUncStatEffb/correction*100.);
1611  fhStatUncEffcFD->SetBinContent(ibin,0.); fhStatUncEffcFD->SetBinError(ibin,correctionUncStatEffc/correction*100.);
1612  // cout << " pt "<< fhStatUncEffcFD->GetBinCenter(ibin) <<" bin "<< ibin<<" fc-stat-c ="<< correctionUncStatEffc/correction <<" fc-stat-b ="<< correctionUncStatEffb/correction <<endl;
1613  }
1614  else{
1615  fgFcConservative->SetPoint(ibin,x,0.);
1616  fgFcConservative->SetPointError(ibin,(fPtBinWidths[ibin-1]/2.),(fPtBinWidths[ibin-1]/2.),0.,0.);
1617  }
1618  }
1619 
1620  }
1621 
1622 }
1623 
1624 
1625 //_________________________________________________________________________________________________________
1632  Double_t uncertainty = 0.;
1633  Double_t trackingUnc = fSystematics->GetTrackingEffErr(pt);
1634  Double_t cutVarUnc = fSystematics->GetCutsEffErr(pt);
1635  uncertainty = TMath::Sqrt( trackingUnc*trackingUnc + cutVarUnc*cutVarUnc );
1636  cout<<" cutVar="<<cutVarUnc<<" trackUnc="<<trackingUnc<<endl;
1637  if(useOnlyCutVar) uncertainty = cutVarUnc;
1638  return uncertainty;
1639 }
1640 
1641 //_________________________________________________________________________________________________________
1642 //void AliHFPtSpectrum::ComputeSystUncertainties(AliHFSystErr *systematics, Bool_t combineFeedDown) {
1644 
1650 
1651  // Estimate the feed-down uncertainty in percentage
1652  Int_t nentries = 0;
1653  TGraphAsymmErrors *grErrFeeddown = 0;
1654  Double_t x=0., y=0., errx=0., erryl=0., erryh=0;
1655  if(fFeedDownOption!=0) {
1656  nentries = fgSigmaCorrConservative->GetN();
1657  grErrFeeddown = new TGraphAsymmErrors(nentries);
1658  for(Int_t i=0; i<nentries; i++) {
1659  x=0.; y=0.; errx=0.; erryl=0.; erryh=0.;
1660  fgSigmaCorrConservative->GetPoint(i,x,y);
1661  if(y>0.){
1662  errx = fgSigmaCorrConservative->GetErrorXlow(i) ;
1663  erryl = fgSigmaCorrConservative->GetErrorYlow(i) / y ;
1664  erryh = fgSigmaCorrConservative->GetErrorYhigh(i) / y ;
1665  }
1666  // cout << " x "<< x << " +- "<<errx<<" , y "<<y<<" + "<<erryh<<" - "<<erryl<<endl;
1667  grErrFeeddown->SetPoint(i,x,0.);
1668  grErrFeeddown->SetPointError(i,errx,errx,erryl,erryh); //i, xl, xh, yl, yh
1669  }
1670  }else{
1671  nentries = fgSigmaCorr->GetN();
1672  }
1673 
1674  // Draw all the systematics independently
1675 // systematics->DrawErrors(grErrFeeddown);
1676  fSystematics->DrawErrors(grErrFeeddown);
1677  // Set the sigma systematic uncertainties
1678  // possibly combine with the feed-down uncertainties
1679  Double_t errylcomb=0., erryhcomb=0;
1680  for(Int_t i=0; i<nentries; i++) {
1681  fgSigmaCorr->GetPoint(i,x,y);
1682  if (fFeedDownOption!=0 && combineFeedDown) {
1683 // errylcomb = systematics->GetTotalSystErr(x,erryl) * y ;
1684 // erryhcomb = systematics->GetTotalSystErr(x,erryh) * y ;
1685  errx = grErrFeeddown->GetErrorXlow(i) ;
1686  erryl = grErrFeeddown->GetErrorYlow(i);
1687  erryh = grErrFeeddown->GetErrorYhigh(i);
1688  errylcomb = fSystematics->GetTotalSystErr(x,erryl) * y ;
1689  erryhcomb = fSystematics->GetTotalSystErr(x,erryh) * y ;
1690  } else {
1691 // errylcomb = systematics->GetTotalSystErr(x) * y ;
1692 // erryhcomb = systematics->GetTotalSystErr(x) * y ;
1693  errylcomb = fSystematics->GetTotalSystErr(x) * y ;
1694  erryhcomb = fSystematics->GetTotalSystErr(x) * y ;
1695  }
1696  fgSigmaCorr->SetPointError(i,errx,errx,errylcomb,erryhcomb);
1697  //
1698  fhSigmaCorrDataSyst->SetBinContent(i,y);
1699 // erryl = systematics->GetTotalSystErr(x) * y ;
1700  erryl = fSystematics->GetTotalSystErr(x) * y ;
1701  fhSigmaCorrDataSyst->SetBinError(i,erryl);
1702  }
1703 
1704 }
1705 
1706 
1707 //_________________________________________________________________________________________________________
1712 
1713  TCanvas *csigma = new TCanvas("csigma","Draw the corrected cross-section & the prediction");
1714  csigma->SetFillColor(0);
1715  gPrediction->GetXaxis()->SetTitleSize(0.05);
1716  gPrediction->GetXaxis()->SetTitleOffset(0.95);
1717  gPrediction->GetYaxis()->SetTitleSize(0.05);
1718  gPrediction->GetYaxis()->SetTitleOffset(0.95);
1719  gPrediction->GetXaxis()->SetTitle("p_{T} [GeV]");
1720  gPrediction->GetYaxis()->SetTitle("BR #times #frac{d#sigma}{dp_{T}} |_{|y|<0.5} [pb/GeV]");
1721  gPrediction->SetLineColor(kGreen+2);
1722  gPrediction->SetLineWidth(3);
1723  gPrediction->SetFillColor(kGreen+1);
1724  gPrediction->Draw("3CA");
1725  fgSigmaCorr->SetLineColor(kRed);
1726  fgSigmaCorr->SetLineWidth(1);
1727  fgSigmaCorr->SetFillColor(kRed);
1728  fgSigmaCorr->SetFillStyle(0);
1729  fgSigmaCorr->Draw("2");
1730  fhSigmaCorr->SetMarkerColor(kRed);
1731  fhSigmaCorr->Draw("esame");
1732  csigma->SetLogy();
1733  TLegend * leg = new TLegend(0.7,0.75,0.87,0.5);
1734  leg->SetBorderSize(0);
1735  leg->SetLineColor(0);
1736  leg->SetFillColor(0);
1737  leg->SetTextFont(42);
1738  leg->AddEntry(gPrediction,"FONLL ","fl");
1739  leg->AddEntry(fhSigmaCorr,"data stat. unc.","pl");
1740  leg->AddEntry(fgSigmaCorr,"data syst. unc.","f");
1741  leg->Draw();
1742  csigma->Draw();
1743 
1744 }
1745 
1746 //_________________________________________________________________________________________________________
1747 TH1D * AliHFPtSpectrum::ReweightHisto(TH1D *hToReweight, TH1D *hReference){
1753 
1754  // check histograms consistency
1755  Bool_t areconsistent=kTRUE;
1756  areconsistent &= CheckHistosConsistency(hToReweight,hReference);
1757  if (!areconsistent) {
1758  AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
1759  return NULL;
1760  }
1761 
1762  // define a new empty histogram
1763  TH1D *hReweighted = (TH1D*)hToReweight->Clone("hReweighted");
1764  hReweighted->Reset();
1765  Double_t weight=1.0;
1766  Double_t yvalue=1.0;
1767  Double_t integralRef = hReference->Integral();
1768  Double_t integralH = hToReweight->Integral();
1769 
1770  // now reweight the spectra
1771  //
1772  // the weight is the relative probability of the given pt bin in the reference histo
1773  // divided by its relative probability (to normalize it) on the histo to re-weight
1774  for (Int_t i=0; i<=hToReweight->GetNbinsX(); i++) {
1775  weight = (hReference->GetBinContent(i)/integralRef) / (hToReweight->GetBinContent(i)/integralH) ;
1776  yvalue = hToReweight->GetBinContent(i);
1777  hReweighted->SetBinContent(i,yvalue*weight);
1778  }
1779 
1780  return (TH1D*)hReweighted;
1781 }
1782 
1783 //_________________________________________________________________________________________________________
1784 TH1D * AliHFPtSpectrum::ReweightRecHisto(TH1D *hRecToReweight, TH1D *hMCToReweight, TH1D *hMCReference){
1790 
1791  // check histograms consistency
1792  Bool_t areconsistent=kTRUE;
1793  areconsistent &= CheckHistosConsistency(hMCToReweight,hMCReference);
1794  areconsistent &= CheckHistosConsistency(hRecToReweight,hMCReference);
1795  if (!areconsistent) {
1796  AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
1797  return NULL;
1798  }
1799 
1800  // define a new empty histogram
1801  TH1D *hReweighted = (TH1D*)hMCToReweight->Clone("hReweighted");
1802  hReweighted->Reset();
1803  TH1D *hRecReweighted = (TH1D*)hRecToReweight->Clone("hRecReweighted");
1804  hRecReweighted->Reset();
1805  Double_t weight=1.0;
1806  Double_t yvalue=1.0, yrecvalue=1.0;
1807  Double_t integralRef = hMCReference->Integral();
1808  Double_t integralH = hMCToReweight->Integral();
1809 
1810  // now reweight the spectra
1811  //
1812  // the weight is the relative probability of the given pt bin
1813  // that should be applied in the MC histo to get the reference histo shape
1814  // Probabilities are properly normalized.
1815  for (Int_t i=0; i<=hMCToReweight->GetNbinsX(); i++) {
1816  weight = (hMCReference->GetBinContent(i)/integralRef) / (hMCToReweight->GetBinContent(i)/integralH) ;
1817  yvalue = hMCToReweight->GetBinContent(i);
1818  hReweighted->SetBinContent(i,yvalue*weight);
1819  yrecvalue = hRecToReweight->GetBinContent(i);
1820  hRecReweighted->SetBinContent(i,yrecvalue*weight);
1821  }
1822 
1823  return (TH1D*)hRecReweighted;
1824 }
1825 
1826 
1827 
1828 //_________________________________________________________________________________________________________
1833 
1834  Int_t nbins = histo->GetNbinsY();
1835  Int_t ybin=0;
1836  for (int j=0; j<=nbins; j++) {
1837  Float_t value = histo->GetYaxis()->GetBinCenter(j);
1838  Float_t width = histo->GetYaxis()->GetBinWidth(j);
1839  // if( TMath::Abs(yvalue-value)<= width/2. ) {
1840  if( TMath::Abs(yvalue-value)<= width ) {
1841  ybin =j;
1842  // cout <<" value "<<value << ", yval "<< yvalue<<", bin width "<<width/2.<< " y ="<<ybin<<endl;
1843  break;
1844  }
1845  }
1846 
1847  return ybin;
1848 }
1849 
1850 //_________________________________________________________________________________________________________
1852 
1853  Int_t entries = fhDirectEffpt->GetNbinsX();
1854  for(Int_t i=0; i<=entries; i++){
1855  fhDirectEffpt->SetBinError(i,0.);
1856  fhFeedDownEffpt->SetBinError(i,0.);
1857  }
1858 
1859 }
TH2D * fhSigmaCorrRcb
Corrected cross-section (syst. unc. from data only)
TGraphAsymmErrors * fgYieldCorr
Corrected yield (stat unc. only) vs the Ratio(c/b eloss)
void SetMCptDistributionsBounds(TH1D *hDirectMax, TH1D *hDirectMin, TH1D *hFeedDownMax, TH1D *hFeedDownMin)
Set the theoretical direct & feeddown pt spectrum upper and lower bounds.
TH1D * fhStatUncEffbFD
Uncertainty on the feed-down correction due to the prompt efficiency statistical uncertainty.
TH1D * fhStatUncEffcFD
Uncertainty on the cross-section due to the feed-down efficiency statistical uncertainty.
void CalculateFeedDownCorrectedSpectrumFc()
Correct the yield for feed-down correction via fc-method.
double Double_t
Definition: External.C:58
TGraphAsymmErrors * fgYieldCorrExtreme
Corrected yield as TGraphAsymmErrors (syst but feed-down)
const char * title
Definition: MakeQAPdf.C:27
void DrawSpectrum(TGraphAsymmErrors *gPrediction)
Drawing the corrected spectrum comparing to theoretical prediction.
Double_t * fPtBinWidths
Int_t fNevts
all reconstructed D Systematic uncertainties
TGraphAsymmErrors * fgYieldCorrConservative
Extreme corrected yield as TGraphAsymmErrors (syst from feed-down)
void CalculateFeedDownCorrectionFc()
Compute the feed-down correction via fc-method.
Double_t CalculateEfficiencyPtDepedentUncertainty(Double_t pt, Bool_t useOnlyCutVar)
Calculate the efficiency pt-dependent uncertainty.
TNtuple * fnHypothesis
Ntuple of the calculation vs the Ratio(c/b eloss)
void SetDirectAccEffCorrection(TH1D *hDirectEff)
Set the acceptance and efficiency corrections for direct.
void EstimateAndSetDirectEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco)
TH1D * fhDirectMCptMin
Input MC maximum c–>D spectra.
TH1D * fhYieldCorrMax
Corrected yield (stat unc. only)
void EstimateAndSetFeedDownEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco)
TH1D * fhSigmaCorrMin
Maximum corrected cross-section.
Double_t GetCutsEffErr(Double_t pt) const
TH2D * fhYieldCorrRcb
Minimum corrected yield.
void CalculateFeedDownCorrectedSpectrumNb(Double_t deltaY, Double_t branchingRatioBintoFinalDecay)
Correct the yield for feed-down correction via Nb-method.
void SetFeedDownMCptDistributionsBounds(TH1D *hFeedDownMax, TH1D *hFeedDownMin)
Set the theoretical feeddown pt spectrum upper and lower bounds.
Bool_t fPbPbElossHypothesis
flag: asymmetric uncertainties are (1) or not (0) considered
Bool_t fAsymUncertainties
feed-down correction flag: 0=none, 1=fc, 2=Nb
AliHFSystErr * fSystematics
Tab parameter and its uncertainty.
TH1D * ReweightRecHisto(TH1D *hRecToReweight, TH1D *hMCToReweight, TH1D *hMCReference)
to reweight the reco-histos: hRecToReweight is reweighted as hReference/hMCToReweight ...
TH1D * fhDirectMCptMax
Input MC b–>D spectra.
int Int_t
Definition: External.C:63
Int_t fnPtBins
flag : when the analysis is done for In/Out of plane, divide the B-prediction by two ...
TH2D * fhFcRcb
Minimum fc histo.
TH1D * fhFeedDownEffpt
c–>D Acceptance and efficiency correction
float Float_t
Definition: External.C:68
TH1D * fhDirectEffpt
Input MC minimum b–>D spectra.
void ResetStatUncEff()
Reset stat unc on the efficiencies.
Bool_t Initialize()
Initialization.
Double_t GetTrackingEffErr(Double_t pt) const
TH1D * fhFeedDownMCptMin
Input MC maximum b–>D spectra.
Definition: External.C:228
Definition: External.C:212
TGraphAsymmErrors * fgFcExtreme
Correction histo fc vs the Ratio(c/b eloss)
TH1D * fhSigmaCorrDataSyst
Minimum corrected cross-section.
TNtuple * fnSigma
Conservative corrected cross-section as TGraphAsymmErrors (syst from feed-down)
Int_t fFeedDownOption
0=pp, 1=Pb-Pb, 2=p-Pb
TH1D * fhFeedDownMCpt
Input MC c–>D spectra.
void SetReconstructedSpectrum(TH1D *hRec)
Set the reconstructed spectrum.
TGraphAsymmErrors * fgSigmaCorrConservative
Extreme corrected cross-section as TGraphAsymmErrors (syst from feed-down)
Double_t GetTotalSystErr(Double_t pt, Double_t feeddownErr=0) const
TH1D * fhFcMin
Maximum fc histo.
Int_t FindTH2YBin(TH2D *histo, Float_t yvalue)
Functionality to find the y-axis bin of a TH2 for a given y-value.
TH1D * fhFeedDownMCptMax
Input MC minimum c–>D spectra.
void SetReconstructedSpectrumSystematics(TGraphAsymmErrors *gRec)
Bool_t fGlobalEfficiencyPtDependent
uncertainties on the efficiency [0]=c, b, [1]=b/c
TH1D * fhSigmaCorr
Conservative corrected yield as TGraphAsymmErrors (syst from feed-down)
Double_t fTab[2]
use a pt-dependent efficiency uncertainty (Nb method unc.)
TGraphAsymmErrors * fgSigmaCorr
Corrected cross-section (stat unc. only) vs the Ratio(c/b eloss)
AliHFPtSpectrum & operator=(const AliHFPtSpectrum &source)
Assignment operator.
void ComputeSystUncertainties(Bool_t combineFeedDown)
TH1D * fhYieldCorrMin
Maximum corrected yield.
Int_t fCollisionType
Ntuple of the calculation vs the Ratio(c/b eloss)
void CalculateCorrectedSpectrumNoFeedDown()
TH1D * fhYieldCorr
Extreme correction as TGraphAsymmErrors.
void SetFeedDownMCptSpectra(TH1D *hFeedDown)
Set the theoretical feeddown pt spectrum.
virtual ~AliHFPtSpectrum()
Destructor.
TH1D * RebinTheoreticalSpectra(TH1D *hTheory, const char *name)
Function to rebin the theoretical spectra in the data-reconstructed spectra binning.
TGraphAsymmErrors * fgRECSystematics
all reconstructed D
Bool_t CheckHistosConsistency(TH1D *h1, TH1D *h2)
Check histograms consistency function.
TH1D * fhStatUncEffbSigma
Uncertainty on the cross-section due to the prompt efficiency statistical uncertainty.
void ComputeHFPtSpectrum(Double_t deltaY=1.0, Double_t branchingRatioC=1.0, Double_t branchingRatioBintoFinalDecay=1.0)
TGraphAsymmErrors * fgSigmaCorrExtreme
Corrected cross-section as TGraphAsymmErrors (syst but feed-down)
Double_t fTrigEfficiency[2]
analyzed luminosity & uncertainty
void SetAccEffCorrection(TH1D *hDirectEff, TH1D *hFeedDownEff)
Set the acceptance and efficiency corrections for direct & feeddown.
TH1D * fhRECpt
b–>D Acceptance and efficiency correction
TH1D * ReweightHisto(TH1D *hToReweight, TH1D *hReference)
TH1D * fhFcMax
Correction histo fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
Double_t fGlobalEfficiencyUncertainties[2]
trigger efficiency & uncertainty
const Int_t nbins
bool Bool_t
Definition: External.C:53
void DrawErrors(TGraphAsymmErrors *grErrFeeddown=0) const
TGraphAsymmErrors * fgFcConservative
Extreme correction as TGraphAsymmErrors.
Int_t fParticleAntiParticle
flag : consider (1) or not (0) the stat unc on the efficiencies
AliHFPtSpectrum(const char *name="AliHFPtSpectrum", const char *title="HF feed down correction class", Int_t option=1)
Constructor.
TH1D * EstimateEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco, const char *name)
Function to estimate the efficiency in the data-reconstructed spectra binning.
Double_t fLuminosity[2]
nb of analyzed events
void SetMCptSpectra(TH1D *hDirect, TH1D *hFeedDown)
TH1D * fhFc
Systematic uncertainy on the raw yields.
Double_t * fPtBinLimits
number of pt bins
Bool_t fIsStatUncEff
flag: whether to do estimates vs Ratio(c/b eloss) hypothesis
TH1D * fhSigmaCorrMax
Corrected cross-section (stat unc. only)
Bool_t fIsEventPlane
1: only one sign, 2: yield is for particle+anti-particle