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 fIsEventPlane(kFALSE),
104 fhStatUncEffcSigma(NULL),
105 fhStatUncEffbSigma(NULL),
106 fhStatUncEffcFD(NULL),
107 fhStatUncEffbFD(NULL)
110 // Default constructor
113 fLuminosity[0]=1.; fLuminosity[1]=0.;
114 fTrigEfficiency[0]=1.; fTrigEfficiency[1]=0.;
115 fGlobalEfficiencyUncertainties[0]=0.; fGlobalEfficiencyUncertainties[1]=0.;
116 fTab[0]=1.; fTab[1]=0.;
120 //_________________________________________________________________________________________________________
121 AliHFPtSpectrum::AliHFPtSpectrum(const AliHFPtSpectrum &rhs):
123 fhDirectMCpt(rhs.fhDirectMCpt),
124 fhFeedDownMCpt(rhs.fhFeedDownMCpt),
125 fhDirectMCptMax(rhs.fhDirectMCptMax),
126 fhDirectMCptMin(rhs.fhDirectMCptMin),
127 fhFeedDownMCptMax(rhs.fhFeedDownMCptMax),
128 fhFeedDownMCptMin(rhs.fhFeedDownMCptMin),
129 fhDirectEffpt(rhs.fhDirectEffpt),
130 fhFeedDownEffpt(rhs.fhFeedDownEffpt),
131 fhRECpt(rhs.fhRECpt),
132 fgRECSystematics(rhs.fgRECSystematics),
136 fGlobalEfficiencyUncertainties(),
139 fhFcMax(rhs.fhFcMax),
140 fhFcMin(rhs.fhFcMin),
141 fhFcRcb(rhs.fhFcRcb),
142 fgFcExtreme(rhs.fgFcExtreme),
143 fgFcConservative(rhs.fgFcConservative),
144 fhYieldCorr(rhs.fhYieldCorr),
145 fhYieldCorrMax(rhs.fhYieldCorrMax),
146 fhYieldCorrMin(rhs.fhYieldCorrMin),
147 fhYieldCorrRcb(rhs.fhYieldCorrRcb),
148 fgYieldCorr(rhs.fgYieldCorr),
149 fgYieldCorrExtreme(rhs.fgYieldCorrExtreme),
150 fgYieldCorrConservative(rhs.fgYieldCorrConservative),
151 fhSigmaCorr(rhs.fhSigmaCorr),
152 fhSigmaCorrMax(rhs.fhSigmaCorrMax),
153 fhSigmaCorrMin(rhs.fhSigmaCorrMin),
154 fhSigmaCorrDataSyst(rhs.fhSigmaCorrDataSyst),
155 fhSigmaCorrRcb(rhs.fhSigmaCorrRcb),
156 fgSigmaCorr(rhs.fgSigmaCorr),
157 fgSigmaCorrExtreme(rhs.fgSigmaCorrExtreme),
158 fgSigmaCorrConservative(rhs.fgSigmaCorrConservative),
159 fnSigma(rhs.fnSigma),
160 fnHypothesis(rhs.fnHypothesis),
161 fFeedDownOption(rhs.fFeedDownOption),
162 fAsymUncertainties(rhs.fAsymUncertainties),
163 fPbPbElossHypothesis(rhs.fPbPbElossHypothesis),
164 fIsStatUncEff(rhs.fIsStatUncEff),
165 fParticleAntiParticle(rhs.fParticleAntiParticle),
166 fIsEventPlane(rhs.fIsEventPlane),
167 fhStatUncEffcSigma(NULL),
168 fhStatUncEffbSigma(NULL),
169 fhStatUncEffcFD(NULL),
170 fhStatUncEffbFD(NULL)
176 for(Int_t i=0; i<2; i++){
177 fLuminosity[i] = rhs.fLuminosity[i];
178 fTrigEfficiency[i] = rhs.fTrigEfficiency[i];
179 fGlobalEfficiencyUncertainties[i] = rhs.fGlobalEfficiencyUncertainties[i];
180 fTab[i] = rhs.fTab[i];
185 //_________________________________________________________________________________________________________
186 AliHFPtSpectrum &AliHFPtSpectrum::operator=(const AliHFPtSpectrum &source){
188 // Assignment operator
191 if (&source == this) return *this;
193 fhDirectMCpt = source.fhDirectMCpt;
194 fhFeedDownMCpt = source.fhFeedDownMCpt;
195 fhDirectMCptMax = source.fhDirectMCptMax;
196 fhDirectMCptMin = source.fhDirectMCptMin;
197 fhFeedDownMCptMax = source.fhFeedDownMCptMax;
198 fhFeedDownMCptMin = source.fhFeedDownMCptMin;
199 fhDirectEffpt = source.fhDirectEffpt;
200 fhFeedDownEffpt = source.fhFeedDownEffpt;
201 fhRECpt = source.fhRECpt;
202 fgRECSystematics = source.fgRECSystematics;
203 fNevts = source.fNevts;
205 fhFcMax = source.fhFcMax;
206 fhFcMin = source.fhFcMin;
207 fhFcRcb = source.fhFcRcb;
208 fgFcExtreme = source.fgFcExtreme;
209 fgFcConservative = source.fgFcConservative;
210 fhYieldCorr = source.fhYieldCorr;
211 fhYieldCorrMax = source.fhYieldCorrMax;
212 fhYieldCorrMin = source.fhYieldCorrMin;
213 fhYieldCorrRcb = source.fhYieldCorrRcb;
214 fgYieldCorr = source.fgYieldCorr;
215 fgYieldCorrExtreme = source.fgYieldCorrExtreme;
216 fgYieldCorrConservative = source.fgYieldCorrConservative;
217 fhSigmaCorr = source.fhSigmaCorr;
218 fhSigmaCorrMax = source.fhSigmaCorrMax;
219 fhSigmaCorrMin = source.fhSigmaCorrMin;
220 fhSigmaCorrDataSyst = source.fhSigmaCorrDataSyst;
221 fhSigmaCorrRcb = source.fhSigmaCorrRcb;
222 fgSigmaCorr = source.fgSigmaCorr;
223 fgSigmaCorrExtreme = source.fgSigmaCorrExtreme;
224 fgSigmaCorrConservative = source.fgSigmaCorrConservative;
225 fnSigma = source.fnSigma;
226 fnHypothesis = source.fnHypothesis;
227 fFeedDownOption = source.fFeedDownOption;
228 fAsymUncertainties = source.fAsymUncertainties;
229 fPbPbElossHypothesis = source.fPbPbElossHypothesis;
230 fIsStatUncEff = source.fIsStatUncEff;
231 fParticleAntiParticle = source.fParticleAntiParticle;
232 fIsEventPlane = source.fIsEventPlane;
234 for(Int_t i=0; i<2; i++){
235 fLuminosity[i] = source.fLuminosity[i];
236 fTrigEfficiency[i] = source.fTrigEfficiency[i];
237 fGlobalEfficiencyUncertainties[i] = source.fGlobalEfficiencyUncertainties[i];
238 fTab[i] = source.fTab[i];
244 //_________________________________________________________________________________________________________
245 AliHFPtSpectrum::~AliHFPtSpectrum(){
249 if (fhDirectMCpt) delete fhDirectMCpt;
250 if (fhFeedDownMCpt) delete fhFeedDownMCpt;
251 if (fhDirectMCptMax) delete fhDirectMCptMax;
252 if (fhDirectMCptMin) delete fhDirectMCptMin;
253 if (fhFeedDownMCptMax) delete fhFeedDownMCptMax;
254 if (fhFeedDownMCptMin) delete fhFeedDownMCptMin;
255 if (fhDirectEffpt) delete fhDirectEffpt;
256 if (fhFeedDownEffpt) delete fhFeedDownEffpt;
257 if (fhRECpt) delete fhRECpt;
258 if (fgRECSystematics) delete fgRECSystematics;
259 if (fhFc) delete fhFc;
260 if (fhFcMax) delete fhFcMax;
261 if (fhFcMin) delete fhFcMin;
262 if (fhFcRcb) delete fhFcRcb;
263 if (fgFcExtreme) delete fgFcExtreme;
264 if (fgFcConservative) delete fgFcConservative;
265 if (fhYieldCorr) delete fhYieldCorr;
266 if (fhYieldCorrMax) delete fhYieldCorrMax;
267 if (fhYieldCorrMin) delete fhYieldCorrMin;
268 if (fhYieldCorrRcb) delete fhYieldCorrRcb;
269 if (fgYieldCorr) delete fgYieldCorr;
270 if (fgYieldCorrExtreme) delete fgYieldCorrExtreme;
271 if (fgYieldCorrConservative) delete fgYieldCorrConservative;
272 if (fhSigmaCorr) delete fhSigmaCorr;
273 if (fhSigmaCorrMax) delete fhSigmaCorrMax;
274 if (fhSigmaCorrMin) delete fhSigmaCorrMin;
275 if (fhSigmaCorrDataSyst) delete fhSigmaCorrDataSyst;
276 if (fgSigmaCorr) delete fgSigmaCorr;
277 if (fgSigmaCorrExtreme) delete fgSigmaCorrExtreme;
278 if (fgSigmaCorrConservative) delete fgSigmaCorrConservative;
279 if (fnSigma) delete fnSigma;
280 if (fnHypothesis) delete fnHypothesis;
284 //_________________________________________________________________________________________________________
285 TH1D * AliHFPtSpectrum::RebinTheoreticalSpectra(TH1D *hTheory, const char *name) {
287 // Function to rebin the theoretical spectrum
288 // with respect to the real-data reconstructed spectrum binning
291 if (!hTheory || !fhRECpt) {
292 AliError("Feed-down or reconstructed spectra don't exist");
297 // Get the reconstructed spectra bins & limits
298 Int_t nbins = fhRECpt->GetNbinsX();
299 Int_t nbinsMC = hTheory->GetNbinsX();
300 Double_t *limits = new Double_t[nbins+1];
301 Double_t xlow=0., binwidth=0.;
302 for (Int_t i=1; i<=nbins; i++) {
303 binwidth = fhRECpt->GetBinWidth(i);
304 xlow = fhRECpt->GetBinLowEdge(i);
307 limits[nbins] = xlow + binwidth;
309 // Check that the reconstructed spectra binning
310 // is larger than the theoretical one
311 Double_t thbinwidth = hTheory->GetBinWidth(1);
312 for (Int_t i=1; i<=nbins; i++) {
313 binwidth = fhRECpt->GetBinWidth(i);
314 if ( thbinwidth > binwidth ) {
315 AliInfo(" Beware it seems that the reconstructed spectra has a smaller binning than the theoretical predictions !! ");
320 // Define a new histogram with the real-data reconstructed spectrum binning
321 TH1D * hTheoryRebin = new TH1D(name," theoretical rebinned prediction",nbins,limits);
323 Double_t sum[nbins], items[nbins];
324 for (Int_t ibin=0; ibin<nbins; ibin++) {
325 sum[ibin]=0.; items[ibin]=0.;
327 for (Int_t ibin=0; ibin<=nbinsMC; ibin++){
329 for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++){
330 if (hTheory->GetBinCenter(ibin)>limits[ibinrec] &&
331 hTheory->GetBinCenter(ibin)<limits[ibinrec+1]){
332 sum[ibinrec]+=hTheory->GetBinContent(ibin);
339 // set the theoretical rebinned spectra to ( sum-bins / n-bins ) per new bin
340 for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++) {
341 hTheoryRebin->SetBinContent(ibinrec+1,sum[ibinrec]/items[ibinrec]);
344 return (TH1D*)hTheoryRebin;
347 //_________________________________________________________________________________________________________
348 void AliHFPtSpectrum::SetMCptSpectra(TH1D *hDirect, TH1D *hFeedDown){
350 // Set the MonteCarlo or Theoretical spectra
351 // both for direct and feed-down contributions
354 if (!hDirect || !hFeedDown || !fhRECpt) {
355 AliError("One or both (direct, feed-down) spectra or the reconstructed spectra don't exist");
359 Bool_t areconsistent = kTRUE;
360 areconsistent = CheckHistosConsistency(hDirect,hFeedDown);
361 if (!areconsistent) {
362 AliInfo("Histograms are not consistent (bin width, bounds)");
367 // Rebin the theoretical predictions to the reconstructed spectra binning
369 fhDirectMCpt = RebinTheoreticalSpectra(hDirect,"fhDirectMCpt");
370 fhDirectMCpt->SetNameTitle("fhDirectMCpt"," direct theoretical prediction");
371 fhFeedDownMCpt = RebinTheoreticalSpectra(hFeedDown,"fhFeedDownMCpt");
372 fhFeedDownMCpt->SetNameTitle("fhFeedDownMCpt"," feed-down theoretical prediction");
376 //_________________________________________________________________________________________________________
377 void AliHFPtSpectrum::SetFeedDownMCptSpectra(TH1D *hFeedDown){
379 // Set the MonteCarlo or Theoretical spectra
380 // for feed-down contribution
383 if (!hFeedDown || !fhRECpt) {
384 AliError("Feed-down or reconstructed spectra don't exist");
389 // Rebin the theoretical predictions to the reconstructed spectra binning
391 fhFeedDownMCpt = RebinTheoreticalSpectra(hFeedDown,"fhFeedDownMCpt");
392 fhFeedDownMCpt->SetNameTitle("fhFeedDownMCpt"," feed-down theoretical prediction");
396 //_________________________________________________________________________________________________________
397 void AliHFPtSpectrum::SetMCptDistributionsBounds(TH1D *hDirectMax, TH1D *hDirectMin, TH1D *hFeedDownMax, TH1D *hFeedDownMin){
399 // Set the maximum and minimum MonteCarlo or Theoretical spectra
400 // both for direct and feed-down contributions
401 // used in case uncertainties are asymmetric and ca not be on the "basic histograms"
404 if (!hDirectMax || !hDirectMin || !hFeedDownMax|| !hFeedDownMin || !fhRECpt) {
405 AliError("One or all of the max/min direct/feed-down or the reconstructed spectra don't exist");
409 Bool_t areconsistent = kTRUE;
410 areconsistent &= CheckHistosConsistency(hDirectMax,hDirectMin);
411 areconsistent &= CheckHistosConsistency(hFeedDownMax,hFeedDownMin);
412 areconsistent &= CheckHistosConsistency(hDirectMax,hFeedDownMax);
413 if (!areconsistent) {
414 AliInfo("Histograms are not consistent (bin width, bounds)");
420 // Rebin the theoretical predictions to the reconstructed spectra binning
422 fhDirectMCptMax = RebinTheoreticalSpectra(hDirectMax,"fhDirectMCptMax");
423 fhDirectMCptMax->SetNameTitle("fhDirectMCptMax"," maximum direct theoretical prediction");
424 fhDirectMCptMin = RebinTheoreticalSpectra(hDirectMin,"fhDirectMCptMin");
425 fhDirectMCptMin->SetNameTitle("fhDirectMCptMin"," minimum direct theoretical prediction");
426 fhFeedDownMCptMax = RebinTheoreticalSpectra(hFeedDownMax,"fhFeedDownMCptMax");
427 fhFeedDownMCptMax->SetNameTitle("fhFeedDownMCptMax"," maximum feed-down theoretical prediction");
428 fhFeedDownMCptMin = RebinTheoreticalSpectra(hFeedDownMin,"fhFeedDownMCptMin");
429 fhFeedDownMCptMin->SetNameTitle("fhFeedDownMCptMin"," minimum feed-down theoretical prediction");
433 //_________________________________________________________________________________________________________
434 void AliHFPtSpectrum::SetFeedDownMCptDistributionsBounds(TH1D *hFeedDownMax, TH1D *hFeedDownMin){
436 // Set the maximum and minimum MonteCarlo or Theoretical spectra
437 // for feed-down contributions
438 // used in case uncertainties are asymmetric and can not be on the "basic histogram"
441 if (!hFeedDownMax || !hFeedDownMin || !fhRECpt) {
442 AliError("One or all of the max/min direct/feed-down spectra don't exist");
446 Bool_t areconsistent = kTRUE;
447 areconsistent &= CheckHistosConsistency(hFeedDownMax,hFeedDownMin);
448 if (!areconsistent) {
449 AliInfo("Histograms are not consistent (bin width, bounds)");
455 // Rebin the theoretical predictions to the reconstructed spectra binning
457 fhFeedDownMCptMax = RebinTheoreticalSpectra(hFeedDownMax,"fhFeedDownMCptMax");
458 fhFeedDownMCptMax->SetNameTitle("fhFeedDownMCptMax"," maximum feed-down theoretical prediction");
459 fhFeedDownMCptMin = RebinTheoreticalSpectra(hFeedDownMin,"fhFeedDownMCptMin");
460 fhFeedDownMCptMin->SetNameTitle("fhFeedDownMCptMin"," minimum feed-down theoretical prediction");
464 //_________________________________________________________________________________________________________
465 void AliHFPtSpectrum::SetDirectAccEffCorrection(TH1D *hDirectEff){
467 // Set the Acceptance and Efficiency corrections
468 // for the direct contribution
472 AliError("The direct acceptance and efficiency corrections doesn't exist");
476 fhDirectEffpt = (TH1D*)hDirectEff->Clone();
477 fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance x efficiency correction");
480 //_________________________________________________________________________________________________________
481 void AliHFPtSpectrum::SetAccEffCorrection(TH1D *hDirectEff, TH1D *hFeedDownEff){
483 // Set the Acceptance and Efficiency corrections
484 // both for direct and feed-down contributions
487 if (!hDirectEff || !hFeedDownEff) {
488 AliError("One or both (direct, feed-down) acceptance and efficiency corrections don't exist");
492 Bool_t areconsistent=kTRUE;
493 areconsistent = CheckHistosConsistency(hDirectEff,hFeedDownEff);
494 if (!areconsistent) {
495 AliInfo("Histograms are not consistent (bin width, bounds)");
499 fhDirectEffpt = (TH1D*)hDirectEff->Clone();
500 fhFeedDownEffpt = (TH1D*)hFeedDownEff->Clone();
501 fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance x efficiency correction");
502 fhFeedDownEffpt->SetNameTitle("fhFeedDownEffpt"," feed-down acceptance x efficiency correction");
505 //_________________________________________________________________________________________________________
506 void AliHFPtSpectrum::SetReconstructedSpectrum(TH1D *hRec) {
508 // Set the reconstructed spectrum
512 AliError("The reconstructed spectrum doesn't exist");
516 fhRECpt = (TH1D*)hRec->Clone();
517 fhRECpt->SetNameTitle("fhRECpt"," reconstructed spectrum");
520 //_________________________________________________________________________________________________________
521 void AliHFPtSpectrum::SetReconstructedSpectrumSystematics(TGraphAsymmErrors *gRec) {
523 // Set the reconstructed spectrum (uncorrected yield) systematic uncertainties
526 // Check the compatibility with the reconstructed spectrum
527 Double_t gbinwidth = gRec->GetErrorXlow(1) + gRec->GetErrorXhigh(1) ;
528 Double_t hbinwidth = fhRECpt->GetBinWidth(1);
529 Double_t gxbincenter=0., gybincenter=0.;
530 gRec->GetPoint(1,gxbincenter,gybincenter);
531 Double_t hbincenter = fhRECpt->GetBinCenter(1);
532 if ( (gbinwidth != hbinwidth) || (gxbincenter!=hbincenter) ) {
533 AliError(" The reconstructed spectrum and its systematics don't seem compatible");
537 fgRECSystematics = gRec;
540 //_________________________________________________________________________________________________________
541 void AliHFPtSpectrum::ComputeHFPtSpectrum(Double_t deltaY, Double_t branchingRatioC, Double_t branchingRatioBintoFinalDecay) {
543 // Main function to compute the corrected cross-section:
544 // variables : analysed delta_y, BR for the final correction,
545 // BR b --> D --> decay (relative to the input theoretical prediction)
547 // Sigma = ( 1. / (lumi * delta_y * BR_c * ParticleAntiPartFactor * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
549 // Uncertainties: (stat) delta_sigma = sigma * sqrt ( (delta_spectra/spectra)^2 )
550 // (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 )
551 // (feed-down syst) delta_sigma = sigma * sqrt ( (delta_spectra_fd/spectra_fd)^2 )
553 // In HIC the feed-down correction varies with an energy loss hypothesis:
554 // Raa(c-->D) / Raa(b-->D) for the "fc" method, Raa(b-->D) for the "Nb" method (see exact formulas in the functions)
558 // First: Initialization
560 Bool_t areHistosOk = Initialize();
562 AliInfo(" Histos not properly initialized. Check : inconsistent binning ? missing histos ?");
565 // Reset the statistical uncertainties on the efficiencies if needed
566 if(!fIsStatUncEff) ResetStatUncEff();
569 // Second: Correct for feed-down
571 Int_t nbins = fhRECpt->GetNbinsX();
572 if (fFeedDownOption==1) {
573 // Compute the feed-down correction via fc-method
574 CalculateFeedDownCorrectionFc();
575 // Correct the yield for feed-down correction via fc-method
576 CalculateFeedDownCorrectedSpectrumFc();
578 else if (fFeedDownOption==2) {
579 // Correct the yield for feed-down correction via Nb-method
580 CalculateFeedDownCorrectedSpectrumNb(deltaY,branchingRatioBintoFinalDecay);
582 else if (fFeedDownOption==0) {
583 // If there is no need for feed-down correction,
584 // the "corrected" yield is equal to the raw yield
585 fhYieldCorr = (TH1D*)fhRECpt->Clone();
586 fhYieldCorr->SetNameTitle("fhYieldCorr","un-corrected yield");
587 fhYieldCorrMax = (TH1D*)fhRECpt->Clone();
588 fhYieldCorrMin = (TH1D*)fhRECpt->Clone();
589 fhYieldCorrMax->SetNameTitle("fhYieldCorrMax","un-corrected yield");
590 fhYieldCorrMin->SetNameTitle("fhYieldCorrMin","un-corrected yield");
591 fgYieldCorr = new TGraphAsymmErrors(nbins+1);
592 fAsymUncertainties=kFALSE;
595 AliInfo(" Are you sure the feed-down correction option is right ?");
599 // Print out information
600 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]);
601 if (fPbPbElossHypothesis) printf("\n\n The considered Tab is %4.2e +- %2.2e \n\n",fTab[0],fTab[1]);
604 // Finally: Correct from yields to cross-section
606 Double_t binwidth = fhRECpt->GetBinWidth(1);
607 Double_t *limits = new Double_t[nbins+1];
608 Double_t *binwidths = new Double_t[nbins];
610 for (Int_t i=1; i<=nbins; i++) {
611 binwidth = fhRECpt->GetBinWidth(i);
612 xlow = fhRECpt->GetBinLowEdge(i);
614 binwidths[i-1] = binwidth;
616 limits[nbins] = xlow + binwidth;
619 // declare the output histograms
620 fhSigmaCorr = new TH1D("fhSigmaCorr","corrected sigma",nbins,limits);
621 fgSigmaCorr = new TGraphAsymmErrors(nbins+1);
623 fhSigmaCorrMax = new TH1D("fhSigmaCorrMax","max corrected sigma",nbins,limits);
624 fhSigmaCorrMin = new TH1D("fhSigmaCorrMin","min corrected sigma",nbins,limits);
625 fhSigmaCorrDataSyst = new TH1D("fhSigmaCorrDataSyst","data syst uncertainties on the corrected sigma",nbins,limits);
627 if (fPbPbElossHypothesis && fFeedDownOption==1) {
628 fhSigmaCorrRcb = new TH2D("fhSigmaCorrRcb","corrected sigma vs Rcb Eloss hypothesis; p_{T} [GeV/c] ; Rcb Eloss hypothesis ; #sigma",nbins,limits,800,0.,4.);
629 fnSigma = new TNtuple("fnSigma"," Sigma ntuple calculation","pt:Signal:Rcb:fc:Yield:Sigma:SigmaStatUnc:SigmaMax:SigmaMin");
631 if (fPbPbElossHypothesis && fFeedDownOption==2) {
632 fhSigmaCorrRcb = new TH2D("fhSigmaCorrRcb","corrected sigma vs Rb Eloss hypothesis; p_{T} [GeV/c] ; Rb Eloss hypothesis ; #sigma",nbins,limits,800,0.,4.);
633 fnSigma = new TNtuple("fnSigma"," Sigma ntuple calculation","pt:Signal:Rb:fc:Yield:Sigma:SigmaStatUnc:SigmaMax:SigmaMin");
635 // and the output TGraphAsymmErrors
636 if (fAsymUncertainties){
637 fgSigmaCorrExtreme = new TGraphAsymmErrors(nbins+1);
638 fgSigmaCorrConservative = new TGraphAsymmErrors(nbins+1);
640 fhStatUncEffcSigma = new TH1D("fhStatUncEffcSigma","direct charm stat unc on the cross section",nbins,limits);
641 fhStatUncEffbSigma = new TH1D("fhStatUncEffbSigma","secondary charm stat unc on the cross section",nbins,limits);
644 // protect against null denominator
645 if (deltaY==0. || fLuminosity[0]==0. || fTrigEfficiency[0]==0. || branchingRatioC==0.) {
646 AliError(" Hey you ! Why luminosity or trigger-efficiency or the c-BR or delta_y are set to zero ?! ");
652 Double_t value=0, errValue=0, errvalueMax=0., errvalueMin=0.;
653 Double_t errvalueExtremeMax=0., errvalueExtremeMin=0.;
654 Double_t errvalueConservativeMax=0., errvalueConservativeMin=0.;
655 Double_t errvalueStatUncEffc=0.;
656 for(Int_t ibin=1; ibin<=nbins; ibin++){
658 // Variables initialization
659 value=0.; errValue=0.; errvalueMax=0.; errvalueMin=0.;
660 errvalueExtremeMax=0.; errvalueExtremeMin=0.;
661 errvalueConservativeMax=0.; errvalueConservativeMin=0.;
662 errvalueStatUncEffc=0.;
665 // Sigma = ( 1. / (lumi * delta_y * BR_c * ParticleAntiPartFactor * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
666 value = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0. && fhRECpt->GetBinContent(ibin)>0.) ?
667 ( fhYieldCorr->GetBinContent(ibin) / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
670 // Sigma statistical uncertainty:
671 // delta_sigma = sigma * sqrt ( (delta_spectra/spectra)^2 )
672 errValue = (value!=0.) ? value * (fhYieldCorr->GetBinError(ibin)/fhYieldCorr->GetBinContent(ibin)) : 0. ;
674 // cout<< " x "<< fhRECpt->GetBinCenter(ibin) << " sigma " << value << " +- "<< errValue << " (stat)"<<endl;
677 // Sigma systematic uncertainties
679 if (fAsymUncertainties && value>0.) {
681 // (syst but feed-down) delta_sigma = sigma * sqrt ( (delta_spectra_syst/spectra)^2 +
682 // (delta_lumi/lumi)^2 + (delta_eff_trig/eff_trig)^2 + (delta_eff/eff)^2 + (global_eff)^2 )
683 errvalueMax = value * TMath::Sqrt( (fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin))*(fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin)) +
684 (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
685 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) +
686 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
687 fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] );
688 errvalueMin = value * TMath::Sqrt( (fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin))*(fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin)) +
689 (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
690 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) +
691 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
692 fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] );
694 // Uncertainties from feed-down
695 // (feed-down syst) delta_sigma = sigma * sqrt ( (delta_spectra_fd/spectra_fd)^2 )
697 errvalueExtremeMax = value * (fgYieldCorrExtreme->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin));
698 errvalueExtremeMin = value * (fgYieldCorrExtreme->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin));
701 errvalueConservativeMax = value * (fgYieldCorrConservative->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin));
702 errvalueConservativeMin = value * (fgYieldCorrConservative->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin));
705 // stat unc of the efficiencies, separately
706 errvalueStatUncEffc = value * (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) ;
710 // protect against null denominator
711 errvalueMax = (value!=0.) ?
712 value * TMath::Sqrt( (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
713 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) +
714 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
715 fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] )
717 errvalueMin = errvalueMax;
721 // Fill the histograms
723 fhSigmaCorr->SetBinContent(ibin,value);
724 fhSigmaCorr->SetBinError(ibin,errValue);
726 Double_t x = fhYieldCorr->GetBinCenter(ibin);
727 fgSigmaCorr->SetPoint(ibin,x,value); // i,x,y
728 fgSigmaCorr->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
729 fhSigmaCorrMax->SetBinContent(ibin,value+errvalueMax);
730 fhSigmaCorrMin->SetBinContent(ibin,value-errvalueMin);
733 // Fill the histos and ntuple vs the Eloss hypothesis
735 if (fPbPbElossHypothesis) {
737 // Loop over the Eloss hypothesis
739 AliError("Error reading the fc hypothesis no ntuple, please check !!");
744 Int_t entriesHypo = fnHypothesis->GetEntries();
745 Float_t pt=0., Rhy=0., fc=0., fcMin=0., fcMax=0.;
746 fnHypothesis->SetBranchAddress("pt",&pt);
747 if (fFeedDownOption==2) fnHypothesis->SetBranchAddress("Rb",&Rhy);
748 else if (fFeedDownOption==1) fnHypothesis->SetBranchAddress("Rcb",&Rhy);
749 fnHypothesis->SetBranchAddress("fc",&fc);
750 fnHypothesis->SetBranchAddress("fcMax",&fcMax);
751 fnHypothesis->SetBranchAddress("fcMin",&fcMin);
753 for (Int_t item=0; item<entriesHypo; item++) {
755 fnHypothesis->GetEntry(item);
756 if ( TMath::Abs( pt - fhDirectEffpt->GetBinCenter(ibin) ) > 0.15 ) continue;
758 Double_t yieldRcbvalue = (fhRECpt->GetBinContent(ibin) ) ? fhRECpt->GetBinContent(ibin) * fc : 0. ;
759 yieldRcbvalue /= fhRECpt->GetBinWidth(ibin) ;
760 Double_t yieldRcbvalueMax = (fhRECpt->GetBinContent(ibin) ) ? fhRECpt->GetBinContent(ibin) * fcMax : 0. ;
761 yieldRcbvalueMax /= fhRECpt->GetBinWidth(ibin) ;
762 Double_t yieldRcbvalueMin = (fhRECpt->GetBinContent(ibin) ) ? fhRECpt->GetBinContent(ibin) * fcMin : 0. ;
763 yieldRcbvalueMin /= fhRECpt->GetBinWidth(ibin) ;
766 // Sigma = ( 1. / (lumi * delta_y * BR_c * ParticleAntiPartFactor * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
767 Double_t sigmaRcbvalue = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)>0.) ?
768 ( yieldRcbvalue / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
770 Double_t sigmaRcbvalueMax = (sigmaRcbvalue!=0.) ?
771 ( yieldRcbvalueMax / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
773 Double_t sigmaRcbvalueMin = (sigmaRcbvalue!=0.) ?
774 ( yieldRcbvalueMin / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
776 // Sigma statistical uncertainty:
777 // delta_sigma = sigma * sqrt ( (delta_spectra/spectra)^2 ) = sigma * ( delta_spectra / (spectra-corr * binwidth) )
778 Double_t sigmaRcbvalueStatUnc = (sigmaRcbvalue!=0.) ?
779 sigmaRcbvalue * ( fhRECpt->GetBinError(ibin) / ( yieldRcbvalue * fhRECpt->GetBinWidth(ibin) ) ) : 0. ;
781 fhSigmaCorrRcb->Fill( fhSigmaCorr->GetBinCenter(ibin) , Rhy, sigmaRcbvalue );
783 // cout << " pt "<< fhRECpt->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<" fc-value="<< fhFcRcb->GetBinContent(ibin,rbin) <<", yield-fcRbvalue="<<yieldRcbvalue<<", sigma-fcRbvalue="<<sigmaRcbvalue<<endl;
784 fnSigma->Fill(fhRECpt->GetBinCenter(ibin), fhRECpt->GetBinContent(ibin),
785 Rhy, fc, yieldRcbvalue, sigmaRcbvalue, sigmaRcbvalueStatUnc, sigmaRcbvalueMax, sigmaRcbvalueMin );
789 // Fill the TGraphAsymmErrors
790 if (fAsymUncertainties) {
791 fgSigmaCorrExtreme->SetPoint(ibin,x,value); // i,x,y
792 fgSigmaCorrExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
793 fgSigmaCorrConservative->SetPoint(ibin,x,value); // i,x,y
794 fgSigmaCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueConservativeMin,errvalueConservativeMax); // i,xl,xh,yl,yh
796 fhStatUncEffcSigma->SetBinContent(ibin,0.);
797 if(value>0.) fhStatUncEffcSigma->SetBinError(ibin,((errvalueStatUncEffc/value)*100.));
798 fhStatUncEffbSigma->SetBinContent(ibin,0.); fhStatUncEffbSigma->SetBinError(ibin,0.);
799 // cout << " pt "<< fhRECpt->GetBinCenter(ibin) <<" bin "<< ibin<<" stat-unc-c-sigma "<< errvalueStatUncEffc/value << endl;
808 //_________________________________________________________________________________________________________
809 TH1D * AliHFPtSpectrum::EstimateEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco, const char *name) {
811 // Function that computes the acceptance and efficiency correction
812 // based on the simulated and reconstructed spectra
813 // and using the reconstructed spectra bin width
815 // eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
819 AliInfo("Hey, the reconstructed histogram was not set yet !");
823 Int_t nbins = fhRECpt->GetNbinsX();
824 Double_t *limits = new Double_t[nbins+1];
825 Double_t xlow=0.,binwidth=0.;
826 for (Int_t i=1; i<=nbins; i++) {
827 binwidth = fhRECpt->GetBinWidth(i);
828 xlow = fhRECpt->GetBinLowEdge(i);
831 limits[nbins] = xlow + binwidth;
833 TH1D * hEfficiency = new TH1D(name," acceptance #times efficiency",nbins,limits);
835 Double_t *sumSimu=new Double_t[nbins];
836 Double_t *sumReco=new Double_t[nbins];
837 for (Int_t ibin=0; ibin<nbins; ibin++){
838 sumSimu[ibin]=0.; sumReco[ibin]=0.;
840 for (Int_t ibin=0; ibin<=hSimu->GetNbinsX(); ibin++){
842 for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++){
843 if ( hSimu->GetBinCenter(ibin)>limits[ibinrec] &&
844 hSimu->GetBinCenter(ibin)<limits[ibinrec+1] ) {
845 sumSimu[ibinrec]+=hSimu->GetBinContent(ibin);
847 if ( hReco->GetBinCenter(ibin)>limits[ibinrec] &&
848 hReco->GetBinCenter(ibin)<limits[ibinrec+1] ) {
849 sumReco[ibinrec]+=hReco->GetBinContent(ibin);
856 // the efficiency is computed as reco/sim (in each bin)
857 // its uncertainty is err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
858 Double_t eff=0., erreff=0.;
859 for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++) {
860 if (sumSimu[ibinrec]!= 0. && sumReco[ibinrec]!=0.) {
861 eff = sumReco[ibinrec] / sumSimu[ibinrec] ;
862 // protection in case eff > 1.0
863 // test calculation (make the argument of the sqrt positive)
864 erreff = TMath::Sqrt( eff * TMath::Abs(1.0 - eff) ) / TMath::Sqrt( sumSimu[ibinrec] );
866 else { eff=0.0; erreff=0.; }
867 hEfficiency->SetBinContent(ibinrec+1,eff);
868 hEfficiency->SetBinError(ibinrec+1,erreff);
874 return (TH1D*)hEfficiency;
877 //_________________________________________________________________________________________________________
878 void AliHFPtSpectrum::EstimateAndSetDirectEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco) {
880 // Function that computes the Direct acceptance and efficiency correction
881 // based on the simulated and reconstructed spectra
882 // and using the reconstructed spectra bin width
884 // eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
887 if(!fhRECpt || !hSimu || !hReco){
888 AliError("Hey, the reconstructed histogram was not set yet !");
892 fhDirectEffpt = EstimateEfficiencyRecoBin(hSimu,hReco,"fhDirectEffpt");
893 fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance #times efficiency");
897 //_________________________________________________________________________________________________________
898 void AliHFPtSpectrum::EstimateAndSetFeedDownEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco) {
900 // Function that computes the Feed-Down acceptance and efficiency correction
901 // based on the simulated and reconstructed spectra
902 // and using the reconstructed spectra bin width
904 // eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
907 if(!fhRECpt || !hSimu || !hReco){
908 AliError("Hey, the reconstructed histogram was not set yet !");
912 fhFeedDownEffpt = EstimateEfficiencyRecoBin(hSimu,hReco,"fhFeedDownEffpt");
913 fhFeedDownEffpt->SetNameTitle("fhFeedDownEffpt"," feed-down acceptance #times efficiency");
917 //_________________________________________________________________________________________________________
918 Bool_t AliHFPtSpectrum::Initialize(){
920 // Initialization of the variables (histograms)
923 if (fFeedDownOption==0) {
924 AliInfo("Getting ready for the corrections without feed-down consideration");
925 } else if (fFeedDownOption==1) {
926 AliInfo("Getting ready for the fc feed-down correction calculation");
927 } else if (fFeedDownOption==2) {
928 AliInfo("Getting ready for the Nb feed-down correction calculation");
929 } else { AliError("The calculation option must be <=2"); return kFALSE; }
931 // Start checking the input histograms consistency
932 Bool_t areconsistent=kTRUE;
935 if (!fhDirectEffpt || !fhRECpt) {
936 AliError(" Reconstructed spectra and/or the Nc efficiency distributions are not defined");
939 areconsistent &= CheckHistosConsistency(fhRECpt,fhDirectEffpt);
940 if (!areconsistent) {
941 AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
944 if (fFeedDownOption==0) return kTRUE;
947 // Common checks for options 1 (fc) & 2(Nb)
948 if (!fhFeedDownMCpt || !fhFeedDownEffpt) {
949 AliError(" Theoretical Nb and/or the Nb efficiency distributions are not defined");
952 areconsistent &= CheckHistosConsistency(fhRECpt,fhFeedDownMCpt);
953 areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownEffpt);
954 if (fAsymUncertainties) {
955 if (!fhFeedDownMCptMax || !fhFeedDownMCptMin) {
956 AliError(" Max/Min theoretical Nb distributions are not defined");
959 areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownMCptMax);
961 if (!areconsistent) {
962 AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
965 if (fFeedDownOption>1) return kTRUE;
968 // Now checks for option 1 (fc correction)
970 AliError("Theoretical Nc distributions is not defined");
973 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhFeedDownMCpt);
974 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectEffpt);
975 if (fAsymUncertainties) {
976 if (!fhDirectMCptMax || !fhDirectMCptMin) {
977 AliError(" Max/Min theoretical Nc distributions are not defined");
980 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectMCptMax);
982 if (!areconsistent) {
983 AliInfo("Histograms required for fc correction are not consistent (bin width, bounds)");
990 //_________________________________________________________________________________________________________
991 Bool_t AliHFPtSpectrum::CheckHistosConsistency(TH1D *h1, TH1D *h2){
993 // Check the histograms consistency (bins, limits)
997 AliError("One or both histograms don't exist");
1001 Double_t binwidth1 = h1->GetBinWidth(1);
1002 Double_t binwidth2 = h2->GetBinWidth(1);
1003 Double_t min1 = h1->GetBinCenter(1) - (binwidth1/2.) ;
1004 // Double_t max1 = h1->GetBinCenter(nbins1) + (binwidth1/2.) ;
1005 Double_t min2 = h2->GetBinCenter(1) - (binwidth2/2.) ;
1006 // Double_t max2 = h2->GetBinCenter(nbins2) + (binwidth2/2.) ;
1008 if (binwidth1!=binwidth2) {
1009 AliInfo(" histograms with different bin width");
1013 AliInfo(" histograms with different minimum");
1016 // if (max1!=max2) {
1017 // AliInfo(" histograms with different maximum");
1024 //_________________________________________________________________________________________________________
1025 void AliHFPtSpectrum::CalculateFeedDownCorrectionFc(){
1027 // Compute fc factor and its uncertainties bin by bin
1028 // fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
1030 // uncertainties: (conservative) combine the upper/lower N_b & N_c predictions together
1031 // (extreme) combine the upper N_b predictions with the lower N_c predictions & viceversa
1032 // systematic uncertainty on the acceptance x efficiency b/c ratio are included
1034 // In addition, in HIC the feed-down correction varies with an energy loss hypothesis: Raa(c-->D) / Raa(b-->D) = Rcb
1035 // fc (Rcb) = ( 1. / ( 1 + (eff_b/eff_c)*(N_b/N_c)* (1/Rcb) ) );
1037 AliInfo("Calculating the feed-down correction factor (fc method)");
1039 // define the variables
1040 Int_t nbins = fhRECpt->GetNbinsX();
1041 Double_t binwidth = fhRECpt->GetBinWidth(1);
1042 Double_t *limits = new Double_t[nbins+1];
1043 Double_t *binwidths = new Double_t[nbins];
1045 for (Int_t i=1; i<=nbins; i++) {
1046 binwidth = fhRECpt->GetBinWidth(i);
1047 xlow = fhRECpt->GetBinLowEdge(i);
1049 binwidths[i-1] = binwidth;
1051 limits[nbins] = xlow + binwidth;
1053 Double_t correction=1.;
1054 Double_t theoryRatio=1.;
1055 Double_t effRatio=1.;
1056 Double_t correctionExtremeA=1., correctionExtremeB=1.;
1057 Double_t theoryRatioExtremeA=1., theoryRatioExtremeB=1.;
1058 Double_t correctionConservativeA=1., correctionConservativeB=1.;
1059 Double_t theoryRatioConservativeA=1., theoryRatioConservativeB=1.;
1060 Double_t correctionUnc=0.;
1061 Double_t correctionExtremeAUnc=0., correctionExtremeBUnc=0.;
1062 Double_t correctionConservativeAUnc=0., correctionConservativeBUnc=0.;
1064 // declare the output histograms
1065 fhFc = new TH1D("fhFc","fc correction factor",nbins,limits);
1066 fhFcMax = new TH1D("fhFcMax","max fc correction factor",nbins,limits);
1067 fhFcMin = new TH1D("fhFcMin","min fc correction factor",nbins,limits);
1068 if(fPbPbElossHypothesis) {
1069 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.);
1070 fnHypothesis = new TNtuple("fnHypothesis"," Feed-down correction vs hypothesis (fc)","pt:Rcb:fc:fcMax:fcMin");
1072 // two local control histograms
1073 TH1D *hTheoryRatio = new TH1D("hTheoryRatio","Theoretical B-->D over c-->D (feed-down/direct) ratio",nbins,limits);
1074 TH1D *hEffRatio = new TH1D("hEffRatio","Efficiency B-->D over c-->D (feed-down/direct) ratio",nbins,limits);
1075 // and the output TGraphAsymmErrors
1076 if (fAsymUncertainties) {
1077 fgFcExtreme = new TGraphAsymmErrors(nbins+1);
1078 fgFcExtreme->SetNameTitle("fgFcExtreme","fgFcExtreme");
1079 fgFcConservative = new TGraphAsymmErrors(nbins+1);
1080 fgFcConservative->SetNameTitle("fgFcConservative","fgFcConservative");
1083 fhStatUncEffcFD = new TH1D("fhStatUncEffcFD","direct charm stat unc on the feed-down correction",nbins,limits);
1084 fhStatUncEffbFD = new TH1D("fhStatUncEffbFD","secondary charm stat unc on the feed-down correction",nbins,limits);
1085 Double_t correctionConservativeAUncStatEffc=0., correctionConservativeBUncStatEffc=0.;
1086 Double_t correctionConservativeAUncStatEffb=0., correctionConservativeBUncStatEffb=0.;
1091 for (Int_t ibin=1; ibin<=nbins; ibin++) {
1093 // Variables initialization
1094 correction=1.; theoryRatio=1.; effRatio=1.;
1095 correctionExtremeA=1.; correctionExtremeB=1.;
1096 theoryRatioExtremeA=1.; theoryRatioExtremeB=1.;
1097 correctionConservativeA=1.; correctionConservativeB=1.;
1098 theoryRatioConservativeA=1.; theoryRatioConservativeB=1.;
1100 correctionExtremeAUnc=0.; correctionExtremeBUnc=0.;
1101 correctionConservativeAUnc=0.; correctionConservativeBUnc=0.;
1102 correctionConservativeAUncStatEffc=0.; correctionConservativeBUncStatEffc=0.;
1103 correctionConservativeAUncStatEffb=0.; correctionConservativeBUncStatEffb=0.;
1105 // theory_ratio = (N_b/N_c)
1106 theoryRatio = (fhDirectMCpt->GetBinContent(ibin)>0. && fhFeedDownMCpt->GetBinContent(ibin)>0.) ?
1107 fhFeedDownMCpt->GetBinContent(ibin) / fhDirectMCpt->GetBinContent(ibin) : 1.0 ;
1110 // Calculate the uncertainty [ considering only the theoretical uncertainties on Nb & Nc for now !!! ]
1112 // extreme A = direct-max, feed-down-min
1113 theoryRatioExtremeA = (fhDirectMCptMax->GetBinContent(ibin)>0. && fhFeedDownMCptMin->GetBinContent(ibin)>0.) ?
1114 fhFeedDownMCptMin->GetBinContent(ibin) / fhDirectMCptMax->GetBinContent(ibin) : 1.0 ;
1115 // extreme B = direct-min, feed-down-max
1116 theoryRatioExtremeB = (fhDirectMCptMin->GetBinContent(ibin)>0. && fhDirectMCptMax->GetBinContent(ibin)>0.) ?
1117 fhFeedDownMCptMax->GetBinContent(ibin) / fhDirectMCptMin->GetBinContent(ibin) : 1.0 ;
1118 // conservative A = direct-max, feed-down-max
1119 theoryRatioConservativeA = (fhDirectMCptMax->GetBinContent(ibin)>0. && fhFeedDownMCptMin->GetBinContent(ibin)>0.) ?
1120 fhFeedDownMCptMax->GetBinContent(ibin) / fhDirectMCptMax->GetBinContent(ibin) : 1.0 ;
1121 // conservative B = direct-min, feed-down-min
1122 theoryRatioConservativeB = (fhDirectMCptMin->GetBinContent(ibin)>0. && fhDirectMCptMax->GetBinContent(ibin)>0.) ?
1123 fhFeedDownMCptMin->GetBinContent(ibin) / fhDirectMCptMin->GetBinContent(ibin) : 1.0 ;
1125 // eff_ratio = (eff_b/eff_c)
1126 effRatio = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ?
1127 fhFeedDownEffpt->GetBinContent(ibin) / fhDirectEffpt->GetBinContent(ibin) : 1.0 ;
1129 // fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
1130 if( TMath::Abs(effRatio - 1.0)<0.0001 || TMath::Abs(theoryRatio - 1.0)<0.0001 ) {
1132 correctionExtremeA = 1.0;
1133 correctionExtremeB = 1.0;
1134 correctionConservativeA = 1.0;
1135 correctionConservativeB = 1.0;
1138 correction = ( 1. / ( 1 + ( effRatio * theoryRatio ) ) );
1139 correctionExtremeA = ( 1. / ( 1 + ( effRatio * theoryRatioExtremeA ) ) );
1140 correctionExtremeB = ( 1. / ( 1 + ( effRatio * theoryRatioExtremeB ) ) );
1141 correctionConservativeA = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeA ) ) );
1142 correctionConservativeB = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeB ) ) );
1146 // fc uncertainty from (eff_b/eff_c) = fc^2 * (N_b/N_c) * delta(eff_b/eff_c)
1147 // delta(eff_b/eff_c) is a percentage = effRatio * sqrt( fGlobalEfficiencyUncertainties[1]^2 + unc_eff_c ^2 + unc_eff_b ^2 )
1148 Double_t relEffUnc = TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] +
1149 (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
1150 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))
1153 correctionUnc = correction*correction * theoryRatio * effRatio * relEffUnc;
1154 correctionExtremeAUnc = correctionExtremeA*correctionExtremeA * theoryRatioExtremeA * effRatio * relEffUnc;
1155 correctionExtremeBUnc = correctionExtremeB*correctionExtremeB * theoryRatioExtremeB * effRatio * relEffUnc;
1157 correctionConservativeAUnc = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA *effRatio * relEffUnc;
1159 correctionConservativeAUncStatEffc = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA *effRatio *
1160 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin));
1161 correctionConservativeAUncStatEffb = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA *effRatio *
1162 (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin));
1164 correctionConservativeBUnc = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB *effRatio * relEffUnc;
1166 correctionConservativeBUncStatEffb = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB *effRatio *
1167 (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin));
1168 correctionConservativeBUncStatEffc = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB *effRatio *
1169 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin));
1172 // Fill in the histograms
1173 hTheoryRatio->SetBinContent(ibin,theoryRatio);
1174 hEffRatio->SetBinContent(ibin,effRatio);
1175 fhFc->SetBinContent(ibin,correction);
1177 // Estimate how the result varies vs charm/beauty Eloss hypothesis
1179 if ( correction>1.0e-16 && fPbPbElossHypothesis){
1180 // Loop over the Eloss hypothesis
1182 for (Float_t rval=0.0025; rval<4.0; rval+=0.005){
1183 // Central fc with Eloss hypothesis
1184 Double_t correctionRcb = ( 1. / ( 1 + ( effRatio * theoryRatio * (1/rval) ) ) );
1185 fhFcRcb->Fill( fhFc->GetBinCenter(ibin) , rval, correctionRcb );
1187 // cout << " pt "<< fhFc->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<", fc-Rcb-value="<<correctionRcb<<endl;
1190 // Upper / lower fc with up / low FONLL bands and Eloss hypothesis
1191 Double_t correctionConservativeARcb = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeA * (1/rval) ) ) );
1192 Double_t correctionConservativeBRcb = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeB * (1/rval) ) ) );
1193 Double_t correctionConservativeARcbUnc = correctionConservativeARcb*correctionConservativeARcb * theoryRatioConservativeA * (1/rval) *effRatio * relEffUnc;
1194 Double_t correctionConservativeBRcbUnc = correctionConservativeBRcb*correctionConservativeBRcb * theoryRatioConservativeB * (1/rval) *effRatio * relEffUnc;
1195 Double_t consvalRcb[4] = { correctionConservativeARcb - correctionConservativeARcbUnc, correctionConservativeARcb + correctionConservativeARcbUnc,
1196 correctionConservativeBRcb - correctionConservativeBRcbUnc, correctionConservativeBRcb + correctionConservativeBRcbUnc};
1197 Double_t uncConservativeRcbMin = correctionRcb - TMath::MinElement(4,consvalRcb);
1198 Double_t uncConservativeRcbMax = TMath::MaxElement(4,consvalRcb) - correctionRcb;
1200 // cout << " pt="<<fhDirectEffpt->GetBinCenter(ibin)<<", hypo="<<rval<<", fc="<<correctionRcb<<" +"<<uncConservativeRcbMax<<" -"<<uncConservativeRcbMin<<endl;
1201 fnHypothesis->Fill( fhDirectEffpt->GetBinCenter(ibin), rval, correctionRcb, correctionRcb+uncConservativeRcbMax, correctionRcb-uncConservativeRcbMin);
1205 // Fill the rest of (asymmetric) histograms
1207 if (fAsymUncertainties) {
1208 Double_t x = fhDirectMCpt->GetBinCenter(ibin);
1209 Double_t val[4] = { correctionExtremeA + correctionExtremeAUnc, correctionExtremeA - correctionExtremeAUnc,
1210 correctionExtremeB + correctionExtremeBUnc, correctionExtremeB - correctionExtremeBUnc };
1211 Double_t uncExtremeMin = correction - TMath::MinElement(4,val);
1212 Double_t uncExtremeMax = TMath::MaxElement(4,val) - correction;
1213 fgFcExtreme->SetPoint(ibin,x,correction); // i,x,y
1214 fgFcExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),uncExtremeMin,uncExtremeMax); // i,xl,xh,yl,yh
1215 fhFcMax->SetBinContent(ibin,correction+uncExtremeMax);
1216 fhFcMin->SetBinContent(ibin,correction-uncExtremeMin);
1217 Double_t consval[4] = { correctionConservativeA - correctionConservativeAUnc, correctionConservativeA + correctionConservativeAUnc,
1218 correctionConservativeB - correctionConservativeBUnc, correctionConservativeB + correctionConservativeBUnc};
1219 Double_t uncConservativeMin = correction - TMath::MinElement(4,consval);
1220 Double_t uncConservativeMax = TMath::MaxElement(4,consval) - correction;
1221 fgFcConservative->SetPoint(ibin,x,correction); // i,x,y
1222 fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),uncConservativeMin,uncConservativeMax); // i,xl,xh,yl,yh
1223 if( !(correction>0.) ){
1224 fgFcExtreme->SetPoint(ibin,x,0.); // i,x,y
1225 fgFcExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.); // i,xl,xh,yl,yh
1226 fgFcConservative->SetPoint(ibin,x,0.); // i,x,y
1227 fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.); // i,xl,xh,yl,yh
1230 Double_t valStatEffc[2] = { correctionConservativeAUncStatEffc/correctionConservativeA,
1231 correctionConservativeBUncStatEffc/correctionConservativeB };
1232 Double_t valStatEffb[2] = { correctionConservativeAUncStatEffb/correctionConservativeA,
1233 correctionConservativeBUncStatEffb/correctionConservativeB };
1234 Double_t uncConservativeStatEffc = TMath::MaxElement(2,valStatEffc);
1235 Double_t uncConservativeStatEffb = TMath::MaxElement(2,valStatEffb);
1236 fhStatUncEffcFD->SetBinContent(ibin,0.); fhStatUncEffcFD->SetBinError(ibin,uncConservativeStatEffc*100.);
1237 fhStatUncEffbFD->SetBinContent(ibin,0.); fhStatUncEffbFD->SetBinError(ibin,uncConservativeStatEffb*100.);
1238 // cout << " pt "<< fhStatUncEffcFD->GetBinCenter(ibin) <<" bin "<< ibin<<" fc-stat-c ="<<uncConservativeStatEffc<<" fc-stat-b ="<<uncConservativeStatEffb<<endl;
1242 delete [] binwidths;
1247 //_________________________________________________________________________________________________________
1248 void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumFc(){
1250 // Compute the feed-down corrected spectrum if feed-down correction is done via fc factor (bin by bin)
1251 // physics = reco * fc / bin-width
1253 // uncertainty: (stat) delta_physics = physics * sqrt ( (delta_reco/reco)^2 )
1254 // (syst but feed-down) delta_physics = physics * sqrt ( (delta_reco_syst/reco)^2 )
1255 // (feed-down syst) delta_physics = physics * sqrt ( (delta_fc/fc)^2 )
1257 // ( Calculation done bin by bin )
1259 // In addition, in HIC the feed-down correction varies with an energy loss hypothesis: Raa(c-->D) / Raa(b-->D) = Rcb
1261 AliInfo(" Calculating the feed-down corrected spectrum (fc method)");
1263 if (!fhFc || !fhRECpt) {
1264 AliError(" Reconstructed or fc distributions are not defined");
1268 Int_t nbins = fhRECpt->GetNbinsX();
1269 Double_t value = 0., errvalue = 0., errvalueMax= 0., errvalueMin= 0.;
1270 Double_t valueExtremeMax= 0., valueExtremeMin= 0.;
1271 Double_t valueConservativeMax= 0., valueConservativeMin= 0.;
1272 Double_t binwidth = fhRECpt->GetBinWidth(1);
1273 Double_t *limits = new Double_t[nbins+1];
1274 Double_t *binwidths = new Double_t[nbins];
1276 for (Int_t i=1; i<=nbins; i++) {
1277 binwidth = fhRECpt->GetBinWidth(i);
1278 xlow = fhRECpt->GetBinLowEdge(i);
1280 binwidths[i-1] = binwidth;
1282 limits[nbins] = xlow + binwidth;
1284 // declare the output histograms
1285 fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by fc)",nbins,limits);
1286 fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by fc)",nbins,limits);
1287 fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by fc)",nbins,limits);
1288 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.);
1289 // and the output TGraphAsymmErrors
1290 fgYieldCorr = new TGraphAsymmErrors(nbins+1);
1291 if (fAsymUncertainties){
1292 fgYieldCorrExtreme = new TGraphAsymmErrors(nbins+1);
1293 fgYieldCorrConservative = new TGraphAsymmErrors(nbins+1);
1297 // Do the calculation
1299 for (Int_t ibin=1; ibin<=nbins; ibin++) {
1301 // Variables initialization
1302 value = 0.; errvalue = 0.; errvalueMax= 0.; errvalueMin= 0.;
1303 valueExtremeMax= 0.; valueExtremeMin= 0.;
1304 valueConservativeMax= 0.; valueConservativeMin= 0.;
1307 // calculate the value
1308 // physics = reco * fc / bin-width
1309 value = (fhRECpt->GetBinContent(ibin) && fhFc->GetBinContent(ibin)) ?
1310 fhRECpt->GetBinContent(ibin) * fhFc->GetBinContent(ibin) : 0. ;
1311 value /= fhRECpt->GetBinWidth(ibin) ;
1313 // Statistical uncertainty
1314 // (stat) delta_physics = physics * sqrt ( (delta_reco/reco)^2 )
1315 errvalue = (value!=0. && fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) ?
1316 value * (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) : 0. ;
1318 // Calculate the systematic uncertainties
1319 // (syst but feed-down) delta_physics = physics * sqrt ( (delta_reco_syst/reco)^2 )
1320 // (feed-down syst) delta_physics = physics * sqrt ( (delta_fc/fc)^2 )
1322 // Protect against null denominator. If so, define uncertainty as null
1323 if (fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) {
1325 if (fAsymUncertainties) {
1327 // Systematics but feed-down
1328 if (fgRECSystematics) {
1329 errvalueMax = value * ( fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinContent(ibin) );
1330 errvalueMin = value * ( fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinContent(ibin) );
1332 else { errvalueMax = 0.; errvalueMin = 0.; }
1334 // Extreme feed-down systematics
1335 valueExtremeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcExtreme->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1336 valueExtremeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcExtreme->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1338 // Conservative feed-down systematics
1339 valueConservativeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcConservative->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1340 valueConservativeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcConservative->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1345 else { errvalueMax = 0.; errvalueMin = 0.; }
1348 // Fill in the histograms
1350 fhYieldCorr->SetBinContent(ibin,value);
1351 fhYieldCorr->SetBinError(ibin,errvalue);
1353 // Fill the histos and ntuple vs the Eloss hypothesis
1355 if (fPbPbElossHypothesis) {
1356 // Loop over the Eloss hypothesis
1357 for (Float_t rval=0.0025; rval<4.0; rval+=0.005){
1358 Int_t rbin = FindTH2YBin(fhYieldCorrRcb,rval);
1359 Double_t fcRcbvalue = fhFcRcb->GetBinContent(ibin,rbin);
1360 // physics = reco * fcRcb / bin-width
1361 Double_t Rcbvalue = (fhRECpt->GetBinContent(ibin) && fcRcbvalue) ?
1362 fhRECpt->GetBinContent(ibin) * fcRcbvalue : 0. ;
1363 Rcbvalue /= fhRECpt->GetBinWidth(ibin) ;
1364 fhYieldCorrRcb->Fill( fhYieldCorr->GetBinCenter(ibin) , rval, Rcbvalue );
1365 // cout << " pt "<< fhRECpt->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<" fc-fcRbvalue="<<fcRcbvalue<<", yield="<<Rcbvalue<<endl;
1368 if (fAsymUncertainties) {
1369 Double_t center = fhYieldCorr->GetBinCenter(ibin);
1370 fgYieldCorr->SetPoint(ibin,center,value); // i,x,y
1371 fgYieldCorr->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
1372 fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax);
1373 fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
1374 fgYieldCorrExtreme->SetPoint(ibin,center,value);
1375 fgYieldCorrExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),value-valueExtremeMin,valueExtremeMax-value);
1376 fgYieldCorrConservative->SetPoint(ibin,center,value);
1377 fgYieldCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),value-valueConservativeMin,valueConservativeMax-value);
1381 delete [] binwidths;
1386 //_________________________________________________________________________________________________________
1387 void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumNb(Double_t deltaY, Double_t branchingRatioBintoFinalDecay) {
1389 // Compute the feed-down corrected spectrum if feed-down correction is done via Nb (bin by bin)
1390 // physics = [ reco - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) ] / bin-width
1392 // uncertainty: (stat) delta_physics = sqrt ( (delta_reco)^2 ) / bin-width
1393 // (syst but feed-down) delta_physics = sqrt ( (delta_reco_syst)^2 ) / bin-width
1394 // (feed-down syst) delta_physics = sqrt ( (k*delta_lumi/lumi)^2 + (k*delta_eff_trig/eff_trig)^2
1395 // + (k*delta_Nb/Nb)^2 + (k*delta_eff/eff)^2 + (k*global_eff_ratio)^2 ) / bin-width
1396 // where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th
1398 // In addition, in HIC the feed-down correction varies with an energy loss hypothesis: Raa(b-->D) = Rb
1399 // physics = [ reco - ( Tab * Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th * Rb ) ] / bin-width
1401 AliInfo("Calculating the feed-down correction factor and spectrum (Nb method)");
1403 Int_t nbins = fhRECpt->GetNbinsX();
1404 Double_t binwidth = fhRECpt->GetBinWidth(1);
1405 Double_t value = 0., errvalue = 0., errvalueMax = 0., errvalueMin = 0., kfactor = 0.;
1406 Double_t errvalueExtremeMax = 0., errvalueExtremeMin = 0.;
1407 Double_t *limits = new Double_t[nbins+1];
1408 Double_t *binwidths = new Double_t[nbins];
1410 for (Int_t i=1; i<=nbins; i++) {
1411 binwidth = fhRECpt->GetBinWidth(i);
1412 xlow = fhRECpt->GetBinLowEdge(i);
1414 binwidths[i-1] = binwidth;
1416 limits[nbins] = xlow + binwidth;
1418 // declare the output histograms
1419 fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by Nb)",nbins,limits);
1420 fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by Nb)",nbins,limits);
1421 fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by Nb)",nbins,limits);
1422 if(fPbPbElossHypothesis) {
1423 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.);
1424 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.);
1425 fnHypothesis = new TNtuple("fnHypothesis"," Feed-down correction vs hypothesis (Nb)","pt:Rb:fc:fcMax:fcMin");
1427 // and the output TGraphAsymmErrors
1428 fgYieldCorr = new TGraphAsymmErrors(nbins+1);
1429 if (fAsymUncertainties){
1430 fgYieldCorrExtreme = new TGraphAsymmErrors(nbins+1);
1431 fgYieldCorrConservative = new TGraphAsymmErrors(nbins+1);
1432 // Define fc-conservative
1433 fgFcConservative = new TGraphAsymmErrors(nbins+1);
1434 AliInfo(" Beware the conservative & extreme uncertainties are equal by definition !");
1437 // variables to define fc-conservative
1438 double correction=0, correctionMax=0., correctionMin=0.;
1440 fhStatUncEffcFD = new TH1D("fhStatUncEffcFD","direct charm stat unc on the feed-down correction",nbins,limits);
1441 fhStatUncEffbFD = new TH1D("fhStatUncEffbFD","secondary charm stat unc on the feed-down correction",nbins,limits);
1442 Double_t correctionUncStatEffc=0.;
1443 Double_t correctionUncStatEffb=0.;
1447 // Do the calculation
1449 for (Int_t ibin=1; ibin<=nbins; ibin++) {
1451 // Calculate the value
1452 // physics = [ reco - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) ] / bin-width
1453 // In HIC : physics = [ reco - ( Tab * Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th * Rb ) ] / bin-width
1456 Double_t frac = 1.0, errfrac =0.;
1458 // Variables initialization
1459 value = 0.; errvalue = 0.; errvalueMax = 0.; errvalueMin = 0.; kfactor = 0.;
1460 errvalueExtremeMax = 0.; errvalueExtremeMin = 0.;
1461 correction=0; correctionMax=0.; correctionMin=0.;
1462 correctionUncStatEffc=0.; correctionUncStatEffb=0.;
1464 if(fPbPbElossHypothesis) {
1465 frac = fTab[0]*fNevts;
1466 if(fIsEventPlane) frac = frac/2.0;
1467 errfrac = frac * TMath::Sqrt( (fTab[1]/fTab[0])*(fTab[1]/fTab[0]) + (1/fNevts) );
1469 frac = fLuminosity[0];
1470 errfrac = fLuminosity[1];
1473 value = ( fhRECpt->GetBinContent(ibin)>0. && fhRECpt->GetBinContent(ibin)!=0. &&
1474 fhFeedDownMCpt->GetBinContent(ibin)>0. && fhFeedDownEffpt->GetBinContent(ibin)>0. ) ?
1475 fhRECpt->GetBinContent(ibin) - frac*(deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) * fhRECpt->GetBinWidth(ibin) )
1477 value /= fhRECpt->GetBinWidth(ibin);
1478 if (value<0.) value =0.;
1480 // Statistical uncertainty: delta_physics = sqrt ( (delta_reco)^2 ) / bin-width
1481 errvalue = (value!=0. && fhRECpt->GetBinError(ibin) && fhRECpt->GetBinError(ibin)!=0.) ?
1482 fhRECpt->GetBinError(ibin) : 0.;
1483 errvalue /= fhRECpt->GetBinWidth(ibin);
1485 // Correction (fc) : Estimate of the relative amount feed-down subtracted
1486 // correction = [ 1 - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th)/reco ]
1487 // in HIC: correction = [ 1 - ( Tab * Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th)/reco ]
1488 correction = (value>0.) ?
1489 1 - (frac*deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) * fhRECpt->GetBinWidth(ibin) ) / fhRECpt->GetBinContent(ibin) : 0. ;
1490 if (correction<0.) correction = 0.;
1492 // 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;
1494 // Systematic uncertainties
1495 // (syst but feed-down) delta_physics = sqrt ( (delta_reco_syst)^2 ) / bin-width
1496 // (feed-down syst) delta_physics = sqrt ( (k*delta_lumi/lumi)^2 + (k*delta_eff_trig/eff_trig)^2
1497 // + (k*delta_Nb/Nb)^2 + (k*delta_eff/eff)^2 + (k*global_eff_ratio)^2 ) / bin-width
1498 // where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th * bin-width
1499 kfactor = frac*deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) * fhRECpt->GetBinWidth(ibin) ;
1501 if (fAsymUncertainties && value>0.) {
1502 Double_t nb = fhFeedDownMCpt->GetBinContent(ibin);
1503 Double_t nbDmax = fhFeedDownMCptMax->GetBinContent(ibin) - fhFeedDownMCpt->GetBinContent(ibin);
1504 Double_t nbDmin = fhFeedDownMCpt->GetBinContent(ibin) - fhFeedDownMCptMin->GetBinContent(ibin);
1506 // Systematics but feed-down
1507 if (fgRECSystematics){
1508 errvalueMax = fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinWidth(ibin) ;
1509 errvalueMin = fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinWidth(ibin);
1511 else { errvalueMax = 0.; errvalueMin = 0.; }
1513 // Feed-down systematics
1514 // min value with the maximum Nb
1515 Double_t errCom = ( (kfactor*errfrac/frac)*(kfactor*errfrac/frac) ) +
1516 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1517 ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
1518 ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) ) ;
1519 errvalueExtremeMin = TMath::Sqrt( errCom + ( (kfactor*nbDmax/nb)*(kfactor*nbDmax/nb) ) ) / fhRECpt->GetBinWidth(ibin);
1520 // max value with the minimum Nb
1521 errvalueExtremeMax = TMath::Sqrt( errCom + ( (kfactor*nbDmin/nb)*(kfactor*nbDmin/nb) ) ) / fhRECpt->GetBinWidth(ibin);
1523 // Correction systematics (fc)
1524 // min value with the maximum Nb
1525 correctionMin = TMath::Sqrt( errCom + ( (kfactor*nbDmax/nb)*(kfactor*nbDmax/nb) ) ) / fhRECpt->GetBinContent(ibin) ;
1526 // max value with the minimum Nb
1527 correctionMax = TMath::Sqrt( errCom + ( (kfactor*nbDmin/nb)*(kfactor*nbDmin/nb) ) ) / fhRECpt->GetBinContent(ibin) ;
1529 correctionUncStatEffb = TMath::Sqrt( ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) )
1530 ) / fhRECpt->GetBinContent(ibin) ;
1531 correctionUncStatEffc = 0.;
1533 else{ // Don't consider Nb uncertainty in this case [ to be tested!!! ]
1534 errvalueExtremeMax = TMath::Sqrt( ( (kfactor*errfrac/frac)*(kfactor*errfrac/frac) ) +
1535 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1536 ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
1537 ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1538 ) / fhRECpt->GetBinWidth(ibin);
1539 errvalueExtremeMin = errvalueExtremeMax ;
1543 // fill in histograms
1544 fhYieldCorr->SetBinContent(ibin,value);
1545 fhYieldCorr->SetBinError(ibin,errvalue);
1547 // Estimate how the result varies vs charm/beauty Eloss hypothesis
1549 if ( correction>1.0e-16 && fPbPbElossHypothesis){
1550 // Loop over the Eloss hypothesis
1552 for (Float_t rval=0.0025; rval<4.0; rval+=0.005){
1553 // correction = [ 1 - (Tab *Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th *binwidth )* (rval) /reco ]
1554 Double_t fcRcbvalue = 1 - (frac*deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin)*fhRECpt->GetBinWidth(ibin) * rval ) / fhRECpt->GetBinContent(ibin) ;
1555 // cout << " rval="<<rval <<" , fc="<<fcRcbvalue<<" ";
1556 if(fcRcbvalue<1.e-16) fcRcbvalue=0.;
1557 fhFcRcb->Fill( fhRECpt->GetBinCenter(ibin) , rval, fcRcbvalue );
1558 // physics = reco * fcRcb / bin-width
1559 Double_t Rcbvalue = (fhRECpt->GetBinContent(ibin) && fcRcbvalue) ?
1560 fhRECpt->GetBinContent(ibin) * fcRcbvalue : 0. ;
1561 Rcbvalue /= fhRECpt->GetBinWidth(ibin) ;
1562 fhYieldCorrRcb->Fill( fhYieldCorr->GetBinCenter(ibin) , rval, Rcbvalue );
1564 // cout << " pt "<< fhFcRcb->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<", fc-Rb-value="<< fcRcbvalue << ", yield-Rb-value="<< Rcbvalue <<endl;
1565 // cout << " pt "<< fhFcRcb->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", fc-Rb-value="<< fcRcbvalue << ", yield-Rb-value="<< Rcbvalue <<endl;
1568 Double_t correctionMaxRcb = correctionMax*rval;
1569 Double_t correctionMinRcb = correctionMin*rval;
1570 fnHypothesis->Fill( fhYieldCorr->GetBinCenter(ibin), rval, fcRcbvalue, fcRcbvalue + correctionMaxRcb, fcRcbvalue - correctionMinRcb);
1572 // cout << " pt="<< fhFcRcb->GetBinCenter(ibin) <<", rval="<<rval<<", fc="<<fcRcbvalue<<" +"<<correctionMaxRcb<<" -"<<correctionMinRcb<<endl;}
1576 // Fill the rest of (asymmetric) histograms
1578 if (fAsymUncertainties) {
1579 Double_t x = fhYieldCorr->GetBinCenter(ibin);
1580 fgYieldCorr->SetPoint(ibin,x,value); // i,x,y
1581 fgYieldCorr->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
1582 fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax);
1583 fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
1584 fgYieldCorrExtreme->SetPoint(ibin,x,value); // i,x,y
1585 fgYieldCorrExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
1586 fgYieldCorrConservative->SetPoint(ibin,x,value); // i,x,y
1587 fgYieldCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
1588 // cout << " bin " << ibin << ", correction " << correction << ", min correction unc " << correctionMin << ", max correction unc " << correctionMax << endl;
1590 fgFcConservative->SetPoint(ibin,x,correction);
1591 fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),correctionMin,correctionMax);
1593 fhStatUncEffbFD->SetBinContent(ibin,0.); fhStatUncEffbFD->SetBinError(ibin,correctionUncStatEffb/correction*100.);
1594 fhStatUncEffcFD->SetBinContent(ibin,0.); fhStatUncEffcFD->SetBinError(ibin,correctionUncStatEffc/correction*100.);
1595 // cout << " pt "<< fhStatUncEffcFD->GetBinCenter(ibin) <<" bin "<< ibin<<" fc-stat-c ="<< correctionUncStatEffc/correction <<" fc-stat-b ="<< correctionUncStatEffb/correction <<endl;
1598 fgFcConservative->SetPoint(ibin,x,0.);
1599 fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.);
1604 delete [] binwidths;
1610 //_________________________________________________________________________________________________________
1611 void AliHFPtSpectrum::ComputeSystUncertainties(AliHFSystErr *systematics, Bool_t combineFeedDown) {
1613 // Function that re-calculates the global systematic uncertainties
1614 // by calling the class AliHFSystErr and combining those
1615 // (in quadrature) with the feed-down subtraction uncertainties
1618 // Estimate the feed-down uncertainty in percentage
1620 TGraphAsymmErrors *grErrFeeddown = 0;
1621 Double_t x=0., y=0., errx=0., erryl=0., erryh=0;
1622 if(fFeedDownOption!=0) {
1623 nentries = fgSigmaCorrConservative->GetN();
1624 grErrFeeddown = new TGraphAsymmErrors(nentries);
1625 for(Int_t i=0; i<nentries; i++) {
1626 x=0.; y=0.; errx=0.; erryl=0.; erryh=0.;
1627 fgSigmaCorrConservative->GetPoint(i,x,y);
1629 errx = fgSigmaCorrConservative->GetErrorXlow(i) ;
1630 erryl = fgSigmaCorrConservative->GetErrorYlow(i) / y ;
1631 erryh = fgSigmaCorrConservative->GetErrorYhigh(i) / y ;
1633 // cout << " x "<< x << " +- "<<errx<<" , y "<<y<<" + "<<erryh<<" - "<<erryl<<endl;
1634 grErrFeeddown->SetPoint(i,x,0.);
1635 grErrFeeddown->SetPointError(i,errx,errx,erryl,erryh); //i, xl, xh, yl, yh
1639 // Draw all the systematics independently
1640 systematics->DrawErrors(grErrFeeddown);
1642 // Set the sigma systematic uncertainties
1643 // possibly combine with the feed-down uncertainties
1644 Double_t errylcomb=0., erryhcomb=0;
1645 for(Int_t i=1; i<nentries; i++) {
1646 fgSigmaCorr->GetPoint(i,x,y);
1647 errx = grErrFeeddown->GetErrorXlow(i) ;
1648 erryl = grErrFeeddown->GetErrorYlow(i);
1649 erryh = grErrFeeddown->GetErrorYhigh(i);
1650 if (combineFeedDown) {
1651 errylcomb = systematics->GetTotalSystErr(x,erryl) * y ;
1652 erryhcomb = systematics->GetTotalSystErr(x,erryh) * y ;
1654 errylcomb = systematics->GetTotalSystErr(x) * y ;
1655 erryhcomb = systematics->GetTotalSystErr(x) * y ;
1657 fgSigmaCorr->SetPointError(i,errx,errx,errylcomb,erryhcomb);
1659 fhSigmaCorrDataSyst->SetBinContent(i,y);
1660 erryl = systematics->GetTotalSystErr(x) * y ;
1661 fhSigmaCorrDataSyst->SetBinError(i,erryl);
1667 //_________________________________________________________________________________________________________
1668 void AliHFPtSpectrum::DrawSpectrum(TGraphAsymmErrors *gPrediction) {
1670 // Example method to draw the corrected spectrum & the theoretical prediction
1673 TCanvas *csigma = new TCanvas("csigma","Draw the corrected cross-section & the prediction");
1674 csigma->SetFillColor(0);
1675 gPrediction->GetXaxis()->SetTitleSize(0.05);
1676 gPrediction->GetXaxis()->SetTitleOffset(0.95);
1677 gPrediction->GetYaxis()->SetTitleSize(0.05);
1678 gPrediction->GetYaxis()->SetTitleOffset(0.95);
1679 gPrediction->GetXaxis()->SetTitle("p_{T} [GeV]");
1680 gPrediction->GetYaxis()->SetTitle("BR #times #frac{d#sigma}{dp_{T}} |_{|y|<0.5} [pb/GeV]");
1681 gPrediction->SetLineColor(kGreen+2);
1682 gPrediction->SetLineWidth(3);
1683 gPrediction->SetFillColor(kGreen+1);
1684 gPrediction->Draw("3CA");
1685 fgSigmaCorr->SetLineColor(kRed);
1686 fgSigmaCorr->SetLineWidth(1);
1687 fgSigmaCorr->SetFillColor(kRed);
1688 fgSigmaCorr->SetFillStyle(0);
1689 fgSigmaCorr->Draw("2");
1690 fhSigmaCorr->SetMarkerColor(kRed);
1691 fhSigmaCorr->Draw("esame");
1693 TLegend * leg = new TLegend(0.7,0.75,0.87,0.5);
1694 leg->SetBorderSize(0);
1695 leg->SetLineColor(0);
1696 leg->SetFillColor(0);
1697 leg->SetTextFont(42);
1698 leg->AddEntry(gPrediction,"FONLL ","fl");
1699 leg->AddEntry(fhSigmaCorr,"data stat. unc.","pl");
1700 leg->AddEntry(fgSigmaCorr,"data syst. unc.","f");
1706 //_________________________________________________________________________________________________________
1707 TH1D * AliHFPtSpectrum::ReweightHisto(TH1D *hToReweight, TH1D *hReference){
1709 // Function to reweight histograms for testing purposes:
1710 // This function takes the histo hToReweight and reweights
1711 // it (its pt shape) with respect to hReference
1714 // check histograms consistency
1715 Bool_t areconsistent=kTRUE;
1716 areconsistent &= CheckHistosConsistency(hToReweight,hReference);
1717 if (!areconsistent) {
1718 AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
1722 // define a new empty histogram
1723 TH1D *hReweighted = (TH1D*)hToReweight->Clone("hReweighted");
1724 hReweighted->Reset();
1725 Double_t weight=1.0;
1726 Double_t yvalue=1.0;
1727 Double_t integralRef = hReference->Integral();
1728 Double_t integralH = hToReweight->Integral();
1730 // now reweight the spectra
1732 // the weight is the relative probability of the given pt bin in the reference histo
1733 // divided by its relative probability (to normalize it) on the histo to re-weight
1734 for (Int_t i=0; i<=hToReweight->GetNbinsX(); i++) {
1735 weight = (hReference->GetBinContent(i)/integralRef) / (hToReweight->GetBinContent(i)/integralH) ;
1736 yvalue = hToReweight->GetBinContent(i);
1737 hReweighted->SetBinContent(i,yvalue*weight);
1740 return (TH1D*)hReweighted;
1743 //_________________________________________________________________________________________________________
1744 TH1D * AliHFPtSpectrum::ReweightRecHisto(TH1D *hRecToReweight, TH1D *hMCToReweight, TH1D *hMCReference){
1746 // Function to reweight histograms for testing purposes:
1747 // This function takes the histo hToReweight and reweights
1748 // it (its pt shape) with respect to hReference /hMCToReweight
1751 // check histograms consistency
1752 Bool_t areconsistent=kTRUE;
1753 areconsistent &= CheckHistosConsistency(hMCToReweight,hMCReference);
1754 areconsistent &= CheckHistosConsistency(hRecToReweight,hMCReference);
1755 if (!areconsistent) {
1756 AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
1760 // define a new empty histogram
1761 TH1D *hReweighted = (TH1D*)hMCToReweight->Clone("hReweighted");
1762 hReweighted->Reset();
1763 TH1D *hRecReweighted = (TH1D*)hRecToReweight->Clone("hRecReweighted");
1764 hRecReweighted->Reset();
1765 Double_t weight=1.0;
1766 Double_t yvalue=1.0, yrecvalue=1.0;
1767 Double_t integralRef = hMCReference->Integral();
1768 Double_t integralH = hMCToReweight->Integral();
1770 // now reweight the spectra
1772 // the weight is the relative probability of the given pt bin
1773 // that should be applied in the MC histo to get the reference histo shape
1774 // Probabilities are properly normalized.
1775 for (Int_t i=0; i<=hMCToReweight->GetNbinsX(); i++) {
1776 weight = (hMCReference->GetBinContent(i)/integralRef) / (hMCToReweight->GetBinContent(i)/integralH) ;
1777 yvalue = hMCToReweight->GetBinContent(i);
1778 hReweighted->SetBinContent(i,yvalue*weight);
1779 yrecvalue = hRecToReweight->GetBinContent(i);
1780 hRecReweighted->SetBinContent(i,yrecvalue*weight);
1783 return (TH1D*)hRecReweighted;
1788 //_________________________________________________________________________________________________________
1789 Int_t AliHFPtSpectrum::FindTH2YBin(TH2D *histo, Float_t yvalue){
1791 // Function to find the y-axis bin of a TH2 for a given y-value
1794 Int_t nbins = histo->GetNbinsY();
1796 for (int j=0; j<=nbins; j++) {
1797 Float_t value = histo->GetYaxis()->GetBinCenter(j);
1798 Float_t width = histo->GetYaxis()->GetBinWidth(j);
1799 // if( TMath::Abs(yvalue-value)<= width/2. ) {
1800 if( TMath::Abs(yvalue-value)<= width ) {
1802 // cout <<" value "<<value << ", yval "<< yvalue<<", bin width "<<width/2.<< " y ="<<ybin<<endl;
1810 //_________________________________________________________________________________________________________
1811 void AliHFPtSpectrum::ResetStatUncEff(){
1813 Int_t entries = fhDirectEffpt->GetNbinsX();
1814 for(Int_t i=0; i<=entries; i++){
1815 fhDirectEffpt->SetBinError(i,0.);
1816 fhFeedDownEffpt->SetBinError(i,0.);