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