1 /**************************************************************************
2 * Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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)
28 // (the corrected yields per bin are divided by the bin-width)
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
35 // Author: Z.Conesa, zconesa@in2p3.fr
36 //***********************************************************************
38 #include <Riostream.h>
46 #include "TGraphAsymmErrors.h"
52 #include "AliHFSystErr.h"
53 #include "AliHFPtSpectrum.h"
55 ClassImp(AliHFPtSpectrum)
57 //_________________________________________________________________________________________________________
58 AliHFPtSpectrum::AliHFPtSpectrum(const char* name, const char* title, Int_t option):
62 fhDirectMCptMax(NULL),
63 fhDirectMCptMin(NULL),
64 fhFeedDownMCptMax(NULL),
65 fhFeedDownMCptMin(NULL),
67 fhFeedDownEffpt(NULL),
69 fgRECSystematics(NULL),
73 fGlobalEfficiencyUncertainties(),
80 fgFcConservative(NULL),
86 fgYieldCorrExtreme(NULL),
87 fgYieldCorrConservative(NULL),
91 fhSigmaCorrDataSyst(NULL),
94 fgSigmaCorrExtreme(NULL),
95 fgSigmaCorrConservative(NULL),
98 fFeedDownOption(option),
99 fAsymUncertainties(kTRUE),
100 fPbPbElossHypothesis(kFALSE),
101 fIsStatUncEff(kTRUE),
102 fParticleAntiParticle(2),
103 fhStatUncEffcSigma(NULL),
104 fhStatUncEffbSigma(NULL),
105 fhStatUncEffcFD(NULL),
106 fhStatUncEffbFD(NULL)
109 // Default constructor
112 fLuminosity[0]=1.; fLuminosity[1]=0.;
113 fTrigEfficiency[0]=1.; fTrigEfficiency[1]=0.;
114 fGlobalEfficiencyUncertainties[0]=0.; fGlobalEfficiencyUncertainties[1]=0.;
115 fTab[0]=1.; fTab[1]=0.;
119 //_________________________________________________________________________________________________________
120 AliHFPtSpectrum::AliHFPtSpectrum(const AliHFPtSpectrum &rhs):
122 fhDirectMCpt(rhs.fhDirectMCpt),
123 fhFeedDownMCpt(rhs.fhFeedDownMCpt),
124 fhDirectMCptMax(rhs.fhDirectMCptMax),
125 fhDirectMCptMin(rhs.fhDirectMCptMin),
126 fhFeedDownMCptMax(rhs.fhFeedDownMCptMax),
127 fhFeedDownMCptMin(rhs.fhFeedDownMCptMin),
128 fhDirectEffpt(rhs.fhDirectEffpt),
129 fhFeedDownEffpt(rhs.fhFeedDownEffpt),
130 fhRECpt(rhs.fhRECpt),
131 fgRECSystematics(rhs.fgRECSystematics),
135 fGlobalEfficiencyUncertainties(),
138 fhFcMax(rhs.fhFcMax),
139 fhFcMin(rhs.fhFcMin),
140 fhFcRcb(rhs.fhFcRcb),
141 fgFcExtreme(rhs.fgFcExtreme),
142 fgFcConservative(rhs.fgFcConservative),
143 fhYieldCorr(rhs.fhYieldCorr),
144 fhYieldCorrMax(rhs.fhYieldCorrMax),
145 fhYieldCorrMin(rhs.fhYieldCorrMin),
146 fhYieldCorrRcb(rhs.fhYieldCorrRcb),
147 fgYieldCorr(rhs.fgYieldCorr),
148 fgYieldCorrExtreme(rhs.fgYieldCorrExtreme),
149 fgYieldCorrConservative(rhs.fgYieldCorrConservative),
150 fhSigmaCorr(rhs.fhSigmaCorr),
151 fhSigmaCorrMax(rhs.fhSigmaCorrMax),
152 fhSigmaCorrMin(rhs.fhSigmaCorrMin),
153 fhSigmaCorrDataSyst(rhs.fhSigmaCorrDataSyst),
154 fhSigmaCorrRcb(rhs.fhSigmaCorrRcb),
155 fgSigmaCorr(rhs.fgSigmaCorr),
156 fgSigmaCorrExtreme(rhs.fgSigmaCorrExtreme),
157 fgSigmaCorrConservative(rhs.fgSigmaCorrConservative),
158 fnSigma(rhs.fnSigma),
159 fnHypothesis(rhs.fnHypothesis),
160 fFeedDownOption(rhs.fFeedDownOption),
161 fAsymUncertainties(rhs.fAsymUncertainties),
162 fPbPbElossHypothesis(rhs.fPbPbElossHypothesis),
163 fIsStatUncEff(rhs.fIsStatUncEff),
164 fParticleAntiParticle(rhs.fParticleAntiParticle),
165 fhStatUncEffcSigma(NULL),
166 fhStatUncEffbSigma(NULL),
167 fhStatUncEffcFD(NULL),
168 fhStatUncEffbFD(NULL)
174 for(Int_t i=0; i<2; i++){
175 fLuminosity[i] = rhs.fLuminosity[i];
176 fTrigEfficiency[i] = rhs.fTrigEfficiency[i];
177 fGlobalEfficiencyUncertainties[i] = rhs.fGlobalEfficiencyUncertainties[i];
178 fTab[i] = rhs.fTab[i];
183 //_________________________________________________________________________________________________________
184 AliHFPtSpectrum &AliHFPtSpectrum::operator=(const AliHFPtSpectrum &source){
186 // Assignment operator
189 if (&source == this) return *this;
191 fhDirectMCpt = source.fhDirectMCpt;
192 fhFeedDownMCpt = source.fhFeedDownMCpt;
193 fhDirectMCptMax = source.fhDirectMCptMax;
194 fhDirectMCptMin = source.fhDirectMCptMin;
195 fhFeedDownMCptMax = source.fhFeedDownMCptMax;
196 fhFeedDownMCptMin = source.fhFeedDownMCptMin;
197 fhDirectEffpt = source.fhDirectEffpt;
198 fhFeedDownEffpt = source.fhFeedDownEffpt;
199 fhRECpt = source.fhRECpt;
200 fgRECSystematics = source.fgRECSystematics;
201 fNevts = source.fNevts;
203 fhFcMax = source.fhFcMax;
204 fhFcMin = source.fhFcMin;
205 fhFcRcb = source.fhFcRcb;
206 fgFcExtreme = source.fgFcExtreme;
207 fgFcConservative = source.fgFcConservative;
208 fhYieldCorr = source.fhYieldCorr;
209 fhYieldCorrMax = source.fhYieldCorrMax;
210 fhYieldCorrMin = source.fhYieldCorrMin;
211 fhYieldCorrRcb = source.fhYieldCorrRcb;
212 fgYieldCorr = source.fgYieldCorr;
213 fgYieldCorrExtreme = source.fgYieldCorrExtreme;
214 fgYieldCorrConservative = source.fgYieldCorrConservative;
215 fhSigmaCorr = source.fhSigmaCorr;
216 fhSigmaCorrMax = source.fhSigmaCorrMax;
217 fhSigmaCorrMin = source.fhSigmaCorrMin;
218 fhSigmaCorrDataSyst = source.fhSigmaCorrDataSyst;
219 fhSigmaCorrRcb = source.fhSigmaCorrRcb;
220 fgSigmaCorr = source.fgSigmaCorr;
221 fgSigmaCorrExtreme = source.fgSigmaCorrExtreme;
222 fgSigmaCorrConservative = source.fgSigmaCorrConservative;
223 fnSigma = source.fnSigma;
224 fnHypothesis = source.fnHypothesis;
225 fFeedDownOption = source.fFeedDownOption;
226 fAsymUncertainties = source.fAsymUncertainties;
227 fPbPbElossHypothesis = source.fPbPbElossHypothesis;
228 fIsStatUncEff = source.fIsStatUncEff;
229 fParticleAntiParticle = source.fParticleAntiParticle;
231 for(Int_t i=0; i<2; i++){
232 fLuminosity[i] = source.fLuminosity[i];
233 fTrigEfficiency[i] = source.fTrigEfficiency[i];
234 fGlobalEfficiencyUncertainties[i] = source.fGlobalEfficiencyUncertainties[i];
235 fTab[i] = source.fTab[i];
241 //_________________________________________________________________________________________________________
242 AliHFPtSpectrum::~AliHFPtSpectrum(){
246 if (fhDirectMCpt) delete fhDirectMCpt;
247 if (fhFeedDownMCpt) delete fhFeedDownMCpt;
248 if (fhDirectMCptMax) delete fhDirectMCptMax;
249 if (fhDirectMCptMin) delete fhDirectMCptMin;
250 if (fhFeedDownMCptMax) delete fhFeedDownMCptMax;
251 if (fhFeedDownMCptMin) delete fhFeedDownMCptMin;
252 if (fhDirectEffpt) delete fhDirectEffpt;
253 if (fhFeedDownEffpt) delete fhFeedDownEffpt;
254 if (fhRECpt) delete fhRECpt;
255 if (fgRECSystematics) delete fgRECSystematics;
256 if (fhFc) delete fhFc;
257 if (fhFcMax) delete fhFcMax;
258 if (fhFcMin) delete fhFcMin;
259 if (fhFcRcb) delete fhFcRcb;
260 if (fgFcExtreme) delete fgFcExtreme;
261 if (fgFcConservative) delete fgFcConservative;
262 if (fhYieldCorr) delete fhYieldCorr;
263 if (fhYieldCorrMax) delete fhYieldCorrMax;
264 if (fhYieldCorrMin) delete fhYieldCorrMin;
265 if (fhYieldCorrRcb) delete fhYieldCorrRcb;
266 if (fgYieldCorr) delete fgYieldCorr;
267 if (fgYieldCorrExtreme) delete fgYieldCorrExtreme;
268 if (fgYieldCorrConservative) delete fgYieldCorrConservative;
269 if (fhSigmaCorr) delete fhSigmaCorr;
270 if (fhSigmaCorrMax) delete fhSigmaCorrMax;
271 if (fhSigmaCorrMin) delete fhSigmaCorrMin;
272 if (fhSigmaCorrDataSyst) delete fhSigmaCorrDataSyst;
273 if (fgSigmaCorr) delete fgSigmaCorr;
274 if (fgSigmaCorrExtreme) delete fgSigmaCorrExtreme;
275 if (fgSigmaCorrConservative) delete fgSigmaCorrConservative;
276 if (fnSigma) delete fnSigma;
277 if (fnHypothesis) delete fnHypothesis;
281 //_________________________________________________________________________________________________________
282 TH1D * AliHFPtSpectrum::RebinTheoreticalSpectra(TH1D *hTheory, const char *name) {
284 // Function to rebin the theoretical spectrum
285 // with respect to the real-data reconstructed spectrum binning
288 if (!hTheory || !fhRECpt) {
289 AliError("Feed-down or reconstructed spectra don't exist");
294 // Get the reconstructed spectra bins & limits
295 Int_t nbins = fhRECpt->GetNbinsX();
296 Int_t nbinsMC = hTheory->GetNbinsX();
297 Double_t *limits = new Double_t[nbins+1];
298 Double_t xlow=0., binwidth=0.;
299 for (Int_t i=1; i<=nbins; i++) {
300 binwidth = fhRECpt->GetBinWidth(i);
301 xlow = fhRECpt->GetBinLowEdge(i);
304 limits[nbins] = xlow + binwidth;
306 // Check that the reconstructed spectra binning
307 // is larger than the theoretical one
308 Double_t thbinwidth = hTheory->GetBinWidth(1);
309 for (Int_t i=1; i<=nbins; i++) {
310 binwidth = fhRECpt->GetBinWidth(i);
311 if ( thbinwidth > binwidth ) {
312 AliInfo(" Beware it seems that the reconstructed spectra has a smaller binning than the theoretical predictions !! ");
317 // Define a new histogram with the real-data reconstructed spectrum binning
318 TH1D * hTheoryRebin = new TH1D(name," theoretical rebinned prediction",nbins,limits);
320 Double_t sum[nbins], items[nbins];
321 for (Int_t ibin=0; ibin<nbins; ibin++) {
322 sum[ibin]=0.; items[ibin]=0.;
324 for (Int_t ibin=0; ibin<=nbinsMC; ibin++){
326 for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++){
327 if (hTheory->GetBinCenter(ibin)>limits[ibinrec] &&
328 hTheory->GetBinCenter(ibin)<limits[ibinrec+1]){
329 sum[ibinrec]+=hTheory->GetBinContent(ibin);
336 // set the theoretical rebinned spectra to ( sum-bins / n-bins ) per new bin
337 for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++) {
338 hTheoryRebin->SetBinContent(ibinrec+1,sum[ibinrec]/items[ibinrec]);
341 return (TH1D*)hTheoryRebin;
344 //_________________________________________________________________________________________________________
345 void AliHFPtSpectrum::SetMCptSpectra(TH1D *hDirect, TH1D *hFeedDown){
347 // Set the MonteCarlo or Theoretical spectra
348 // both for direct and feed-down contributions
351 if (!hDirect || !hFeedDown || !fhRECpt) {
352 AliError("One or both (direct, feed-down) spectra or the reconstructed spectra don't exist");
356 Bool_t areconsistent = kTRUE;
357 areconsistent = CheckHistosConsistency(hDirect,hFeedDown);
358 if (!areconsistent) {
359 AliInfo("Histograms are not consistent (bin width, bounds)");
364 // Rebin the theoretical predictions to the reconstructed spectra binning
366 fhDirectMCpt = RebinTheoreticalSpectra(hDirect,"fhDirectMCpt");
367 fhDirectMCpt->SetNameTitle("fhDirectMCpt"," direct theoretical prediction");
368 fhFeedDownMCpt = RebinTheoreticalSpectra(hFeedDown,"fhFeedDownMCpt");
369 fhFeedDownMCpt->SetNameTitle("fhFeedDownMCpt"," feed-down theoretical prediction");
373 //_________________________________________________________________________________________________________
374 void AliHFPtSpectrum::SetFeedDownMCptSpectra(TH1D *hFeedDown){
376 // Set the MonteCarlo or Theoretical spectra
377 // for feed-down contribution
380 if (!hFeedDown || !fhRECpt) {
381 AliError("Feed-down or reconstructed spectra don't exist");
386 // Rebin the theoretical predictions to the reconstructed spectra binning
388 fhFeedDownMCpt = RebinTheoreticalSpectra(hFeedDown,"fhFeedDownMCpt");
389 fhFeedDownMCpt->SetNameTitle("fhFeedDownMCpt"," feed-down theoretical prediction");
393 //_________________________________________________________________________________________________________
394 void AliHFPtSpectrum::SetMCptDistributionsBounds(TH1D *hDirectMax, TH1D *hDirectMin, TH1D *hFeedDownMax, TH1D *hFeedDownMin){
396 // Set the maximum and minimum MonteCarlo or Theoretical spectra
397 // both for direct and feed-down contributions
398 // used in case uncertainties are asymmetric and ca not be on the "basic histograms"
401 if (!hDirectMax || !hDirectMin || !hFeedDownMax|| !hFeedDownMin || !fhRECpt) {
402 AliError("One or all of the max/min direct/feed-down or the reconstructed spectra don't exist");
406 Bool_t areconsistent = kTRUE;
407 areconsistent &= CheckHistosConsistency(hDirectMax,hDirectMin);
408 areconsistent &= CheckHistosConsistency(hFeedDownMax,hFeedDownMin);
409 areconsistent &= CheckHistosConsistency(hDirectMax,hFeedDownMax);
410 if (!areconsistent) {
411 AliInfo("Histograms are not consistent (bin width, bounds)");
417 // Rebin the theoretical predictions to the reconstructed spectra binning
419 fhDirectMCptMax = RebinTheoreticalSpectra(hDirectMax,"fhDirectMCptMax");
420 fhDirectMCptMax->SetNameTitle("fhDirectMCptMax"," maximum direct theoretical prediction");
421 fhDirectMCptMin = RebinTheoreticalSpectra(hDirectMin,"fhDirectMCptMin");
422 fhDirectMCptMin->SetNameTitle("fhDirectMCptMin"," minimum direct theoretical prediction");
423 fhFeedDownMCptMax = RebinTheoreticalSpectra(hFeedDownMax,"fhFeedDownMCptMax");
424 fhFeedDownMCptMax->SetNameTitle("fhFeedDownMCptMax"," maximum feed-down theoretical prediction");
425 fhFeedDownMCptMin = RebinTheoreticalSpectra(hFeedDownMin,"fhFeedDownMCptMin");
426 fhFeedDownMCptMin->SetNameTitle("fhFeedDownMCptMin"," minimum feed-down theoretical prediction");
430 //_________________________________________________________________________________________________________
431 void AliHFPtSpectrum::SetFeedDownMCptDistributionsBounds(TH1D *hFeedDownMax, TH1D *hFeedDownMin){
433 // Set the maximum and minimum MonteCarlo or Theoretical spectra
434 // for feed-down contributions
435 // used in case uncertainties are asymmetric and can not be on the "basic histogram"
438 if (!hFeedDownMax || !hFeedDownMin || !fhRECpt) {
439 AliError("One or all of the max/min direct/feed-down spectra don't exist");
443 Bool_t areconsistent = kTRUE;
444 areconsistent &= CheckHistosConsistency(hFeedDownMax,hFeedDownMin);
445 if (!areconsistent) {
446 AliInfo("Histograms are not consistent (bin width, bounds)");
452 // Rebin the theoretical predictions to the reconstructed spectra binning
454 fhFeedDownMCptMax = RebinTheoreticalSpectra(hFeedDownMax,"fhFeedDownMCptMax");
455 fhFeedDownMCptMax->SetNameTitle("fhFeedDownMCptMax"," maximum feed-down theoretical prediction");
456 fhFeedDownMCptMin = RebinTheoreticalSpectra(hFeedDownMin,"fhFeedDownMCptMin");
457 fhFeedDownMCptMin->SetNameTitle("fhFeedDownMCptMin"," minimum feed-down theoretical prediction");
461 //_________________________________________________________________________________________________________
462 void AliHFPtSpectrum::SetDirectAccEffCorrection(TH1D *hDirectEff){
464 // Set the Acceptance and Efficiency corrections
465 // for the direct contribution
469 AliError("The direct acceptance and efficiency corrections doesn't exist");
473 fhDirectEffpt = (TH1D*)hDirectEff->Clone();
474 fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance x efficiency correction");
477 //_________________________________________________________________________________________________________
478 void AliHFPtSpectrum::SetAccEffCorrection(TH1D *hDirectEff, TH1D *hFeedDownEff){
480 // Set the Acceptance and Efficiency corrections
481 // both for direct and feed-down contributions
484 if (!hDirectEff || !hFeedDownEff) {
485 AliError("One or both (direct, feed-down) acceptance and efficiency corrections don't exist");
489 Bool_t areconsistent=kTRUE;
490 areconsistent = CheckHistosConsistency(hDirectEff,hFeedDownEff);
491 if (!areconsistent) {
492 AliInfo("Histograms are not consistent (bin width, bounds)");
496 fhDirectEffpt = (TH1D*)hDirectEff->Clone();
497 fhFeedDownEffpt = (TH1D*)hFeedDownEff->Clone();
498 fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance x efficiency correction");
499 fhFeedDownEffpt->SetNameTitle("fhFeedDownEffpt"," feed-down acceptance x efficiency correction");
502 //_________________________________________________________________________________________________________
503 void AliHFPtSpectrum::SetReconstructedSpectrum(TH1D *hRec) {
505 // Set the reconstructed spectrum
509 AliError("The reconstructed spectrum doesn't exist");
513 fhRECpt = (TH1D*)hRec->Clone();
514 fhRECpt->SetNameTitle("fhRECpt"," reconstructed spectrum");
517 //_________________________________________________________________________________________________________
518 void AliHFPtSpectrum::SetReconstructedSpectrumSystematics(TGraphAsymmErrors *gRec) {
520 // Set the reconstructed spectrum (uncorrected yield) systematic uncertainties
523 // Check the compatibility with the reconstructed spectrum
524 Double_t gbinwidth = gRec->GetErrorXlow(1) + gRec->GetErrorXhigh(1) ;
525 Double_t hbinwidth = fhRECpt->GetBinWidth(1);
526 Double_t gxbincenter=0., gybincenter=0.;
527 gRec->GetPoint(1,gxbincenter,gybincenter);
528 Double_t hbincenter = fhRECpt->GetBinCenter(1);
529 if ( (gbinwidth != hbinwidth) || (gxbincenter!=hbincenter) ) {
530 AliError(" The reconstructed spectrum and its systematics don't seem compatible");
534 fgRECSystematics = gRec;
537 //_________________________________________________________________________________________________________
538 void AliHFPtSpectrum::ComputeHFPtSpectrum(Double_t deltaY, Double_t branchingRatioC, Double_t branchingRatioBintoFinalDecay) {
540 // Main function to compute the corrected cross-section:
541 // variables : analysed delta_y, BR for the final correction,
542 // BR b --> D --> decay (relative to the input theoretical prediction)
544 // Sigma = ( 1. / (lumi * delta_y * BR_c * ParticleAntiPartFactor * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
546 // Uncertainties: (stat) delta_sigma = sigma * sqrt ( (delta_spectra/spectra)^2 )
547 // (syst but feed-down) delta_sigma = sigma * sqrt ( (delta_spectra_syst/spectra)^2 + (delta_lumi/lumi)^2 + (delta_eff_trig/eff_trig)^2 + (delta_eff/eff)^2 )
548 // (feed-down syst) delta_sigma = sigma * sqrt ( (delta_spectra_fd/spectra_fd)^2 )
550 // In HIC the feed-down correction varies with an energy loss hypothesis:
551 // Raa(c-->D) / Raa(b-->D) for the "fc" method, Raa(b-->D) for the "Nb" method (see exact formulas in the functions)
555 // First: Initialization
557 Bool_t areHistosOk = Initialize();
559 AliInfo(" Histos not properly initialized. Check : inconsistent binning ? missing histos ?");
562 // Reset the statistical uncertainties on the efficiencies if needed
563 if(!fIsStatUncEff) ResetStatUncEff();
566 // Second: Correct for feed-down
568 if (fFeedDownOption==1) {
569 // Compute the feed-down correction via fc-method
570 CalculateFeedDownCorrectionFc();
571 // Correct the yield for feed-down correction via fc-method
572 CalculateFeedDownCorrectedSpectrumFc();
574 else if (fFeedDownOption==2) {
575 // Correct the yield for feed-down correction via Nb-method
576 CalculateFeedDownCorrectedSpectrumNb(deltaY,branchingRatioBintoFinalDecay);
578 else if (fFeedDownOption==0) {
579 // If there is no need for feed-down correction,
580 // the "corrected" yield is equal to the raw yield
581 fhYieldCorr = (TH1D*)fhRECpt->Clone();
582 fhYieldCorr->SetNameTitle("fhYieldCorr","un-corrected yield");
583 fhYieldCorrMax = (TH1D*)fhRECpt->Clone();
584 fhYieldCorrMin = (TH1D*)fhRECpt->Clone();
585 fhYieldCorrMax->SetNameTitle("fhYieldCorrMax","un-corrected yield");
586 fhYieldCorrMin->SetNameTitle("fhYieldCorrMin","un-corrected yield");
587 fAsymUncertainties=kFALSE;
590 AliInfo(" Are you sure the feed-down correction option is right ?");
594 // Print out information
595 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]);
596 if (fPbPbElossHypothesis) printf("\n\n The considered Tab is %4.2e +- %2.2e \n\n",fTab[0],fTab[1]);
599 // Finally: Correct from yields to cross-section
601 Int_t nbins = fhRECpt->GetNbinsX();
602 Double_t binwidth = fhRECpt->GetBinWidth(1);
603 Double_t *limits = new Double_t[nbins+1];
604 Double_t *binwidths = new Double_t[nbins];
606 for (Int_t i=1; i<=nbins; i++) {
607 binwidth = fhRECpt->GetBinWidth(i);
608 xlow = fhRECpt->GetBinLowEdge(i);
610 binwidths[i-1] = binwidth;
612 limits[nbins] = xlow + binwidth;
615 // declare the output histograms
616 fhSigmaCorr = new TH1D("fhSigmaCorr","corrected sigma",nbins,limits);
617 fhSigmaCorrMax = new TH1D("fhSigmaCorrMax","max corrected sigma",nbins,limits);
618 fhSigmaCorrMin = new TH1D("fhSigmaCorrMin","min corrected sigma",nbins,limits);
619 fhSigmaCorrDataSyst = new TH1D("fhSigmaCorrDataSyst","data syst uncertainties on the corrected sigma",nbins,limits);
620 if (fPbPbElossHypothesis && fFeedDownOption==1) {
621 fhSigmaCorrRcb = new TH2D("fhSigmaCorrRcb","corrected sigma vs Rcb Eloss hypothesis; p_{T} [GeV/c] ; Rcb Eloss hypothesis ; #sigma",nbins,limits,800,0.,4.);
622 fnSigma = new TNtuple("fnSigma"," Sigma ntuple calculation","pt:Signal:Rcb:fc:Yield:Sigma:SigmaStatUnc:SigmaMax:SigmaMin");
624 if (fPbPbElossHypothesis && fFeedDownOption==2) {
625 fhSigmaCorrRcb = new TH2D("fhSigmaCorrRcb","corrected sigma vs Rb Eloss hypothesis; p_{T} [GeV/c] ; Rb Eloss hypothesis ; #sigma",nbins,limits,800,0.,4.);
626 fnSigma = new TNtuple("fnSigma"," Sigma ntuple calculation","pt:Signal:Rb:fc:Yield:Sigma:SigmaStatUnc:SigmaMax:SigmaMin");
628 // and the output TGraphAsymmErrors
629 if (fAsymUncertainties){
630 fgSigmaCorr = new TGraphAsymmErrors(nbins+1);
631 fgSigmaCorrExtreme = new TGraphAsymmErrors(nbins+1);
632 fgSigmaCorrConservative = new TGraphAsymmErrors(nbins+1);
634 fhStatUncEffcSigma = new TH1D("fhStatUncEffcSigma","direct charm stat unc on the cross section",nbins,limits);
635 fhStatUncEffbSigma = new TH1D("fhStatUncEffbSigma","secondary charm stat unc on the cross section",nbins,limits);
638 // protect against null denominator
639 if (deltaY==0. || fLuminosity[0]==0. || fTrigEfficiency[0]==0. || branchingRatioC==0.) {
640 AliError(" Hey you ! Why luminosity or trigger-efficiency or the c-BR or delta_y are set to zero ?! ");
646 Double_t value=0, errValue=0, errvalueMax=0., errvalueMin=0.;
647 Double_t errvalueExtremeMax=0., errvalueExtremeMin=0.;
648 Double_t errvalueConservativeMax=0., errvalueConservativeMin=0.;
649 Double_t errvalueStatUncEffc=0., errvalueStatUncEffb=0.;
650 for(Int_t ibin=1; ibin<=nbins; ibin++){
652 // Variables initialization
653 value=0.; errValue=0.; errvalueMax=0.; errvalueMin=0.;
654 errvalueExtremeMax=0.; errvalueExtremeMin=0.;
655 errvalueConservativeMax=0.; errvalueConservativeMin=0.;
656 errvalueStatUncEffc=0.; errvalueStatUncEffb=0.;
659 // Sigma = ( 1. / (lumi * delta_y * BR_c * ParticleAntiPartFactor * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
660 value = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0. && fhRECpt->GetBinContent(ibin)>0.) ?
661 ( fhYieldCorr->GetBinContent(ibin) / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
664 // Sigma statistical uncertainty:
665 // delta_sigma = sigma * sqrt ( (delta_spectra/spectra)^2 )
666 errValue = (value!=0.) ? value * (fhYieldCorr->GetBinError(ibin)/fhYieldCorr->GetBinContent(ibin)) : 0. ;
668 // cout<< " x "<< fhRECpt->GetBinCenter(ibin) << " sigma " << value << " +- "<< errValue << " (stat)"<<endl;
671 // Sigma systematic uncertainties
673 if (fAsymUncertainties && value>0.) {
675 // (syst but feed-down) delta_sigma = sigma * sqrt ( (delta_spectra_syst/spectra)^2 +
676 // (delta_lumi/lumi)^2 + (delta_eff_trig/eff_trig)^2 + (delta_eff/eff)^2 + (global_eff)^2 )
677 errvalueMax = value * TMath::Sqrt( (fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin))*(fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin)) +
678 (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
679 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) +
680 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
681 fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] );
682 errvalueMin = value * TMath::Sqrt( (fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin))*(fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin)) +
683 (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
684 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) +
685 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
686 fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] );
688 // Uncertainties from feed-down
689 // (feed-down syst) delta_sigma = sigma * sqrt ( (delta_spectra_fd/spectra_fd)^2 )
691 errvalueExtremeMax = value * (fgYieldCorrExtreme->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin));
692 errvalueExtremeMin = value * (fgYieldCorrExtreme->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin));
695 errvalueConservativeMax = value * (fgYieldCorrConservative->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin));
696 errvalueConservativeMin = value * (fgYieldCorrConservative->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin));
699 // stat unc of the efficiencies, separately
700 errvalueStatUncEffc = value * (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) ;
701 errvalueStatUncEffb = 0.;
705 // protect against null denominator
706 errvalueMax = (value!=0.) ?
707 value * TMath::Sqrt( (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
708 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) +
709 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
710 fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] )
712 errvalueMin = errvalueMax;
716 // Fill the histograms
718 fhSigmaCorr->SetBinContent(ibin,value);
719 fhSigmaCorr->SetBinError(ibin,errValue);
722 // Fill the histos and ntuple vs the Eloss hypothesis
724 if (fPbPbElossHypothesis) {
726 // Loop over the Eloss hypothesis
728 AliError("Error reading the fc hypothesis no ntuple, please check !!");
733 Int_t entriesHypo = fnHypothesis->GetEntries();
734 Float_t pt=0., Rhy=0., fc=0., fcMin=0., fcMax=0.;
735 fnHypothesis->SetBranchAddress("pt",&pt);
736 if (fFeedDownOption==2) fnHypothesis->SetBranchAddress("Rb",&Rhy);
737 else if (fFeedDownOption==1) fnHypothesis->SetBranchAddress("Rcb",&Rhy);
738 fnHypothesis->SetBranchAddress("fc",&fc);
739 fnHypothesis->SetBranchAddress("fcMax",&fcMax);
740 fnHypothesis->SetBranchAddress("fcMin",&fcMin);
742 for (Int_t item=0; item<entriesHypo; item++) {
744 fnHypothesis->GetEntry(item);
745 if ( TMath::Abs( pt - fhDirectEffpt->GetBinCenter(ibin) ) > 0.15 ) continue;
747 Double_t yieldRcbvalue = (fhRECpt->GetBinContent(ibin) ) ? fhRECpt->GetBinContent(ibin) * fc : 0. ;
748 yieldRcbvalue /= fhRECpt->GetBinWidth(ibin) ;
749 Double_t yieldRcbvalueMax = (fhRECpt->GetBinContent(ibin) ) ? fhRECpt->GetBinContent(ibin) * fcMax : 0. ;
750 yieldRcbvalueMax /= fhRECpt->GetBinWidth(ibin) ;
751 Double_t yieldRcbvalueMin = (fhRECpt->GetBinContent(ibin) ) ? fhRECpt->GetBinContent(ibin) * fcMin : 0. ;
752 yieldRcbvalueMin /= fhRECpt->GetBinWidth(ibin) ;
755 // Sigma = ( 1. / (lumi * delta_y * BR_c * ParticleAntiPartFactor * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
756 Double_t sigmaRcbvalue = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)>0.) ?
757 ( yieldRcbvalue / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
759 Double_t sigmaRcbvalueMax = (sigmaRcbvalue!=0.) ?
760 ( yieldRcbvalueMax / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
762 Double_t sigmaRcbvalueMin = (sigmaRcbvalue!=0.) ?
763 ( yieldRcbvalueMin / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
765 // Sigma statistical uncertainty:
766 // delta_sigma = sigma * sqrt ( (delta_spectra/spectra)^2 ) = sigma * ( delta_spectra / (spectra-corr * binwidth) )
767 Double_t sigmaRcbvalueStatUnc = (sigmaRcbvalue!=0.) ?
768 sigmaRcbvalue * ( fhRECpt->GetBinError(ibin) / ( yieldRcbvalue * fhRECpt->GetBinWidth(ibin) ) ) : 0. ;
770 fhSigmaCorrRcb->Fill( fhSigmaCorr->GetBinCenter(ibin) , Rhy, sigmaRcbvalue );
772 // cout << " pt "<< fhRECpt->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<" fc-value="<< fhFcRcb->GetBinContent(ibin,rbin) <<", yield-fcRbvalue="<<yieldRcbvalue<<", sigma-fcRbvalue="<<sigmaRcbvalue<<endl;
773 fnSigma->Fill(fhRECpt->GetBinCenter(ibin), fhRECpt->GetBinContent(ibin),
774 Rhy, fc, yieldRcbvalue, sigmaRcbvalue, sigmaRcbvalueStatUnc, sigmaRcbvalueMax, sigmaRcbvalueMin );
778 // Fill the TGraphAsymmErrors
779 if (fAsymUncertainties) {
780 Double_t x = fhYieldCorr->GetBinCenter(ibin);
781 fgSigmaCorr->SetPoint(ibin,x,value); // i,x,y
782 fgSigmaCorr->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
783 fhSigmaCorrMax->SetBinContent(ibin,value+errvalueMax);
784 fhSigmaCorrMin->SetBinContent(ibin,value-errvalueMin);
785 fgSigmaCorrExtreme->SetPoint(ibin,x,value); // i,x,y
786 fgSigmaCorrExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
787 fgSigmaCorrConservative->SetPoint(ibin,x,value); // i,x,y
788 fgSigmaCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueConservativeMin,errvalueConservativeMax); // i,xl,xh,yl,yh
790 fhStatUncEffcSigma->SetBinContent(ibin,0.);
791 if(value>0.) fhStatUncEffcSigma->SetBinError(ibin,((errvalueStatUncEffc/value)*100.));
792 fhStatUncEffbSigma->SetBinContent(ibin,0.); fhStatUncEffbSigma->SetBinError(ibin,0.);
793 // cout << " pt "<< fhRECpt->GetBinCenter(ibin) <<" bin "<< ibin<<" stat-unc-c-sigma "<< errvalueStatUncEffc/value << endl;
802 //_________________________________________________________________________________________________________
803 TH1D * AliHFPtSpectrum::EstimateEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco, const char *name) {
805 // Function that computes the acceptance and efficiency correction
806 // based on the simulated and reconstructed spectra
807 // and using the reconstructed spectra bin width
809 // eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
813 AliInfo("Hey, the reconstructed histogram was not set yet !");
817 Int_t nbins = fhRECpt->GetNbinsX();
818 Double_t *limits = new Double_t[nbins+1];
819 Double_t xlow=0.,binwidth=0.;
820 for (Int_t i=1; i<=nbins; i++) {
821 binwidth = fhRECpt->GetBinWidth(i);
822 xlow = fhRECpt->GetBinLowEdge(i);
825 limits[nbins] = xlow + binwidth;
827 TH1D * hEfficiency = new TH1D(name," acceptance #times efficiency",nbins,limits);
829 Double_t *sumSimu=new Double_t[nbins];
830 Double_t *sumReco=new Double_t[nbins];
831 for (Int_t ibin=0; ibin<nbins; ibin++){
832 sumSimu[ibin]=0.; sumReco[ibin]=0.;
834 for (Int_t ibin=0; ibin<=hSimu->GetNbinsX(); ibin++){
836 for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++){
837 if ( hSimu->GetBinCenter(ibin)>limits[ibinrec] &&
838 hSimu->GetBinCenter(ibin)<limits[ibinrec+1] ) {
839 sumSimu[ibinrec]+=hSimu->GetBinContent(ibin);
841 if ( hReco->GetBinCenter(ibin)>limits[ibinrec] &&
842 hReco->GetBinCenter(ibin)<limits[ibinrec+1] ) {
843 sumReco[ibinrec]+=hReco->GetBinContent(ibin);
850 // the efficiency is computed as reco/sim (in each bin)
851 // its uncertainty is err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
852 Double_t eff=0., erreff=0.;
853 for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++) {
854 if (sumSimu[ibinrec]!= 0. && sumReco[ibinrec]!=0.) {
855 eff = sumReco[ibinrec] / sumSimu[ibinrec] ;
856 // protection in case eff > 1.0
857 // test calculation (make the argument of the sqrt positive)
858 erreff = TMath::Sqrt( eff * TMath::Abs(1.0 - eff) ) / TMath::Sqrt( sumSimu[ibinrec] );
860 else { eff=0.0; erreff=0.; }
861 hEfficiency->SetBinContent(ibinrec+1,eff);
862 hEfficiency->SetBinError(ibinrec+1,erreff);
868 return (TH1D*)hEfficiency;
871 //_________________________________________________________________________________________________________
872 void AliHFPtSpectrum::EstimateAndSetDirectEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco) {
874 // Function that computes the Direct acceptance and efficiency correction
875 // based on the simulated and reconstructed spectra
876 // and using the reconstructed spectra bin width
878 // eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
881 if(!fhRECpt || !hSimu || !hReco){
882 AliError("Hey, the reconstructed histogram was not set yet !");
886 fhDirectEffpt = EstimateEfficiencyRecoBin(hSimu,hReco,"fhDirectEffpt");
887 fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance #times efficiency");
891 //_________________________________________________________________________________________________________
892 void AliHFPtSpectrum::EstimateAndSetFeedDownEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco) {
894 // Function that computes the Feed-Down acceptance and efficiency correction
895 // based on the simulated and reconstructed spectra
896 // and using the reconstructed spectra bin width
898 // eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
901 if(!fhRECpt || !hSimu || !hReco){
902 AliError("Hey, the reconstructed histogram was not set yet !");
906 fhFeedDownEffpt = EstimateEfficiencyRecoBin(hSimu,hReco,"fhFeedDownEffpt");
907 fhFeedDownEffpt->SetNameTitle("fhFeedDownEffpt"," feed-down acceptance #times efficiency");
911 //_________________________________________________________________________________________________________
912 Bool_t AliHFPtSpectrum::Initialize(){
914 // Initialization of the variables (histograms)
917 if (fFeedDownOption==0) {
918 AliInfo("Getting ready for the corrections without feed-down consideration");
919 } else if (fFeedDownOption==1) {
920 AliInfo("Getting ready for the fc feed-down correction calculation");
921 } else if (fFeedDownOption==2) {
922 AliInfo("Getting ready for the Nb feed-down correction calculation");
923 } else { AliError("The calculation option must be <=2"); return kFALSE; }
925 // Start checking the input histograms consistency
926 Bool_t areconsistent=kTRUE;
929 if (!fhDirectEffpt || !fhRECpt) {
930 AliError(" Reconstructed spectra and/or the Nc efficiency distributions are not defined");
933 areconsistent &= CheckHistosConsistency(fhRECpt,fhDirectEffpt);
934 if (!areconsistent) {
935 AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
938 if (fFeedDownOption==0) return kTRUE;
941 // Common checks for options 1 (fc) & 2(Nb)
942 if (!fhFeedDownMCpt || !fhFeedDownEffpt) {
943 AliError(" Theoretical Nb and/or the Nb efficiency distributions are not defined");
946 areconsistent &= CheckHistosConsistency(fhRECpt,fhFeedDownMCpt);
947 areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownEffpt);
948 if (fAsymUncertainties) {
949 if (!fhFeedDownMCptMax || !fhFeedDownMCptMin) {
950 AliError(" Max/Min theoretical Nb distributions are not defined");
953 areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownMCptMax);
955 if (!areconsistent) {
956 AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
959 if (fFeedDownOption>1) return kTRUE;
962 // Now checks for option 1 (fc correction)
964 AliError("Theoretical Nc distributions is not defined");
967 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhFeedDownMCpt);
968 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectEffpt);
969 if (fAsymUncertainties) {
970 if (!fhDirectMCptMax || !fhDirectMCptMin) {
971 AliError(" Max/Min theoretical Nc distributions are not defined");
974 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectMCptMax);
976 if (!areconsistent) {
977 AliInfo("Histograms required for fc correction are not consistent (bin width, bounds)");
984 //_________________________________________________________________________________________________________
985 Bool_t AliHFPtSpectrum::CheckHistosConsistency(TH1D *h1, TH1D *h2){
987 // Check the histograms consistency (bins, limits)
991 AliError("One or both histograms don't exist");
995 Double_t binwidth1 = h1->GetBinWidth(1);
996 Double_t binwidth2 = h2->GetBinWidth(1);
997 Double_t min1 = h1->GetBinCenter(1) - (binwidth1/2.) ;
998 // Double_t max1 = h1->GetBinCenter(nbins1) + (binwidth1/2.) ;
999 Double_t min2 = h2->GetBinCenter(1) - (binwidth2/2.) ;
1000 // Double_t max2 = h2->GetBinCenter(nbins2) + (binwidth2/2.) ;
1002 if (binwidth1!=binwidth2) {
1003 AliInfo(" histograms with different bin width");
1007 AliInfo(" histograms with different minimum");
1010 // if (max1!=max2) {
1011 // AliInfo(" histograms with different maximum");
1018 //_________________________________________________________________________________________________________
1019 void AliHFPtSpectrum::CalculateFeedDownCorrectionFc(){
1021 // Compute fc factor and its uncertainties bin by bin
1022 // fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
1024 // uncertainties: (conservative) combine the upper/lower N_b & N_c predictions together
1025 // (extreme) combine the upper N_b predictions with the lower N_c predictions & viceversa
1026 // systematic uncertainty on the acceptance x efficiency b/c ratio are included
1028 // In addition, in HIC the feed-down correction varies with an energy loss hypothesis: Raa(c-->D) / Raa(b-->D) = Rcb
1029 // fc (Rcb) = ( 1. / ( 1 + (eff_b/eff_c)*(N_b/N_c)* (1/Rcb) ) );
1031 AliInfo("Calculating the feed-down correction factor (fc method)");
1033 // define the variables
1034 Int_t nbins = fhRECpt->GetNbinsX();
1035 Double_t binwidth = fhRECpt->GetBinWidth(1);
1036 Double_t *limits = new Double_t[nbins+1];
1037 Double_t *binwidths = new Double_t[nbins];
1039 for (Int_t i=1; i<=nbins; i++) {
1040 binwidth = fhRECpt->GetBinWidth(i);
1041 xlow = fhRECpt->GetBinLowEdge(i);
1043 binwidths[i-1] = binwidth;
1045 limits[nbins] = xlow + binwidth;
1047 Double_t correction=1.;
1048 Double_t theoryRatio=1.;
1049 Double_t effRatio=1.;
1050 Double_t correctionExtremeA=1., correctionExtremeB=1.;
1051 Double_t theoryRatioExtremeA=1., theoryRatioExtremeB=1.;
1052 Double_t correctionConservativeA=1., correctionConservativeB=1.;
1053 Double_t theoryRatioConservativeA=1., theoryRatioConservativeB=1.;
1054 Double_t correctionUnc=0.;
1055 Double_t correctionExtremeAUnc=0., correctionExtremeBUnc=0.;
1056 Double_t correctionConservativeAUnc=0., correctionConservativeBUnc=0.;
1058 // declare the output histograms
1059 fhFc = new TH1D("fhFc","fc correction factor",nbins,limits);
1060 fhFcMax = new TH1D("fhFcMax","max fc correction factor",nbins,limits);
1061 fhFcMin = new TH1D("fhFcMin","min fc correction factor",nbins,limits);
1062 if(fPbPbElossHypothesis) {
1063 fhFcRcb = new TH2D("fhFcRcb","fc correction factor vs Rcb Eloss hypothesis; p_{T} [GeV/c] ; Rcb Eloss hypothesis ; fc correction",nbins,limits,800,0.,4.);
1064 fnHypothesis = new TNtuple("fnHypothesis"," Feed-down correction vs hypothesis (fc)","pt:Rcb:fc:fcMax:fcMin");
1066 // two local control histograms
1067 TH1D *hTheoryRatio = new TH1D("hTheoryRatio","Theoretical B-->D over c-->D (feed-down/direct) ratio",nbins,limits);
1068 TH1D *hEffRatio = new TH1D("hEffRatio","Efficiency B-->D over c-->D (feed-down/direct) ratio",nbins,limits);
1069 // and the output TGraphAsymmErrors
1070 if (fAsymUncertainties) {
1071 fgFcExtreme = new TGraphAsymmErrors(nbins+1);
1072 fgFcExtreme->SetNameTitle("fgFcExtreme","fgFcExtreme");
1073 fgFcConservative = new TGraphAsymmErrors(nbins+1);
1074 fgFcConservative->SetNameTitle("fgFcConservative","fgFcConservative");
1077 fhStatUncEffcFD = new TH1D("fhStatUncEffcFD","direct charm stat unc on the feed-down correction",nbins,limits);
1078 fhStatUncEffbFD = new TH1D("fhStatUncEffbFD","secondary charm stat unc on the feed-down correction",nbins,limits);
1079 Double_t correctionConservativeAUncStatEffc=0., correctionConservativeBUncStatEffc=0.;
1080 Double_t correctionConservativeAUncStatEffb=0., correctionConservativeBUncStatEffb=0.;
1085 for (Int_t ibin=1; ibin<=nbins; ibin++) {
1087 // theory_ratio = (N_b/N_c)
1088 theoryRatio = (fhDirectMCpt->GetBinContent(ibin)>0. && fhFeedDownMCpt->GetBinContent(ibin)>0.) ?
1089 fhFeedDownMCpt->GetBinContent(ibin) / fhDirectMCpt->GetBinContent(ibin) : 1.0 ;
1092 // Calculate the uncertainty [ considering only the theoretical uncertainties on Nb & Nc for now !!! ]
1094 // extreme A = direct-max, feed-down-min
1095 theoryRatioExtremeA = (fhDirectMCptMax->GetBinContent(ibin)>0. && fhFeedDownMCptMin->GetBinContent(ibin)>0.) ?
1096 fhFeedDownMCptMin->GetBinContent(ibin) / fhDirectMCptMax->GetBinContent(ibin) : 1.0 ;
1097 // extreme B = direct-min, feed-down-max
1098 theoryRatioExtremeB = (fhDirectMCptMin->GetBinContent(ibin)>0. && fhDirectMCptMax->GetBinContent(ibin)>0.) ?
1099 fhFeedDownMCptMax->GetBinContent(ibin) / fhDirectMCptMin->GetBinContent(ibin) : 1.0 ;
1100 // conservative A = direct-max, feed-down-max
1101 theoryRatioConservativeA = (fhDirectMCptMax->GetBinContent(ibin)>0. && fhFeedDownMCptMin->GetBinContent(ibin)>0.) ?
1102 fhFeedDownMCptMax->GetBinContent(ibin) / fhDirectMCptMax->GetBinContent(ibin) : 1.0 ;
1103 // conservative B = direct-min, feed-down-min
1104 theoryRatioConservativeB = (fhDirectMCptMin->GetBinContent(ibin)>0. && fhDirectMCptMax->GetBinContent(ibin)>0.) ?
1105 fhFeedDownMCptMin->GetBinContent(ibin) / fhDirectMCptMin->GetBinContent(ibin) : 1.0 ;
1107 // eff_ratio = (eff_b/eff_c)
1108 effRatio = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ?
1109 fhFeedDownEffpt->GetBinContent(ibin) / fhDirectEffpt->GetBinContent(ibin) : 1.0 ;
1111 // fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
1112 if( TMath::Abs(effRatio - 1.0)<0.01 || TMath::Abs(theoryRatio - 1.0)<0.01 ) {
1114 correctionExtremeA = 1.0;
1115 correctionExtremeB = 1.0;
1116 correctionConservativeA = 1.0;
1117 correctionConservativeB = 1.0;
1120 correction = ( 1. / ( 1 + ( effRatio * theoryRatio ) ) );
1121 correctionExtremeA = ( 1. / ( 1 + ( effRatio * theoryRatioExtremeA ) ) );
1122 correctionExtremeB = ( 1. / ( 1 + ( effRatio * theoryRatioExtremeB ) ) );
1123 correctionConservativeA = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeA ) ) );
1124 correctionConservativeB = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeB ) ) );
1128 // fc uncertainty from (eff_b/eff_c) = fc^2 * (N_b/N_c) * delta(eff_b/eff_c)
1129 // delta(eff_b/eff_c) is a percentage = effRatio * sqrt( fGlobalEfficiencyUncertainties[1]^2 + unc_eff_c ^2 + unc_eff_b ^2 )
1130 Double_t relEffUnc = TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] +
1131 (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
1132 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))
1135 correctionUnc = correction*correction * theoryRatio * effRatio * relEffUnc;
1136 correctionExtremeAUnc = correctionExtremeA*correctionExtremeA * theoryRatioExtremeA * effRatio * relEffUnc;
1137 correctionExtremeBUnc = correctionExtremeB*correctionExtremeB * theoryRatioExtremeB * effRatio * relEffUnc;
1139 correctionConservativeAUnc = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA *effRatio * relEffUnc;
1141 correctionConservativeAUncStatEffc = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA *effRatio *
1142 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin));
1143 correctionConservativeAUncStatEffb = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA *effRatio *
1144 (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin));
1146 correctionConservativeBUnc = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB *effRatio * relEffUnc;
1148 correctionConservativeBUncStatEffb = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB *effRatio *
1149 (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin));
1150 correctionConservativeBUncStatEffc = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB *effRatio *
1151 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin));
1154 // Fill in the histograms
1155 hTheoryRatio->SetBinContent(ibin,theoryRatio);
1156 hEffRatio->SetBinContent(ibin,effRatio);
1157 fhFc->SetBinContent(ibin,correction);
1159 // Estimate how the result varies vs charm/beauty Eloss hypothesis
1161 if ( TMath::Abs(correction-1.0)>0.0001 && fPbPbElossHypothesis){
1162 // Loop over the Eloss hypothesis
1164 for (Float_t rval=0.0025; rval<4.0; rval+=0.005){
1165 // Central fc with Eloss hypothesis
1166 Double_t correctionRcb = ( 1. / ( 1 + ( effRatio * theoryRatio * (1/rval) ) ) );
1167 fhFcRcb->Fill( fhFc->GetBinCenter(ibin) , rval, correctionRcb );
1169 // cout << " pt "<< fhFc->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<", fc-Rcb-value="<<correctionRcb<<endl;
1172 // Upper / lower fc with up / low FONLL bands and Eloss hypothesis
1173 Double_t correctionConservativeARcb = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeA * (1/rval) ) ) );
1174 Double_t correctionConservativeBRcb = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeB * (1/rval) ) ) );
1175 Double_t correctionConservativeARcbUnc = correctionConservativeARcb*correctionConservativeARcb * theoryRatioConservativeA * (1/rval) *effRatio * relEffUnc;
1176 Double_t correctionConservativeBRcbUnc = correctionConservativeBRcb*correctionConservativeBRcb * theoryRatioConservativeB * (1/rval) *effRatio * relEffUnc;
1177 Double_t consvalRcb[4] = { correctionConservativeARcb - correctionConservativeARcbUnc, correctionConservativeARcb + correctionConservativeARcbUnc,
1178 correctionConservativeBRcb - correctionConservativeBRcbUnc, correctionConservativeBRcb + correctionConservativeBRcbUnc};
1179 Double_t uncConservativeRcbMin = correctionRcb - TMath::MinElement(4,consvalRcb);
1180 Double_t uncConservativeRcbMax = TMath::MaxElement(4,consvalRcb) - correctionRcb;
1182 // cout << " pt="<<fhDirectEffpt->GetBinCenter(ibin)<<", hypo="<<rval<<", fc="<<correctionRcb<<" +"<<uncConservativeRcbMax<<" -"<<uncConservativeRcbMin<<endl;
1183 fnHypothesis->Fill( fhDirectEffpt->GetBinCenter(ibin), rval, correctionRcb, correctionRcb+uncConservativeRcbMax, correctionRcb-uncConservativeRcbMin);
1187 // Fill the rest of (asymmetric) histograms
1189 if (fAsymUncertainties) {
1190 Double_t x = fhDirectMCpt->GetBinCenter(ibin);
1191 Double_t val[4] = { correctionExtremeA + correctionExtremeAUnc, correctionExtremeA - correctionExtremeAUnc,
1192 correctionExtremeB + correctionExtremeBUnc, correctionExtremeB - correctionExtremeBUnc };
1193 Double_t uncExtremeMin = correction - TMath::MinElement(4,val);
1194 Double_t uncExtremeMax = TMath::MaxElement(4,val) - correction;
1195 fgFcExtreme->SetPoint(ibin,x,correction); // i,x,y
1196 fgFcExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),uncExtremeMin,uncExtremeMax); // i,xl,xh,yl,yh
1197 fhFcMax->SetBinContent(ibin,correction+uncExtremeMax);
1198 fhFcMin->SetBinContent(ibin,correction-uncExtremeMin);
1199 Double_t consval[4] = { correctionConservativeA - correctionConservativeAUnc, correctionConservativeA + correctionConservativeAUnc,
1200 correctionConservativeB - correctionConservativeBUnc, correctionConservativeB + correctionConservativeBUnc};
1201 Double_t uncConservativeMin = correction - TMath::MinElement(4,consval);
1202 Double_t uncConservativeMax = TMath::MaxElement(4,consval) - correction;
1203 fgFcConservative->SetPoint(ibin,x,correction); // i,x,y
1204 fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),uncConservativeMin,uncConservativeMax); // i,xl,xh,yl,yh
1205 if( !(correction>0.) ){
1206 fgFcExtreme->SetPoint(ibin,x,0.); // i,x,y
1207 fgFcExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.); // i,xl,xh,yl,yh
1208 fgFcConservative->SetPoint(ibin,x,0.); // i,x,y
1209 fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.); // i,xl,xh,yl,yh
1212 Double_t valStatEffc[2] = { correctionConservativeAUncStatEffc/correctionConservativeA,
1213 correctionConservativeBUncStatEffc/correctionConservativeB };
1214 Double_t valStatEffb[2] = { correctionConservativeAUncStatEffb/correctionConservativeA,
1215 correctionConservativeBUncStatEffb/correctionConservativeB };
1216 Double_t uncConservativeStatEffc = TMath::MaxElement(2,valStatEffc);
1217 Double_t uncConservativeStatEffb = TMath::MaxElement(2,valStatEffb);
1218 fhStatUncEffcFD->SetBinContent(ibin,0.); fhStatUncEffcFD->SetBinError(ibin,uncConservativeStatEffc*100.);
1219 fhStatUncEffbFD->SetBinContent(ibin,0.); fhStatUncEffbFD->SetBinError(ibin,uncConservativeStatEffb*100.);
1220 // cout << " pt "<< fhStatUncEffcFD->GetBinCenter(ibin) <<" bin "<< ibin<<" fc-stat-c ="<<uncConservativeStatEffc<<" fc-stat-b ="<<uncConservativeStatEffb<<endl;
1224 delete [] binwidths;
1229 //_________________________________________________________________________________________________________
1230 void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumFc(){
1232 // Compute the feed-down corrected spectrum if feed-down correction is done via fc factor (bin by bin)
1233 // physics = reco * fc / bin-width
1235 // uncertainty: (stat) delta_physics = physics * sqrt ( (delta_reco/reco)^2 )
1236 // (syst but feed-down) delta_physics = physics * sqrt ( (delta_reco_syst/reco)^2 )
1237 // (feed-down syst) delta_physics = physics * sqrt ( (delta_fc/fc)^2 )
1239 // ( Calculation done bin by bin )
1241 // In addition, in HIC the feed-down correction varies with an energy loss hypothesis: Raa(c-->D) / Raa(b-->D) = Rcb
1243 AliInfo(" Calculating the feed-down corrected spectrum (fc method)");
1245 if (!fhFc || !fhRECpt) {
1246 AliError(" Reconstructed or fc distributions are not defined");
1250 Int_t nbins = fhRECpt->GetNbinsX();
1251 Double_t value = 0., errvalue = 0., errvalueMax= 0., errvalueMin= 0.;
1252 Double_t valueExtremeMax= 0., valueExtremeMin= 0.;
1253 Double_t valueConservativeMax= 0., valueConservativeMin= 0.;
1254 Double_t binwidth = fhRECpt->GetBinWidth(1);
1255 Double_t *limits = new Double_t[nbins+1];
1256 Double_t *binwidths = new Double_t[nbins];
1258 for (Int_t i=1; i<=nbins; i++) {
1259 binwidth = fhRECpt->GetBinWidth(i);
1260 xlow = fhRECpt->GetBinLowEdge(i);
1262 binwidths[i-1] = binwidth;
1264 limits[nbins] = xlow + binwidth;
1266 // declare the output histograms
1267 fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by fc)",nbins,limits);
1268 fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by fc)",nbins,limits);
1269 fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by fc)",nbins,limits);
1270 if(fPbPbElossHypothesis) fhYieldCorrRcb = new TH2D("fhYieldCorrRcb","corrected yield (by fc) vs Rcb Eloss hypothesis; p_{T} [GeV/c] ; Rcb Eloss hypothesis ; corrected yield",nbins,limits,800,0.,4.);
1271 // and the output TGraphAsymmErrors
1272 if (fAsymUncertainties){
1273 fgYieldCorr = new TGraphAsymmErrors(nbins+1);
1274 fgYieldCorrExtreme = new TGraphAsymmErrors(nbins+1);
1275 fgYieldCorrConservative = new TGraphAsymmErrors(nbins+1);
1279 // Do the calculation
1281 for (Int_t ibin=1; ibin<=nbins; ibin++) {
1283 // calculate the value
1284 // physics = reco * fc / bin-width
1285 value = (fhRECpt->GetBinContent(ibin) && fhFc->GetBinContent(ibin)) ?
1286 fhRECpt->GetBinContent(ibin) * fhFc->GetBinContent(ibin) : 0. ;
1287 value /= fhRECpt->GetBinWidth(ibin) ;
1289 // Statistical uncertainty
1290 // (stat) delta_physics = physics * sqrt ( (delta_reco/reco)^2 )
1291 errvalue = (value!=0. && fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) ?
1292 value * (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) : 0. ;
1294 // Calculate the systematic uncertainties
1295 // (syst but feed-down) delta_physics = physics * sqrt ( (delta_reco_syst/reco)^2 )
1296 // (feed-down syst) delta_physics = physics * sqrt ( (delta_fc/fc)^2 )
1298 // Protect against null denominator. If so, define uncertainty as null
1299 if (fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) {
1301 if (fAsymUncertainties) {
1303 // Systematics but feed-down
1304 if (fgRECSystematics) {
1305 errvalueMax = value * ( fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinContent(ibin) );
1306 errvalueMin = value * ( fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinContent(ibin) );
1308 else { errvalueMax = 0.; errvalueMin = 0.; }
1310 // Extreme feed-down systematics
1311 valueExtremeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcExtreme->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1312 valueExtremeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcExtreme->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1314 // Conservative feed-down systematics
1315 valueConservativeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcConservative->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1316 valueConservativeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcConservative->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1321 else { errvalueMax = 0.; errvalueMin = 0.; }
1324 // Fill in the histograms
1326 fhYieldCorr->SetBinContent(ibin,value);
1327 fhYieldCorr->SetBinError(ibin,errvalue);
1329 // Fill the histos and ntuple vs the Eloss hypothesis
1331 if (fPbPbElossHypothesis) {
1332 // Loop over the Eloss hypothesis
1333 for (Float_t rval=0.0025; rval<4.0; rval+=0.005){
1334 Int_t rbin = FindTH2YBin(fhYieldCorrRcb,rval);
1335 Double_t fcRcbvalue = fhFcRcb->GetBinContent(ibin,rbin);
1336 // physics = reco * fcRcb / bin-width
1337 Double_t Rcbvalue = (fhRECpt->GetBinContent(ibin) && fcRcbvalue) ?
1338 fhRECpt->GetBinContent(ibin) * fcRcbvalue : 0. ;
1339 Rcbvalue /= fhRECpt->GetBinWidth(ibin) ;
1340 fhYieldCorrRcb->Fill( fhYieldCorr->GetBinCenter(ibin) , rval, Rcbvalue );
1341 // cout << " pt "<< fhRECpt->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<" fc-fcRbvalue="<<fcRcbvalue<<", yield="<<Rcbvalue<<endl;
1344 if (fAsymUncertainties) {
1345 Double_t center = fhYieldCorr->GetBinCenter(ibin);
1346 fgYieldCorr->SetPoint(ibin,center,value); // i,x,y
1347 fgYieldCorr->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
1348 fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax);
1349 fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
1350 fgYieldCorrExtreme->SetPoint(ibin,center,value);
1351 fgYieldCorrExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),value-valueExtremeMin,valueExtremeMax-value);
1352 fgYieldCorrConservative->SetPoint(ibin,center,value);
1353 fgYieldCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),value-valueConservativeMin,valueConservativeMax-value);
1357 delete [] binwidths;
1362 //_________________________________________________________________________________________________________
1363 void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumNb(Double_t deltaY, Double_t branchingRatioBintoFinalDecay) {
1365 // Compute the feed-down corrected spectrum if feed-down correction is done via Nb (bin by bin)
1366 // physics = [ reco - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) ] / bin-width
1368 // uncertainty: (stat) delta_physics = sqrt ( (delta_reco)^2 ) / bin-width
1369 // (syst but feed-down) delta_physics = sqrt ( (delta_reco_syst)^2 ) / bin-width
1370 // (feed-down syst) delta_physics = sqrt ( (k*delta_lumi/lumi)^2 + (k*delta_eff_trig/eff_trig)^2
1371 // + (k*delta_Nb/Nb)^2 + (k*delta_eff/eff)^2 + (k*global_eff_ratio)^2 ) / bin-width
1372 // where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th
1374 // In addition, in HIC the feed-down correction varies with an energy loss hypothesis: Raa(b-->D) = Rb
1375 // physics = [ reco - ( Tab * Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th * Rb ) ] / bin-width
1377 AliInfo("Calculating the feed-down correction factor and spectrum (Nb method)");
1379 Int_t nbins = fhRECpt->GetNbinsX();
1380 Double_t binwidth = fhRECpt->GetBinWidth(1);
1381 Double_t value = 0., errvalue = 0., errvalueMax = 0., errvalueMin = 0., kfactor = 0.;
1382 Double_t errvalueExtremeMax = 0., errvalueExtremeMin = 0.;
1383 Double_t *limits = new Double_t[nbins+1];
1384 Double_t *binwidths = new Double_t[nbins];
1386 for (Int_t i=1; i<=nbins; i++) {
1387 binwidth = fhRECpt->GetBinWidth(i);
1388 xlow = fhRECpt->GetBinLowEdge(i);
1390 binwidths[i-1] = binwidth;
1392 limits[nbins] = xlow + binwidth;
1394 // declare the output histograms
1395 fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by Nb)",nbins,limits);
1396 fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by Nb)",nbins,limits);
1397 fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by Nb)",nbins,limits);
1398 if(fPbPbElossHypothesis) {
1399 fhFcRcb = new TH2D("fhFcRcb","fc correction factor (Nb method) vs Rb Eloss hypothesis; p_{T} [GeV/c] ; Rb Eloss hypothesis ; fc correction",nbins,limits,800,0.,4.);
1400 fhYieldCorrRcb = new TH2D("fhYieldCorrRcb","corrected yield (by Nb) vs Rb Eloss hypothesis; p_{T} [GeV/c] ; Rb Eloss hypothesis ; corrected yield",nbins,limits,800,0.,4.);
1401 fnHypothesis = new TNtuple("fnHypothesis"," Feed-down correction vs hypothesis (Nb)","pt:Rb:fc:fcMax:fcMin");
1403 // and the output TGraphAsymmErrors
1404 if (fAsymUncertainties){
1405 fgYieldCorr = new TGraphAsymmErrors(nbins+1);
1406 fgYieldCorrExtreme = new TGraphAsymmErrors(nbins+1);
1407 fgYieldCorrConservative = new TGraphAsymmErrors(nbins+1);
1408 // Define fc-conservative
1409 fgFcConservative = new TGraphAsymmErrors(nbins+1);
1410 AliInfo(" Beware the conservative & extreme uncertainties are equal by definition !");
1413 // variables to define fc-conservative
1414 double correction=0, correctionMax=0., correctionMin=0.;
1416 fhStatUncEffcFD = new TH1D("fhStatUncEffcFD","direct charm stat unc on the feed-down correction",nbins,limits);
1417 fhStatUncEffbFD = new TH1D("fhStatUncEffbFD","secondary charm stat unc on the feed-down correction",nbins,limits);
1418 Double_t correctionUncStatEffc=0.;
1419 Double_t correctionUncStatEffb=0.;
1423 // Do the calculation
1425 for (Int_t ibin=1; ibin<=nbins; ibin++) {
1427 // Calculate the value
1428 // physics = [ reco - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) ] / bin-width
1429 // In HIC : physics = [ reco - ( Tab * Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th * Rb ) ] / bin-width
1432 Double_t frac = 1.0, errfrac =0.;
1433 if(fPbPbElossHypothesis) {
1434 frac = fTab[0]*fNevts;
1435 errfrac = frac * TMath::Sqrt( (fTab[1]/fTab[0])*(fTab[1]/fTab[0]) + (1/fNevts) );
1437 frac = fLuminosity[0];
1438 errfrac = fLuminosity[1];
1441 value = ( fhRECpt->GetBinContent(ibin)>0. && fhRECpt->GetBinContent(ibin)!=0. &&
1442 fhFeedDownMCpt->GetBinContent(ibin)>0. && fhFeedDownEffpt->GetBinContent(ibin)>0. ) ?
1443 fhRECpt->GetBinContent(ibin) - frac*(deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) * fhRECpt->GetBinWidth(ibin) )
1445 value /= fhRECpt->GetBinWidth(ibin);
1446 if (value<0.) value =0.;
1448 // Statistical uncertainty: delta_physics = sqrt ( (delta_reco)^2 ) / bin-width
1449 errvalue = (value!=0. && fhRECpt->GetBinError(ibin) && fhRECpt->GetBinError(ibin)!=0.) ?
1450 fhRECpt->GetBinError(ibin) : 0.;
1451 errvalue /= fhRECpt->GetBinWidth(ibin);
1453 // Correction (fc) : Estimate of the relative amount feed-down subtracted
1454 // correction = [ 1 - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th)/reco ]
1455 // in HIC: correction = [ 1 - ( Tab * Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th)/reco ]
1456 correction = (value>0.) ?
1457 1 - (frac*deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) * fhRECpt->GetBinWidth(ibin) ) / fhRECpt->GetBinContent(ibin) : 0. ;
1458 if (correction<0.) correction = 0.;
1460 // Systematic uncertainties
1461 // (syst but feed-down) delta_physics = sqrt ( (delta_reco_syst)^2 ) / bin-width
1462 // (feed-down syst) delta_physics = sqrt ( (k*delta_lumi/lumi)^2 + (k*delta_eff_trig/eff_trig)^2
1463 // + (k*delta_Nb/Nb)^2 + (k*delta_eff/eff)^2 + (k*global_eff_ratio)^2 ) / bin-width
1464 // where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th * bin-width
1465 kfactor = frac*deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) * fhRECpt->GetBinWidth(ibin) ;
1467 if (fAsymUncertainties && value>0.) {
1468 Double_t nb = fhFeedDownMCpt->GetBinContent(ibin);
1469 Double_t nbDmax = fhFeedDownMCptMax->GetBinContent(ibin) - fhFeedDownMCpt->GetBinContent(ibin);
1470 Double_t nbDmin = fhFeedDownMCpt->GetBinContent(ibin) - fhFeedDownMCptMin->GetBinContent(ibin);
1472 // Systematics but feed-down
1473 if (fgRECSystematics){
1474 errvalueMax = fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinWidth(ibin) ;
1475 errvalueMin = fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinWidth(ibin);
1477 else { errvalueMax = 0.; errvalueMin = 0.; }
1479 // Feed-down systematics
1480 // min value with the maximum Nb
1481 Double_t errCom = ( (kfactor*errfrac/frac)*(kfactor*errfrac/frac) ) +
1482 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1483 ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
1484 ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) ) ;
1485 errvalueExtremeMin = TMath::Sqrt( errCom + ( (kfactor*nbDmax/nb)*(kfactor*nbDmax/nb) ) ) / fhRECpt->GetBinWidth(ibin);
1486 // max value with the minimum Nb
1487 errvalueExtremeMax = TMath::Sqrt( errCom + ( (kfactor*nbDmin/nb)*(kfactor*nbDmin/nb) ) ) / fhRECpt->GetBinWidth(ibin);
1489 // Correction systematics (fc)
1490 // min value with the maximum Nb
1491 correctionMin = TMath::Sqrt( errCom + ( (kfactor*nbDmax/nb)*(kfactor*nbDmax/nb) ) ) / fhRECpt->GetBinContent(ibin) ;
1492 // max value with the minimum Nb
1493 correctionMax = TMath::Sqrt( errCom + ( (kfactor*nbDmin/nb)*(kfactor*nbDmin/nb) ) ) / fhRECpt->GetBinContent(ibin) ;
1495 correctionUncStatEffb = TMath::Sqrt( ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) )
1496 ) / fhRECpt->GetBinContent(ibin) ;
1497 correctionUncStatEffc = 0.;
1499 else{ // Don't consider Nb uncertainty in this case [ to be tested!!! ]
1500 errvalueExtremeMax = TMath::Sqrt( ( (kfactor*errfrac/frac)*(kfactor*errfrac/frac) ) +
1501 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1502 ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
1503 ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1504 ) / fhRECpt->GetBinWidth(ibin);
1505 errvalueExtremeMin = errvalueExtremeMax ;
1509 // fill in histograms
1510 fhYieldCorr->SetBinContent(ibin,value);
1511 fhYieldCorr->SetBinError(ibin,errvalue);
1513 // Estimate how the result varies vs charm/beauty Eloss hypothesis
1515 if ( correction>0.0001 && fPbPbElossHypothesis){
1516 // Loop over the Eloss hypothesis
1518 for (Float_t rval=0.0025; rval<4.0; rval+=0.005){
1519 // correction = [ 1 - (Tab *Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th *binwidth )* (rval) /reco ]
1520 Double_t fcRcbvalue = 1 - (fTab[0]*fNevts*deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin)*fhRECpt->GetBinWidth(ibin) * rval ) / fhRECpt->GetBinContent(ibin) ;
1521 if(fcRcbvalue<0.) fcRcbvalue=0.;
1522 fhFcRcb->Fill( fhRECpt->GetBinCenter(ibin) , rval, fcRcbvalue );
1523 // physics = reco * fcRcb / bin-width
1524 Double_t Rcbvalue = (fhRECpt->GetBinContent(ibin) && fcRcbvalue) ?
1525 fhRECpt->GetBinContent(ibin) * fcRcbvalue : 0. ;
1526 Rcbvalue /= fhRECpt->GetBinWidth(ibin) ;
1527 fhYieldCorrRcb->Fill( fhYieldCorr->GetBinCenter(ibin) , rval, Rcbvalue );
1529 // cout << " pt "<< fhFcRcb->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<", fc-Rb-value="<< fcRcbvalue << ", yield-Rb-value="<< Rcbvalue <<endl;
1530 // cout << " pt "<< fhFcRcb->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", fc-Rb-value="<< fcRcbvalue << ", yield-Rb-value="<< Rcbvalue <<endl;
1533 Double_t correctionMaxRcb = correctionMax*rval;
1534 Double_t correctionMinRcb = correctionMin*rval;
1535 fnHypothesis->Fill( fhYieldCorr->GetBinCenter(ibin), rval, fcRcbvalue, fcRcbvalue + correctionMaxRcb, fcRcbvalue - correctionMinRcb);
1537 // cout << " pt="<< fhFcRcb->GetBinCenter(ibin) <<", rval="<<rval<<", fc="<<fcRcbvalue<<" +"<<correctionMaxRcb<<" -"<<correctionMinRcb<<endl;}
1541 // Fill the rest of (asymmetric) histograms
1543 if (fAsymUncertainties) {
1544 Double_t x = fhYieldCorr->GetBinCenter(ibin);
1545 fgYieldCorr->SetPoint(ibin,x,value); // i,x,y
1546 fgYieldCorr->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
1547 fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax);
1548 fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
1549 fgYieldCorrExtreme->SetPoint(ibin,x,value); // i,x,y
1550 fgYieldCorrExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
1551 fgYieldCorrConservative->SetPoint(ibin,x,value); // i,x,y
1552 fgYieldCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
1553 // cout << " bin " << ibin << ", correction " << correction << ", min correction unc " << correctionMin << ", max correction unc " << correctionMax << endl;
1555 fgFcConservative->SetPoint(ibin,x,correction);
1556 fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),correctionMin,correctionMax);
1558 fhStatUncEffbFD->SetBinContent(ibin,0.); fhStatUncEffbFD->SetBinError(ibin,correctionUncStatEffb/correction*100.);
1559 fhStatUncEffcFD->SetBinContent(ibin,0.); fhStatUncEffcFD->SetBinError(ibin,correctionUncStatEffc/correction*100.);
1560 // cout << " pt "<< fhStatUncEffcFD->GetBinCenter(ibin) <<" bin "<< ibin<<" fc-stat-c ="<< correctionUncStatEffc/correction <<" fc-stat-b ="<< correctionUncStatEffb/correction <<endl;
1563 fgFcConservative->SetPoint(ibin,x,0.);
1564 fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.);
1569 delete [] binwidths;
1575 //_________________________________________________________________________________________________________
1576 void AliHFPtSpectrum::ComputeSystUncertainties(AliHFSystErr *systematics, Bool_t combineFeedDown) {
1578 // Function that re-calculates the global systematic uncertainties
1579 // by calling the class AliHFSystErr and combining those
1580 // (in quadrature) with the feed-down subtraction uncertainties
1583 // Estimate the feed-down uncertainty in percentage
1584 Int_t nentries = fgSigmaCorrConservative->GetN();
1585 TGraphAsymmErrors *grErrFeeddown = new TGraphAsymmErrors(nentries);
1586 Double_t x=0., y=0., errx=0., erryl=0., erryh=0;
1587 for(Int_t i=0; i<nentries; i++) {
1588 x=0.; y=0.; errx=0.; erryl=0.; erryh=0.;
1589 fgSigmaCorrConservative->GetPoint(i,x,y);
1591 errx = fgSigmaCorrConservative->GetErrorXlow(i) ;
1592 erryl = fgSigmaCorrConservative->GetErrorYlow(i) / y ;
1593 erryh = fgSigmaCorrConservative->GetErrorYhigh(i) / y ;
1595 // cout << " x "<< x << " +- "<<errx<<" , y "<<y<<" + "<<erryh<<" - "<<erryl<<endl;
1596 grErrFeeddown->SetPoint(i,x,0.);
1597 grErrFeeddown->SetPointError(i,errx,errx,erryl,erryh); //i, xl, xh, yl, yh
1600 // Draw all the systematics independently
1601 systematics->DrawErrors(grErrFeeddown);
1603 // Set the sigma systematic uncertainties
1604 // possibly combine with the feed-down uncertainties
1605 Double_t errylcomb=0., erryhcomb=0;
1606 for(Int_t i=1; i<nentries; i++) {
1607 fgSigmaCorr->GetPoint(i,x,y);
1608 errx = grErrFeeddown->GetErrorXlow(i) ;
1609 erryl = grErrFeeddown->GetErrorYlow(i);
1610 erryh = grErrFeeddown->GetErrorYhigh(i);
1611 if (combineFeedDown) {
1612 errylcomb = systematics->GetTotalSystErr(x,erryl) * y ;
1613 erryhcomb = systematics->GetTotalSystErr(x,erryh) * y ;
1615 errylcomb = systematics->GetTotalSystErr(x) * y ;
1616 erryhcomb = systematics->GetTotalSystErr(x) * y ;
1618 fgSigmaCorr->SetPointError(i,errx,errx,errylcomb,erryhcomb);
1620 fhSigmaCorrDataSyst->SetBinContent(i,y);
1621 erryl = systematics->GetTotalSystErr(x) * y ;
1622 fhSigmaCorrDataSyst->SetBinError(i,erryl);
1628 //_________________________________________________________________________________________________________
1629 void AliHFPtSpectrum::DrawSpectrum(TGraphAsymmErrors *gPrediction) {
1631 // Example method to draw the corrected spectrum & the theoretical prediction
1634 TCanvas *csigma = new TCanvas("csigma","Draw the corrected cross-section & the prediction");
1635 csigma->SetFillColor(0);
1636 gPrediction->GetXaxis()->SetTitleSize(0.05);
1637 gPrediction->GetXaxis()->SetTitleOffset(0.95);
1638 gPrediction->GetYaxis()->SetTitleSize(0.05);
1639 gPrediction->GetYaxis()->SetTitleOffset(0.95);
1640 gPrediction->GetXaxis()->SetTitle("p_{T} [GeV]");
1641 gPrediction->GetYaxis()->SetTitle("BR #times #frac{d#sigma}{dp_{T}} |_{|y|<0.5} [pb/GeV]");
1642 gPrediction->SetLineColor(kGreen+2);
1643 gPrediction->SetLineWidth(3);
1644 gPrediction->SetFillColor(kGreen+1);
1645 gPrediction->Draw("3CA");
1646 fgSigmaCorr->SetLineColor(kRed);
1647 fgSigmaCorr->SetLineWidth(1);
1648 fgSigmaCorr->SetFillColor(kRed);
1649 fgSigmaCorr->SetFillStyle(0);
1650 fgSigmaCorr->Draw("2");
1651 fhSigmaCorr->SetMarkerColor(kRed);
1652 fhSigmaCorr->Draw("esame");
1654 TLegend * leg = new TLegend(0.7,0.75,0.87,0.5);
1655 leg->SetBorderSize(0);
1656 leg->SetLineColor(0);
1657 leg->SetFillColor(0);
1658 leg->SetTextFont(42);
1659 leg->AddEntry(gPrediction,"FONLL ","fl");
1660 leg->AddEntry(fhSigmaCorr,"data stat. unc.","pl");
1661 leg->AddEntry(fgSigmaCorr,"data syst. unc.","f");
1667 //_________________________________________________________________________________________________________
1668 TH1D * AliHFPtSpectrum::ReweightHisto(TH1D *hToReweight, TH1D *hReference){
1670 // Function to reweight histograms for testing purposes:
1671 // This function takes the histo hToReweight and reweights
1672 // it (its pt shape) with respect to hReference
1675 // check histograms consistency
1676 Bool_t areconsistent=kTRUE;
1677 areconsistent &= CheckHistosConsistency(hToReweight,hReference);
1678 if (!areconsistent) {
1679 AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
1683 // define a new empty histogram
1684 TH1D *hReweighted = (TH1D*)hToReweight->Clone("hReweighted");
1685 hReweighted->Reset();
1686 Double_t weight=1.0;
1687 Double_t yvalue=1.0;
1688 Double_t integralRef = hReference->Integral();
1689 Double_t integralH = hToReweight->Integral();
1691 // now reweight the spectra
1693 // the weight is the relative probability of the given pt bin in the reference histo
1694 // divided by its relative probability (to normalize it) on the histo to re-weight
1695 for (Int_t i=0; i<=hToReweight->GetNbinsX(); i++) {
1696 weight = (hReference->GetBinContent(i)/integralRef) / (hToReweight->GetBinContent(i)/integralH) ;
1697 yvalue = hToReweight->GetBinContent(i);
1698 hReweighted->SetBinContent(i,yvalue*weight);
1701 return (TH1D*)hReweighted;
1704 //_________________________________________________________________________________________________________
1705 TH1D * AliHFPtSpectrum::ReweightRecHisto(TH1D *hRecToReweight, TH1D *hMCToReweight, TH1D *hMCReference){
1707 // Function to reweight histograms for testing purposes:
1708 // This function takes the histo hToReweight and reweights
1709 // it (its pt shape) with respect to hReference /hMCToReweight
1712 // check histograms consistency
1713 Bool_t areconsistent=kTRUE;
1714 areconsistent &= CheckHistosConsistency(hMCToReweight,hMCReference);
1715 areconsistent &= CheckHistosConsistency(hRecToReweight,hMCReference);
1716 if (!areconsistent) {
1717 AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
1721 // define a new empty histogram
1722 TH1D *hReweighted = (TH1D*)hMCToReweight->Clone("hReweighted");
1723 hReweighted->Reset();
1724 TH1D *hRecReweighted = (TH1D*)hRecToReweight->Clone("hRecReweighted");
1725 hRecReweighted->Reset();
1726 Double_t weight=1.0;
1727 Double_t yvalue=1.0, yrecvalue=1.0;
1728 Double_t integralRef = hMCReference->Integral();
1729 Double_t integralH = hMCToReweight->Integral();
1731 // now reweight the spectra
1733 // the weight is the relative probability of the given pt bin
1734 // that should be applied in the MC histo to get the reference histo shape
1735 // Probabilities are properly normalized.
1736 for (Int_t i=0; i<=hMCToReweight->GetNbinsX(); i++) {
1737 weight = (hMCReference->GetBinContent(i)/integralRef) / (hMCToReweight->GetBinContent(i)/integralH) ;
1738 yvalue = hMCToReweight->GetBinContent(i);
1739 hReweighted->SetBinContent(i,yvalue*weight);
1740 yrecvalue = hRecToReweight->GetBinContent(i);
1741 hRecReweighted->SetBinContent(i,yrecvalue*weight);
1744 return (TH1D*)hRecReweighted;
1749 //_________________________________________________________________________________________________________
1750 Int_t AliHFPtSpectrum::FindTH2YBin(TH2D *histo, Float_t yvalue){
1752 // Function to find the y-axis bin of a TH2 for a given y-value
1755 Int_t nbins = histo->GetNbinsY();
1757 for (int j=0; j<=nbins; j++) {
1758 Float_t value = histo->GetYaxis()->GetBinCenter(j);
1759 Float_t width = histo->GetYaxis()->GetBinWidth(j);
1760 // if( TMath::Abs(yvalue-value)<= width/2. ) {
1761 if( TMath::Abs(yvalue-value)<= width ) {
1763 // cout <<" value "<<value << ", yval "<< yvalue<<", bin width "<<width/2.<< " y ="<<ybin<<endl;
1771 //_________________________________________________________________________________________________________
1772 void AliHFPtSpectrum::ResetStatUncEff(){
1774 Int_t entries = fhDirectEffpt->GetNbinsX();
1775 for(Int_t i=0; i<=entries; i++){
1776 fhDirectEffpt->SetBinError(i,0.);
1777 fhFeedDownEffpt->SetBinError(i,0.);