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)
30 // Author: Z.Conesa, zconesa@in2p3.fr
31 //***********************************************************************
33 #include <Riostream.h>
38 #include "TGraphAsymmErrors.h"
44 #include "AliHFSystErr.h"
45 #include "AliHFPtSpectrum.h"
47 ClassImp(AliHFPtSpectrum)
49 //_________________________________________________________________________________________________________
50 AliHFPtSpectrum::AliHFPtSpectrum(const char* name, const char* title, Int_t option):
54 fhDirectMCptMax(NULL),
55 fhDirectMCptMin(NULL),
56 fhFeedDownMCptMax(NULL),
57 fhFeedDownMCptMin(NULL),
59 fhFeedDownEffpt(NULL),
61 fgRECSystematics(NULL),
64 fGlobalEfficiencyUncertainties(),
69 fgFcConservative(NULL),
74 fgYieldCorrExtreme(NULL),
75 fgYieldCorrConservative(NULL),
80 fgSigmaCorrExtreme(NULL),
81 fgSigmaCorrConservative(NULL),
82 fFeedDownOption(option),
83 fAsymUncertainties(kTRUE)
86 // Default constructor
89 fLuminosity[0]=1.; fLuminosity[1]=0.;
90 fTrigEfficiency[0]=1.; fTrigEfficiency[1]=0.;
91 fGlobalEfficiencyUncertainties[0]=0.; fGlobalEfficiencyUncertainties[1]=0.;
95 //_________________________________________________________________________________________________________
96 AliHFPtSpectrum::AliHFPtSpectrum(const AliHFPtSpectrum &rhs):
98 fhDirectMCpt(rhs.fhDirectMCpt),
99 fhFeedDownMCpt(rhs.fhFeedDownMCpt),
100 fhDirectMCptMax(rhs.fhDirectMCptMax),
101 fhDirectMCptMin(rhs.fhDirectMCptMin),
102 fhFeedDownMCptMax(rhs.fhFeedDownMCptMax),
103 fhFeedDownMCptMin(rhs.fhFeedDownMCptMin),
104 fhDirectEffpt(rhs.fhDirectEffpt),
105 fhFeedDownEffpt(rhs.fhFeedDownEffpt),
106 fhRECpt(rhs.fhRECpt),
107 fgRECSystematics(rhs.fgRECSystematics),
110 fGlobalEfficiencyUncertainties(),
112 fhFcMax(rhs.fhFcMax),
113 fhFcMin(rhs.fhFcMin),
114 fgFcExtreme(rhs.fgFcExtreme),
115 fgFcConservative(rhs.fgFcConservative),
116 fhYieldCorr(rhs.fhYieldCorr),
117 fhYieldCorrMax(rhs.fhYieldCorrMax),
118 fhYieldCorrMin(rhs.fhYieldCorrMin),
119 fgYieldCorr(rhs.fgYieldCorr),
120 fgYieldCorrExtreme(rhs.fgYieldCorrExtreme),
121 fgYieldCorrConservative(rhs.fgYieldCorrConservative),
122 fhSigmaCorr(rhs.fhSigmaCorr),
123 fhSigmaCorrMax(rhs.fhSigmaCorrMax),
124 fhSigmaCorrMin(rhs.fhSigmaCorrMin),
125 fgSigmaCorr(rhs.fgSigmaCorr),
126 fgSigmaCorrExtreme(rhs.fgSigmaCorrExtreme),
127 fgSigmaCorrConservative(rhs.fgSigmaCorrConservative),
128 fFeedDownOption(rhs.fFeedDownOption),
129 fAsymUncertainties(rhs.fAsymUncertainties)
135 for(Int_t i=0; i<2; i++){
136 fLuminosity[i] = rhs.fLuminosity[i];
137 fTrigEfficiency[i] = rhs.fTrigEfficiency[i];
138 fGlobalEfficiencyUncertainties[i] = rhs.fGlobalEfficiencyUncertainties[i];
143 //_________________________________________________________________________________________________________
144 AliHFPtSpectrum &AliHFPtSpectrum::operator=(const AliHFPtSpectrum &source){
146 // Assignment operator
149 if (&source == this) return *this;
151 fhDirectMCpt = source.fhDirectMCpt;
152 fhFeedDownMCpt = source.fhFeedDownMCpt;
153 fhDirectMCptMax = source.fhDirectMCptMax;
154 fhDirectMCptMin = source.fhDirectMCptMin;
155 fhFeedDownMCptMax = source.fhFeedDownMCptMax;
156 fhFeedDownMCptMin = source.fhFeedDownMCptMin;
157 fhDirectEffpt = source.fhDirectEffpt;
158 fhFeedDownEffpt = source.fhFeedDownEffpt;
159 fhRECpt = source.fhRECpt;
160 fgRECSystematics = source.fgRECSystematics;
162 fhFcMax = source.fhFcMax;
163 fhFcMin = source.fhFcMin;
164 fgFcExtreme = source.fgFcExtreme;
165 fgFcConservative = source.fgFcConservative;
166 fhYieldCorr = source.fhYieldCorr;
167 fhYieldCorrMax = source.fhYieldCorrMax;
168 fhYieldCorrMin = source.fhYieldCorrMin;
169 fgYieldCorr = source.fgYieldCorr;
170 fgYieldCorrExtreme = source.fgYieldCorrExtreme;
171 fgYieldCorrConservative = source.fgYieldCorrConservative;
172 fhSigmaCorr = source.fhSigmaCorr;
173 fhSigmaCorrMax = source.fhSigmaCorrMax;
174 fhSigmaCorrMin = source.fhSigmaCorrMin;
175 fgSigmaCorr = source.fgSigmaCorr;
176 fgSigmaCorrExtreme = source.fgSigmaCorrExtreme;
177 fgSigmaCorrConservative = source.fgSigmaCorrConservative;
178 fFeedDownOption = source.fFeedDownOption;
179 fAsymUncertainties = source.fAsymUncertainties;
181 for(Int_t i=0; i<2; i++){
182 fLuminosity[i] = source.fLuminosity[i];
183 fTrigEfficiency[i] = source.fTrigEfficiency[i];
184 fGlobalEfficiencyUncertainties[i] = source.fGlobalEfficiencyUncertainties[i];
190 //_________________________________________________________________________________________________________
191 AliHFPtSpectrum::~AliHFPtSpectrum(){
195 if (fhDirectMCpt) delete fhDirectMCpt;
196 if (fhFeedDownMCpt) delete fhFeedDownMCpt;
197 if (fhDirectMCptMax) delete fhDirectMCptMax;
198 if (fhDirectMCptMin) delete fhDirectMCptMin;
199 if (fhFeedDownMCptMax) delete fhFeedDownMCptMax;
200 if (fhFeedDownMCptMin) delete fhFeedDownMCptMin;
201 if (fhDirectEffpt) delete fhDirectEffpt;
202 if (fhFeedDownEffpt) delete fhFeedDownEffpt;
203 if (fhRECpt) delete fhRECpt;
204 if (fgRECSystematics) delete fgRECSystematics;
205 if (fhFc) delete fhFc;
206 if (fhFcMax) delete fhFcMax;
207 if (fhFcMin) delete fhFcMin;
208 if (fgFcExtreme) delete fgFcExtreme;
209 if (fgFcConservative) delete fgFcConservative;
210 if (fhYieldCorr) delete fhYieldCorr;
211 if (fhYieldCorrMax) delete fhYieldCorrMax;
212 if (fhYieldCorrMin) delete fhYieldCorrMin;
213 if (fgYieldCorr) delete fgYieldCorr;
214 if (fgYieldCorrExtreme) delete fgYieldCorrExtreme;
215 if (fgYieldCorrConservative) delete fgYieldCorrConservative;
216 if (fhSigmaCorr) delete fhSigmaCorr;
217 if (fhSigmaCorrMax) delete fhSigmaCorrMax;
218 if (fhSigmaCorrMin) delete fhSigmaCorrMin;
219 if (fgSigmaCorr) delete fgSigmaCorr;
220 if (fgSigmaCorrExtreme) delete fgSigmaCorrExtreme;
221 if (fgSigmaCorrConservative) delete fgSigmaCorrConservative;
225 //_________________________________________________________________________________________________________
226 TH1D * AliHFPtSpectrum::RebinTheoreticalSpectra(TH1D *hTheory, const char *name) {
228 // Function to rebin the theoretical spectrum
229 // with respect to the real-data reconstructed spectrum binning
232 if (!hTheory || !fhRECpt) {
233 AliError("Feed-down or reconstructed spectra don't exist");
238 // Get the reconstructed spectra bins & limits
239 Int_t nbins = fhRECpt->GetNbinsX();
240 Int_t nbinsMC = hTheory->GetNbinsX();
241 Double_t *limits = new Double_t[nbins+1];
242 Double_t xlow=0., binwidth=0.;
243 for (Int_t i=1; i<=nbins; i++) {
244 binwidth = fhRECpt->GetBinWidth(i);
245 xlow = fhRECpt->GetBinLowEdge(i);
248 limits[nbins] = xlow + binwidth;
250 // Check that the reconstructed spectra binning
251 // is larger than the theoretical one
252 Double_t thbinwidth = hTheory->GetBinWidth(1);
253 for (Int_t i=1; i<=nbins; i++) {
254 binwidth = fhRECpt->GetBinWidth(i);
255 if ( thbinwidth > binwidth ) {
256 AliInfo(" Beware it seems that the reconstructed spectra has a smaller binning than the theoretical predictions !! ");
261 // Define a new histogram with the real-data reconstructed spectrum binning
262 TH1D * hTheoryRebin = new TH1D(name," theoretical rebinned prediction",nbins,limits);
264 Double_t sum[nbins], items[nbins];
265 for (Int_t ibin=0; ibin<nbins; ibin++) {
266 sum[ibin]=0.; items[ibin]=0.;
268 for (Int_t ibin=0; ibin<=nbinsMC; ibin++){
270 for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++){
271 if (hTheory->GetBinCenter(ibin)>limits[ibinrec] &&
272 hTheory->GetBinCenter(ibin)<limits[ibinrec+1]){
273 sum[ibinrec]+=hTheory->GetBinContent(ibin);
280 // set the theoretical rebinned spectra to ( sum-bins / n-bins ) per new bin
281 for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++) {
282 hTheoryRebin->SetBinContent(ibinrec+1,sum[ibinrec]/items[ibinrec]);
285 return (TH1D*)hTheoryRebin;
288 //_________________________________________________________________________________________________________
289 void AliHFPtSpectrum::SetMCptSpectra(TH1D *hDirect, TH1D *hFeedDown){
291 // Set the MonteCarlo or Theoretical spectra
292 // both for direct and feed-down contributions
295 if (!hDirect || !hFeedDown || !fhRECpt) {
296 AliError("One or both (direct, feed-down) spectra or the reconstructed spectra don't exist");
300 Bool_t areconsistent = kTRUE;
301 areconsistent = CheckHistosConsistency(hDirect,hFeedDown);
302 if (!areconsistent) {
303 AliInfo("Histograms are not consistent (bin width, bounds)");
308 // Rebin the theoretical predictions to the reconstructed spectra binning
310 fhDirectMCpt = RebinTheoreticalSpectra(hDirect,"fhDirectMCpt");
311 fhDirectMCpt->SetNameTitle("fhDirectMCpt"," direct theoretical prediction");
312 fhFeedDownMCpt = RebinTheoreticalSpectra(hFeedDown,"fhFeedDownMCpt");
313 fhFeedDownMCpt->SetNameTitle("fhFeedDownMCpt"," feed-down theoretical prediction");
317 //_________________________________________________________________________________________________________
318 void AliHFPtSpectrum::SetFeedDownMCptSpectra(TH1D *hFeedDown){
320 // Set the MonteCarlo or Theoretical spectra
321 // for feed-down contribution
324 if (!hFeedDown || !fhRECpt) {
325 AliError("Feed-down or reconstructed spectra don't exist");
330 // Rebin the theoretical predictions to the reconstructed spectra binning
332 fhFeedDownMCpt = RebinTheoreticalSpectra(hFeedDown,"fhFeedDownMCpt");
333 fhFeedDownMCpt->SetNameTitle("fhFeedDownMCpt"," feed-down theoretical prediction");
337 //_________________________________________________________________________________________________________
338 void AliHFPtSpectrum::SetMCptDistributionsBounds(TH1D *hDirectMax, TH1D *hDirectMin, TH1D *hFeedDownMax, TH1D *hFeedDownMin){
340 // Set the maximum and minimum MonteCarlo or Theoretical spectra
341 // both for direct and feed-down contributions
342 // used in case uncertainties are asymmetric and ca not be on the "basic histograms"
345 if (!hDirectMax || !hDirectMin || !hFeedDownMax|| !hFeedDownMin || !fhRECpt) {
346 AliError("One or all of the max/min direct/feed-down or the reconstructed spectra don't exist");
350 Bool_t areconsistent = kTRUE;
351 areconsistent &= CheckHistosConsistency(hDirectMax,hDirectMin);
352 areconsistent &= CheckHistosConsistency(hFeedDownMax,hFeedDownMin);
353 areconsistent &= CheckHistosConsistency(hDirectMax,hFeedDownMax);
354 if (!areconsistent) {
355 AliInfo("Histograms are not consistent (bin width, bounds)");
361 // Rebin the theoretical predictions to the reconstructed spectra binning
363 fhDirectMCptMax = RebinTheoreticalSpectra(hDirectMax,"fhDirectMCptMax");
364 fhDirectMCptMax->SetNameTitle("fhDirectMCptMax"," maximum direct theoretical prediction");
365 fhDirectMCptMin = RebinTheoreticalSpectra(hDirectMin,"fhDirectMCptMin");
366 fhDirectMCptMin->SetNameTitle("fhDirectMCptMin"," minimum direct theoretical prediction");
367 fhFeedDownMCptMax = RebinTheoreticalSpectra(hFeedDownMax,"fhFeedDownMCptMax");
368 fhFeedDownMCptMax->SetNameTitle("fhFeedDownMCptMax"," maximum feed-down theoretical prediction");
369 fhFeedDownMCptMin = RebinTheoreticalSpectra(hFeedDownMin,"fhFeedDownMCptMin");
370 fhFeedDownMCptMin->SetNameTitle("fhFeedDownMCptMin"," minimum feed-down theoretical prediction");
374 //_________________________________________________________________________________________________________
375 void AliHFPtSpectrum::SetFeedDownMCptDistributionsBounds(TH1D *hFeedDownMax, TH1D *hFeedDownMin){
377 // Set the maximum and minimum MonteCarlo or Theoretical spectra
378 // for feed-down contributions
379 // used in case uncertainties are asymmetric and can not be on the "basic histogram"
382 if (!hFeedDownMax || !hFeedDownMin || !fhRECpt) {
383 AliError("One or all of the max/min direct/feed-down spectra don't exist");
387 Bool_t areconsistent = kTRUE;
388 areconsistent &= CheckHistosConsistency(hFeedDownMax,hFeedDownMin);
389 if (!areconsistent) {
390 AliInfo("Histograms are not consistent (bin width, bounds)");
396 // Rebin the theoretical predictions to the reconstructed spectra binning
398 fhFeedDownMCptMax = RebinTheoreticalSpectra(hFeedDownMax,"fhFeedDownMCptMax");
399 fhFeedDownMCptMax->SetNameTitle("fhFeedDownMCptMax"," maximum feed-down theoretical prediction");
400 fhFeedDownMCptMin = RebinTheoreticalSpectra(hFeedDownMin,"fhFeedDownMCptMin");
401 fhFeedDownMCptMin->SetNameTitle("fhFeedDownMCptMin"," minimum feed-down theoretical prediction");
405 //_________________________________________________________________________________________________________
406 void AliHFPtSpectrum::SetDirectAccEffCorrection(TH1D *hDirectEff){
408 // Set the Acceptance and Efficiency corrections
409 // for the direct contribution
413 AliError("The direct acceptance and efficiency corrections doesn't exist");
417 fhDirectEffpt = (TH1D*)hDirectEff->Clone();
418 fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance x efficiency correction");
421 //_________________________________________________________________________________________________________
422 void AliHFPtSpectrum::SetAccEffCorrection(TH1D *hDirectEff, TH1D *hFeedDownEff){
424 // Set the Acceptance and Efficiency corrections
425 // both for direct and feed-down contributions
428 if (!hDirectEff || !hFeedDownEff) {
429 AliError("One or both (direct, feed-down) acceptance and efficiency corrections don't exist");
433 Bool_t areconsistent=kTRUE;
434 areconsistent = CheckHistosConsistency(hDirectEff,hFeedDownEff);
435 if (!areconsistent) {
436 AliInfo("Histograms are not consistent (bin width, bounds)");
440 fhDirectEffpt = (TH1D*)hDirectEff->Clone();
441 fhFeedDownEffpt = (TH1D*)hFeedDownEff->Clone();
442 fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance x efficiency correction");
443 fhFeedDownEffpt->SetNameTitle("fhFeedDownEffpt"," feed-down acceptance x efficiency correction");
446 //_________________________________________________________________________________________________________
447 void AliHFPtSpectrum::SetReconstructedSpectrum(TH1D *hRec) {
449 // Set the reconstructed spectrum
453 AliError("The reconstructed spectrum doesn't exist");
457 fhRECpt = (TH1D*)hRec->Clone();
458 fhRECpt->SetNameTitle("fhRECpt"," reconstructed spectrum");
461 //_________________________________________________________________________________________________________
462 void AliHFPtSpectrum::SetReconstructedSpectrumSystematics(TGraphAsymmErrors *gRec) {
464 // Set the reconstructed spectrum (uncorrected yield) systematic uncertainties
467 // Check the compatibility with the reconstructed spectrum
468 Double_t gbinwidth = gRec->GetErrorXlow(1) + gRec->GetErrorXhigh(1) ;
469 Double_t hbinwidth = fhRECpt->GetBinWidth(1);
470 Double_t gxbincenter=0., gybincenter=0.;
471 gRec->GetPoint(1,gxbincenter,gybincenter);
472 Double_t hbincenter = fhRECpt->GetBinCenter(1);
473 if ( (gbinwidth != hbinwidth) || (gxbincenter!=hbincenter) ) {
474 AliError(" The reconstructed spectrum and its systematics don't seem compatible");
478 fgRECSystematics = gRec;
481 //_________________________________________________________________________________________________________
482 void AliHFPtSpectrum::ComputeHFPtSpectrum(Double_t deltaY, Double_t branchingRatioC, Double_t branchingRatioBintoFinalDecay) {
484 // Main function to compute the corrected cross-section:
485 // variables : analysed delta_y, BR for the final correction,
486 // BR b --> D --> decay (relative to the input theoretical prediction)
488 // Sigma = ( 1. / (lumi * delta_y * BR_c * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
490 // Uncertainties: (stat) delta_sigma = sigma * sqrt ( (delta_spectra/spectra)^2 )
491 // (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 )
492 // (feed-down syst) delta_sigma = sigma * sqrt ( (delta_spectra_fd/spectra_fd)^2 )
495 // First: Initialization
497 Bool_t areHistosOk = Initialize();
499 AliInfo(" Histos not properly initialized. Check : inconsistent binning ? missing histos ?");
504 // Second: Correct for feed-down
506 if (fFeedDownOption==1) {
507 // Compute the feed-down correction via fc-method
508 CalculateFeedDownCorrectionFc();
509 // Correct the yield for feed-down correction via fc-method
510 CalculateFeedDownCorrectedSpectrumFc();
512 else if (fFeedDownOption==2) {
513 // Correct the yield for feed-down correction via Nb-method
514 CalculateFeedDownCorrectedSpectrumNb(deltaY,branchingRatioBintoFinalDecay);
516 else if (fFeedDownOption==0) {
517 // If there is no need for feed-down correction,
518 // the "corrected" yield is equal to the raw yield
519 fhYieldCorr = (TH1D*)fhRECpt->Clone();
520 fhYieldCorr->SetNameTitle("fhYieldCorr","un-corrected yield");
521 fhYieldCorrMax = (TH1D*)fhRECpt->Clone();
522 fhYieldCorrMin = (TH1D*)fhRECpt->Clone();
523 fhYieldCorrMax->SetNameTitle("fhYieldCorrMax","un-corrected yield");
524 fhYieldCorrMin->SetNameTitle("fhYieldCorrMin","un-corrected yield");
525 fAsymUncertainties=kFALSE;
528 AliInfo(" Are you sure the feed-down correction option is right ?");
531 // Print out information
532 // 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\n",fLuminosity[0],fLuminosity[1],fTrigEfficiency[0],fTrigEfficiency[1],deltaY,branchingRatioC,branchingRatioBintoFinalDecay);
533 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]);
536 // Finally: Correct from yields to cross-section
538 Int_t nbins = fhRECpt->GetNbinsX();
539 Double_t binwidth = fhRECpt->GetBinWidth(1);
540 Double_t *limits = new Double_t[nbins+1];
541 Double_t *binwidths = new Double_t[nbins];
543 for (Int_t i=1; i<=nbins; i++) {
544 binwidth = fhRECpt->GetBinWidth(i);
545 xlow = fhRECpt->GetBinLowEdge(i);
547 binwidths[i-1] = binwidth;
549 limits[nbins] = xlow + binwidth;
552 // declare the output histograms
553 fhSigmaCorr = new TH1D("fhSigmaCorr","corrected sigma",nbins,limits);
554 fhSigmaCorrMax = new TH1D("fhSigmaCorrMax","max corrected sigma",nbins,limits);
555 fhSigmaCorrMin = new TH1D("fhSigmaCorrMin","min corrected sigma",nbins,limits);
556 // and the output TGraphAsymmErrors
557 if (fAsymUncertainties){
558 fgSigmaCorr = new TGraphAsymmErrors(nbins+1);
559 fgSigmaCorrExtreme = new TGraphAsymmErrors(nbins+1);
560 fgSigmaCorrConservative = new TGraphAsymmErrors(nbins+1);
563 // protect against null denominator
564 if (deltaY==0. || fLuminosity[0]==0. || fTrigEfficiency[0]==0. || branchingRatioC==0.) {
565 AliError(" Hey you ! Why luminosity or trigger-efficiency or the c-BR or delta_y are set to zero ?! ");
571 Double_t value=0, errValue=0, errvalueMax=0., errvalueMin=0.;
572 Double_t errvalueExtremeMax=0., errvalueExtremeMin=0.;
573 Double_t errvalueConservativeMax=0., errvalueConservativeMin=0.;
574 for(Int_t ibin=1; ibin<=nbins; ibin++){
577 // Sigma = ( 1. / (lumi * delta_y * BR_c * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
578 value = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ?
579 ( fhYieldCorr->GetBinContent(ibin) / ( deltaY * branchingRatioC * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
582 // Sigma statistical uncertainty:
583 // delta_sigma = sigma * sqrt ( (delta_spectra/spectra)^2 )
584 errValue = (value!=0.) ? value * (fhYieldCorr->GetBinError(ibin)/fhYieldCorr->GetBinContent(ibin)) : 0. ;
588 // Sigma systematic uncertainties
590 if (fAsymUncertainties) {
592 // (syst but feed-down) delta_sigma = sigma * sqrt ( (delta_spectra_syst/spectra)^2 +
593 // (delta_lumi/lumi)^2 + (delta_eff_trig/eff_trig)^2 + (delta_eff/eff)^2 + (global_eff)^2 )
594 errvalueMax = value * TMath::Sqrt( (fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin))*(fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin)) +
595 (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
596 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) +
597 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
598 fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] );
599 errvalueMin = value * TMath::Sqrt( (fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin))*(fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin)) +
600 (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
601 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) +
602 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
603 fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] );
605 // Uncertainties from feed-down
606 // (feed-down syst) delta_sigma = sigma * sqrt ( (delta_spectra_fd/spectra_fd)^2 )
608 errvalueExtremeMax = value * (fgYieldCorrExtreme->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin));
609 errvalueExtremeMin = value * (fgYieldCorrExtreme->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin));
612 errvalueConservativeMax = value * (fgYieldCorrConservative->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin));
613 errvalueConservativeMin = value * (fgYieldCorrConservative->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin));
617 // protect against null denominator
618 errvalueMax = (value!=0.) ?
619 value * TMath::Sqrt( (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
620 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) +
621 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
622 fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] )
624 errvalueMin = errvalueMax;
627 // Fill the histograms
628 fhSigmaCorr->SetBinContent(ibin,value);
629 fhSigmaCorr->SetBinError(ibin,errValue);
630 // Fill the TGraphAsymmErrors
631 if (fAsymUncertainties) {
632 Double_t x = fhYieldCorr->GetBinCenter(ibin);
633 fgSigmaCorr->SetPoint(ibin,x,value); // i,x,y
634 fgSigmaCorr->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
635 fhSigmaCorrMax->SetBinContent(ibin,value+errvalueMax);
636 fhSigmaCorrMin->SetBinContent(ibin,value-errvalueMin);
637 fgSigmaCorrExtreme->SetPoint(ibin,x,value); // i,x,y
638 fgSigmaCorrExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
639 fgSigmaCorrConservative->SetPoint(ibin,x,value); // i,x,y
640 fgSigmaCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueConservativeMin,errvalueConservativeMax); // i,xl,xh,yl,yh
650 //_________________________________________________________________________________________________________
651 TH1D * AliHFPtSpectrum::EstimateEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco, const char *name) {
653 // Function that computes the acceptance and efficiency correction
654 // based on the simulated and reconstructed spectra
655 // and using the reconstructed spectra bin width
657 // eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
661 AliInfo("Hey, the reconstructed histogram was not set yet !");
665 Int_t nbins = fhRECpt->GetNbinsX();
666 Double_t *limits = new Double_t[nbins+1];
667 Double_t xlow=0.,binwidth=0.;
668 for (Int_t i=1; i<=nbins; i++) {
669 binwidth = fhRECpt->GetBinWidth(i);
670 xlow = fhRECpt->GetBinLowEdge(i);
673 limits[nbins] = xlow + binwidth;
675 TH1D * hEfficiency = new TH1D(name," acceptance #times efficiency",nbins,limits);
677 Double_t *sumSimu=new Double_t[nbins];
678 Double_t *sumReco=new Double_t[nbins];
679 for (Int_t ibin=0; ibin<nbins; ibin++){
680 sumSimu[ibin]=0.; sumReco[ibin]=0.;
682 for (Int_t ibin=0; ibin<=hSimu->GetNbinsX(); ibin++){
684 for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++){
685 if ( hSimu->GetBinCenter(ibin)>limits[ibinrec] &&
686 hSimu->GetBinCenter(ibin)<limits[ibinrec+1] ) {
687 sumSimu[ibinrec]+=hSimu->GetBinContent(ibin);
689 if ( hReco->GetBinCenter(ibin)>limits[ibinrec] &&
690 hReco->GetBinCenter(ibin)<limits[ibinrec+1] ) {
691 sumReco[ibinrec]+=hReco->GetBinContent(ibin);
698 // the efficiency is computed as reco/sim (in each bin)
699 // its uncertainty is err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
700 Double_t eff=0., erreff=0.;
701 for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++) {
702 if (sumSimu[ibinrec]!= 0. && sumReco[ibinrec]!=0.) {
703 eff = sumReco[ibinrec] / sumSimu[ibinrec] ;
704 // protection in case eff > 1.0
705 // test calculation (make the argument of the sqrt positive)
706 erreff = TMath::Sqrt( eff * TMath::Abs(1.0 - eff) ) / TMath::Sqrt( sumSimu[ibinrec] );
708 else { eff=0.0; erreff=0.; }
709 hEfficiency->SetBinContent(ibinrec+1,eff);
710 hEfficiency->SetBinError(ibinrec+1,erreff);
716 return (TH1D*)hEfficiency;
719 //_________________________________________________________________________________________________________
720 void AliHFPtSpectrum::EstimateAndSetDirectEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco) {
722 // Function that computes the Direct acceptance and efficiency correction
723 // based on the simulated and reconstructed spectra
724 // and using the reconstructed spectra bin width
726 // eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
729 if(!fhRECpt || !hSimu || !hReco){
730 AliError("Hey, the reconstructed histogram was not set yet !");
734 fhDirectEffpt = EstimateEfficiencyRecoBin(hSimu,hReco,"fhDirectEffpt");
735 fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance #times efficiency");
739 //_________________________________________________________________________________________________________
740 void AliHFPtSpectrum::EstimateAndSetFeedDownEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco) {
742 // Function that computes the Feed-Down acceptance and efficiency correction
743 // based on the simulated and reconstructed spectra
744 // and using the reconstructed spectra bin width
746 // eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
749 if(!fhRECpt || !hSimu || !hReco){
750 AliError("Hey, the reconstructed histogram was not set yet !");
754 fhFeedDownEffpt = EstimateEfficiencyRecoBin(hSimu,hReco,"fhFeedDownEffpt");
755 fhFeedDownEffpt->SetNameTitle("fhFeedDownEffpt"," feed-down acceptance #times efficiency");
759 //_________________________________________________________________________________________________________
760 Bool_t AliHFPtSpectrum::Initialize(){
762 // Initialization of the variables (histograms)
765 if (fFeedDownOption==0) {
766 AliInfo("Getting ready for the corrections without feed-down consideration");
767 } else if (fFeedDownOption==1) {
768 AliInfo("Getting ready for the fc feed-down correction calculation");
769 } else if (fFeedDownOption==2) {
770 AliInfo("Getting ready for the Nb feed-down correction calculation");
771 } else { AliError("The calculation option must be <=2"); return kFALSE; }
773 // Start checking the input histograms consistency
774 Bool_t areconsistent=kTRUE;
777 if (!fhDirectEffpt || !fhRECpt) {
778 AliError(" Reconstructed spectra and/or the Nc efficiency distributions are not defined");
781 areconsistent &= CheckHistosConsistency(fhRECpt,fhDirectEffpt);
782 if (!areconsistent) {
783 AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
786 if (fFeedDownOption==0) return kTRUE;
789 // Common checks for options 1 (fc) & 2(Nb)
790 if (!fhFeedDownMCpt || !fhFeedDownEffpt) {
791 AliError(" Theoretical Nb and/or the Nb efficiency distributions are not defined");
794 areconsistent &= CheckHistosConsistency(fhRECpt,fhFeedDownMCpt);
795 areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownEffpt);
796 if (fAsymUncertainties) {
797 if (!fhFeedDownMCptMax || !fhFeedDownMCptMin) {
798 AliError(" Max/Min theoretical Nb distributions are not defined");
801 areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownMCptMax);
803 if (!areconsistent) {
804 AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
807 if (fFeedDownOption>1) return kTRUE;
810 // Now checks for option 1 (fc correction)
812 AliError("Theoretical Nc distributions is not defined");
815 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhFeedDownMCpt);
816 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectEffpt);
817 if (fAsymUncertainties) {
818 if (!fhDirectMCptMax || !fhDirectMCptMin) {
819 AliError(" Max/Min theoretical Nc distributions are not defined");
822 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectMCptMax);
824 if (!areconsistent) {
825 AliInfo("Histograms required for fc correction are not consistent (bin width, bounds)");
832 //_________________________________________________________________________________________________________
833 Bool_t AliHFPtSpectrum::CheckHistosConsistency(TH1D *h1, TH1D *h2){
835 // Check the histograms consistency (bins, limits)
839 AliError("One or both histograms don't exist");
843 Double_t binwidth1 = h1->GetBinWidth(1);
844 Double_t binwidth2 = h2->GetBinWidth(1);
845 Double_t min1 = h1->GetBinCenter(1) - (binwidth1/2.) ;
846 // Double_t max1 = h1->GetBinCenter(nbins1) + (binwidth1/2.) ;
847 Double_t min2 = h2->GetBinCenter(1) - (binwidth2/2.) ;
848 // Double_t max2 = h2->GetBinCenter(nbins2) + (binwidth2/2.) ;
850 if (binwidth1!=binwidth2) {
851 AliInfo(" histograms with different bin width");
855 AliInfo(" histograms with different minimum");
859 // AliInfo(" histograms with different maximum");
866 //_________________________________________________________________________________________________________
867 void AliHFPtSpectrum::CalculateFeedDownCorrectionFc(){
869 // Compute fc factor and its uncertainties bin by bin
870 // fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
872 // uncertainties: (conservative) combine the upper/lower N_b & N_c predictions together
873 // (extreme) combine the upper N_b predictions with the lower N_c predictions & viceversa
874 // systematic uncertainty on the acceptance x efficiency b/c ratio are included
877 // define the variables
878 Int_t nbins = fhRECpt->GetNbinsX();
879 Double_t binwidth = fhRECpt->GetBinWidth(1);
880 Double_t *limits = new Double_t[nbins+1];
881 Double_t *binwidths = new Double_t[nbins];
883 for (Int_t i=1; i<=nbins; i++) {
884 binwidth = fhRECpt->GetBinWidth(i);
885 xlow = fhRECpt->GetBinLowEdge(i);
887 binwidths[i-1] = binwidth;
889 limits[nbins] = xlow + binwidth;
891 Double_t correction=1.;
892 Double_t theoryRatio=1.;
893 Double_t effRatio=1.;
894 Double_t correctionExtremeA=1., correctionExtremeB=1.;
895 Double_t theoryRatioExtremeA=1., theoryRatioExtremeB=1.;
896 Double_t correctionConservativeA=1., correctionConservativeB=1.;
897 Double_t theoryRatioConservativeA=1., theoryRatioConservativeB=1.;
898 Double_t correctionUnc=0.;
899 Double_t correctionExtremeAUnc=0., correctionExtremeBUnc=0.;
900 Double_t correctionConservativeAUnc=0., correctionConservativeBUnc=0.;
902 // declare the output histograms
903 fhFc = new TH1D("fhFc","fc correction factor",nbins,limits);
904 fhFcMax = new TH1D("fhFcMax","max fc correction factor",nbins,limits);
905 fhFcMin = new TH1D("fhFcMin","min fc correction factor",nbins,limits);
906 // two local control histograms
907 TH1D *hTheoryRatio = new TH1D("hTheoryRatio","Theoretical B-->D over c-->D (feed-down/direct) ratio",nbins,limits);
908 TH1D *hEffRatio = new TH1D("hEffRatio","Efficiency B-->D over c-->D (feed-down/direct) ratio",nbins,limits);
909 // and the output TGraphAsymmErrors
910 if (fAsymUncertainties) {
911 fgFcExtreme = new TGraphAsymmErrors(nbins+1);
912 fgFcExtreme->SetNameTitle("fgFcExtreme","fgFcExtreme");
913 fgFcConservative = new TGraphAsymmErrors(nbins+1);
914 fgFcConservative->SetNameTitle("fgFcConservative","fgFcConservative");
921 for (Int_t ibin=1; ibin<=nbins; ibin++) {
923 // theory_ratio = (N_b/N_c)
924 theoryRatio = (fhDirectMCpt->GetBinContent(ibin)>0. && fhFeedDownMCpt->GetBinContent(ibin)>0.) ?
925 fhFeedDownMCpt->GetBinContent(ibin) / fhDirectMCpt->GetBinContent(ibin) : 1.0 ;
928 // Calculate the uncertainty [ considering only the theoretical uncertainties on Nb & Nc for now !!! ]
930 // extreme A = direct-max, feed-down-min
931 theoryRatioExtremeA = (fhDirectMCptMax->GetBinContent(ibin)>0. && fhFeedDownMCptMin->GetBinContent(ibin)>0.) ?
932 fhFeedDownMCptMin->GetBinContent(ibin) / fhDirectMCptMax->GetBinContent(ibin) : 1.0 ;
933 // extreme B = direct-min, feed-down-max
934 theoryRatioExtremeB = (fhDirectMCptMin->GetBinContent(ibin)>0. && fhDirectMCptMax->GetBinContent(ibin)>0.) ?
935 fhFeedDownMCptMax->GetBinContent(ibin) / fhDirectMCptMin->GetBinContent(ibin) : 1.0 ;
936 // conservative A = direct-max, feed-down-max
937 theoryRatioConservativeA = (fhDirectMCptMax->GetBinContent(ibin)>0. && fhFeedDownMCptMin->GetBinContent(ibin)>0.) ?
938 fhFeedDownMCptMax->GetBinContent(ibin) / fhDirectMCptMax->GetBinContent(ibin) : 1.0 ;
939 // conservative B = direct-min, feed-down-min
940 theoryRatioConservativeB = (fhDirectMCptMin->GetBinContent(ibin)>0. && fhDirectMCptMax->GetBinContent(ibin)>0.) ?
941 fhFeedDownMCptMin->GetBinContent(ibin) / fhDirectMCptMin->GetBinContent(ibin) : 1.0 ;
943 // eff_ratio = (eff_b/eff_c)
944 effRatio = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ?
945 fhFeedDownEffpt->GetBinContent(ibin) / fhDirectEffpt->GetBinContent(ibin) : 1.0 ;
947 // fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
948 if( TMath::Abs(effRatio - 1.0)<0.01 || TMath::Abs(theoryRatio - 1.0)<0.01 ) {
950 correctionExtremeA = 1.0;
951 correctionExtremeB = 1.0;
952 correctionConservativeA = 1.0;
953 correctionConservativeB = 1.0;
956 correction = ( 1. / ( 1 + ( effRatio * theoryRatio ) ) );
957 correctionExtremeA = ( 1. / ( 1 + ( effRatio * theoryRatioExtremeA ) ) );
958 correctionExtremeB = ( 1. / ( 1 + ( effRatio * theoryRatioExtremeB ) ) );
959 correctionConservativeA = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeA ) ) );
960 correctionConservativeB = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeB ) ) );
964 // fc uncertainty from (eff_b/eff_c) = fc^2 * (N_b/N_c) * delta(eff_b/eff_c)
965 // delta(eff_b/eff_c) is a percentage = fGlobalEfficiencyUncertainties[1]*effRatio
966 // correctionUnc = correction*correction * theoryRatio * fGlobalEfficiencyUncertainties[1]*effRatio;
967 // correctionExtremeAUnc = correctionExtremeA*correctionExtremeA * theoryRatioExtremeA * fGlobalEfficiencyUncertainties[1]*effRatio;
968 // correctionExtremeBUnc = correctionExtremeB*correctionExtremeB * theoryRatioExtremeB * fGlobalEfficiencyUncertainties[1]*effRatio;
969 // correctionConservativeAUnc = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA * fGlobalEfficiencyUncertainties[1]*effRatio;
970 // correctionConservativeBUnc = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB * fGlobalEfficiencyUncertainties[1]*effRatio;
971 correctionUnc = correction*correction * theoryRatio * effRatio *
972 TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] +
973 (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
974 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))
976 correctionExtremeAUnc = correctionExtremeA*correctionExtremeA * theoryRatioExtremeA * effRatio *
977 TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] +
978 (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
979 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))
981 correctionExtremeBUnc = correctionExtremeB*correctionExtremeB * theoryRatioExtremeB * effRatio *
982 TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] +
983 (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
984 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))
986 correctionConservativeAUnc = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA *effRatio *
987 TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] +
988 (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
989 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))
991 correctionConservativeBUnc = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB *effRatio *
992 TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] +
993 (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
994 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))
998 // Fill in the histograms
999 hTheoryRatio->SetBinContent(ibin,theoryRatio);
1000 hEffRatio->SetBinContent(ibin,effRatio);
1001 fhFc->SetBinContent(ibin,correction);
1002 if (fAsymUncertainties) {
1003 Double_t x = fhDirectMCpt->GetBinCenter(ibin);
1004 Double_t val[4] = { correctionExtremeA + correctionExtremeAUnc, correctionExtremeA - correctionExtremeAUnc,
1005 correctionExtremeB + correctionExtremeBUnc, correctionExtremeB - correctionExtremeBUnc };
1006 Double_t uncExtremeMin = correction - TMath::MinElement(4,val);
1007 Double_t uncExtremeMax = TMath::MaxElement(4,val) - correction;
1008 fgFcExtreme->SetPoint(ibin,x,correction); // i,x,y
1009 fgFcExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),uncExtremeMin,uncExtremeMax); // i,xl,xh,yl,yh
1010 fhFcMax->SetBinContent(ibin,correction+uncExtremeMax);
1011 fhFcMin->SetBinContent(ibin,correction-uncExtremeMin);
1012 Double_t consval[4] = { correctionConservativeA - correctionConservativeAUnc, correctionConservativeA + correctionConservativeAUnc,
1013 correctionConservativeB - correctionConservativeBUnc, correctionConservativeB + correctionConservativeBUnc};
1014 Double_t uncConservativeMin = correction - TMath::MinElement(4,consval);
1015 Double_t uncConservativeMax = TMath::MaxElement(4,consval) - correction;
1016 fgFcConservative->SetPoint(ibin,x,correction); // i,x,y
1017 fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),uncConservativeMin,uncConservativeMax); // i,xl,xh,yl,yh
1018 if( !(correction>0.) ){
1019 fgFcExtreme->SetPoint(ibin,x,0.); // i,x,y
1020 fgFcExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.); // i,xl,xh,yl,yh
1021 fgFcConservative->SetPoint(ibin,x,0.); // i,x,y
1022 fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.); // i,xl,xh,yl,yh
1027 delete [] binwidths;
1032 //_________________________________________________________________________________________________________
1033 void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumFc(){
1035 // Compute the feed-down corrected spectrum if feed-down correction is done via fc factor (bin by bin)
1036 // physics = reco * fc / bin-width
1038 // uncertainty: (stat) delta_physics = physics * sqrt ( (delta_reco/reco)^2 )
1039 // (syst but feed-down) delta_physics = physics * sqrt ( (delta_reco_syst/reco)^2 )
1040 // (feed-down syst) delta_physics = physics * sqrt ( (delta_fc/fc)^2 )
1042 // ( Calculation done bin by bin )
1044 if (!fhFc || !fhRECpt) {
1045 AliError(" Reconstructed or fc distributions are not defined");
1049 Int_t nbins = fhRECpt->GetNbinsX();
1050 Double_t value = 0., errvalue = 0., errvalueMax= 0., errvalueMin= 0.;
1051 Double_t valueExtremeMax= 0., valueExtremeMin= 0.;
1052 Double_t valueConservativeMax= 0., valueConservativeMin= 0.;
1053 Double_t binwidth = fhRECpt->GetBinWidth(1);
1054 Double_t *limits = new Double_t[nbins+1];
1055 Double_t *binwidths = new Double_t[nbins];
1057 for (Int_t i=1; i<=nbins; i++) {
1058 binwidth = fhRECpt->GetBinWidth(i);
1059 xlow = fhRECpt->GetBinLowEdge(i);
1061 binwidths[i-1] = binwidth;
1063 limits[nbins] = xlow + binwidth;
1065 // declare the output histograms
1066 fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by fc)",nbins,limits);
1067 fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by fc)",nbins,limits);
1068 fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by fc)",nbins,limits);
1069 // and the output TGraphAsymmErrors
1070 if (fAsymUncertainties){
1071 fgYieldCorr = new TGraphAsymmErrors(nbins+1);
1072 fgYieldCorrExtreme = new TGraphAsymmErrors(nbins+1);
1073 fgYieldCorrConservative = new TGraphAsymmErrors(nbins+1);
1077 // Do the calculation
1079 for (Int_t ibin=1; ibin<=nbins; ibin++) {
1081 // calculate the value
1082 // physics = reco * fc / bin-width
1083 value = (fhRECpt->GetBinContent(ibin) && fhFc->GetBinContent(ibin)) ?
1084 fhRECpt->GetBinContent(ibin) * fhFc->GetBinContent(ibin) : 0. ;
1085 value /= fhRECpt->GetBinWidth(ibin) ;
1087 // Statistical uncertainty
1088 // (stat) delta_physics = physics * sqrt ( (delta_reco/reco)^2 )
1089 errvalue = (value!=0. && fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) ?
1090 value * (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) : 0. ;
1092 // Calculate the systematic uncertainties
1093 // (syst but feed-down) delta_physics = physics * sqrt ( (delta_reco_syst/reco)^2 )
1094 // (feed-down syst) delta_physics = physics * sqrt ( (delta_fc/fc)^2 )
1096 // Protect against null denominator. If so, define uncertainty as null
1097 if (fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) {
1099 if (fAsymUncertainties) {
1101 // Systematics but feed-down
1102 if (fgRECSystematics) {
1103 errvalueMax = value * ( fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinContent(ibin) );
1104 errvalueMin = value * ( fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinContent(ibin) );
1106 else { errvalueMax = 0.; errvalueMin = 0.; }
1108 // Extreme feed-down systematics
1109 valueExtremeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcExtreme->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1110 valueExtremeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcExtreme->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1112 // Conservative feed-down systematics
1113 valueConservativeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcConservative->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1114 valueConservativeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcConservative->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1119 else { errvalueMax = 0.; errvalueMin = 0.; }
1121 // fill in the histograms
1122 fhYieldCorr->SetBinContent(ibin,value);
1123 fhYieldCorr->SetBinError(ibin,errvalue);
1124 if (fAsymUncertainties) {
1125 Double_t center = fhYieldCorr->GetBinCenter(ibin);
1126 fgYieldCorr->SetPoint(ibin,center,value); // i,x,y
1127 fgYieldCorr->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
1128 fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax);
1129 fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
1130 fgYieldCorrExtreme->SetPoint(ibin,center,value);
1131 fgYieldCorrExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),value-valueExtremeMin,valueExtremeMax-value);
1132 fgYieldCorrConservative->SetPoint(ibin,center,value);
1133 fgYieldCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),value-valueConservativeMin,valueConservativeMax-value);
1137 delete [] binwidths;
1142 //_________________________________________________________________________________________________________
1143 void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumNb(Double_t deltaY, Double_t branchingRatioBintoFinalDecay) {
1145 // Compute the feed-down corrected spectrum if feed-down correction is done via Nb (bin by bin)
1146 // physics = [ reco - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) ] / bin-width
1148 // uncertainty: (stat) delta_physics = sqrt ( (delta_reco)^2 ) / bin-width
1149 // (syst but feed-down) delta_physics = sqrt ( (delta_reco_syst)^2 ) / bin-width
1150 // (feed-down syst) delta_physics = sqrt ( (k*delta_lumi/lumi)^2 + (k*delta_eff_trig/eff_trig)^2
1151 // + (k*delta_Nb/Nb)^2 + (k*delta_eff/eff)^2 + (k*global_eff_ratio)^2 ) / bin-width
1152 // where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th
1155 Int_t nbins = fhRECpt->GetNbinsX();
1156 Double_t binwidth = fhRECpt->GetBinWidth(1);
1157 Double_t value = 0., errvalue = 0., errvalueMax = 0., errvalueMin = 0., kfactor = 0.;
1158 Double_t errvalueExtremeMax = 0., errvalueExtremeMin = 0.;
1159 Double_t *limits = new Double_t[nbins+1];
1160 Double_t *binwidths = new Double_t[nbins];
1162 for (Int_t i=1; i<=nbins; i++) {
1163 binwidth = fhRECpt->GetBinWidth(i);
1164 xlow = fhRECpt->GetBinLowEdge(i);
1166 binwidths[i-1] = binwidth;
1168 limits[nbins] = xlow + binwidth;
1170 // declare the output histograms
1171 fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by Nb)",nbins,limits);
1172 fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by Nb)",nbins,limits);
1173 fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by Nb)",nbins,limits);
1174 // and the output TGraphAsymmErrors
1175 if (fAsymUncertainties){
1176 fgYieldCorr = new TGraphAsymmErrors(nbins+1);
1177 fgYieldCorrExtreme = new TGraphAsymmErrors(nbins+1);
1178 fgYieldCorrConservative = new TGraphAsymmErrors(nbins+1);
1179 // Define fc-conservative
1180 fgFcConservative = new TGraphAsymmErrors(nbins+1);
1181 AliInfo(" Beware the conservative & extreme uncertainties are equal by definition !");
1184 // variables to define fc-conservative
1185 double correction=0, correctionMax=0., correctionMin=0.;
1188 // Do the calculation
1190 for (Int_t ibin=1; ibin<=nbins; ibin++) {
1192 // Calculate the value
1193 // physics = [ reco - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) ] / bin-width
1194 value = ( fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0. &&
1195 fhFeedDownMCpt->GetBinContent(ibin)>0. && fhFeedDownEffpt->GetBinContent(ibin)>0. ) ?
1196 fhRECpt->GetBinContent(ibin) - (deltaY*branchingRatioBintoFinalDecay*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) )
1198 value /= fhRECpt->GetBinWidth(ibin);
1200 // Statistical uncertainty: delta_physics = sqrt ( (delta_reco)^2 ) / bin-width
1201 errvalue = (value!=0. && fhRECpt->GetBinError(ibin) && fhRECpt->GetBinError(ibin)!=0.) ?
1202 fhRECpt->GetBinError(ibin) : 0.;
1203 errvalue /= fhRECpt->GetBinWidth(ibin);
1205 // Correction (fc) : Estimate of the relative amount feed-down subtracted
1206 // correction = [ 1 - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th)/reco ]
1207 correction = 1 - (deltaY*branchingRatioBintoFinalDecay*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) ) / fhRECpt->GetBinContent(ibin) ;
1209 // Systematic uncertainties
1210 // (syst but feed-down) delta_physics = sqrt ( (delta_reco_syst)^2 ) / bin-width
1211 // (feed-down syst) delta_physics = sqrt ( (k*delta_lumi/lumi)^2 + (k*delta_eff_trig/eff_trig)^2
1212 // + (k*delta_Nb/Nb)^2 + (k*delta_eff/eff)^2 + (k*global_eff_ratio)^2 ) / bin-width
1213 // where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th
1214 kfactor = deltaY*branchingRatioBintoFinalDecay*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) ;
1217 if (fAsymUncertainties) {
1218 Double_t nb = fhFeedDownMCpt->GetBinContent(ibin);
1219 Double_t nbDmax = fhFeedDownMCptMax->GetBinContent(ibin) - fhFeedDownMCpt->GetBinContent(ibin);
1220 Double_t nbDmin = fhFeedDownMCpt->GetBinContent(ibin) - fhFeedDownMCptMin->GetBinContent(ibin);
1222 // Systematics but feed-down
1223 if (fgRECSystematics){
1224 errvalueMax = fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinWidth(ibin) ;
1225 errvalueMin = fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinWidth(ibin);
1227 else { errvalueMax = 0.; errvalueMin = 0.; }
1229 // Feed-down systematics
1230 // min value with the maximum Nb
1231 errvalueExtremeMin = TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
1232 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1233 ( (kfactor*nbDmax/nb)*(kfactor*nbDmax/nb) ) +
1234 ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
1235 ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1236 ) / fhRECpt->GetBinWidth(ibin);
1237 // max value with the minimum Nb
1238 errvalueExtremeMax = TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
1239 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1240 ( (kfactor*nbDmin/nb)*(kfactor*nbDmin/nb) ) +
1241 ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
1242 ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1243 ) / fhRECpt->GetBinWidth(ibin);
1245 // Correction systematics (fc)
1246 // min value with the maximum Nb
1247 correctionMin = TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
1248 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1249 ( (kfactor*nbDmax/nb)*(kfactor*nbDmax/nb) ) +
1250 ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
1251 ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1252 ) / fhRECpt->GetBinContent(ibin) ;
1253 // max value with the minimum Nb
1254 correctionMax = TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
1255 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1256 ( (kfactor*nbDmin/nb)*(kfactor*nbDmin/nb) ) +
1257 ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
1258 ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1259 ) / fhRECpt->GetBinContent(ibin) ;
1261 else{ // Don't consider Nb uncertainty in this case [ to be tested!!! ]
1262 errvalueExtremeMax = TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
1263 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1264 ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
1265 ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1266 ) / fhRECpt->GetBinWidth(ibin);
1267 errvalueExtremeMin = errvalueExtremeMax ;
1271 // fill in histograms
1272 fhYieldCorr->SetBinContent(ibin,value);
1273 fhYieldCorr->SetBinError(ibin,errvalue);
1274 if (fAsymUncertainties) {
1275 Double_t x = fhYieldCorr->GetBinCenter(ibin);
1276 fgYieldCorr->SetPoint(ibin,x,value); // i,x,y
1277 fgYieldCorr->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
1278 fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax);
1279 fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
1280 fgYieldCorrExtreme->SetPoint(ibin,x,value); // i,x,y
1281 fgYieldCorrExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
1282 fgYieldCorrConservative->SetPoint(ibin,x,value); // i,x,y
1283 fgYieldCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
1284 // cout << " bin " << ibin << ", correction " << correction << ", min correction unc " << correctionMin << ", max correction unc " << correctionMax << endl;
1286 fgFcConservative->SetPoint(ibin,x,correction);
1287 fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),correctionMin,correctionMax);
1290 fgFcConservative->SetPoint(ibin,x,0.);
1291 fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.);
1296 delete [] binwidths;
1302 //_________________________________________________________________________________________________________
1303 void AliHFPtSpectrum::ComputeSystUncertainties(AliHFSystErr *systematics, Bool_t combineFeedDown) {
1305 // Function that re-calculates the global systematic uncertainties
1306 // by calling the class AliHFSystErr and combining those
1307 // (in quadrature) with the feed-down subtraction uncertainties
1311 // Estimate the feed-down uncertainty in percentage
1312 Int_t nentries = fgSigmaCorrConservative->GetN();
1313 TGraphAsymmErrors *grErrFeeddown = new TGraphAsymmErrors(nentries);
1314 Double_t x=0., y=0., errx=0., erryl=0., erryh=0;
1315 for(Int_t i=0; i<nentries; i++) {
1316 fgSigmaCorrConservative->GetPoint(i,x,y);
1317 errx = fgSigmaCorrConservative->GetErrorXlow(i) ;
1318 erryl = fgSigmaCorrConservative->GetErrorYlow(i) / y ;
1319 erryh = fgSigmaCorrConservative->GetErrorYhigh(i) / y ;
1320 grErrFeeddown->SetPoint(i,x,0.);
1321 grErrFeeddown->SetPointError(i,errx,errx,erryl,erryh); //i, xl, xh, yl, yh
1324 // Draw all the systematics independently
1325 systematics->DrawErrors(grErrFeeddown);
1327 // Set the sigma systematic uncertainties
1328 // possibly combine with the feed-down uncertainties
1329 Double_t errylcomb=0., erryhcomb=0;
1330 for(Int_t i=1; i<nentries; i++) {
1331 fgSigmaCorr->GetPoint(i,x,y);
1332 errx = grErrFeeddown->GetErrorXlow(i) ;
1333 erryl = grErrFeeddown->GetErrorYlow(i);
1334 erryh = grErrFeeddown->GetErrorYhigh(i);
1335 if (combineFeedDown) {
1336 errylcomb = systematics->GetTotalSystErr(x,erryl) * y ;
1337 erryhcomb = systematics->GetTotalSystErr(x,erryh) * y ;
1339 errylcomb = systematics->GetTotalSystErr(x) * y ;
1340 erryhcomb = systematics->GetTotalSystErr(x) * y ;
1342 fgSigmaCorr->SetPointError(i,errx,errx,errylcomb,erryhcomb);
1348 //_________________________________________________________________________________________________________
1349 void AliHFPtSpectrum::DrawSpectrum(TGraphAsymmErrors *gPrediction) {
1351 // Example method to draw the corrected spectrum & the theoretical prediction
1354 TCanvas *csigma = new TCanvas("csigma","Draw the corrected cross-section & the prediction");
1355 csigma->SetFillColor(0);
1356 gPrediction->GetXaxis()->SetTitleSize(0.05);
1357 gPrediction->GetXaxis()->SetTitleOffset(0.95);
1358 gPrediction->GetYaxis()->SetTitleSize(0.05);
1359 gPrediction->GetYaxis()->SetTitleOffset(0.95);
1360 gPrediction->GetXaxis()->SetTitle("p_{T} [GeV]");
1361 gPrediction->GetYaxis()->SetTitle("BR #times #frac{d#sigma}{dp_{T}} |_{|y|<0.5} [pb/GeV]");
1362 gPrediction->SetLineColor(kGreen+2);
1363 gPrediction->SetLineWidth(3);
1364 gPrediction->SetFillColor(kGreen+1);
1365 gPrediction->Draw("3CA");
1366 fgSigmaCorr->SetLineColor(kRed);
1367 fgSigmaCorr->SetLineWidth(1);
1368 fgSigmaCorr->SetFillColor(kRed);
1369 fgSigmaCorr->SetFillStyle(0);
1370 fgSigmaCorr->Draw("2");
1371 fhSigmaCorr->SetMarkerColor(kRed);
1372 fhSigmaCorr->Draw("esame");
1374 TLegend * leg = new TLegend(0.7,0.75,0.87,0.5);
1375 leg->SetBorderSize(0);
1376 leg->SetLineColor(0);
1377 leg->SetFillColor(0);
1378 leg->SetTextFont(42);
1379 leg->AddEntry(gPrediction,"FONLL ","fl");
1380 leg->AddEntry(fhSigmaCorr,"data stat. unc.","pl");
1381 leg->AddEntry(fgSigmaCorr,"data syst. unc.","f");
1387 //_________________________________________________________________________________________________________
1388 TH1D * AliHFPtSpectrum::ReweightHisto(TH1D *hToReweight, TH1D *hReference){
1390 // Function to reweight histograms for testing purposes:
1391 // This function takes the histo hToReweight and reweights
1392 // it (its pt shape) with respect to hReference
1395 // check histograms consistency
1396 Bool_t areconsistent=kTRUE;
1397 areconsistent &= CheckHistosConsistency(hToReweight,hReference);
1398 if (!areconsistent) {
1399 AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
1403 // define a new empty histogram
1404 TH1D *hReweighted = (TH1D*)hToReweight->Clone("hReweighted");
1405 hReweighted->Reset();
1406 Double_t weight=1.0;
1407 Double_t yvalue=1.0;
1408 Double_t integralRef = hReference->Integral();
1409 Double_t integralH = hToReweight->Integral();
1411 // now reweight the spectra
1413 // the weight is the relative probability of the given pt bin in the reference histo
1414 // divided by its relative probability (to normalize it) on the histo to re-weight
1415 for (Int_t i=0; i<=hToReweight->GetNbinsX(); i++) {
1416 weight = (hReference->GetBinContent(i)/integralRef) / (hToReweight->GetBinContent(i)/integralH) ;
1417 yvalue = hToReweight->GetBinContent(i);
1418 hReweighted->SetBinContent(i,yvalue*weight);
1421 return (TH1D*)hReweighted;
1424 //_________________________________________________________________________________________________________
1425 TH1D * AliHFPtSpectrum::ReweightRecHisto(TH1D *hRecToReweight, TH1D *hMCToReweight, TH1D *hMCReference){
1427 // Function to reweight histograms for testing purposes:
1428 // This function takes the histo hToReweight and reweights
1429 // it (its pt shape) with respect to hReference /hMCToReweight
1432 // check histograms consistency
1433 Bool_t areconsistent=kTRUE;
1434 areconsistent &= CheckHistosConsistency(hMCToReweight,hMCReference);
1435 areconsistent &= CheckHistosConsistency(hRecToReweight,hMCReference);
1436 if (!areconsistent) {
1437 AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
1441 // define a new empty histogram
1442 TH1D *hReweighted = (TH1D*)hMCToReweight->Clone("hReweighted");
1443 hReweighted->Reset();
1444 TH1D *hRecReweighted = (TH1D*)hRecToReweight->Clone("hRecReweighted");
1445 hRecReweighted->Reset();
1446 Double_t weight=1.0;
1447 Double_t yvalue=1.0, yrecvalue=1.0;
1448 Double_t integralRef = hMCReference->Integral();
1449 Double_t integralH = hMCToReweight->Integral();
1451 // now reweight the spectra
1453 // the weight is the relative probability of the given pt bin
1454 // that should be applied in the MC histo to get the reference histo shape
1455 // Probabilities are properly normalized.
1456 for (Int_t i=0; i<=hMCToReweight->GetNbinsX(); i++) {
1457 weight = (hMCReference->GetBinContent(i)/integralRef) / (hMCToReweight->GetBinContent(i)/integralH) ;
1458 yvalue = hMCToReweight->GetBinContent(i);
1459 hReweighted->SetBinContent(i,yvalue*weight);
1460 yrecvalue = hRecToReweight->GetBinContent(i);
1461 hRecReweighted->SetBinContent(i,yrecvalue*weight);
1464 return (TH1D*)hRecReweighted;