Added booking of event selection (Renu)
[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>
bb427707 32#include "TNamed.h"
65e55bbd 33
34#include "AliLog.h"
35#include "AliHFPtSpectrum.h"
36
37ClassImp(AliHFPtSpectrum)
38
39//_________________________________________________________________________________________________________
40AliHFPtSpectrum::AliHFPtSpectrum(const char* name, const char* title, Int_t option):
41 TNamed(name,title),
86bdcd8c 42 fhDirectMCpt(NULL),
43 fhFeedDownMCpt(NULL),
44 fhDirectMCptMax(NULL),
45 fhDirectMCptMin(NULL),
46 fhFeedDownMCptMax(NULL),
47 fhFeedDownMCptMin(NULL),
48 fhDirectEffpt(NULL),
49 fhFeedDownEffpt(NULL),
50 fhRECpt(NULL),
51 fgRECSystematics(NULL),
65e55bbd 52 fLuminosity(),
53 fTrigEfficiency(),
86bdcd8c 54 fhFc(NULL),
55 fhFcMax(NULL),
56 fhFcMin(NULL),
57 fgFcExtreme(NULL),
58 fgFcConservative(NULL),
59 fhYieldCorr(NULL),
60 fhYieldCorrMax(NULL),
61 fhYieldCorrMin(NULL),
62 fgYieldCorr(NULL),
63 fgYieldCorrExtreme(NULL),
64 fgYieldCorrConservative(NULL),
65 fhSigmaCorr(NULL),
66 fhSigmaCorrMax(NULL),
67 fhSigmaCorrMin(NULL),
68 fgSigmaCorr(NULL),
69 fgSigmaCorrExtreme(NULL),
70 fgSigmaCorrConservative(NULL),
65e55bbd 71 fFeedDownOption(option),
72 fAsymUncertainties(kTRUE)
73{
74 //
75 // Default constructor
76 //
77
78 fLuminosity[0]=1.; fLuminosity[1]=0.;
79 fTrigEfficiency[0]=1.; fTrigEfficiency[1]=0.;
80
81}
82
83//_________________________________________________________________________________________________________
84AliHFPtSpectrum::AliHFPtSpectrum(const AliHFPtSpectrum &rhs):
85 TNamed(rhs),
86 fhDirectMCpt(rhs.fhDirectMCpt),
87 fhFeedDownMCpt(rhs.fhFeedDownMCpt),
a17b17dd 88 fhDirectMCptMax(rhs.fhDirectMCptMax),
89 fhDirectMCptMin(rhs.fhDirectMCptMin),
90 fhFeedDownMCptMax(rhs.fhFeedDownMCptMax),
91 fhFeedDownMCptMin(rhs.fhFeedDownMCptMin),
65e55bbd 92 fhDirectEffpt(rhs.fhDirectEffpt),
93 fhFeedDownEffpt(rhs.fhFeedDownEffpt),
94 fhRECpt(rhs.fhRECpt),
86bdcd8c 95 fgRECSystematics(rhs.fgRECSystematics),
65e55bbd 96 fLuminosity(),
97 fTrigEfficiency(),
98 fhFc(rhs.fhFc),
a17b17dd 99 fhFcMax(rhs.fhFcMax),
100 fhFcMin(rhs.fhFcMin),
86bdcd8c 101 fgFcExtreme(rhs.fgFcExtreme),
102 fgFcConservative(rhs.fgFcConservative),
65e55bbd 103 fhYieldCorr(rhs.fhYieldCorr),
a17b17dd 104 fhYieldCorrMax(rhs.fhYieldCorrMax),
105 fhYieldCorrMin(rhs.fhYieldCorrMin),
65e55bbd 106 fgYieldCorr(rhs.fgYieldCorr),
86bdcd8c 107 fgYieldCorrExtreme(rhs.fgYieldCorrExtreme),
108 fgYieldCorrConservative(rhs.fgYieldCorrConservative),
65e55bbd 109 fhSigmaCorr(rhs.fhSigmaCorr),
a17b17dd 110 fhSigmaCorrMax(rhs.fhSigmaCorrMax),
111 fhSigmaCorrMin(rhs.fhSigmaCorrMin),
65e55bbd 112 fgSigmaCorr(rhs.fgSigmaCorr),
86bdcd8c 113 fgSigmaCorrExtreme(rhs.fgSigmaCorrExtreme),
114 fgSigmaCorrConservative(rhs.fgSigmaCorrConservative),
65e55bbd 115 fFeedDownOption(rhs.fFeedDownOption),
116 fAsymUncertainties(rhs.fAsymUncertainties)
117{
118 //
119 // Copy constructor
120 //
121
122 for(Int_t i=0; i<2; i++){
123 fLuminosity[i] = rhs.fLuminosity[i];
124 fTrigEfficiency[i] = rhs.fTrigEfficiency[i];
125 }
126
127}
128
129//_________________________________________________________________________________________________________
130AliHFPtSpectrum &AliHFPtSpectrum::operator=(const AliHFPtSpectrum &source){
131 //
132 // Assignment operator
133 //
134
135 if (&source == this) return *this;
136
137 fhDirectMCpt = source.fhDirectMCpt;
138 fhFeedDownMCpt = source.fhFeedDownMCpt;
a17b17dd 139 fhDirectMCptMax = source.fhDirectMCptMax;
140 fhDirectMCptMin = source.fhDirectMCptMin;
141 fhFeedDownMCptMax = source.fhFeedDownMCptMax;
142 fhFeedDownMCptMin = source.fhFeedDownMCptMin;
65e55bbd 143 fhDirectEffpt = source.fhDirectEffpt;
144 fhFeedDownEffpt = source.fhFeedDownEffpt;
145 fhRECpt = source.fhRECpt;
86bdcd8c 146 fgRECSystematics = source.fgRECSystematics;
65e55bbd 147 fhFc = source.fhFc;
a17b17dd 148 fhFcMax = source.fhFcMax;
149 fhFcMin = source.fhFcMin;
86bdcd8c 150 fgFcExtreme = source.fgFcExtreme;
151 fgFcConservative = source.fgFcConservative;
65e55bbd 152 fhYieldCorr = source.fhYieldCorr;
a17b17dd 153 fhYieldCorrMax = source.fhYieldCorrMax;
154 fhYieldCorrMin = source.fhYieldCorrMin;
65e55bbd 155 fgYieldCorr = source.fgYieldCorr;
86bdcd8c 156 fgYieldCorrExtreme = source.fgYieldCorrExtreme;
157 fgYieldCorrConservative = source.fgYieldCorrConservative;
65e55bbd 158 fhSigmaCorr = source.fhSigmaCorr;
a17b17dd 159 fhSigmaCorrMax = source.fhSigmaCorrMax;
160 fhSigmaCorrMin = source.fhSigmaCorrMin;
65e55bbd 161 fgSigmaCorr = source.fgSigmaCorr;
86bdcd8c 162 fgSigmaCorrExtreme = source.fgSigmaCorrExtreme;
163 fgSigmaCorrConservative = source.fgSigmaCorrConservative;
65e55bbd 164 fFeedDownOption = source.fFeedDownOption;
165 fAsymUncertainties = source.fAsymUncertainties;
166
167 for(Int_t i=0; i<2; i++){
168 fLuminosity[i] = source.fLuminosity[i];
169 fTrigEfficiency[i] = source.fTrigEfficiency[i];
170 }
171
172 return *this;
173}
174
175//_________________________________________________________________________________________________________
176AliHFPtSpectrum::~AliHFPtSpectrum(){
177 //
178 // Destructor
179 //
86bdcd8c 180 if (fhDirectMCpt) delete fhDirectMCpt;
181 if (fhFeedDownMCpt) delete fhFeedDownMCpt;
182 if (fhDirectMCptMax) delete fhDirectMCptMax;
183 if (fhDirectMCptMin) delete fhDirectMCptMin;
184 if (fhFeedDownMCptMax) delete fhFeedDownMCptMax;
185 if (fhFeedDownMCptMin) delete fhFeedDownMCptMin;
186 if (fhDirectEffpt) delete fhDirectEffpt;
187 if (fhFeedDownEffpt) delete fhFeedDownEffpt;
188 if (fhRECpt) delete fhRECpt;
189 if (fgRECSystematics) delete fgRECSystematics;
190 if (fhFc) delete fhFc;
191 if (fhFcMax) delete fhFcMax;
192 if (fhFcMin) delete fhFcMin;
193 if (fgFcExtreme) delete fgFcExtreme;
194 if (fgFcConservative) delete fgFcConservative;
195 if (fhYieldCorr) delete fhYieldCorr;
196 if (fhYieldCorrMax) delete fhYieldCorrMax;
197 if (fhYieldCorrMin) delete fhYieldCorrMin;
198 if (fgYieldCorr) delete fgYieldCorr;
199 if (fgYieldCorrExtreme) delete fgYieldCorrExtreme;
200 if (fgYieldCorrConservative) delete fgYieldCorrConservative;
201 if (fhSigmaCorr) delete fhSigmaCorr;
202 if (fhSigmaCorrMax) delete fhSigmaCorrMax;
203 if (fhSigmaCorrMin) delete fhSigmaCorrMin;
204 if (fgSigmaCorr) delete fgSigmaCorr;
205 if (fgSigmaCorrExtreme) delete fgSigmaCorrExtreme;
206 if (fgSigmaCorrConservative) delete fgSigmaCorrConservative;
65e55bbd 207}
208
209
210//_________________________________________________________________________________________________________
86bdcd8c 211TH1D * AliHFPtSpectrum::RebinTheoreticalSpectra(TH1D *hTheory, const char *name) {
5f3c1b97 212 //
213 // Function to rebin the theoretical spectrum
214 // with respect to the real-data reconstructed spectrum binning
215 //
216
217 if (!hTheory || !fhRECpt) {
218 AliError("Feed-down or reconstructed spectra don't exist");
219 return NULL;
220 }
221
222 //
223 // Get the reconstructed spectra bins & limits
224 Int_t nbins = fhRECpt->GetNbinsX();
225 Int_t nbinsMC = hTheory->GetNbinsX();
226 Double_t *limits = new Double_t[nbins+1];
227 Double_t xlow=0., binwidth=0.;
86bdcd8c 228 for (Int_t i=1; i<=nbins; i++) {
5f3c1b97 229 binwidth = fhRECpt->GetBinWidth(i);
230 xlow = fhRECpt->GetBinLowEdge(i);
231 limits[i-1] = xlow;
232 }
233 limits[nbins] = xlow + binwidth;
234
235 //
236 // Define a new histogram with the real-data reconstructed spectrum binning
86bdcd8c 237 TH1D * hTheoryRebin = new TH1D(name," theoretical rebinned prediction",nbins,limits);
5f3c1b97 238
239 Double_t sum[nbins], items[nbins];
240 for (Int_t ibin=0; ibin<=nbinsMC; ibin++){
241
242 for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++){
86bdcd8c 243 if (hTheory->GetBinCenter(ibin)>limits[ibinrec] &&
5f3c1b97 244 hTheory->GetBinCenter(ibin)<limits[ibinrec+1]){
245 sum[ibinrec]+=hTheory->GetBinContent(ibin);
246 items[ibinrec]+=1.;
247 }
248 }
249
250 }
251
252 // set the theoretical rebinned spectra to ( sum-bins / n-bins ) per new bin
253 for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++) {
254 hTheoryRebin->SetBinContent(ibinrec+1,sum[ibinrec]/items[ibinrec]);
255 }
256
86bdcd8c 257 return (TH1D*)hTheoryRebin;
5f3c1b97 258}
259
260//_________________________________________________________________________________________________________
86bdcd8c 261void AliHFPtSpectrum::SetMCptSpectra(TH1D *hDirect, TH1D *hFeedDown){
65e55bbd 262 //
263 // Set the MonteCarlo or Theoretical spectra
264 // both for direct and feed-down contributions
265 //
266
bb427707 267 if (!hDirect || !hFeedDown || !fhRECpt) {
268 AliError("One or both (direct, feed-down) spectra or the reconstructed spectra don't exist");
65e55bbd 269 return;
270 }
271
272 Bool_t areconsistent = kTRUE;
273 areconsistent = CheckHistosConsistency(hDirect,hFeedDown);
274 if (!areconsistent) {
275 AliInfo("Histograms are not consistent (bin width, bounds)");
276 return;
277 }
278
bb427707 279 //
bb427707 280 // If the predictions shape do not have the same
281 // bin width (check via number of bins) as the reconstructed spectra, change them
248ae766 282 if ( TMath::Abs(hDirect->GetBinWidth(1) - fhRECpt->GetBinWidth(1))< 0.0001 ) {
bb427707 283
86bdcd8c 284 fhDirectMCpt = (TH1D*)hDirect->Clone();
285 fhDirectMCpt->SetNameTitle("fhDirectMCpt"," direct theoretical prediction");
286 fhFeedDownMCpt = (TH1D*)hFeedDown->Clone();
287 fhFeedDownMCpt->SetNameTitle("fhFeedDownMCpt"," feed-down theoretical prediction");
bb427707 288
86bdcd8c 289 }
bb427707 290 else {
291
86bdcd8c 292 fhDirectMCpt = RebinTheoreticalSpectra(hDirect,"fhDirectMCpt");
5f3c1b97 293 fhDirectMCpt->SetNameTitle("fhDirectMCpt"," direct theoretical prediction");
86bdcd8c 294 fhFeedDownMCpt = RebinTheoreticalSpectra(hFeedDown,"fhFeedDownMCpt");
5f3c1b97 295 fhFeedDownMCpt->SetNameTitle("fhFeedDownMCpt"," feed-down theoretical prediction");
bb427707 296
bb427707 297 }
298
65e55bbd 299}
300
301//_________________________________________________________________________________________________________
86bdcd8c 302void AliHFPtSpectrum::SetFeedDownMCptSpectra(TH1D *hFeedDown){
65e55bbd 303 //
304 // Set the MonteCarlo or Theoretical spectra
305 // for feed-down contribution
306 //
307
bb427707 308 if (!hFeedDown || !fhRECpt) {
309 AliError("Feed-down or reconstructed spectra don't exist");
65e55bbd 310 return;
311 }
bb427707 312
313 //
bb427707 314 // If the predictions shape do not have the same
315 // bin width (check via number of bins) as the reconstructed spectra, change them
248ae766 316 if ( TMath::Abs(hFeedDown->GetBinWidth(1) - fhRECpt->GetBinWidth(1))< 0.0001 ) {
bb427707 317
86bdcd8c 318 fhFeedDownMCpt = (TH1D*)hFeedDown->Clone();
319 fhFeedDownMCpt->SetNameTitle("fhFeedDownMCpt"," feed-down theoretical prediction");
bb427707 320
321 }
322 else {
323
86bdcd8c 324 fhFeedDownMCpt = RebinTheoreticalSpectra(hFeedDown,"fhFeedDownMCpt");
5f3c1b97 325 fhFeedDownMCpt->SetNameTitle("fhFeedDownMCpt"," feed-down theoretical prediction");
bb427707 326
bb427707 327 }
328
65e55bbd 329}
330
331//_________________________________________________________________________________________________________
86bdcd8c 332void AliHFPtSpectrum::SetMCptDistributionsBounds(TH1D *hDirectMax, TH1D *hDirectMin, TH1D *hFeedDownMax, TH1D *hFeedDownMin){
65e55bbd 333 //
334 // Set the maximum and minimum MonteCarlo or Theoretical spectra
335 // both for direct and feed-down contributions
336 // used in case uncertainties are asymmetric and ca not be on the "basic histograms"
337 //
338
bb427707 339 if (!hDirectMax || !hDirectMin || !hFeedDownMax|| !hFeedDownMin || !fhRECpt) {
340 AliError("One or all of the max/min direct/feed-down or the reconstructed spectra don't exist");
65e55bbd 341 return;
342 }
343
344 Bool_t areconsistent = kTRUE;
345 areconsistent &= CheckHistosConsistency(hDirectMax,hDirectMin);
346 areconsistent &= CheckHistosConsistency(hFeedDownMax,hFeedDownMin);
347 areconsistent &= CheckHistosConsistency(hDirectMax,hFeedDownMax);
348 if (!areconsistent) {
349 AliInfo("Histograms are not consistent (bin width, bounds)");
350 return;
351 }
352
bb427707 353
354 //
355 // If the predictions shape do not have the same
356 // bin width (check via number of bins) as the reconstructed spectra, change them
248ae766 357 if ( TMath::Abs(hFeedDownMax->GetBinWidth(1) - fhRECpt->GetBinWidth(1))< 0.0001 ) {
bb427707 358
86bdcd8c 359 fhDirectMCptMax = (TH1D*)hDirectMax->Clone();
360 fhDirectMCptMin = (TH1D*)hDirectMin->Clone();
361 fhFeedDownMCptMax = (TH1D*)hFeedDownMax->Clone();
362 fhFeedDownMCptMin = (TH1D*)hFeedDownMin->Clone();
363 fhDirectMCptMax->SetNameTitle("fhDirectMCptMax"," maximum direct theoretical prediction");
364 fhDirectMCptMin->SetNameTitle("fhDirectMCptMin"," minimum direct theoretical prediction");
365 fhFeedDownMCptMax->SetNameTitle("fhFeedDownMCptMax"," maximum feed-down theoretical prediction");
366 fhFeedDownMCptMin->SetNameTitle("fhFeedDownMCptMin"," minimum feed-down theoretical prediction");
bb427707 367
368 }
369 else {
370
86bdcd8c 371 fhDirectMCptMax = RebinTheoreticalSpectra(hDirectMax,"fhDirectMCptMax");
5f3c1b97 372 fhDirectMCptMax->SetNameTitle("fhDirectMCptMax"," maximum direct theoretical prediction");
86bdcd8c 373 fhDirectMCptMin = RebinTheoreticalSpectra(hDirectMin,"fhDirectMCptMin");
5f3c1b97 374 fhDirectMCptMin->SetNameTitle("fhDirectMCptMin"," minimum direct theoretical prediction");
86bdcd8c 375 fhFeedDownMCptMax = RebinTheoreticalSpectra(hFeedDownMax,"fhFeedDownMCptMax");
5f3c1b97 376 fhFeedDownMCptMax->SetNameTitle("fhFeedDownMCptMax"," maximum feed-down theoretical prediction");
86bdcd8c 377 fhFeedDownMCptMin = RebinTheoreticalSpectra(hFeedDownMin,"fhFeedDownMCptMin");
5f3c1b97 378 fhFeedDownMCptMin->SetNameTitle("fhFeedDownMCptMin"," minimum feed-down theoretical prediction");
bb427707 379
bb427707 380 }
381
65e55bbd 382}
383
384//_________________________________________________________________________________________________________
86bdcd8c 385void AliHFPtSpectrum::SetFeedDownMCptDistributionsBounds(TH1D *hFeedDownMax, TH1D *hFeedDownMin){
65e55bbd 386 //
387 // Set the maximum and minimum MonteCarlo or Theoretical spectra
388 // for feed-down contributions
389 // used in case uncertainties are asymmetric and can not be on the "basic histogram"
390 //
391
bb427707 392 if (!hFeedDownMax || !hFeedDownMin || !fhRECpt) {
65e55bbd 393 AliError("One or all of the max/min direct/feed-down spectra don't exist");
394 return;
395 }
396
397 Bool_t areconsistent = kTRUE;
398 areconsistent &= CheckHistosConsistency(hFeedDownMax,hFeedDownMin);
399 if (!areconsistent) {
400 AliInfo("Histograms are not consistent (bin width, bounds)");
401 return;
402 }
403
bb427707 404
405 //
406 // If the predictions shape do not have the same
407 // bin width (check via number of bins) as the reconstructed spectra, change them
248ae766 408 if ( TMath::Abs(hFeedDownMax->GetBinWidth(1) - fhRECpt->GetBinWidth(1))< 0.0001 ) {
bb427707 409
86bdcd8c 410 fhFeedDownMCptMax = (TH1D*)hFeedDownMax->Clone();
411 fhFeedDownMCptMin = (TH1D*)hFeedDownMin->Clone();
412 fhFeedDownMCptMax->SetNameTitle("fhFeedDownMCptMax"," maximum feed-down theoretical prediction");
413 fhFeedDownMCptMin->SetNameTitle("fhFeedDownMCptMin"," minimum feed-down theoretical prediction");
bb427707 414
415 }
416 else {
417
86bdcd8c 418 fhFeedDownMCptMax = RebinTheoreticalSpectra(hFeedDownMax,"fhFeedDownMCptMax");
5f3c1b97 419 fhFeedDownMCptMax->SetNameTitle("fhFeedDownMCptMax"," maximum feed-down theoretical prediction");
86bdcd8c 420 fhFeedDownMCptMin = RebinTheoreticalSpectra(hFeedDownMin,"fhFeedDownMCptMin");
5f3c1b97 421 fhFeedDownMCptMin->SetNameTitle("fhFeedDownMCptMin"," minimum feed-down theoretical prediction");
bb427707 422
bb427707 423 }
424
65e55bbd 425}
426
427//_________________________________________________________________________________________________________
86bdcd8c 428void AliHFPtSpectrum::SetDirectAccEffCorrection(TH1D *hDirectEff){
65e55bbd 429 //
430 // Set the Acceptance and Efficiency corrections
431 // for the direct contribution
432 //
433
434 if (!hDirectEff) {
435 AliError("The direct acceptance and efficiency corrections doesn't exist");
436 return;
437 }
438
86bdcd8c 439 fhDirectEffpt = (TH1D*)hDirectEff->Clone();
440 fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance x efficiency correction");
65e55bbd 441}
442
443//_________________________________________________________________________________________________________
86bdcd8c 444void AliHFPtSpectrum::SetAccEffCorrection(TH1D *hDirectEff, TH1D *hFeedDownEff){
65e55bbd 445 //
446 // Set the Acceptance and Efficiency corrections
447 // both for direct and feed-down contributions
448 //
449
450 if (!hDirectEff || !hFeedDownEff) {
451 AliError("One or both (direct, feed-down) acceptance and efficiency corrections don't exist");
452 return;
453 }
454
455 Bool_t areconsistent=kTRUE;
456 areconsistent = CheckHistosConsistency(hDirectEff,hFeedDownEff);
457 if (!areconsistent) {
458 AliInfo("Histograms are not consistent (bin width, bounds)");
459 return;
460 }
461
86bdcd8c 462 fhDirectEffpt = (TH1D*)hDirectEff->Clone();
463 fhFeedDownEffpt = (TH1D*)hFeedDownEff->Clone();
464 fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance x efficiency correction");
465 fhFeedDownEffpt->SetNameTitle("fhFeedDownEffpt"," feed-down acceptance x efficiency correction");
65e55bbd 466}
467
468//_________________________________________________________________________________________________________
86bdcd8c 469void AliHFPtSpectrum::SetReconstructedSpectrum(TH1D *hRec) {
65e55bbd 470 //
471 // Set the reconstructed spectrum
472 //
473
474 if (!hRec) {
475 AliError("The reconstructed spectrum doesn't exist");
476 return;
477 }
478
86bdcd8c 479 fhRECpt = (TH1D*)hRec->Clone();
480 fhRECpt->SetNameTitle("fhRECpt"," reconstructed spectrum");
481}
482
483//_________________________________________________________________________________________________________
484void AliHFPtSpectrum::SetReconstructedSpectrumSystematics(TGraphAsymmErrors *gRec) {
485 //
486 // Set the reconstructed spectrum systematic uncertainties
487 //
488
489 // Check the compatibility with the reconstructed spectrum
490 Double_t gbinwidth = fgRECSystematics->GetErrorXlow(1) + fgRECSystematics->GetErrorXhigh(1) ;
491 Double_t hbinwidth = fhRECpt->GetBinWidth(1);
492 Double_t gxbincenter=0., gybincenter=0.;
493 fgRECSystematics->GetPoint(1,gxbincenter,gybincenter);
494 Double_t hbincenter = fhRECpt->GetBinCenter(1);
248ae766 495 if ( TMath::Abs(gbinwidth - hbinwidth)>0.0001 || TMath::Abs(gxbincenter - hbincenter)>0.0001 ) {
86bdcd8c 496 AliError(" The reconstructed spectrum and its systematics don't seem compatible");
497 return;
498 }
499
500 fgRECSystematics = gRec;
65e55bbd 501}
502
503//_________________________________________________________________________________________________________
a17b17dd 504void AliHFPtSpectrum::ComputeHFPtSpectrum(Double_t deltaY, Double_t branchingRatioC, Double_t branchingRatioBintoFinalDecay) {
65e55bbd 505 //
506 // Main function to compute the corrected cross-section:
a17b17dd 507 // variables : analysed delta_y, BR for the final correction,
508 // BR b --> D --> decay (relative to the input theoretical prediction)
65e55bbd 509 //
510 // Sigma = ( 1. / (lumi * delta_y * BR_c * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
511 //
86bdcd8c 512 // Uncertainties: (stat) delta_sigma = sigma * sqrt ( (delta_spectra/spectra)^2 )
513 // (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 )
514 // (feed-down syst) delta_sigma = sigma * sqrt ( (delta_spectra_fd/spectra_fd)^2 )
65e55bbd 515
516 //
517 // First: Initialization
518 //
519 Bool_t areHistosOk = Initialize();
520 if (!areHistosOk) {
521 AliInfo(" Histos not properly initialized. Check : inconsistent binning ? missing histos ?");
522 return;
523 }
524
525 //
526 // Second: Correct for feed-down
527 //
528 if (fFeedDownOption==1) {
529 // Compute the feed-down correction via fc-method
a17b17dd 530 CalculateFeedDownCorrectionFc();
65e55bbd 531 // Correct the yield for feed-down correction via fc-method
a17b17dd 532 CalculateFeedDownCorrectedSpectrumFc();
65e55bbd 533 }
534 else if (fFeedDownOption==2) {
535 // Correct the yield for feed-down correction via Nb-method
a17b17dd 536 CalculateFeedDownCorrectedSpectrumNb(deltaY,branchingRatioBintoFinalDecay);
65e55bbd 537 }
538 else if (fFeedDownOption==0) {
539 // If there is no need for feed-down correction,
540 // the "corrected" yield is equal to the raw yield
86bdcd8c 541 fhYieldCorr = (TH1D*)fhRECpt->Clone();
65e55bbd 542 fhYieldCorr->SetNameTitle("fhYieldCorr","un-corrected yield");
86bdcd8c 543 fhYieldCorrMax = (TH1D*)fhRECpt->Clone();
544 fhYieldCorrMin = (TH1D*)fhRECpt->Clone();
a17b17dd 545 fhYieldCorrMax->SetNameTitle("fhYieldCorrMax","un-corrected yield");
546 fhYieldCorrMin->SetNameTitle("fhYieldCorrMin","un-corrected yield");
65e55bbd 547 fAsymUncertainties=kFALSE;
548 }
549 else {
550 AliInfo(" Are you sure the feed-down correction option is right ?");
551 }
552
553 // Print out information
a17b17dd 554 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 555
556 //
557 // Finally: Correct from yields to cross-section
558 //
559 Int_t nbins = fhRECpt->GetNbinsX();
560 Double_t binwidth = fhRECpt->GetBinWidth(1);
bb427707 561 Double_t *limits = new Double_t[nbins+1];
de97013c 562 Double_t *binwidths = new Double_t[nbins+1];
bb427707 563 Double_t xlow=0.;
86bdcd8c 564 for (Int_t i=1; i<=nbins; i++) {
bb427707 565 binwidth = fhRECpt->GetBinWidth(i);
566 xlow = fhRECpt->GetBinLowEdge(i);
567 limits[i-1] = xlow;
de97013c 568 binwidths[i] = binwidth;
bb427707 569 }
570 limits[nbins] = xlow + binwidth;
571
65e55bbd 572
573 // declare the output histograms
86bdcd8c 574 if (!fhSigmaCorr) fhSigmaCorr = new TH1D("fhSigmaCorr","corrected sigma",nbins,limits);
575 if (!fhSigmaCorrMax) fhSigmaCorrMax = new TH1D("fhSigmaCorrMax","max corrected sigma",nbins,limits);
576 if (!fhSigmaCorrMin) fhSigmaCorrMin = new TH1D("fhSigmaCorrMin","min corrected sigma",nbins,limits);
65e55bbd 577 // and the output TGraphAsymmErrors
86bdcd8c 578 if (fAsymUncertainties & !fgSigmaCorr) fgSigmaCorr = new TGraphAsymmErrors(nbins+1);
579 if (fAsymUncertainties & !fgSigmaCorrExtreme) fgSigmaCorrExtreme = new TGraphAsymmErrors(nbins+1);
580 if (fAsymUncertainties & !fgSigmaCorrConservative) fgSigmaCorrConservative = new TGraphAsymmErrors(nbins+1);
65e55bbd 581
582 // protect against null denominator
248ae766 583 if (deltaY<0.01 || fLuminosity[0]<0.0000001 || fTrigEfficiency[0]<0.0000001 || branchingRatioC<0.000000001) {
65e55bbd 584 AliError(" Hey you ! Why luminosity or trigger-efficiency or the c-BR or delta_y are set to zero ?! ");
585 return ;
586 }
587
86bdcd8c 588 Double_t value=0, errValue=0, errvalueMax=0., errvalueMin=0.;
589 Double_t errvalueExtremeMax=0., errvalueExtremeMin=0.;
590 Double_t errvalueConservativeMax=0., errvalueConservativeMin=0.;
65e55bbd 591 for(Int_t ibin=0; ibin<=nbins; ibin++){
592
86bdcd8c 593 // Sigma calculation
65e55bbd 594 // Sigma = ( 1. / (lumi * delta_y * BR_c * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
595 value = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ?
a17b17dd 596 ( fhYieldCorr->GetBinContent(ibin) / ( deltaY * branchingRatioC * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
65e55bbd 597 : 0. ;
86bdcd8c 598
599 // Sigma statistical uncertainty:
600 // delta_sigma = sigma * sqrt ( (delta_spectra/spectra)^2 )
601 errValue = (value!=0.) ? value * (fhYieldCorr->GetBinError(ibin)/fhYieldCorr->GetBinContent(ibin)) : 0. ;
602
603 //
604 // Sigma systematic uncertainties
605 //
65e55bbd 606 if (fAsymUncertainties) {
86bdcd8c 607
608 // (syst but feed-down) delta_sigma = sigma * sqrt ( (delta_spectra_syst/spectra)^2 +
609 // (delta_lumi/lumi)^2 + (delta_eff_trig/eff_trig)^2 + (delta_eff/eff)^2 )
610 errvalueMax = value * TMath::Sqrt( (fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin))*(fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin)) +
611 (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
612 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) +
613 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) );
614 errvalueMin = value * TMath::Sqrt( (fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin))*(fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin)) +
615 (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
616 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) +
617 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) );
618
619 // Uncertainties from feed-down
620 // (feed-down syst) delta_sigma = sigma * sqrt ( (delta_spectra_fd/spectra_fd)^2 )
621 // extreme case
622 errvalueExtremeMax = value * (fgYieldCorrExtreme->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin));
623 errvalueExtremeMin = value * (fgYieldCorrExtreme->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin));
624 //
625 // conservative case
626 errvalueConservativeMax = value * (fgYieldCorrConservative->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin));
627 errvalueConservativeMin = value * (fgYieldCorrConservative->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin));
628
65e55bbd 629 }
630 else {
631 // protect against null denominator
86bdcd8c 632 errvalueMax = (value!=0.) ?
633 value * TMath::Sqrt( (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
634 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) +
635 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) )
65e55bbd 636 : 0. ;
86bdcd8c 637 errvalueMin = errvalueMax;
65e55bbd 638 }
639
640 // Fill the histograms
86bdcd8c 641 fhSigmaCorr->SetBinContent(ibin,value);
642 fhSigmaCorr->SetBinError(ibin,errValue);
65e55bbd 643 // Fill the TGraphAsymmErrors
644 if (fAsymUncertainties) {
645 Double_t x = fhYieldCorr->GetBinCenter(ibin);
646 fgSigmaCorr->SetPoint(ibin,x,value); // i,x,y
86bdcd8c 647 fgSigmaCorr->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
648 fhSigmaCorrMax->SetBinContent(ibin,value+errvalueMax);
649 fhSigmaCorrMin->SetBinContent(ibin,value-errvalueMin);
650 fgSigmaCorrExtreme->SetPoint(ibin,x,value); // i,x,y
651 fgSigmaCorrExtreme->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
652 fgSigmaCorrConservative->SetPoint(ibin,x,value); // i,x,y
653 fgSigmaCorrConservative->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),errvalueConservativeMin,errvalueConservativeMax); // i,xl,xh,yl,yh
654
65e55bbd 655 }
656
657 }
658
65e55bbd 659}
660
661//_________________________________________________________________________________________________________
86bdcd8c 662TH1D * AliHFPtSpectrum::EstimateEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco, const char *name) {
bb427707 663 //
5f3c1b97 664 // Function that computes the acceptance and efficiency correction
bb427707 665 // based on the simulated and reconstructed spectra
666 // and using the reconstructed spectra bin width
667 //
668 // eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
669 //
670
671 if(!fhRECpt){
672 AliInfo("Hey, the reconstructed histogram was not set yet !");
673 return NULL;
674 }
675
676 Int_t nbins = fhRECpt->GetNbinsX();
677 Double_t *limits = new Double_t[nbins+1];
678 Double_t xlow=0.,binwidth=0.;
86bdcd8c 679 for (Int_t i=1; i<=nbins; i++) {
bb427707 680 binwidth = fhRECpt->GetBinWidth(i);
681 xlow = fhRECpt->GetBinLowEdge(i);
682 limits[i-1] = xlow;
683 }
684 limits[nbins] = xlow + binwidth;
685
86bdcd8c 686 TH1D * hEfficiency = new TH1D(name," acceptance #times efficiency",nbins,limits);
bb427707 687
688 Double_t sumSimu[nbins], sumReco[nbins];
689 for (Int_t ibin=0; ibin<=hSimu->GetNbinsX(); ibin++){
690
691 for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++){
692 if ( hSimu->GetBinCenter(ibin)>limits[ibinrec] &&
693 hSimu->GetBinCenter(ibin)<limits[ibinrec+1] ) {
694 sumSimu[ibinrec]+=hSimu->GetBinContent(ibin);
695 }
696 if ( hReco->GetBinCenter(ibin)>limits[ibinrec] &&
697 hReco->GetBinCenter(ibin)<limits[ibinrec+1] ) {
698 sumReco[ibinrec]+=hReco->GetBinContent(ibin);
699 }
700 }
701
702 }
703
86bdcd8c 704
bb427707 705 // the efficiency is computed as reco/sim (in each bin)
706 // its uncertainty is err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
707 Double_t eff=0., erreff=0.;
86bdcd8c 708 for (Int_t ibinrec=0; ibinrec<=nbins; ibinrec++) {
709 if (sumSimu[ibinrec]!= 0. && sumReco[ibinrec]!=0.) {
710 eff = sumReco[ibinrec] / sumSimu[ibinrec] ;
711 erreff = TMath::Sqrt( eff * (1.0 - eff) ) / TMath::Sqrt( sumSimu[ibinrec] );
712 }
713 else { eff=0.0; erreff=0.; }
5f3c1b97 714 hEfficiency->SetBinContent(ibinrec+1,eff);
715 hEfficiency->SetBinError(ibinrec+1,erreff);
bb427707 716 }
717
86bdcd8c 718 return (TH1D*)hEfficiency;
bb427707 719}
720
721//_________________________________________________________________________________________________________
86bdcd8c 722void AliHFPtSpectrum::EstimateAndSetDirectEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco) {
bb427707 723 //
5f3c1b97 724 // Function that computes the Direct acceptance and efficiency correction
bb427707 725 // based on the simulated and reconstructed spectra
726 // and using the reconstructed spectra bin width
727 //
728 // eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
729 //
5f3c1b97 730
86bdcd8c 731 if(!fhRECpt || !hSimu || !hReco){
732 AliError("Hey, the reconstructed histogram was not set yet !");
733 return;
bb427707 734 }
735
86bdcd8c 736 fhDirectEffpt = EstimateEfficiencyRecoBin(hSimu,hReco,"fhDirectEffpt");
5f3c1b97 737 fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance #times efficiency");
bb427707 738
5f3c1b97 739}
bb427707 740
5f3c1b97 741//_________________________________________________________________________________________________________
86bdcd8c 742void AliHFPtSpectrum::EstimateAndSetFeedDownEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco) {
5f3c1b97 743 //
744 // Function that computes the Feed-Down acceptance and efficiency correction
745 // based on the simulated and reconstructed spectra
746 // and using the reconstructed spectra bin width
747 //
748 // eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
749 //
bb427707 750
86bdcd8c 751 if(!fhRECpt || !hSimu || !hReco){
752 AliError("Hey, the reconstructed histogram was not set yet !");
753 return;
bb427707 754 }
755
86bdcd8c 756 fhFeedDownEffpt = EstimateEfficiencyRecoBin(hSimu,hReco,"fhFeedDownEffpt");
5f3c1b97 757 fhFeedDownEffpt->SetNameTitle("fhFeedDownEffpt"," feed-down acceptance #times efficiency");
758
bb427707 759}
760
761//_________________________________________________________________________________________________________
65e55bbd 762Bool_t AliHFPtSpectrum::Initialize(){
763 //
764 // Initialization of the variables (histograms)
765 //
766
767 if (fFeedDownOption==0) {
768 AliInfo("Getting ready for the corrections without feed-down consideration");
769 } else if (fFeedDownOption==1) {
770 AliInfo("Getting ready for the fc feed-down correction calculation");
771 } else if (fFeedDownOption==2) {
772 AliInfo("Getting ready for the Nb feed-down correction calculation");
773 } else { AliError("The calculation option must be <=2"); return kFALSE; }
774
775 // Start checking the input histograms consistency
776 Bool_t areconsistent=kTRUE;
777
778 // General checks
779 if (!fhDirectEffpt || !fhRECpt) {
780 AliError(" Reconstructed spectra and/or the Nc efficiency distributions are not defined");
781 return kFALSE;
782 }
783 areconsistent &= CheckHistosConsistency(fhRECpt,fhDirectEffpt);
784 if (!areconsistent) {
785 AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
786 return kFALSE;
787 }
788 if (fFeedDownOption==0) return kTRUE;
789
790 //
791 // Common checks for options 1 (fc) & 2(Nb)
792 if (!fhFeedDownMCpt || !fhFeedDownEffpt) {
793 AliError(" Theoretical Nb and/or the Nb efficiency distributions are not defined");
794 return kFALSE;
795 }
796 areconsistent &= CheckHistosConsistency(fhRECpt,fhFeedDownMCpt);
797 areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownEffpt);
798 if (fAsymUncertainties) {
a17b17dd 799 if (!fhFeedDownMCptMax || !fhFeedDownMCptMin) {
65e55bbd 800 AliError(" Max/Min theoretical Nb distributions are not defined");
801 return kFALSE;
802 }
a17b17dd 803 areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownMCptMax);
65e55bbd 804 }
805 if (!areconsistent) {
806 AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
807 return kFALSE;
808 }
809 if (fFeedDownOption>1) return kTRUE;
810
811 //
812 // Now checks for option 1 (fc correction)
813 if (!fhDirectMCpt) {
814 AliError("Theoretical Nc distributions is not defined");
815 return kFALSE;
816 }
817 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhFeedDownMCpt);
818 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectEffpt);
819 if (fAsymUncertainties) {
a17b17dd 820 if (!fhDirectMCptMax || !fhDirectMCptMin) {
65e55bbd 821 AliError(" Max/Min theoretical Nc distributions are not defined");
822 return kFALSE;
823 }
a17b17dd 824 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectMCptMax);
65e55bbd 825 }
826 if (!areconsistent) {
827 AliInfo("Histograms required for fc correction are not consistent (bin width, bounds)");
828 return kFALSE;
829 }
830
831 return kTRUE;
832}
833
834//_________________________________________________________________________________________________________
86bdcd8c 835Bool_t AliHFPtSpectrum::CheckHistosConsistency(TH1D *h1, TH1D *h2){
65e55bbd 836 //
837 // Check the histograms consistency (bins, limits)
838 //
839
840 if (!h1 || !h2) {
841 AliError("One or both histograms don't exist");
842 return kFALSE;
843 }
844
845 Double_t binwidth1 = h1->GetBinWidth(1);
846 Double_t binwidth2 = h2->GetBinWidth(1);
847 Double_t min1 = h1->GetBinCenter(1) - (binwidth1/2.) ;
848// Double_t max1 = h1->GetBinCenter(nbins1) + (binwidth1/2.) ;
849 Double_t min2 = h2->GetBinCenter(1) - (binwidth2/2.) ;
850// Double_t max2 = h2->GetBinCenter(nbins2) + (binwidth2/2.) ;
851
852 if (binwidth1!=binwidth2) {
853 AliInfo(" histograms with different bin width");
854 return kFALSE;
855 }
856 if (min1!=min2) {
857 AliInfo(" histograms with different minimum");
858 return kFALSE;
859 }
860// if (max1!=max2) {
861// AliInfo(" histograms with different maximum");
862// return kFALSE;
863// }
864
865 return kTRUE;
866}
867
868//_________________________________________________________________________________________________________
a17b17dd 869void AliHFPtSpectrum::CalculateFeedDownCorrectionFc(){
65e55bbd 870 //
871 // Compute fc factor and its uncertainties bin by bin
872 // fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
873 //
86bdcd8c 874 // uncertainties: (conservative) combine the upper/lower N_b & N_c predictions together
875 // (extreme) combine the upper N_b predictions with the lower N_c predictions & viceversa
876 // the uncertainties on the acceptance x efficiency corrections are not included here
877 //
65e55bbd 878
879 // define the variables
880 Int_t nbins = fhRECpt->GetNbinsX();
881 Double_t binwidth = fhRECpt->GetBinWidth(1);
bb427707 882 Double_t *limits = new Double_t[nbins+1];
662607e3 883 Double_t *binwidths = new Double_t[nbins+1];
bb427707 884 Double_t xlow=0.;
86bdcd8c 885 for (Int_t i=1; i<=nbins; i++) {
bb427707 886 binwidth = fhRECpt->GetBinWidth(i);
887 xlow = fhRECpt->GetBinLowEdge(i);
888 limits[i-1] = xlow;
de97013c 889 binwidths[i] = binwidth;
bb427707 890 }
891 limits[nbins] = xlow + binwidth;
892
65e55bbd 893 Double_t correction=1.;
a17b17dd 894 Double_t theoryRatio=1.;
895 Double_t effRatio=1.;
86bdcd8c 896 Double_t correctionExtremeA=1., correctionExtremeB=1.;
897 Double_t theoryRatioExtremeA=1., theoryRatioExtremeB=1.;
898 Double_t correctionConservativeA=1., correctionConservativeB=1.;
899 Double_t theoryRatioConservativeA=1., theoryRatioConservativeB=1.;
900
65e55bbd 901 // declare the output histograms
86bdcd8c 902 if (!fhFc) fhFc = new TH1D("fhFc","fc correction factor",nbins,limits);
903 if (!fhFcMax) fhFcMax = new TH1D("fhFcMax","max fc correction factor",nbins,limits);
904 if (!fhFcMin) fhFcMin = new TH1D("fhFcMin","min fc correction factor",nbins,limits);
65e55bbd 905 // two local control histograms
bb427707 906 TH1D *hTheoryRatio = new TH1D("hTheoryRatio","Theoretical B-->D over c-->D (feed-down/direct) ratio",nbins,limits);
907 TH1D *hEffRatio = new TH1D("hEffRatio","Efficiency B-->D over c-->D (feed-down/direct) ratio",nbins,limits);
65e55bbd 908 // and the output TGraphAsymmErrors
86bdcd8c 909 if (fAsymUncertainties & !fgFcExtreme) {
910 fgFcExtreme = new TGraphAsymmErrors(nbins+1);
911 fgFcExtreme->SetNameTitle("fgFcExtreme","fgFcExtreme");
912 }
913 if (fAsymUncertainties & !fgFcConservative) {
914 fgFcConservative = new TGraphAsymmErrors(nbins+1);
915 fgFcConservative->SetNameTitle("fgFcConservative","fgFcConservative");
916 }
917
65e55bbd 918
919 //
920 // Compute fc
921 //
922 for (Int_t ibin=0; ibin<=nbins; ibin++) {
923
924 // theory_ratio = (N_b/N_c)
86bdcd8c 925 theoryRatio = (fhDirectMCpt->GetBinContent(ibin) && fhDirectMCpt->GetBinContent(ibin)!=0.) ?
926 fhFeedDownMCpt->GetBinContent(ibin) / fhDirectMCpt->GetBinContent(ibin) : 1.0 ;
927 //
928 // Calculate the uncertainty [ considering only the theoretical uncertainties on Nb & Nc for now !!! ]
929 //
930 // extreme A = direct-max, feed-down-min
931 theoryRatioExtremeA = (fhDirectMCptMax->GetBinContent(ibin) && fhDirectMCptMax->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) && fhDirectMCptMin->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) && fhDirectMCptMax->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) && fhDirectMCptMin->GetBinContent(ibin)!=0.) ?
941 fhFeedDownMCptMin->GetBinContent(ibin) / fhDirectMCptMin->GetBinContent(ibin) : 1.0 ;
bb427707 942
65e55bbd 943 // eff_ratio = (eff_b/eff_c)
86bdcd8c 944 effRatio = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ?
945 fhFeedDownEffpt->GetBinContent(ibin) / fhDirectEffpt->GetBinContent(ibin) : 1.0 ;
bb427707 946
65e55bbd 947 // fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
86bdcd8c 948 correction = ( 1. / ( 1 + ( effRatio * theoryRatio ) ) );
949 correctionExtremeA = ( 1. / ( 1 + ( effRatio * theoryRatioExtremeA ) ) );
950 correctionExtremeB = ( 1. / ( 1 + ( effRatio * theoryRatioExtremeB ) ) );
951 correctionConservativeA = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeA ) ) );
952 correctionConservativeB = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeB ) ) );
65e55bbd 953
954 // Fill in the histograms
a17b17dd 955 hTheoryRatio->SetBinContent(ibin,theoryRatio);
956 hEffRatio->SetBinContent(ibin,effRatio);
86bdcd8c 957 fhFc->SetBinContent(ibin,correction);
65e55bbd 958 if (fAsymUncertainties) {
959 Double_t x = fhDirectMCpt->GetBinCenter(ibin);
86bdcd8c 960 Double_t val[2] = { correctionExtremeA, correctionExtremeB };
961 Double_t uncExtremeMin = correction - TMath::MinElement(2,val);
962 Double_t uncExtremeMax = TMath::MaxElement(2,val) - correction;
963 fgFcExtreme->SetPoint(ibin,x,correction); // i,x,y
964 fgFcExtreme->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),uncExtremeMin,uncExtremeMax); // i,xl,xh,yl,yh
965 fhFcMax->SetBinContent(ibin,correction+uncExtremeMax);
966 fhFcMin->SetBinContent(ibin,correction-uncExtremeMin);
967 Double_t uncConservativeMin = correction - correctionConservativeA;
968 Double_t uncConservativeMax = correctionConservativeB - correction;
969 fgFcConservative->SetPoint(ibin,x,correction); // i,x,y
970 fgFcConservative->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),uncConservativeMin,uncConservativeMax); // i,xl,xh,yl,yh
65e55bbd 971 }
972
973 }
974
65e55bbd 975}
976
977//_________________________________________________________________________________________________________
a17b17dd 978void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumFc(){
65e55bbd 979 //
980 // Compute the feed-down corrected spectrum if feed-down correction is done via fc factor (bin by bin)
bb427707 981 // physics = reco * fc / bin-width
65e55bbd 982 //
86bdcd8c 983 // uncertainty: (stat) delta_physics = physics * sqrt ( (delta_reco/reco)^2 )
984 // (syst but feed-down) delta_physics = physics * sqrt ( (delta_reco_syst/reco)^2 )
985 // (feed-down syst) delta_physics = physics * sqrt ( (delta_fc/fc)^2 )
65e55bbd 986 //
987 // ( Calculation done bin by bin )
988
989 if (!fhFc || !fhRECpt) {
990 AliError(" Reconstructed or fc distributions are not defined");
991 return;
992 }
993
994 Int_t nbins = fhRECpt->GetNbinsX();
86bdcd8c 995 Double_t value = 0., errvalue = 0., errvalueMax= 0., errvalueMin= 0.;
996 Double_t valueExtremeMax= 0., valueExtremeMin= 0.;
997 Double_t valueConservativeMax= 0., valueConservativeMin= 0.;
65e55bbd 998 Double_t binwidth = fhRECpt->GetBinWidth(1);
bb427707 999 Double_t *limits = new Double_t[nbins+1];
662607e3 1000 Double_t *binwidths = new Double_t[nbins+1];
bb427707 1001 Double_t xlow=0.;
86bdcd8c 1002 for (Int_t i=1; i<=nbins; i++) {
bb427707 1003 binwidth = fhRECpt->GetBinWidth(i);
1004 xlow = fhRECpt->GetBinLowEdge(i);
1005 limits[i-1] = xlow;
de97013c 1006 binwidths[i] = binwidth;
bb427707 1007 }
86bdcd8c 1008 limits[nbins] = xlow + binwidth;
65e55bbd 1009
1010 // declare the output histograms
86bdcd8c 1011 if (!fhYieldCorr) fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by fc)",nbins,limits);
1012 if (!fhYieldCorrMax) fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by fc)",nbins,limits);
1013 if (!fhYieldCorrMin) fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by fc)",nbins,limits);
65e55bbd 1014 // and the output TGraphAsymmErrors
86bdcd8c 1015 if (fAsymUncertainties & !fgYieldCorr) fgYieldCorr = new TGraphAsymmErrors(nbins+1);
1016 if (fAsymUncertainties & !fgYieldCorrExtreme) fgYieldCorrExtreme = new TGraphAsymmErrors(nbins+1);
1017 if (fAsymUncertainties & !fgYieldCorrConservative) fgYieldCorrConservative = new TGraphAsymmErrors(nbins+1);
65e55bbd 1018
1019 //
1020 // Do the calculation
1021 //
1022 for (Int_t ibin=0; ibin<=nbins; ibin++) {
1023
1024 // calculate the value
86bdcd8c 1025 // physics = reco * fc / bin-width
65e55bbd 1026 value = fhRECpt->GetBinContent(ibin) * fhFc->GetBinContent(ibin) ;
bb427707 1027 value /= fhRECpt->GetBinWidth(ibin) ;
86bdcd8c 1028
1029 // Statistical uncertainty
1030 // (stat) delta_physics = physics * sqrt ( (delta_reco/reco)^2 )
1031 errvalue = value * (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) ;
1032
1033 // Calculate the systematic uncertainties
1034 // (syst but feed-down) delta_physics = physics * sqrt ( (delta_reco_syst/reco)^2 )
1035 // (feed-down syst) delta_physics = physics * sqrt ( (delta_fc/fc)^2 )
1036 //
65e55bbd 1037 // Protect against null denominator. If so, define uncertainty as null
1038 if (fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) {
1039
1040 if (fAsymUncertainties) {
1041
86bdcd8c 1042 // Systematics but feed-down
1043 if (fgRECSystematics) {
1044 errvalueMax = value * ( fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinContent(ibin) );
1045 errvalueMin = value * ( fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinContent(ibin) );
65e55bbd 1046 }
86bdcd8c 1047
1048 // Extreme feed-down systematics
1049 valueExtremeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcExtreme->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1050 valueExtremeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcExtreme->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1051
1052 // Conservative feed-down systematics
1053 valueConservativeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcConservative->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1054 valueConservativeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcConservative->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
65e55bbd 1055
1056 }
65e55bbd 1057
1058 }
86bdcd8c 1059 else { errvalueMax = 0.; errvalueMin = 0.; }
65e55bbd 1060
1061 // fill in the histograms
86bdcd8c 1062 fhYieldCorr->SetBinContent(ibin,value);
1063 fhYieldCorr->SetBinError(ibin,errvalue);
65e55bbd 1064 if (fAsymUncertainties) {
86bdcd8c 1065 Double_t center = fhYieldCorr->GetBinCenter(ibin);
65e55bbd 1066 fgYieldCorr->SetPoint(ibin,center,value); // i,x,y
86bdcd8c 1067 fgYieldCorr->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
1068 fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax);
1069 fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
1070 fgYieldCorrExtreme->SetPoint(ibin,center,value);
1071 fgYieldCorrExtreme->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),value-valueExtremeMin,valueExtremeMax-value);
1072 fgYieldCorrConservative->SetPoint(ibin,center,value);
1073 fgYieldCorrConservative->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),value-valueConservativeMin,valueConservativeMax-value);
65e55bbd 1074 }
1075
1076 }
1077
65e55bbd 1078}
1079
1080//_________________________________________________________________________________________________________
a17b17dd 1081void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumNb(Double_t deltaY, Double_t branchingRatioBintoFinalDecay) {
65e55bbd 1082 //
1083 // Compute the feed-down corrected spectrum if feed-down correction is done via Nb (bin by bin)
bb427707 1084 // physics = [ reco - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) ] / bin-width
65e55bbd 1085 //
86bdcd8c 1086 // uncertainty: (stat) delta_physics = sqrt ( (delta_reco)^2 ) / bin-width
1087 // (syst but feed-down) delta_physics = sqrt ( (delta_reco_syst)^2 ) / bin-width
1088 // (feed-down syst) delta_physics = sqrt ( (k*delta_lumi/lumi)^2 + (k*delta_eff_trig/eff_trig)^2
1089 // + (k*delta_Nb/Nb)^2 + (k*delta_eff/eff)^2 ) / bin-width
65e55bbd 1090 // where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th
1091 //
1092
1093 Int_t nbins = fhRECpt->GetNbinsX();
1094 Double_t binwidth = fhRECpt->GetBinWidth(1);
86bdcd8c 1095 Double_t value = 0., errvalue = 0., errvalueMax = 0., errvalueMin = 0., kfactor = 0.;
1096 Double_t errvalueExtremeMax = 0., errvalueExtremeMin = 0.;
bb427707 1097 Double_t *limits = new Double_t[nbins+1];
de97013c 1098 Double_t *binwidths = new Double_t[nbins+1];
bb427707 1099 Double_t xlow=0.;
86bdcd8c 1100 for (Int_t i=1; i<=nbins; i++) {
bb427707 1101 binwidth = fhRECpt->GetBinWidth(i);
1102 xlow = fhRECpt->GetBinLowEdge(i);
1103 limits[i-1] = xlow;
de97013c 1104 binwidths[i] = binwidth;
bb427707 1105 }
86bdcd8c 1106 limits[nbins] = xlow + binwidth;
65e55bbd 1107
1108 // declare the output histograms
86bdcd8c 1109 if (!fhYieldCorr) fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by Nb)",nbins,limits);
1110 if (!fhYieldCorrMax) fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by Nb)",nbins,limits);
1111 if (!fhYieldCorrMin) fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by Nb)",nbins,limits);
65e55bbd 1112 // and the output TGraphAsymmErrors
86bdcd8c 1113 if (fAsymUncertainties & !fgYieldCorr) fgYieldCorr = new TGraphAsymmErrors(nbins+1);
1114 if (fAsymUncertainties & !fgYieldCorrExtreme) fgYieldCorrExtreme = new TGraphAsymmErrors(nbins+1);
1115 if (fAsymUncertainties & !fgYieldCorrConservative) {
1116 fgYieldCorrConservative = new TGraphAsymmErrors(nbins+1);
1117 AliInfo(" Beware the conservative & extreme uncertainties are equal by definition !");
1118 }
65e55bbd 1119
1120 //
1121 // Do the calculation
1122 //
1123 for (Int_t ibin=0; ibin<=nbins; ibin++) {
1124
86bdcd8c 1125 // Calculate the value
1126 // physics = [ reco - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) ] / bin-width
a17b17dd 1127 value = fhRECpt->GetBinContent(ibin) - (deltaY*branchingRatioBintoFinalDecay*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) );
bb427707 1128 value /= fhRECpt->GetBinWidth(ibin);
65e55bbd 1129
86bdcd8c 1130 // Statistical uncertainty: delta_physics = sqrt ( (delta_reco)^2 ) / bin-width
1131 errvalue = fhRECpt->GetBinError(ibin);
1132 errvalue /= fhRECpt->GetBinWidth(ibin);
65e55bbd 1133
86bdcd8c 1134 // Systematic uncertainties
1135 // (syst but feed-down) delta_physics = sqrt ( (delta_reco_syst)^2 ) / bin-width
1136 // (feed-down syst) delta_physics = sqrt ( (k*delta_lumi/lumi)^2 + (k*delta_eff_trig/eff_trig)^2
1137 // + (k*delta_Nb/Nb)^2 + (k*delta_eff/eff)^2 ) / bin-width
1138 // where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th
1139 kfactor = deltaY*branchingRatioBintoFinalDecay*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) ;
1140 //
65e55bbd 1141 if (fAsymUncertainties) {
248ae766 1142 Double_t nb = fhFeedDownMCpt->GetBinContent(ibin);
1143 Double_t nbDmax = fhFeedDownMCptMax->GetBinContent(ibin) - fhFeedDownMCpt->GetBinContent(ibin);
1144 Double_t nbDmin = fhFeedDownMCpt->GetBinContent(ibin) - fhFeedDownMCptMin->GetBinContent(ibin);
86bdcd8c 1145
1146 // Systematics but feed-down
1147 if (fgRECSystematics){
1148 errvalueMax = fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinWidth(ibin) ;
1149 errvalueMin = fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinWidth(ibin);
1150 }
1151
1152 // Feed-down systematics
1153 // min value with the maximum Nb
1154 errvalueExtremeMin = TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
1155 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
248ae766 1156 ( (kfactor*nbDmax/nb)*(kfactor*nbDmax/nb) ) +
86bdcd8c 1157 ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) )
1158 / fhRECpt->GetBinWidth(ibin);
1159 // max value with the minimum Nb
1160 errvalueExtremeMax = TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
1161 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
248ae766 1162 ( (kfactor*nbDmin/nb)*(kfactor*nbDmin/nb) ) +
86bdcd8c 1163 ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) )
1164 / fhRECpt->GetBinWidth(ibin);
65e55bbd 1165 }
1166 else{ // Don't consider Nb uncertainty in this case [ to be tested!!! ]
86bdcd8c 1167 errvalueExtremeMax = TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
1168 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1169 ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) )
1170 / fhRECpt->GetBinWidth(ibin);
1171 errvalueExtremeMin = errvalueExtremeMax ;
65e55bbd 1172 }
1173
1174 // fill in histograms
86bdcd8c 1175 fhYieldCorr->SetBinContent(ibin,value);
1176 fhYieldCorr->SetBinError(ibin,errvalue);
65e55bbd 1177 if (fAsymUncertainties) {
86bdcd8c 1178 Double_t x = fhYieldCorr->GetBinCenter(ibin);
65e55bbd 1179 fgYieldCorr->SetPoint(ibin,x,value); // i,x,y
86bdcd8c 1180 fgYieldCorr->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
1181 fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax);
1182 fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
1183 fgYieldCorrExtreme->SetPoint(ibin,x,value); // i,x,y
1184 fgYieldCorrExtreme->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
1185 fgYieldCorrConservative->SetPoint(ibin,x,value); // i,x,y
1186 fgYieldCorrConservative->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
65e55bbd 1187 }
1188
1189 }
1190
65e55bbd 1191}
1192
1193
1194//_________________________________________________________________________________________________________
86bdcd8c 1195TH1D * AliHFPtSpectrum::ReweightHisto(TH1D *hToReweight, TH1D *hReference){
65e55bbd 1196 //
1197 // Function to reweight histograms for testing purposes:
1198 // This function takes the histo hToReweight and reweights
1199 // it (its pt shape) with respect to hReference
1200 //
1201
1202 // check histograms consistency
1203 Bool_t areconsistent=kTRUE;
1204 areconsistent &= CheckHistosConsistency(hToReweight,hReference);
1205 if (!areconsistent) {
1206 AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
1207 return NULL;
1208 }
1209
1210 // define a new empty histogram
86bdcd8c 1211 TH1D *hReweighted = (TH1D*)hToReweight->Clone("hReweighted");
65e55bbd 1212 hReweighted->Reset();
1213 Double_t weight=1.0;
1214 Double_t yvalue=1.0;
a17b17dd 1215 Double_t integralRef = hReference->Integral();
1216 Double_t integralH = hToReweight->Integral();
65e55bbd 1217
1218 // now reweight the spectra
1219 //
1220 // the weight is the relative probability of the given pt bin in the reference histo
1221 // divided by its relative probability (to normalize it) on the histo to re-weight
1222 for (Int_t i=0; i<=hToReweight->GetNbinsX(); i++) {
a17b17dd 1223 weight = (hReference->GetBinContent(i)/integralRef) / (hToReweight->GetBinContent(i)/integralH) ;
65e55bbd 1224 yvalue = hToReweight->GetBinContent(i);
1225 hReweighted->SetBinContent(i,yvalue*weight);
1226 }
1227
86bdcd8c 1228 return (TH1D*)hReweighted;
65e55bbd 1229}
1230
1231//_________________________________________________________________________________________________________
86bdcd8c 1232TH1D * AliHFPtSpectrum::ReweightRecHisto(TH1D *hRecToReweight, TH1D *hMCToReweight, TH1D *hMCReference){
65e55bbd 1233 //
1234 // Function to reweight histograms for testing purposes:
1235 // This function takes the histo hToReweight and reweights
1236 // it (its pt shape) with respect to hReference /hMCToReweight
1237 //
1238
1239 // check histograms consistency
1240 Bool_t areconsistent=kTRUE;
1241 areconsistent &= CheckHistosConsistency(hMCToReweight,hMCReference);
1242 areconsistent &= CheckHistosConsistency(hRecToReweight,hMCReference);
1243 if (!areconsistent) {
1244 AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
1245 return NULL;
1246 }
1247
1248 // define a new empty histogram
86bdcd8c 1249 TH1D *hReweighted = (TH1D*)hMCToReweight->Clone("hReweighted");
65e55bbd 1250 hReweighted->Reset();
86bdcd8c 1251 TH1D *hRecReweighted = (TH1D*)hRecToReweight->Clone("hRecReweighted");
65e55bbd 1252 hRecReweighted->Reset();
1253 Double_t weight=1.0;
1254 Double_t yvalue=1.0, yrecvalue=1.0;
a17b17dd 1255 Double_t integralRef = hMCReference->Integral();
1256 Double_t integralH = hMCToReweight->Integral();
65e55bbd 1257
1258 // now reweight the spectra
1259 //
1260 // the weight is the relative probability of the given pt bin
1261 // that should be applied in the MC histo to get the reference histo shape
1262 // Probabilities are properly normalized.
1263 for (Int_t i=0; i<=hMCToReweight->GetNbinsX(); i++) {
a17b17dd 1264 weight = (hMCReference->GetBinContent(i)/integralRef) / (hMCToReweight->GetBinContent(i)/integralH) ;
65e55bbd 1265 yvalue = hMCToReweight->GetBinContent(i);
1266 hReweighted->SetBinContent(i,yvalue*weight);
1267 yrecvalue = hRecToReweight->GetBinContent(i);
1268 hRecReweighted->SetBinContent(i,yrecvalue*weight);
1269 }
1270
86bdcd8c 1271 return (TH1D*)hRecReweighted;
65e55bbd 1272}
1273