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 **************************************************************************/
16 //***********************************************************************
17 // Class AliHFPtSpectrum
18 // Base class for feed-down corrections on heavy-flavour decays
19 // computes the cross-section via one of the three implemented methods:
20 // 0) Consider no feed-down prediction
21 // 1) Subtract the feed-down with the "fc" method
22 // Yield = Reco * fc; where fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) ) ;
23 // 2) Subtract the feed-down with the "Nb" method
24 // Yield = Reco - Feed-down (exact formula on the function implementation)
26 // Author: Z.Conesa, zconesa@in2p3.fr
27 //***********************************************************************
29 #include <Riostream.h>
34 #include "TGraphAsymmErrors.h"
37 #include "AliHFPtSpectrum.h"
39 ClassImp(AliHFPtSpectrum)
41 //_________________________________________________________________________________________________________
42 AliHFPtSpectrum::AliHFPtSpectrum(const char* name, const char* title, Int_t option):
67 fFeedDownOption(option),
68 fAsymUncertainties(kTRUE)
71 // Default constructor
74 fLuminosity[0]=1.; fLuminosity[1]=0.;
75 fTrigEfficiency[0]=1.; fTrigEfficiency[1]=0.;
79 //_________________________________________________________________________________________________________
80 AliHFPtSpectrum::AliHFPtSpectrum(const AliHFPtSpectrum &rhs):
82 fhDirectMCpt(rhs.fhDirectMCpt),
83 fhFeedDownMCpt(rhs.fhFeedDownMCpt),
84 fhDirectMCptMax(rhs.fhDirectMCptMax),
85 fhDirectMCptMin(rhs.fhDirectMCptMin),
86 fhFeedDownMCptMax(rhs.fhFeedDownMCptMax),
87 fhFeedDownMCptMin(rhs.fhFeedDownMCptMin),
88 fhDirectEffpt(rhs.fhDirectEffpt),
89 fhFeedDownEffpt(rhs.fhFeedDownEffpt),
97 fhYieldCorr(rhs.fhYieldCorr),
98 fhYieldCorrMax(rhs.fhYieldCorrMax),
99 fhYieldCorrMin(rhs.fhYieldCorrMin),
100 fgYieldCorr(rhs.fgYieldCorr),
101 fhSigmaCorr(rhs.fhSigmaCorr),
102 fhSigmaCorrMax(rhs.fhSigmaCorrMax),
103 fhSigmaCorrMin(rhs.fhSigmaCorrMin),
104 fgSigmaCorr(rhs.fgSigmaCorr),
105 fFeedDownOption(rhs.fFeedDownOption),
106 fAsymUncertainties(rhs.fAsymUncertainties)
112 for(Int_t i=0; i<2; i++){
113 fLuminosity[i] = rhs.fLuminosity[i];
114 fTrigEfficiency[i] = rhs.fTrigEfficiency[i];
119 //_________________________________________________________________________________________________________
120 AliHFPtSpectrum &AliHFPtSpectrum::operator=(const AliHFPtSpectrum &source){
122 // Assignment operator
125 if (&source == this) return *this;
127 fhDirectMCpt = source.fhDirectMCpt;
128 fhFeedDownMCpt = source.fhFeedDownMCpt;
129 fhDirectMCptMax = source.fhDirectMCptMax;
130 fhDirectMCptMin = source.fhDirectMCptMin;
131 fhFeedDownMCptMax = source.fhFeedDownMCptMax;
132 fhFeedDownMCptMin = source.fhFeedDownMCptMin;
133 fhDirectEffpt = source.fhDirectEffpt;
134 fhFeedDownEffpt = source.fhFeedDownEffpt;
135 fhRECpt = source.fhRECpt;
137 fhFcMax = source.fhFcMax;
138 fhFcMin = source.fhFcMin;
140 fhYieldCorr = source.fhYieldCorr;
141 fhYieldCorrMax = source.fhYieldCorrMax;
142 fhYieldCorrMin = source.fhYieldCorrMin;
143 fgYieldCorr = source.fgYieldCorr;
144 fhSigmaCorr = source.fhSigmaCorr;
145 fhSigmaCorrMax = source.fhSigmaCorrMax;
146 fhSigmaCorrMin = source.fhSigmaCorrMin;
147 fgSigmaCorr = source.fgSigmaCorr;
148 fFeedDownOption = source.fFeedDownOption;
149 fAsymUncertainties = source.fAsymUncertainties;
151 for(Int_t i=0; i<2; i++){
152 fLuminosity[i] = source.fLuminosity[i];
153 fTrigEfficiency[i] = source.fTrigEfficiency[i];
159 //_________________________________________________________________________________________________________
160 AliHFPtSpectrum::~AliHFPtSpectrum(){
168 //_________________________________________________________________________________________________________
169 void AliHFPtSpectrum::SetMCptSpectra(TH1 *hDirect, TH1 *hFeedDown){
171 // Set the MonteCarlo or Theoretical spectra
172 // both for direct and feed-down contributions
175 if (!hDirect || !hFeedDown) {
176 AliError("One or both (direct, feed-down) spectra don't exist");
180 Bool_t areconsistent = kTRUE;
181 areconsistent = CheckHistosConsistency(hDirect,hFeedDown);
182 if (!areconsistent) {
183 AliInfo("Histograms are not consistent (bin width, bounds)");
187 fhDirectMCpt = hDirect;
188 fhFeedDownMCpt = hFeedDown;
191 //_________________________________________________________________________________________________________
192 void AliHFPtSpectrum::SetFeedDownMCptSpectra(TH1 *hFeedDown){
194 // Set the MonteCarlo or Theoretical spectra
195 // for feed-down contribution
199 AliError("Feed-down spectra don't exist");
202 fhFeedDownMCpt = hFeedDown;
205 //_________________________________________________________________________________________________________
206 void AliHFPtSpectrum::SetMCptDistributionsBounds(TH1 *hDirectMax, TH1 *hDirectMin, TH1 *hFeedDownMax, TH1 *hFeedDownMin){
208 // Set the maximum and minimum MonteCarlo or Theoretical spectra
209 // both for direct and feed-down contributions
210 // used in case uncertainties are asymmetric and ca not be on the "basic histograms"
213 if (!hDirectMax || !hDirectMin || !hFeedDownMax|| !hFeedDownMin) {
214 AliError("One or all of the max/min direct/feed-down spectra don't exist");
218 Bool_t areconsistent = kTRUE;
219 areconsistent &= CheckHistosConsistency(hDirectMax,hDirectMin);
220 areconsistent &= CheckHistosConsistency(hFeedDownMax,hFeedDownMin);
221 areconsistent &= CheckHistosConsistency(hDirectMax,hFeedDownMax);
222 if (!areconsistent) {
223 AliInfo("Histograms are not consistent (bin width, bounds)");
227 fhDirectMCptMax = hDirectMax;
228 fhDirectMCptMin = hDirectMin;
229 fhFeedDownMCptMax = hFeedDownMax;
230 fhFeedDownMCptMin = hFeedDownMin;
233 //_________________________________________________________________________________________________________
234 void AliHFPtSpectrum::SetFeedDownMCptDistributionsBounds(TH1 *hFeedDownMax, TH1 *hFeedDownMin){
236 // Set the maximum and minimum MonteCarlo or Theoretical spectra
237 // for feed-down contributions
238 // used in case uncertainties are asymmetric and can not be on the "basic histogram"
241 if (!hFeedDownMax || !hFeedDownMin) {
242 AliError("One or all of the max/min direct/feed-down spectra don't exist");
246 Bool_t areconsistent = kTRUE;
247 areconsistent &= CheckHistosConsistency(hFeedDownMax,hFeedDownMin);
248 if (!areconsistent) {
249 AliInfo("Histograms are not consistent (bin width, bounds)");
253 fhFeedDownMCptMax = hFeedDownMax;
254 fhFeedDownMCptMin = hFeedDownMin;
257 //_________________________________________________________________________________________________________
258 void AliHFPtSpectrum::SetDirectAccEffCorrection(TH1 *hDirectEff){
260 // Set the Acceptance and Efficiency corrections
261 // for the direct contribution
265 AliError("The direct acceptance and efficiency corrections doesn't exist");
269 fhDirectEffpt = hDirectEff;
272 //_________________________________________________________________________________________________________
273 void AliHFPtSpectrum::SetAccEffCorrection(TH1 *hDirectEff, TH1 *hFeedDownEff){
275 // Set the Acceptance and Efficiency corrections
276 // both for direct and feed-down contributions
279 if (!hDirectEff || !hFeedDownEff) {
280 AliError("One or both (direct, feed-down) acceptance and efficiency corrections don't exist");
284 Bool_t areconsistent=kTRUE;
285 areconsistent = CheckHistosConsistency(hDirectEff,hFeedDownEff);
286 if (!areconsistent) {
287 AliInfo("Histograms are not consistent (bin width, bounds)");
291 fhDirectEffpt = hDirectEff;
292 fhFeedDownEffpt = hFeedDownEff;
295 //_________________________________________________________________________________________________________
296 void AliHFPtSpectrum::SetReconstructedSpectrum(TH1 *hRec) {
298 // Set the reconstructed spectrum
302 AliError("The reconstructed spectrum doesn't exist");
309 //_________________________________________________________________________________________________________
310 void AliHFPtSpectrum::ComputeHFPtSpectrum(Double_t deltaY, Double_t branchingRatioC, Double_t branchingRatioBintoFinalDecay) {
312 // Main function to compute the corrected cross-section:
313 // variables : analysed delta_y, BR for the final correction,
314 // BR b --> D --> decay (relative to the input theoretical prediction)
316 // Sigma = ( 1. / (lumi * delta_y * BR_c * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
318 // Uncertainties: delta_sigma = sigma * sqrt ( (delta_reco/reco)^2 + (delta_lumi/lumi)^2 + (delta_eff_trig/eff_trig)^2 )
321 // First: Initialization
323 Bool_t areHistosOk = Initialize();
325 AliInfo(" Histos not properly initialized. Check : inconsistent binning ? missing histos ?");
330 // Second: Correct for feed-down
332 if (fFeedDownOption==1) {
333 // Compute the feed-down correction via fc-method
334 CalculateFeedDownCorrectionFc();
335 // Correct the yield for feed-down correction via fc-method
336 CalculateFeedDownCorrectedSpectrumFc();
338 else if (fFeedDownOption==2) {
339 // Correct the yield for feed-down correction via Nb-method
340 CalculateFeedDownCorrectedSpectrumNb(deltaY,branchingRatioBintoFinalDecay);
342 else if (fFeedDownOption==0) {
343 // If there is no need for feed-down correction,
344 // the "corrected" yield is equal to the raw yield
345 fhYieldCorr = fhRECpt;
346 fhYieldCorr->SetNameTitle("fhYieldCorr","un-corrected yield");
347 fhYieldCorrMax = fhRECpt;
348 fhYieldCorrMin = fhRECpt;
349 fhYieldCorrMax->SetNameTitle("fhYieldCorrMax","un-corrected yield");
350 fhYieldCorrMin->SetNameTitle("fhYieldCorrMin","un-corrected yield");
351 fAsymUncertainties=kFALSE;
354 AliInfo(" Are you sure the feed-down correction option is right ?");
357 // Print out information
358 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);
361 // Finally: Correct from yields to cross-section
363 Int_t nbins = fhRECpt->GetNbinsX();
364 Double_t binwidth = fhRECpt->GetBinWidth(1);
365 Double_t xmin = fhRECpt->GetBinCenter(1) - (binwidth/2.) ;
366 Double_t xmax = fhRECpt->GetBinCenter(nbins) + (binwidth/2.) ;
368 // declare the output histograms
369 TH1D *hSigmaCorr = new TH1D("hSigmaCorr","corrected sigma",nbins,xmin,xmax);
370 TH1D *hSigmaCorrMax = new TH1D("hSigmaCorrMax","max corrected sigma",nbins,xmin,xmax);
371 TH1D *hSigmaCorrMin = new TH1D("hSigmaCorrMin","min corrected sigma",nbins,xmin,xmax);
372 // and the output TGraphAsymmErrors
373 if (fAsymUncertainties & !fgSigmaCorr) fgSigmaCorr = new TGraphAsymmErrors(nbins);
375 // protect against null denominator
376 if (deltaY==0. || fLuminosity[0]==0. || fTrigEfficiency[0]==0. || branchingRatioC==0.) {
377 AliError(" Hey you ! Why luminosity or trigger-efficiency or the c-BR or delta_y are set to zero ?! ");
381 Double_t value=0, valueMax=0., valueMin=0.;
382 for(Int_t ibin=0; ibin<=nbins; ibin++){
384 // Sigma = ( 1. / (lumi * delta_y * BR_c * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
385 value = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ?
386 ( fhYieldCorr->GetBinContent(ibin) / ( deltaY * branchingRatioC * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
389 // Uncertainties: delta_sigma = sigma * sqrt ( (delta_reco/reco)^2 + (delta_lumi/lumi)^2 + (delta_eff_trig/eff_trig)^2 )
390 if (fAsymUncertainties) {
391 valueMax = value * TMath::Sqrt( (fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin))* (fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin)) +
392 (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
393 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) );
394 valueMin = value * TMath::Sqrt( (fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin))* (fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin)) +
395 (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
396 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) );
399 // protect against null denominator
400 valueMax = (value!=0.) ?
401 value * TMath::Sqrt( (fhYieldCorr->GetBinError(ibin)/fhYieldCorr->GetBinContent(ibin))* (fhYieldCorr->GetBinError(ibin)/fhYieldCorr->GetBinContent(ibin)) +
402 (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
403 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) )
408 // Fill the histograms
409 hSigmaCorr->SetBinContent(ibin,value);
410 hSigmaCorr->SetBinError(ibin,valueMax);
411 hSigmaCorrMax->SetBinContent(ibin,valueMax);
412 hSigmaCorrMin->SetBinContent(ibin,valueMin);
413 // Fill the TGraphAsymmErrors
414 if (fAsymUncertainties) {
415 Double_t x = fhYieldCorr->GetBinCenter(ibin);
416 fgSigmaCorr->SetPoint(ibin,x,value); // i,x,y
417 fgSigmaCorr->SetPointError(ibin,(binwidth/2.),(binwidth/2.),valueMin,valueMax); // i,xl,xh,yl,yh
422 fhSigmaCorr = hSigmaCorr ;
423 fhSigmaCorrMax = hSigmaCorrMax ;
424 fhSigmaCorrMin = hSigmaCorrMin ;
427 //_________________________________________________________________________________________________________
428 Bool_t AliHFPtSpectrum::Initialize(){
430 // Initialization of the variables (histograms)
433 if (fFeedDownOption==0) {
434 AliInfo("Getting ready for the corrections without feed-down consideration");
435 } else if (fFeedDownOption==1) {
436 AliInfo("Getting ready for the fc feed-down correction calculation");
437 } else if (fFeedDownOption==2) {
438 AliInfo("Getting ready for the Nb feed-down correction calculation");
439 } else { AliError("The calculation option must be <=2"); return kFALSE; }
441 // Start checking the input histograms consistency
442 Bool_t areconsistent=kTRUE;
445 if (!fhDirectEffpt || !fhRECpt) {
446 AliError(" Reconstructed spectra and/or the Nc efficiency distributions are not defined");
449 areconsistent &= CheckHistosConsistency(fhRECpt,fhDirectEffpt);
450 if (!areconsistent) {
451 AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
454 if (fFeedDownOption==0) return kTRUE;
457 // Common checks for options 1 (fc) & 2(Nb)
458 if (!fhFeedDownMCpt || !fhFeedDownEffpt) {
459 AliError(" Theoretical Nb and/or the Nb efficiency distributions are not defined");
462 areconsistent &= CheckHistosConsistency(fhRECpt,fhFeedDownMCpt);
463 areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownEffpt);
464 if (fAsymUncertainties) {
465 if (!fhFeedDownMCptMax || !fhFeedDownMCptMin) {
466 AliError(" Max/Min theoretical Nb distributions are not defined");
469 areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownMCptMax);
471 if (!areconsistent) {
472 AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
475 if (fFeedDownOption>1) return kTRUE;
478 // Now checks for option 1 (fc correction)
480 AliError("Theoretical Nc distributions is not defined");
483 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhFeedDownMCpt);
484 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectEffpt);
485 if (fAsymUncertainties) {
486 if (!fhDirectMCptMax || !fhDirectMCptMin) {
487 AliError(" Max/Min theoretical Nc distributions are not defined");
490 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectMCptMax);
492 if (!areconsistent) {
493 AliInfo("Histograms required for fc correction are not consistent (bin width, bounds)");
500 //_________________________________________________________________________________________________________
501 Bool_t AliHFPtSpectrum::CheckHistosConsistency(TH1 *h1, TH1 *h2){
503 // Check the histograms consistency (bins, limits)
507 AliError("One or both histograms don't exist");
511 Double_t binwidth1 = h1->GetBinWidth(1);
512 Double_t binwidth2 = h2->GetBinWidth(1);
513 Double_t min1 = h1->GetBinCenter(1) - (binwidth1/2.) ;
514 // Double_t max1 = h1->GetBinCenter(nbins1) + (binwidth1/2.) ;
515 Double_t min2 = h2->GetBinCenter(1) - (binwidth2/2.) ;
516 // Double_t max2 = h2->GetBinCenter(nbins2) + (binwidth2/2.) ;
518 if (binwidth1!=binwidth2) {
519 AliInfo(" histograms with different bin width");
523 AliInfo(" histograms with different minimum");
527 // AliInfo(" histograms with different maximum");
534 //_________________________________________________________________________________________________________
535 void AliHFPtSpectrum::CalculateFeedDownCorrectionFc(){
537 // Compute fc factor and its uncertainties bin by bin
538 // fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
541 // define the variables
542 Int_t nbins = fhRECpt->GetNbinsX();
543 Double_t binwidth = fhRECpt->GetBinWidth(1);
544 Double_t xmin = fhRECpt->GetBinCenter(1) - (binwidth/2.) ;
545 Double_t xmax = fhRECpt->GetBinCenter(nbins) + (binwidth/2.) ;
546 Double_t correction=1.;
547 Double_t correctionMax=1., correctionMin=1.;
548 Double_t theoryRatio=1.;
549 Double_t effRatio=1.;
551 // declare the output histograms
552 TH1D *hfc = new TH1D("hfc","fc correction factor",nbins,xmin,xmax);
553 TH1D *hfcMax = new TH1D("hfcMax","max fc correction factor",nbins,xmin,xmax);
554 TH1D *hfcMin = new TH1D("hfcMin","min fc correction factor",nbins,xmin,xmax);
555 // two local control histograms
556 TH1D *hTheoryRatio = new TH1D("hTheoryRatio","Theoretical B-->D over c-->D (feed-down/direct) ratio",nbins,xmin,xmax);
557 TH1D *hEffRatio = new TH1D("hEffRatio","Efficiency B-->D over c-->D (feed-down/direct) ratio",nbins,xmin,xmax);
558 // and the output TGraphAsymmErrors
559 if (fAsymUncertainties & !fgFc) fgFc = new TGraphAsymmErrors(nbins);
564 for (Int_t ibin=0; ibin<=nbins; ibin++) {
566 // theory_ratio = (N_b/N_c)
567 theoryRatio = (fhDirectMCpt->GetBinContent(ibin) && fhDirectMCpt->GetBinContent(ibin)!=0.) ? fhFeedDownMCpt->GetBinContent(ibin) / fhDirectMCpt->GetBinContent(ibin) : 1.0 ;
568 // eff_ratio = (eff_b/eff_c)
569 effRatio = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ? fhFeedDownEffpt->GetBinContent(ibin) / fhDirectEffpt->GetBinContent(ibin) : 1.0 ;
570 // fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
571 correction = (effRatio && theoryRatio) ? ( 1. / ( 1 + ( effRatio * theoryRatio ) ) ) : 1.0 ;
573 // Calculate the uncertainty [ considering only the theoretical uncertainties on Nb & Nc for now !!! ]
574 // delta_fc = fc^2 * (Nb/Nc) * sqrt ( (delta_Nb/Nb)^2 + (delta_Nc/Nc)^2 )
575 Double_t deltaNbMax = fhFeedDownMCptMax->GetBinContent(ibin) - fhFeedDownMCpt->GetBinContent(ibin) ;
576 Double_t deltaNbMin = fhFeedDownMCpt->GetBinContent(ibin) - fhFeedDownMCptMin->GetBinContent(ibin) ;
577 Double_t deltaNcMax = fhDirectMCptMax->GetBinContent(ibin) - fhDirectMCpt->GetBinContent(ibin) ;
578 Double_t deltaNcMin = fhDirectMCpt->GetBinContent(ibin) - fhDirectMCptMin->GetBinContent(ibin) ;
580 // Protect against null denominator. If so, define uncertainty as null
581 if (fhFeedDownMCpt->GetBinContent(ibin) && fhFeedDownMCpt->GetBinContent(ibin)!=0. &&
582 fhDirectMCpt->GetBinContent(ibin) && fhDirectMCpt->GetBinContent(ibin)!=0. ) {
583 correctionMax = correction*correction*theoryRatio *
585 (deltaNbMax/fhFeedDownMCpt->GetBinContent(ibin))*(deltaNbMax/fhFeedDownMCpt->GetBinContent(ibin)) +
586 (deltaNcMax/fhDirectMCpt->GetBinContent(ibin))*(deltaNcMax/fhDirectMCpt->GetBinContent(ibin))
588 correctionMin = correction*correction*theoryRatio *
590 (deltaNbMin/fhFeedDownMCpt->GetBinContent(ibin))*(deltaNbMin/fhFeedDownMCpt->GetBinContent(ibin)) +
591 (deltaNcMin/fhDirectMCpt->GetBinContent(ibin))*(deltaNcMin/fhDirectMCpt->GetBinContent(ibin))
594 else { correctionMax = 0.; correctionMin = 0.; }
597 // Fill in the histograms
598 hTheoryRatio->SetBinContent(ibin,theoryRatio);
599 hEffRatio->SetBinContent(ibin,effRatio);
600 hfc->SetBinContent(ibin,correction);
601 hfcMax->SetBinContent(ibin,correction+correctionMax);
602 hfcMin->SetBinContent(ibin,correction-correctionMin);
603 if (fAsymUncertainties) {
604 Double_t x = fhDirectMCpt->GetBinCenter(ibin);
605 fgFc->SetPoint(ibin,x,correction); // i,x,y
606 fgFc->SetPointError(ibin,(binwidth/2.),(binwidth/2.),correctionMin,correctionMax); // i,xl,xh,yl,yh
616 //_________________________________________________________________________________________________________
617 void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumFc(){
619 // Compute the feed-down corrected spectrum if feed-down correction is done via fc factor (bin by bin)
620 // physics = reco * fc
622 // uncertainty: delta_physics = physics * sqrt ( (delta_reco/reco)^2 + (delta_fc/fc)^2 )
624 // ( Calculation done bin by bin )
626 if (!fhFc || !fhRECpt) {
627 AliError(" Reconstructed or fc distributions are not defined");
631 Int_t nbins = fhRECpt->GetNbinsX();
632 Double_t value = 0., valueDmax= 0., valueDmin= 0.;
633 Double_t binwidth = fhRECpt->GetBinWidth(1);
634 Double_t xmin = fhRECpt->GetBinCenter(1) - (binwidth/2.) ;
635 Double_t xmax = fhRECpt->GetBinCenter(nbins) + (binwidth/2.) ;
637 // declare the output histograms
638 TH1D *hYield = new TH1D("hYield","corrected yield (by fc)",nbins,xmin,xmax);
639 TH1D *hYieldMax = new TH1D("hYieldMax","max corrected yield (by fc)",nbins,xmin,xmax);
640 TH1D *hYieldMin = new TH1D("hYieldMin","min corrected yield (by fc)",nbins,xmin,xmax);
641 // and the output TGraphAsymmErrors
642 if (fAsymUncertainties & !fgYieldCorr) fgYieldCorr = new TGraphAsymmErrors(nbins);
645 // Do the calculation
647 for (Int_t ibin=0; ibin<=nbins; ibin++) {
649 // calculate the value
650 value = fhRECpt->GetBinContent(ibin) * fhFc->GetBinContent(ibin) ;
652 // calculate the value uncertainty
653 // Protect against null denominator. If so, define uncertainty as null
654 if (fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) {
656 if (fAsymUncertainties) {
658 if (fhFc->GetBinContent(ibin) && fhFc->GetBinContent(ibin)!=0.) {
659 valueDmax = value * TMath::Sqrt( ( (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin))*(fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) ) + ( (fgFc->GetErrorYhigh(ibin)/fhFc->GetBinContent(ibin))*(fgFc->GetErrorYhigh(ibin)/fhFc->GetBinContent(ibin)) ) );
660 valueDmin = value * TMath::Sqrt( ( (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin))*(fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) ) + ( (fgFc->GetErrorYlow(ibin)/fhFc->GetBinContent(ibin))*(fgFc->GetErrorYlow(ibin)/fhFc->GetBinContent(ibin)) ) );
662 else { valueDmax = 0.; valueDmin = 0.; }
665 else { // Don't consider fc uncertainty in this case [ to be tested!!! ]
666 valueDmax = value * (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) ;
667 valueDmin = value * (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) ;
671 else { valueDmax = 0.; valueDmin = 0.; }
673 // fill in the histograms
674 hYield->SetBinContent(ibin,value);
675 hYield->SetBinError(ibin,valueDmax);
676 hYieldMax->SetBinContent(ibin,value+valueDmax);
677 hYieldMin->SetBinContent(ibin,value-valueDmin);
678 if (fAsymUncertainties) {
679 Double_t center = hYield->GetBinCenter(ibin);
680 fgYieldCorr->SetPoint(ibin,center,value); // i,x,y
681 fgYieldCorr->SetPointError(ibin,(binwidth/2.),(binwidth/2.),valueDmin,valueDmax); // i,xl,xh,yl,yh
686 fhYieldCorr = hYield;
687 fhYieldCorrMax = hYieldMax;
688 fhYieldCorrMin = hYieldMin;
691 //_________________________________________________________________________________________________________
692 void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumNb(Double_t deltaY, Double_t branchingRatioBintoFinalDecay) {
694 // Compute the feed-down corrected spectrum if feed-down correction is done via Nb (bin by bin)
695 // physics = reco - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th)
697 // uncertainty: delta_physics = sqrt ( (delta_reco)^2 + (k*delta_lumi/lumi)^2 +
698 // (k*delta_eff_trig/eff_trig)^2 + (k*delta_Nb/Nb)^2 )
699 // where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th
702 Int_t nbins = fhRECpt->GetNbinsX();
703 Double_t binwidth = fhRECpt->GetBinWidth(1);
704 Double_t value = 0., valueDmax = 0., valueDmin = 0., kfactor = 0.;
705 Double_t xmin = fhRECpt->GetBinCenter(1) - (binwidth/2.) ;
706 Double_t xmax = fhRECpt->GetBinCenter(nbins) + (binwidth/2.) ;
708 // declare the output histograms
709 TH1D *hYield = new TH1D("hYield","corrected yield (by Nb)",nbins,xmin,xmax);
710 TH1D *hYieldMax = new TH1D("hYieldMax","max corrected yield (by Nb)",nbins,xmin,xmax);
711 TH1D *hYieldMin = new TH1D("hYieldMin","min corrected yield (by Nb)",nbins,xmin,xmax);
712 // and the output TGraphAsymmErrors
713 if (fAsymUncertainties & !fgYieldCorr) fgYieldCorr = new TGraphAsymmErrors(nbins);
716 // Do the calculation
718 for (Int_t ibin=0; ibin<=nbins; ibin++) {
720 // calculate the value
721 value = fhRECpt->GetBinContent(ibin) - (deltaY*branchingRatioBintoFinalDecay*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) );
723 kfactor = deltaY*branchingRatioBintoFinalDecay*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) ;
725 // calculate the value uncertainty
726 if (fAsymUncertainties) {
727 Double_t Nb = fhFeedDownMCpt->GetBinContent(ibin);
728 Double_t NbDmax = fhFeedDownMCptMax->GetBinContent(ibin) - fhFeedDownMCpt->GetBinContent(ibin);
729 Double_t NbDmin = fhFeedDownMCpt->GetBinContent(ibin) - fhFeedDownMCptMin->GetBinContent(ibin);
730 valueDmax = TMath::Sqrt( ( fhRECpt->GetBinError(ibin)*fhRECpt->GetBinError(ibin) ) +
731 ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
732 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
733 ( (kfactor*NbDmax/Nb)*(kfactor*NbDmax/Nb) ) );
734 valueDmin = TMath::Sqrt( ( fhRECpt->GetBinError(ibin)*fhRECpt->GetBinError(ibin) ) +
735 ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
736 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
737 ( (kfactor*NbDmin/Nb)*(kfactor*NbDmin/Nb) ) );
739 else{ // Don't consider Nb uncertainty in this case [ to be tested!!! ]
740 valueDmax = TMath::Sqrt( ( fhRECpt->GetBinError(ibin)*fhRECpt->GetBinError(ibin) ) +
741 ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
742 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) );
743 valueDmin = valueDmax ;
746 // fill in histograms
747 hYield->SetBinContent(ibin,value);
748 hYield->SetBinError(ibin,valueDmax);
749 hYieldMax->SetBinContent(ibin,value+valueDmax);
750 hYieldMin->SetBinContent(ibin,value-valueDmin);
751 if (fAsymUncertainties) {
752 Double_t x = hYield->GetBinCenter(ibin);
753 fgYieldCorr->SetPoint(ibin,x,value); // i,x,y
754 fgYieldCorr->SetPointError(ibin,(binwidth/2.),(binwidth/2.),valueDmin,valueDmax); // i,xl,xh,yl,yh
759 fhYieldCorr = hYield;
760 fhYieldCorrMax = hYieldMax;
761 fhYieldCorrMin = hYieldMin;
765 //_________________________________________________________________________________________________________
766 TH1 * AliHFPtSpectrum::ReweightHisto(TH1 *hToReweight, TH1 *hReference){
768 // Function to reweight histograms for testing purposes:
769 // This function takes the histo hToReweight and reweights
770 // it (its pt shape) with respect to hReference
773 // check histograms consistency
774 Bool_t areconsistent=kTRUE;
775 areconsistent &= CheckHistosConsistency(hToReweight,hReference);
776 if (!areconsistent) {
777 AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
781 // define a new empty histogram
782 TH1 *hReweighted = (TH1*)hToReweight->Clone("hReweighted");
783 hReweighted->Reset();
786 Double_t integralRef = hReference->Integral();
787 Double_t integralH = hToReweight->Integral();
789 // now reweight the spectra
791 // the weight is the relative probability of the given pt bin in the reference histo
792 // divided by its relative probability (to normalize it) on the histo to re-weight
793 for (Int_t i=0; i<=hToReweight->GetNbinsX(); i++) {
794 weight = (hReference->GetBinContent(i)/integralRef) / (hToReweight->GetBinContent(i)/integralH) ;
795 yvalue = hToReweight->GetBinContent(i);
796 hReweighted->SetBinContent(i,yvalue*weight);
799 return (TH1*)hReweighted;
802 //_________________________________________________________________________________________________________
803 TH1 * AliHFPtSpectrum::ReweightRecHisto(TH1 *hRecToReweight, TH1 *hMCToReweight, TH1 *hMCReference){
805 // Function to reweight histograms for testing purposes:
806 // This function takes the histo hToReweight and reweights
807 // it (its pt shape) with respect to hReference /hMCToReweight
810 // check histograms consistency
811 Bool_t areconsistent=kTRUE;
812 areconsistent &= CheckHistosConsistency(hMCToReweight,hMCReference);
813 areconsistent &= CheckHistosConsistency(hRecToReweight,hMCReference);
814 if (!areconsistent) {
815 AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
819 // define a new empty histogram
820 TH1 *hReweighted = (TH1*)hMCToReweight->Clone("hReweighted");
821 hReweighted->Reset();
822 TH1 *hRecReweighted = (TH1*)hRecToReweight->Clone("hRecReweighted");
823 hRecReweighted->Reset();
825 Double_t yvalue=1.0, yrecvalue=1.0;
826 Double_t integralRef = hMCReference->Integral();
827 Double_t integralH = hMCToReweight->Integral();
829 // now reweight the spectra
831 // the weight is the relative probability of the given pt bin
832 // that should be applied in the MC histo to get the reference histo shape
833 // Probabilities are properly normalized.
834 for (Int_t i=0; i<=hMCToReweight->GetNbinsX(); i++) {
835 weight = (hMCReference->GetBinContent(i)/integralRef) / (hMCToReweight->GetBinContent(i)/integralH) ;
836 yvalue = hMCToReweight->GetBinContent(i);
837 hReweighted->SetBinContent(i,yvalue*weight);
838 yrecvalue = hRecToReweight->GetBinContent(i);
839 hRecReweighted->SetBinContent(i,yrecvalue*weight);
842 return (TH1*)hRecReweighted;