]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/vertexingHF/AliHFPtSpectrum.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Bug fix for poly fit (ChiaraB)
[u/mrichter/AliRoot.git] / PWG3 / vertexingHF / AliHFPtSpectrum.cxx
CommitLineData
65e55bbd 1/**************************************************************************
2 * Copyright(c) 1998-2010, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
6 * *
7 * Permission to use, copy, modify and distribute this software and its *
8 * documentation strictly for non-commercial purposes is hereby granted *
9 * without fee, provided that the above copyright notice appears in all *
10 * copies and that both the copyright notice and this permission notice *
11 * appear in the supporting documentation. The authors make no claims *
12 * about the suitability of this software for any purpose. It is *
13 * provided "as is" without express or implied warranty. *
14 **************************************************************************/
15
16//***********************************************************************
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)
25//
bb427707 26// (the corrected yields per bin are divided by the bin-width)
27//
65e55bbd 28// Author: Z.Conesa, zconesa@in2p3.fr
29//***********************************************************************
30
31#include <Riostream.h>
32
33#include "TMath.h"
34#include "TH1.h"
35#include "TH1D.h"
36#include "TGraphAsymmErrors.h"
bb427707 37#include "TNamed.h"
65e55bbd 38
39#include "AliLog.h"
40#include "AliHFPtSpectrum.h"
41
42ClassImp(AliHFPtSpectrum)
43
44//_________________________________________________________________________________________________________
45AliHFPtSpectrum::AliHFPtSpectrum(const char* name, const char* title, Int_t option):
46 TNamed(name,title),
47 fhDirectMCpt(),
48 fhFeedDownMCpt(),
a17b17dd 49 fhDirectMCptMax(),
50 fhDirectMCptMin(),
51 fhFeedDownMCptMax(),
52 fhFeedDownMCptMin(),
65e55bbd 53 fhDirectEffpt(),
54 fhFeedDownEffpt(),
55 fhRECpt(),
56 fLuminosity(),
57 fTrigEfficiency(),
58 fhFc(),
a17b17dd 59 fhFcMax(),
60 fhFcMin(),
65e55bbd 61 fgFc(),
62 fhYieldCorr(),
a17b17dd 63 fhYieldCorrMax(),
64 fhYieldCorrMin(),
65e55bbd 65 fgYieldCorr(),
66 fhSigmaCorr(),
a17b17dd 67 fhSigmaCorrMax(),
68 fhSigmaCorrMin(),
65e55bbd 69 fgSigmaCorr(),
70 fFeedDownOption(option),
71 fAsymUncertainties(kTRUE)
72{
73 //
74 // Default constructor
75 //
76
77 fLuminosity[0]=1.; fLuminosity[1]=0.;
78 fTrigEfficiency[0]=1.; fTrigEfficiency[1]=0.;
79
80}
81
82//_________________________________________________________________________________________________________
83AliHFPtSpectrum::AliHFPtSpectrum(const AliHFPtSpectrum &rhs):
84 TNamed(rhs),
85 fhDirectMCpt(rhs.fhDirectMCpt),
86 fhFeedDownMCpt(rhs.fhFeedDownMCpt),
a17b17dd 87 fhDirectMCptMax(rhs.fhDirectMCptMax),
88 fhDirectMCptMin(rhs.fhDirectMCptMin),
89 fhFeedDownMCptMax(rhs.fhFeedDownMCptMax),
90 fhFeedDownMCptMin(rhs.fhFeedDownMCptMin),
65e55bbd 91 fhDirectEffpt(rhs.fhDirectEffpt),
92 fhFeedDownEffpt(rhs.fhFeedDownEffpt),
93 fhRECpt(rhs.fhRECpt),
94 fLuminosity(),
95 fTrigEfficiency(),
96 fhFc(rhs.fhFc),
a17b17dd 97 fhFcMax(rhs.fhFcMax),
98 fhFcMin(rhs.fhFcMin),
65e55bbd 99 fgFc(rhs.fgFc),
100 fhYieldCorr(rhs.fhYieldCorr),
a17b17dd 101 fhYieldCorrMax(rhs.fhYieldCorrMax),
102 fhYieldCorrMin(rhs.fhYieldCorrMin),
65e55bbd 103 fgYieldCorr(rhs.fgYieldCorr),
104 fhSigmaCorr(rhs.fhSigmaCorr),
a17b17dd 105 fhSigmaCorrMax(rhs.fhSigmaCorrMax),
106 fhSigmaCorrMin(rhs.fhSigmaCorrMin),
65e55bbd 107 fgSigmaCorr(rhs.fgSigmaCorr),
108 fFeedDownOption(rhs.fFeedDownOption),
109 fAsymUncertainties(rhs.fAsymUncertainties)
110{
111 //
112 // Copy constructor
113 //
114
115 for(Int_t i=0; i<2; i++){
116 fLuminosity[i] = rhs.fLuminosity[i];
117 fTrigEfficiency[i] = rhs.fTrigEfficiency[i];
118 }
119
120}
121
122//_________________________________________________________________________________________________________
123AliHFPtSpectrum &AliHFPtSpectrum::operator=(const AliHFPtSpectrum &source){
124 //
125 // Assignment operator
126 //
127
128 if (&source == this) return *this;
129
130 fhDirectMCpt = source.fhDirectMCpt;
131 fhFeedDownMCpt = source.fhFeedDownMCpt;
a17b17dd 132 fhDirectMCptMax = source.fhDirectMCptMax;
133 fhDirectMCptMin = source.fhDirectMCptMin;
134 fhFeedDownMCptMax = source.fhFeedDownMCptMax;
135 fhFeedDownMCptMin = source.fhFeedDownMCptMin;
65e55bbd 136 fhDirectEffpt = source.fhDirectEffpt;
137 fhFeedDownEffpt = source.fhFeedDownEffpt;
138 fhRECpt = source.fhRECpt;
139 fhFc = source.fhFc;
a17b17dd 140 fhFcMax = source.fhFcMax;
141 fhFcMin = source.fhFcMin;
65e55bbd 142 fgFc = source.fgFc;
143 fhYieldCorr = source.fhYieldCorr;
a17b17dd 144 fhYieldCorrMax = source.fhYieldCorrMax;
145 fhYieldCorrMin = source.fhYieldCorrMin;
65e55bbd 146 fgYieldCorr = source.fgYieldCorr;
147 fhSigmaCorr = source.fhSigmaCorr;
a17b17dd 148 fhSigmaCorrMax = source.fhSigmaCorrMax;
149 fhSigmaCorrMin = source.fhSigmaCorrMin;
65e55bbd 150 fgSigmaCorr = source.fgSigmaCorr;
151 fFeedDownOption = source.fFeedDownOption;
152 fAsymUncertainties = source.fAsymUncertainties;
153
154 for(Int_t i=0; i<2; i++){
155 fLuminosity[i] = source.fLuminosity[i];
156 fTrigEfficiency[i] = source.fTrigEfficiency[i];
157 }
158
159 return *this;
160}
161
162//_________________________________________________________________________________________________________
163AliHFPtSpectrum::~AliHFPtSpectrum(){
164 //
165 // Destructor
166 //
167 ;
168}
169
170
5f3c1b97 171//_________________________________________________________________________________________________________
172TH1 * AliHFPtSpectrum::RebinTheoreticalSpectra(TH1 *hTheory) {
173 //
174 // Function to rebin the theoretical spectrum
175 // with respect to the real-data reconstructed spectrum binning
176 //
177
178 if (!hTheory || !fhRECpt) {
179 AliError("Feed-down or reconstructed spectra don't exist");
180 return NULL;
181 }
182
183 //
184 // Get the reconstructed spectra bins & limits
185 Int_t nbins = fhRECpt->GetNbinsX();
186 Int_t nbinsMC = hTheory->GetNbinsX();
187 Double_t *limits = new Double_t[nbins+1];
188 Double_t xlow=0., binwidth=0.;
189 for (Int_t i=0; i<=nbins; i++) {
190 binwidth = fhRECpt->GetBinWidth(i);
191 xlow = fhRECpt->GetBinLowEdge(i);
192 limits[i-1] = xlow;
193 }
194 limits[nbins] = xlow + binwidth;
195
196 //
197 // Define a new histogram with the real-data reconstructed spectrum binning
198 TH1D * hTheoryRebin = new TH1D("hTheoryRebin"," theoretical rebinned prediction",nbins,limits);
199
200 Double_t sum[nbins], items[nbins];
201 for (Int_t ibin=0; ibin<=nbinsMC; ibin++){
202
203 for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++){
204 if( hTheory->GetBinCenter(ibin)>limits[ibinrec] &&
205 hTheory->GetBinCenter(ibin)<limits[ibinrec+1]){
206 sum[ibinrec]+=hTheory->GetBinContent(ibin);
207 items[ibinrec]+=1.;
208 }
209 }
210
211 }
212
213 // set the theoretical rebinned spectra to ( sum-bins / n-bins ) per new bin
214 for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++) {
215 hTheoryRebin->SetBinContent(ibinrec+1,sum[ibinrec]/items[ibinrec]);
216 }
217
218 return (TH1*)hTheoryRebin;
219}
220
65e55bbd 221//_________________________________________________________________________________________________________
222void AliHFPtSpectrum::SetMCptSpectra(TH1 *hDirect, TH1 *hFeedDown){
223 //
224 // Set the MonteCarlo or Theoretical spectra
225 // both for direct and feed-down contributions
226 //
227
bb427707 228 if (!hDirect || !hFeedDown || !fhRECpt) {
229 AliError("One or both (direct, feed-down) spectra or the reconstructed spectra don't exist");
65e55bbd 230 return;
231 }
232
233 Bool_t areconsistent = kTRUE;
234 areconsistent = CheckHistosConsistency(hDirect,hFeedDown);
235 if (!areconsistent) {
236 AliInfo("Histograms are not consistent (bin width, bounds)");
237 return;
238 }
239
bb427707 240 //
241 // If the predictions shape do not have the same
242 // bin width (check via number of bins) as the reconstructed spectra, change them
bb427707 243 if (hDirect->GetBinWidth(1) == fhRECpt->GetBinWidth(1)) {
244
245 fhDirectMCpt = hDirect;
246 fhFeedDownMCpt = hFeedDown;
247
248 }
249 else {
250
5f3c1b97 251 fhDirectMCpt = RebinTheoreticalSpectra(hDirect);
252 fhDirectMCpt->SetNameTitle("fhDirectMCpt"," direct theoretical prediction");
253 fhFeedDownMCpt = RebinTheoreticalSpectra(hFeedDown);
254 fhFeedDownMCpt->SetNameTitle("fhFeedDownMCpt"," feed-down theoretical prediction");
bb427707 255
bb427707 256 }
257
65e55bbd 258}
259
260//_________________________________________________________________________________________________________
261void AliHFPtSpectrum::SetFeedDownMCptSpectra(TH1 *hFeedDown){
262 //
263 // Set the MonteCarlo or Theoretical spectra
264 // for feed-down contribution
265 //
266
bb427707 267 if (!hFeedDown || !fhRECpt) {
268 AliError("Feed-down or reconstructed spectra don't exist");
65e55bbd 269 return;
270 }
bb427707 271
bb427707 272 //
273 // If the predictions shape do not have the same
274 // bin width (check via number of bins) as the reconstructed spectra, change them
bb427707 275 if (hFeedDown->GetBinWidth(1) == fhRECpt->GetBinWidth(1)) {
276
277 fhFeedDownMCpt = hFeedDown;
278
279 }
280 else {
281
5f3c1b97 282 fhFeedDownMCpt = RebinTheoreticalSpectra(hFeedDown);
283 fhFeedDownMCpt->SetNameTitle("fhFeedDownMCpt"," feed-down theoretical prediction");
bb427707 284
bb427707 285 }
286
65e55bbd 287}
288
289//_________________________________________________________________________________________________________
290void AliHFPtSpectrum::SetMCptDistributionsBounds(TH1 *hDirectMax, TH1 *hDirectMin, TH1 *hFeedDownMax, TH1 *hFeedDownMin){
291 //
292 // Set the maximum and minimum MonteCarlo or Theoretical spectra
293 // both for direct and feed-down contributions
294 // used in case uncertainties are asymmetric and ca not be on the "basic histograms"
295 //
296
bb427707 297 if (!hDirectMax || !hDirectMin || !hFeedDownMax|| !hFeedDownMin || !fhRECpt) {
298 AliError("One or all of the max/min direct/feed-down or the reconstructed spectra don't exist");
65e55bbd 299 return;
300 }
301
302 Bool_t areconsistent = kTRUE;
303 areconsistent &= CheckHistosConsistency(hDirectMax,hDirectMin);
304 areconsistent &= CheckHistosConsistency(hFeedDownMax,hFeedDownMin);
305 areconsistent &= CheckHistosConsistency(hDirectMax,hFeedDownMax);
306 if (!areconsistent) {
307 AliInfo("Histograms are not consistent (bin width, bounds)");
308 return;
309 }
310
bb427707 311
312 //
313 // If the predictions shape do not have the same
314 // bin width (check via number of bins) as the reconstructed spectra, change them
bb427707 315 if (hFeedDownMax->GetBinWidth(1) == fhRECpt->GetBinWidth(1)) {
316
317 fhDirectMCptMax = hDirectMax;
318 fhDirectMCptMin = hDirectMin;
319 fhFeedDownMCptMax = hFeedDownMax;
320 fhFeedDownMCptMin = hFeedDownMin;
321
322 }
323 else {
324
5f3c1b97 325 fhDirectMCptMax = RebinTheoreticalSpectra(hDirectMax);
326 fhDirectMCptMax->SetNameTitle("fhDirectMCptMax"," maximum direct theoretical prediction");
327 fhDirectMCptMin = RebinTheoreticalSpectra(hDirectMin);
328 fhDirectMCptMin->SetNameTitle("fhDirectMCptMin"," minimum direct theoretical prediction");
329 fhFeedDownMCptMax = RebinTheoreticalSpectra(hFeedDownMax);
330 fhFeedDownMCptMax->SetNameTitle("fhFeedDownMCptMax"," maximum feed-down theoretical prediction");
331 fhFeedDownMCptMin = RebinTheoreticalSpectra(hFeedDownMin);
332 fhFeedDownMCptMin->SetNameTitle("fhFeedDownMCptMin"," minimum feed-down theoretical prediction");
bb427707 333
bb427707 334 }
335
65e55bbd 336}
337
338//_________________________________________________________________________________________________________
339void AliHFPtSpectrum::SetFeedDownMCptDistributionsBounds(TH1 *hFeedDownMax, TH1 *hFeedDownMin){
340 //
341 // Set the maximum and minimum MonteCarlo or Theoretical spectra
342 // for feed-down contributions
343 // used in case uncertainties are asymmetric and can not be on the "basic histogram"
344 //
345
bb427707 346 if (!hFeedDownMax || !hFeedDownMin || !fhRECpt) {
65e55bbd 347 AliError("One or all of the max/min direct/feed-down spectra don't exist");
348 return;
349 }
350
351 Bool_t areconsistent = kTRUE;
352 areconsistent &= CheckHistosConsistency(hFeedDownMax,hFeedDownMin);
353 if (!areconsistent) {
354 AliInfo("Histograms are not consistent (bin width, bounds)");
355 return;
356 }
357
bb427707 358
359 //
360 // If the predictions shape do not have the same
361 // bin width (check via number of bins) as the reconstructed spectra, change them
bb427707 362 if (hFeedDownMax->GetBinWidth(1) == fhRECpt->GetBinWidth(1)) {
363
364 fhFeedDownMCptMax = hFeedDownMax;
365 fhFeedDownMCptMin = hFeedDownMin;
366
367 }
368 else {
369
5f3c1b97 370 fhFeedDownMCptMax = RebinTheoreticalSpectra(hFeedDownMax);
371 fhFeedDownMCptMax->SetNameTitle("fhFeedDownMCptMax"," maximum feed-down theoretical prediction");
372 fhFeedDownMCptMin = RebinTheoreticalSpectra(hFeedDownMin);
373 fhFeedDownMCptMin->SetNameTitle("fhFeedDownMCptMin"," minimum feed-down theoretical prediction");
bb427707 374
bb427707 375 }
376
65e55bbd 377}
378
379//_________________________________________________________________________________________________________
380void AliHFPtSpectrum::SetDirectAccEffCorrection(TH1 *hDirectEff){
381 //
382 // Set the Acceptance and Efficiency corrections
383 // for the direct contribution
384 //
385
386 if (!hDirectEff) {
387 AliError("The direct acceptance and efficiency corrections doesn't exist");
388 return;
389 }
390
391 fhDirectEffpt = hDirectEff;
392}
393
394//_________________________________________________________________________________________________________
395void AliHFPtSpectrum::SetAccEffCorrection(TH1 *hDirectEff, TH1 *hFeedDownEff){
396 //
397 // Set the Acceptance and Efficiency corrections
398 // both for direct and feed-down contributions
399 //
400
401 if (!hDirectEff || !hFeedDownEff) {
402 AliError("One or both (direct, feed-down) acceptance and efficiency corrections don't exist");
403 return;
404 }
405
406 Bool_t areconsistent=kTRUE;
407 areconsistent = CheckHistosConsistency(hDirectEff,hFeedDownEff);
408 if (!areconsistent) {
409 AliInfo("Histograms are not consistent (bin width, bounds)");
410 return;
411 }
412
413 fhDirectEffpt = hDirectEff;
414 fhFeedDownEffpt = hFeedDownEff;
415}
416
417//_________________________________________________________________________________________________________
418void AliHFPtSpectrum::SetReconstructedSpectrum(TH1 *hRec) {
419 //
420 // Set the reconstructed spectrum
421 //
422
423 if (!hRec) {
424 AliError("The reconstructed spectrum doesn't exist");
425 return;
426 }
427
428 fhRECpt = hRec;
429}
430
431//_________________________________________________________________________________________________________
a17b17dd 432void AliHFPtSpectrum::ComputeHFPtSpectrum(Double_t deltaY, Double_t branchingRatioC, Double_t branchingRatioBintoFinalDecay) {
65e55bbd 433 //
434 // Main function to compute the corrected cross-section:
a17b17dd 435 // variables : analysed delta_y, BR for the final correction,
436 // BR b --> D --> decay (relative to the input theoretical prediction)
65e55bbd 437 //
438 // Sigma = ( 1. / (lumi * delta_y * BR_c * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
439 //
440 // Uncertainties: delta_sigma = sigma * sqrt ( (delta_reco/reco)^2 + (delta_lumi/lumi)^2 + (delta_eff_trig/eff_trig)^2 )
441
442 //
443 // First: Initialization
444 //
445 Bool_t areHistosOk = Initialize();
446 if (!areHistosOk) {
447 AliInfo(" Histos not properly initialized. Check : inconsistent binning ? missing histos ?");
448 return;
449 }
450
451 //
452 // Second: Correct for feed-down
453 //
454 if (fFeedDownOption==1) {
455 // Compute the feed-down correction via fc-method
a17b17dd 456 CalculateFeedDownCorrectionFc();
65e55bbd 457 // Correct the yield for feed-down correction via fc-method
a17b17dd 458 CalculateFeedDownCorrectedSpectrumFc();
65e55bbd 459 }
460 else if (fFeedDownOption==2) {
461 // Correct the yield for feed-down correction via Nb-method
a17b17dd 462 CalculateFeedDownCorrectedSpectrumNb(deltaY,branchingRatioBintoFinalDecay);
65e55bbd 463 }
464 else if (fFeedDownOption==0) {
465 // If there is no need for feed-down correction,
466 // the "corrected" yield is equal to the raw yield
467 fhYieldCorr = fhRECpt;
468 fhYieldCorr->SetNameTitle("fhYieldCorr","un-corrected yield");
a17b17dd 469 fhYieldCorrMax = fhRECpt;
470 fhYieldCorrMin = fhRECpt;
471 fhYieldCorrMax->SetNameTitle("fhYieldCorrMax","un-corrected yield");
472 fhYieldCorrMin->SetNameTitle("fhYieldCorrMin","un-corrected yield");
65e55bbd 473 fAsymUncertainties=kFALSE;
474 }
475 else {
476 AliInfo(" Are you sure the feed-down correction option is right ?");
477 }
478
479 // Print out information
a17b17dd 480 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);
65e55bbd 481
482 //
483 // Finally: Correct from yields to cross-section
484 //
485 Int_t nbins = fhRECpt->GetNbinsX();
486 Double_t binwidth = fhRECpt->GetBinWidth(1);
bb427707 487 Double_t *limits = new Double_t[nbins+1];
de97013c 488 Double_t *binwidths = new Double_t[nbins+1];
bb427707 489 Double_t xlow=0.;
490 for (Int_t i=0; i<=nbins; i++) {
491 binwidth = fhRECpt->GetBinWidth(i);
492 xlow = fhRECpt->GetBinLowEdge(i);
493 limits[i-1] = xlow;
de97013c 494 binwidths[i] = binwidth;
bb427707 495 }
496 limits[nbins] = xlow + binwidth;
497
65e55bbd 498
499 // declare the output histograms
bb427707 500 TH1D *hSigmaCorr = new TH1D("hSigmaCorr","corrected sigma",nbins,limits);
501 TH1D *hSigmaCorrMax = new TH1D("hSigmaCorrMax","max corrected sigma",nbins,limits);
502 TH1D *hSigmaCorrMin = new TH1D("hSigmaCorrMin","min corrected sigma",nbins,limits);
65e55bbd 503 // and the output TGraphAsymmErrors
504 if (fAsymUncertainties & !fgSigmaCorr) fgSigmaCorr = new TGraphAsymmErrors(nbins);
505
506 // protect against null denominator
a17b17dd 507 if (deltaY==0. || fLuminosity[0]==0. || fTrigEfficiency[0]==0. || branchingRatioC==0.) {
65e55bbd 508 AliError(" Hey you ! Why luminosity or trigger-efficiency or the c-BR or delta_y are set to zero ?! ");
509 return ;
510 }
511
a17b17dd 512 Double_t value=0, valueMax=0., valueMin=0.;
65e55bbd 513 for(Int_t ibin=0; ibin<=nbins; ibin++){
514
515 // Sigma = ( 1. / (lumi * delta_y * BR_c * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
516 value = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ?
a17b17dd 517 ( fhYieldCorr->GetBinContent(ibin) / ( deltaY * branchingRatioC * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
65e55bbd 518 : 0. ;
bb427707 519 // cout << " bin="<< ibin << " yield-corr="<< fhYieldCorr->GetBinContent(ibin) <<", eff_D="<< fhDirectEffpt->GetBinContent(ibin) <<", value="<<value<<endl;
65e55bbd 520
521 // Uncertainties: delta_sigma = sigma * sqrt ( (delta_reco/reco)^2 + (delta_lumi/lumi)^2 + (delta_eff_trig/eff_trig)^2 )
522 if (fAsymUncertainties) {
a17b17dd 523 valueMax = value * TMath::Sqrt( (fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin))* (fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin)) +
65e55bbd 524 (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
525 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) );
a17b17dd 526 valueMin = value * TMath::Sqrt( (fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin))* (fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin)) +
65e55bbd 527 (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
528 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) );
529 }
530 else {
531 // protect against null denominator
a17b17dd 532 valueMax = (value!=0.) ?
65e55bbd 533 value * TMath::Sqrt( (fhYieldCorr->GetBinError(ibin)/fhYieldCorr->GetBinContent(ibin))* (fhYieldCorr->GetBinError(ibin)/fhYieldCorr->GetBinContent(ibin)) +
534 (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
535 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) )
536 : 0. ;
a17b17dd 537 valueMin = valueMax;
65e55bbd 538 }
539
540 // Fill the histograms
541 hSigmaCorr->SetBinContent(ibin,value);
a17b17dd 542 hSigmaCorr->SetBinError(ibin,valueMax);
543 hSigmaCorrMax->SetBinContent(ibin,valueMax);
544 hSigmaCorrMin->SetBinContent(ibin,valueMin);
65e55bbd 545 // Fill the TGraphAsymmErrors
546 if (fAsymUncertainties) {
547 Double_t x = fhYieldCorr->GetBinCenter(ibin);
548 fgSigmaCorr->SetPoint(ibin,x,value); // i,x,y
de97013c 549 fgSigmaCorr->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),valueMin,valueMax); // i,xl,xh,yl,yh
65e55bbd 550 }
551
552 }
553
554 fhSigmaCorr = hSigmaCorr ;
a17b17dd 555 fhSigmaCorrMax = hSigmaCorrMax ;
556 fhSigmaCorrMin = hSigmaCorrMin ;
65e55bbd 557}
558
bb427707 559//_________________________________________________________________________________________________________
5f3c1b97 560TH1 * AliHFPtSpectrum::EstimateEfficiencyRecoBin(TH1 *hSimu, TH1 *hReco) {
bb427707 561 //
5f3c1b97 562 // Function that computes the acceptance and efficiency correction
bb427707 563 // based on the simulated and reconstructed spectra
564 // and using the reconstructed spectra bin width
565 //
566 // eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
567 //
568
569 if(!fhRECpt){
570 AliInfo("Hey, the reconstructed histogram was not set yet !");
571 return NULL;
572 }
573
574 Int_t nbins = fhRECpt->GetNbinsX();
575 Double_t *limits = new Double_t[nbins+1];
576 Double_t xlow=0.,binwidth=0.;
577 for (Int_t i=0; i<=nbins; i++) {
578 binwidth = fhRECpt->GetBinWidth(i);
579 xlow = fhRECpt->GetBinLowEdge(i);
580 limits[i-1] = xlow;
581 }
582 limits[nbins] = xlow + binwidth;
583
5f3c1b97 584 TH1D * hEfficiency = new TH1D("hEfficiency"," acceptance #times efficiency",nbins,limits);
bb427707 585
586 Double_t sumSimu[nbins], sumReco[nbins];
587 for (Int_t ibin=0; ibin<=hSimu->GetNbinsX(); ibin++){
588
589 for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++){
590 if ( hSimu->GetBinCenter(ibin)>limits[ibinrec] &&
591 hSimu->GetBinCenter(ibin)<limits[ibinrec+1] ) {
592 sumSimu[ibinrec]+=hSimu->GetBinContent(ibin);
593 }
594 if ( hReco->GetBinCenter(ibin)>limits[ibinrec] &&
595 hReco->GetBinCenter(ibin)<limits[ibinrec+1] ) {
596 sumReco[ibinrec]+=hReco->GetBinContent(ibin);
597 }
598 }
599
600 }
601
602 // the efficiency is computed as reco/sim (in each bin)
603 // its uncertainty is err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
604 Double_t eff=0., erreff=0.;
605 for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++) {
606 eff = sumReco[ibinrec] / sumSimu[ibinrec] ;
607 erreff = TMath::Sqrt( eff * (1.0 - eff) ) / TMath::Sqrt( sumSimu[ibinrec] );
5f3c1b97 608 hEfficiency->SetBinContent(ibinrec+1,eff);
609 hEfficiency->SetBinError(ibinrec+1,erreff);
bb427707 610 }
611
5f3c1b97 612 return (TH1*)hEfficiency;
bb427707 613}
614
615//_________________________________________________________________________________________________________
5f3c1b97 616TH1 * AliHFPtSpectrum::EstimateAndSetDirectEfficiencyRecoBin(TH1 *hSimu, TH1 *hReco) {
bb427707 617 //
5f3c1b97 618 // Function that computes the Direct acceptance and efficiency correction
bb427707 619 // based on the simulated and reconstructed spectra
620 // and using the reconstructed spectra bin width
621 //
622 // eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
623 //
5f3c1b97 624
bb427707 625 if(!fhRECpt){
626 AliInfo("Hey, the reconstructed histogram was not set yet !");
627 return NULL;
628 }
629
5f3c1b97 630 fhDirectEffpt = EstimateEfficiencyRecoBin(hSimu,hReco);
631 fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance #times efficiency");
bb427707 632
5f3c1b97 633 return (TH1*)fhDirectEffpt;
634}
bb427707 635
5f3c1b97 636//_________________________________________________________________________________________________________
637TH1 * AliHFPtSpectrum::EstimateAndSetFeedDownEfficiencyRecoBin(TH1 *hSimu, TH1 *hReco) {
638 //
639 // Function that computes the Feed-Down acceptance and efficiency correction
640 // based on the simulated and reconstructed spectra
641 // and using the reconstructed spectra bin width
642 //
643 // eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
644 //
bb427707 645
5f3c1b97 646 if(!fhRECpt){
647 AliInfo("Hey, the reconstructed histogram was not set yet !");
648 return NULL;
bb427707 649 }
650
5f3c1b97 651 fhFeedDownEffpt = EstimateEfficiencyRecoBin(hSimu,hReco);
652 fhFeedDownEffpt->SetNameTitle("fhFeedDownEffpt"," feed-down acceptance #times efficiency");
653
bb427707 654 return (TH1*)fhFeedDownEffpt;
655}
656
65e55bbd 657//_________________________________________________________________________________________________________
658Bool_t AliHFPtSpectrum::Initialize(){
659 //
660 // Initialization of the variables (histograms)
661 //
662
663 if (fFeedDownOption==0) {
664 AliInfo("Getting ready for the corrections without feed-down consideration");
665 } else if (fFeedDownOption==1) {
666 AliInfo("Getting ready for the fc feed-down correction calculation");
667 } else if (fFeedDownOption==2) {
668 AliInfo("Getting ready for the Nb feed-down correction calculation");
669 } else { AliError("The calculation option must be <=2"); return kFALSE; }
670
671 // Start checking the input histograms consistency
672 Bool_t areconsistent=kTRUE;
673
674 // General checks
675 if (!fhDirectEffpt || !fhRECpt) {
676 AliError(" Reconstructed spectra and/or the Nc efficiency distributions are not defined");
677 return kFALSE;
678 }
679 areconsistent &= CheckHistosConsistency(fhRECpt,fhDirectEffpt);
680 if (!areconsistent) {
681 AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
682 return kFALSE;
683 }
684 if (fFeedDownOption==0) return kTRUE;
685
686 //
687 // Common checks for options 1 (fc) & 2(Nb)
688 if (!fhFeedDownMCpt || !fhFeedDownEffpt) {
689 AliError(" Theoretical Nb and/or the Nb efficiency distributions are not defined");
690 return kFALSE;
691 }
692 areconsistent &= CheckHistosConsistency(fhRECpt,fhFeedDownMCpt);
693 areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownEffpt);
694 if (fAsymUncertainties) {
a17b17dd 695 if (!fhFeedDownMCptMax || !fhFeedDownMCptMin) {
65e55bbd 696 AliError(" Max/Min theoretical Nb distributions are not defined");
697 return kFALSE;
698 }
a17b17dd 699 areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownMCptMax);
65e55bbd 700 }
701 if (!areconsistent) {
702 AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
703 return kFALSE;
704 }
705 if (fFeedDownOption>1) return kTRUE;
706
707 //
708 // Now checks for option 1 (fc correction)
709 if (!fhDirectMCpt) {
710 AliError("Theoretical Nc distributions is not defined");
711 return kFALSE;
712 }
713 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhFeedDownMCpt);
714 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectEffpt);
715 if (fAsymUncertainties) {
a17b17dd 716 if (!fhDirectMCptMax || !fhDirectMCptMin) {
65e55bbd 717 AliError(" Max/Min theoretical Nc distributions are not defined");
718 return kFALSE;
719 }
a17b17dd 720 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectMCptMax);
65e55bbd 721 }
722 if (!areconsistent) {
723 AliInfo("Histograms required for fc correction are not consistent (bin width, bounds)");
724 return kFALSE;
725 }
726
727 return kTRUE;
728}
729
730//_________________________________________________________________________________________________________
731Bool_t AliHFPtSpectrum::CheckHistosConsistency(TH1 *h1, TH1 *h2){
732 //
733 // Check the histograms consistency (bins, limits)
734 //
735
736 if (!h1 || !h2) {
737 AliError("One or both histograms don't exist");
738 return kFALSE;
739 }
740
741 Double_t binwidth1 = h1->GetBinWidth(1);
742 Double_t binwidth2 = h2->GetBinWidth(1);
743 Double_t min1 = h1->GetBinCenter(1) - (binwidth1/2.) ;
744// Double_t max1 = h1->GetBinCenter(nbins1) + (binwidth1/2.) ;
745 Double_t min2 = h2->GetBinCenter(1) - (binwidth2/2.) ;
746// Double_t max2 = h2->GetBinCenter(nbins2) + (binwidth2/2.) ;
747
748 if (binwidth1!=binwidth2) {
749 AliInfo(" histograms with different bin width");
750 return kFALSE;
751 }
752 if (min1!=min2) {
753 AliInfo(" histograms with different minimum");
754 return kFALSE;
755 }
756// if (max1!=max2) {
757// AliInfo(" histograms with different maximum");
758// return kFALSE;
759// }
760
761 return kTRUE;
762}
763
764//_________________________________________________________________________________________________________
a17b17dd 765void AliHFPtSpectrum::CalculateFeedDownCorrectionFc(){
65e55bbd 766 //
767 // Compute fc factor and its uncertainties bin by bin
768 // fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
769 //
770
771 // define the variables
772 Int_t nbins = fhRECpt->GetNbinsX();
773 Double_t binwidth = fhRECpt->GetBinWidth(1);
bb427707 774 Double_t *limits = new Double_t[nbins+1];
de97013c 775 Double_t *binwidths = new Double_t[nbins];
bb427707 776 Double_t xlow=0.;
777 for (Int_t i=0; i<=nbins; i++) {
778 binwidth = fhRECpt->GetBinWidth(i);
779 xlow = fhRECpt->GetBinLowEdge(i);
780 limits[i-1] = xlow;
de97013c 781 binwidths[i] = binwidth;
bb427707 782 }
783 limits[nbins] = xlow + binwidth;
784
65e55bbd 785 Double_t correction=1.;
a17b17dd 786 Double_t correctionMax=1., correctionMin=1.;
787 Double_t theoryRatio=1.;
788 Double_t effRatio=1.;
65e55bbd 789
790 // declare the output histograms
bb427707 791 TH1D *hfc = new TH1D("hfc","fc correction factor",nbins,limits);
792 TH1D *hfcMax = new TH1D("hfcMax","max fc correction factor",nbins,limits);
793 TH1D *hfcMin = new TH1D("hfcMin","min fc correction factor",nbins,limits);
65e55bbd 794 // two local control histograms
bb427707 795 TH1D *hTheoryRatio = new TH1D("hTheoryRatio","Theoretical B-->D over c-->D (feed-down/direct) ratio",nbins,limits);
796 TH1D *hEffRatio = new TH1D("hEffRatio","Efficiency B-->D over c-->D (feed-down/direct) ratio",nbins,limits);
65e55bbd 797 // and the output TGraphAsymmErrors
798 if (fAsymUncertainties & !fgFc) fgFc = new TGraphAsymmErrors(nbins);
799
800 //
801 // Compute fc
802 //
803 for (Int_t ibin=0; ibin<=nbins; ibin++) {
804
805 // theory_ratio = (N_b/N_c)
a17b17dd 806 theoryRatio = (fhDirectMCpt->GetBinContent(ibin) && fhDirectMCpt->GetBinContent(ibin)!=0.) ? fhFeedDownMCpt->GetBinContent(ibin) / fhDirectMCpt->GetBinContent(ibin) : 1.0 ;
bb427707 807
65e55bbd 808 // eff_ratio = (eff_b/eff_c)
a17b17dd 809 effRatio = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ? fhFeedDownEffpt->GetBinContent(ibin) / fhDirectEffpt->GetBinContent(ibin) : 1.0 ;
bb427707 810
65e55bbd 811 // fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
a17b17dd 812 correction = (effRatio && theoryRatio) ? ( 1. / ( 1 + ( effRatio * theoryRatio ) ) ) : 1.0 ;
65e55bbd 813
814 // Calculate the uncertainty [ considering only the theoretical uncertainties on Nb & Nc for now !!! ]
815 // delta_fc = fc^2 * (Nb/Nc) * sqrt ( (delta_Nb/Nb)^2 + (delta_Nc/Nc)^2 )
a17b17dd 816 Double_t deltaNbMax = fhFeedDownMCptMax->GetBinContent(ibin) - fhFeedDownMCpt->GetBinContent(ibin) ;
817 Double_t deltaNbMin = fhFeedDownMCpt->GetBinContent(ibin) - fhFeedDownMCptMin->GetBinContent(ibin) ;
818 Double_t deltaNcMax = fhDirectMCptMax->GetBinContent(ibin) - fhDirectMCpt->GetBinContent(ibin) ;
819 Double_t deltaNcMin = fhDirectMCpt->GetBinContent(ibin) - fhDirectMCptMin->GetBinContent(ibin) ;
65e55bbd 820
821 // Protect against null denominator. If so, define uncertainty as null
822 if (fhFeedDownMCpt->GetBinContent(ibin) && fhFeedDownMCpt->GetBinContent(ibin)!=0. &&
823 fhDirectMCpt->GetBinContent(ibin) && fhDirectMCpt->GetBinContent(ibin)!=0. ) {
a17b17dd 824 correctionMax = correction*correction*theoryRatio *
65e55bbd 825 TMath::Sqrt(
a17b17dd 826 (deltaNbMax/fhFeedDownMCpt->GetBinContent(ibin))*(deltaNbMax/fhFeedDownMCpt->GetBinContent(ibin)) +
827 (deltaNcMax/fhDirectMCpt->GetBinContent(ibin))*(deltaNcMax/fhDirectMCpt->GetBinContent(ibin))
65e55bbd 828 );
a17b17dd 829 correctionMin = correction*correction*theoryRatio *
65e55bbd 830 TMath::Sqrt(
a17b17dd 831 (deltaNbMin/fhFeedDownMCpt->GetBinContent(ibin))*(deltaNbMin/fhFeedDownMCpt->GetBinContent(ibin)) +
832 (deltaNcMin/fhDirectMCpt->GetBinContent(ibin))*(deltaNcMin/fhDirectMCpt->GetBinContent(ibin))
65e55bbd 833 );
834 }
a17b17dd 835 else { correctionMax = 0.; correctionMin = 0.; }
65e55bbd 836
837
838 // Fill in the histograms
a17b17dd 839 hTheoryRatio->SetBinContent(ibin,theoryRatio);
840 hEffRatio->SetBinContent(ibin,effRatio);
65e55bbd 841 hfc->SetBinContent(ibin,correction);
a17b17dd 842 hfcMax->SetBinContent(ibin,correction+correctionMax);
843 hfcMin->SetBinContent(ibin,correction-correctionMin);
65e55bbd 844 if (fAsymUncertainties) {
845 Double_t x = fhDirectMCpt->GetBinCenter(ibin);
846 fgFc->SetPoint(ibin,x,correction); // i,x,y
de97013c 847 fgFc->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),correctionMin,correctionMax); // i,xl,xh,yl,yh
65e55bbd 848 }
849
850 }
851
852 fhFc = hfc;
a17b17dd 853 fhFcMax = hfcMax;
854 fhFcMin = hfcMin;
65e55bbd 855}
856
857//_________________________________________________________________________________________________________
a17b17dd 858void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumFc(){
65e55bbd 859 //
860 // Compute the feed-down corrected spectrum if feed-down correction is done via fc factor (bin by bin)
bb427707 861 // physics = reco * fc / bin-width
65e55bbd 862 //
863 // uncertainty: delta_physics = physics * sqrt ( (delta_reco/reco)^2 + (delta_fc/fc)^2 )
864 //
865 // ( Calculation done bin by bin )
866
867 if (!fhFc || !fhRECpt) {
868 AliError(" Reconstructed or fc distributions are not defined");
869 return;
870 }
871
872 Int_t nbins = fhRECpt->GetNbinsX();
a17b17dd 873 Double_t value = 0., valueDmax= 0., valueDmin= 0.;
65e55bbd 874 Double_t binwidth = fhRECpt->GetBinWidth(1);
bb427707 875 Double_t *limits = new Double_t[nbins+1];
de97013c 876 Double_t *binwidths = new Double_t[nbins];
bb427707 877 Double_t xlow=0.;
878 for (Int_t i=0; i<=nbins; i++) {
879 binwidth = fhRECpt->GetBinWidth(i);
880 xlow = fhRECpt->GetBinLowEdge(i);
881 limits[i-1] = xlow;
de97013c 882 binwidths[i] = binwidth;
bb427707 883 }
884 limits[nbins] = xlow + binwidth;
65e55bbd 885
886 // declare the output histograms
bb427707 887 TH1D *hYield = new TH1D("hYield","corrected yield (by fc)",nbins,limits);
888 TH1D *hYieldMax = new TH1D("hYieldMax","max corrected yield (by fc)",nbins,limits);
889 TH1D *hYieldMin = new TH1D("hYieldMin","min corrected yield (by fc)",nbins,limits);
65e55bbd 890 // and the output TGraphAsymmErrors
891 if (fAsymUncertainties & !fgYieldCorr) fgYieldCorr = new TGraphAsymmErrors(nbins);
892
893 //
894 // Do the calculation
895 //
896 for (Int_t ibin=0; ibin<=nbins; ibin++) {
897
898 // calculate the value
899 value = fhRECpt->GetBinContent(ibin) * fhFc->GetBinContent(ibin) ;
bb427707 900 value /= fhRECpt->GetBinWidth(ibin) ;
901 // cout << " bin="<< ibin << " reco="<<fhRECpt->GetBinContent(ibin) <<", fc="<< fhFc->GetBinContent(ibin) <<", value="<<value<<endl;
65e55bbd 902
903 // calculate the value uncertainty
904 // Protect against null denominator. If so, define uncertainty as null
905 if (fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) {
906
907 if (fAsymUncertainties) {
908
909 if (fhFc->GetBinContent(ibin) && fhFc->GetBinContent(ibin)!=0.) {
a17b17dd 910 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)) ) );
911 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)) ) );
65e55bbd 912 }
a17b17dd 913 else { valueDmax = 0.; valueDmin = 0.; }
65e55bbd 914
915 }
916 else { // Don't consider fc uncertainty in this case [ to be tested!!! ]
a17b17dd 917 valueDmax = value * (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) ;
918 valueDmin = value * (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) ;
65e55bbd 919 }
920
921 }
a17b17dd 922 else { valueDmax = 0.; valueDmin = 0.; }
65e55bbd 923
924 // fill in the histograms
925 hYield->SetBinContent(ibin,value);
a17b17dd 926 hYield->SetBinError(ibin,valueDmax);
927 hYieldMax->SetBinContent(ibin,value+valueDmax);
928 hYieldMin->SetBinContent(ibin,value-valueDmin);
65e55bbd 929 if (fAsymUncertainties) {
930 Double_t center = hYield->GetBinCenter(ibin);
931 fgYieldCorr->SetPoint(ibin,center,value); // i,x,y
de97013c 932 fgYieldCorr->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),valueDmin,valueDmax); // i,xl,xh,yl,yh
65e55bbd 933 }
934
935 }
936
937 fhYieldCorr = hYield;
a17b17dd 938 fhYieldCorrMax = hYieldMax;
939 fhYieldCorrMin = hYieldMin;
65e55bbd 940}
941
942//_________________________________________________________________________________________________________
a17b17dd 943void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumNb(Double_t deltaY, Double_t branchingRatioBintoFinalDecay) {
65e55bbd 944 //
945 // Compute the feed-down corrected spectrum if feed-down correction is done via Nb (bin by bin)
bb427707 946 // physics = [ reco - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) ] / bin-width
65e55bbd 947 //
948 // uncertainty: delta_physics = sqrt ( (delta_reco)^2 + (k*delta_lumi/lumi)^2 +
949 // (k*delta_eff_trig/eff_trig)^2 + (k*delta_Nb/Nb)^2 )
950 // where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th
951 //
952
953 Int_t nbins = fhRECpt->GetNbinsX();
954 Double_t binwidth = fhRECpt->GetBinWidth(1);
a17b17dd 955 Double_t value = 0., valueDmax = 0., valueDmin = 0., kfactor = 0.;
bb427707 956 Double_t *limits = new Double_t[nbins+1];
de97013c 957 Double_t *binwidths = new Double_t[nbins+1];
bb427707 958 Double_t xlow=0.;
959 for (Int_t i=0; i<=nbins; i++) {
960 binwidth = fhRECpt->GetBinWidth(i);
961 xlow = fhRECpt->GetBinLowEdge(i);
962 limits[i-1] = xlow;
de97013c 963 binwidths[i] = binwidth;
bb427707 964 }
965 limits[nbins] = xlow + binwidth;
65e55bbd 966
967 // declare the output histograms
bb427707 968 TH1D *hYield = new TH1D("hYield","corrected yield (by Nb)",nbins,limits);
969 TH1D *hYieldMax = new TH1D("hYieldMax","max corrected yield (by Nb)",nbins,limits);
970 TH1D *hYieldMin = new TH1D("hYieldMin","min corrected yield (by Nb)",nbins,limits);
65e55bbd 971 // and the output TGraphAsymmErrors
972 if (fAsymUncertainties & !fgYieldCorr) fgYieldCorr = new TGraphAsymmErrors(nbins);
973
974 //
975 // Do the calculation
976 //
977 for (Int_t ibin=0; ibin<=nbins; ibin++) {
978
979 // calculate the value
a17b17dd 980 value = fhRECpt->GetBinContent(ibin) - (deltaY*branchingRatioBintoFinalDecay*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) );
bb427707 981 value /= fhRECpt->GetBinWidth(ibin);
65e55bbd 982
a17b17dd 983 kfactor = deltaY*branchingRatioBintoFinalDecay*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) ;
65e55bbd 984
985 // calculate the value uncertainty
986 if (fAsymUncertainties) {
987 Double_t Nb = fhFeedDownMCpt->GetBinContent(ibin);
a17b17dd 988 Double_t NbDmax = fhFeedDownMCptMax->GetBinContent(ibin) - fhFeedDownMCpt->GetBinContent(ibin);
989 Double_t NbDmin = fhFeedDownMCpt->GetBinContent(ibin) - fhFeedDownMCptMin->GetBinContent(ibin);
990 valueDmax = TMath::Sqrt( ( fhRECpt->GetBinError(ibin)*fhRECpt->GetBinError(ibin) ) +
991 ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
992 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
993 ( (kfactor*NbDmax/Nb)*(kfactor*NbDmax/Nb) ) );
994 valueDmin = TMath::Sqrt( ( fhRECpt->GetBinError(ibin)*fhRECpt->GetBinError(ibin) ) +
65e55bbd 995 ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
996 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
a17b17dd 997 ( (kfactor*NbDmin/Nb)*(kfactor*NbDmin/Nb) ) );
65e55bbd 998 }
999 else{ // Don't consider Nb uncertainty in this case [ to be tested!!! ]
a17b17dd 1000 valueDmax = TMath::Sqrt( ( fhRECpt->GetBinError(ibin)*fhRECpt->GetBinError(ibin) ) +
65e55bbd 1001 ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
1002 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) );
a17b17dd 1003 valueDmin = valueDmax ;
65e55bbd 1004 }
1005
1006 // fill in histograms
1007 hYield->SetBinContent(ibin,value);
a17b17dd 1008 hYield->SetBinError(ibin,valueDmax);
1009 hYieldMax->SetBinContent(ibin,value+valueDmax);
1010 hYieldMin->SetBinContent(ibin,value-valueDmin);
65e55bbd 1011 if (fAsymUncertainties) {
1012 Double_t x = hYield->GetBinCenter(ibin);
1013 fgYieldCorr->SetPoint(ibin,x,value); // i,x,y
de97013c 1014 fgYieldCorr->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),valueDmin,valueDmax); // i,xl,xh,yl,yh
65e55bbd 1015 }
1016
1017 }
1018
1019 fhYieldCorr = hYield;
a17b17dd 1020 fhYieldCorrMax = hYieldMax;
1021 fhYieldCorrMin = hYieldMin;
65e55bbd 1022}
1023
1024
1025//_________________________________________________________________________________________________________
1026TH1 * AliHFPtSpectrum::ReweightHisto(TH1 *hToReweight, TH1 *hReference){
1027 //
1028 // Function to reweight histograms for testing purposes:
1029 // This function takes the histo hToReweight and reweights
1030 // it (its pt shape) with respect to hReference
1031 //
1032
1033 // check histograms consistency
1034 Bool_t areconsistent=kTRUE;
1035 areconsistent &= CheckHistosConsistency(hToReweight,hReference);
1036 if (!areconsistent) {
1037 AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
1038 return NULL;
1039 }
1040
1041 // define a new empty histogram
1042 TH1 *hReweighted = (TH1*)hToReweight->Clone("hReweighted");
1043 hReweighted->Reset();
1044 Double_t weight=1.0;
1045 Double_t yvalue=1.0;
a17b17dd 1046 Double_t integralRef = hReference->Integral();
1047 Double_t integralH = hToReweight->Integral();
65e55bbd 1048
1049 // now reweight the spectra
1050 //
1051 // the weight is the relative probability of the given pt bin in the reference histo
1052 // divided by its relative probability (to normalize it) on the histo to re-weight
1053 for (Int_t i=0; i<=hToReweight->GetNbinsX(); i++) {
a17b17dd 1054 weight = (hReference->GetBinContent(i)/integralRef) / (hToReweight->GetBinContent(i)/integralH) ;
65e55bbd 1055 yvalue = hToReweight->GetBinContent(i);
1056 hReweighted->SetBinContent(i,yvalue*weight);
1057 }
1058
1059 return (TH1*)hReweighted;
1060}
1061
1062//_________________________________________________________________________________________________________
1063TH1 * AliHFPtSpectrum::ReweightRecHisto(TH1 *hRecToReweight, TH1 *hMCToReweight, TH1 *hMCReference){
1064 //
1065 // Function to reweight histograms for testing purposes:
1066 // This function takes the histo hToReweight and reweights
1067 // it (its pt shape) with respect to hReference /hMCToReweight
1068 //
1069
1070 // check histograms consistency
1071 Bool_t areconsistent=kTRUE;
1072 areconsistent &= CheckHistosConsistency(hMCToReweight,hMCReference);
1073 areconsistent &= CheckHistosConsistency(hRecToReweight,hMCReference);
1074 if (!areconsistent) {
1075 AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
1076 return NULL;
1077 }
1078
1079 // define a new empty histogram
1080 TH1 *hReweighted = (TH1*)hMCToReweight->Clone("hReweighted");
1081 hReweighted->Reset();
1082 TH1 *hRecReweighted = (TH1*)hRecToReweight->Clone("hRecReweighted");
1083 hRecReweighted->Reset();
1084 Double_t weight=1.0;
1085 Double_t yvalue=1.0, yrecvalue=1.0;
a17b17dd 1086 Double_t integralRef = hMCReference->Integral();
1087 Double_t integralH = hMCToReweight->Integral();
65e55bbd 1088
1089 // now reweight the spectra
1090 //
1091 // the weight is the relative probability of the given pt bin
1092 // that should be applied in the MC histo to get the reference histo shape
1093 // Probabilities are properly normalized.
1094 for (Int_t i=0; i<=hMCToReweight->GetNbinsX(); i++) {
a17b17dd 1095 weight = (hMCReference->GetBinContent(i)/integralRef) / (hMCToReweight->GetBinContent(i)/integralH) ;
65e55bbd 1096 yvalue = hMCToReweight->GetBinContent(i);
1097 hReweighted->SetBinContent(i,yvalue*weight);
1098 yrecvalue = hRecToReweight->GetBinContent(i);
1099 hRecReweighted->SetBinContent(i,yrecvalue*weight);
1100 }
1101
1102 return (TH1*)hRecReweighted;
1103}
1104