]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGHF/vertexingHF/AliHFPtSpectrum.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fix typo
[u/mrichter/AliRoot.git] / PWGHF / vertexingHF / AliHFPtSpectrum.cxx
CommitLineData
b188dc47 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/* $Id$ */
17
18//***********************************************************************
19// Class AliHFPtSpectrum
20// Base class for feed-down corrections on heavy-flavour decays
21// computes the cross-section via one of the three implemented methods:
22// 0) Consider no feed-down prediction
23// 1) Subtract the feed-down with the "fc" method
24// Yield = Reco * fc; where fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) ) ;
25// 2) Subtract the feed-down with the "Nb" method
26// Yield = Reco - Feed-down (exact formula on the function implementation)
27//
28// (the corrected yields per bin are divided by the bin-width)
29//
30//
31// In HIC you can also evaluate how the feed-down correction is influenced by an energy loss hypothesis:
32// Raa(c-->D) / Raa(b-->D) defined here as Rcb for the "fc" method
33// Raa(b-->D) defined here as Rb for the "Nb" method
34//
35// Author: Z.Conesa, zconesa@in2p3.fr
36//***********************************************************************
37
38#include <Riostream.h>
39
40#include "TMath.h"
41#include "TH1.h"
42#include "TH1D.h"
43#include "TH2.h"
44#include "TH2D.h"
45#include "TNtuple.h"
46#include "TGraphAsymmErrors.h"
47#include "TNamed.h"
48#include "TCanvas.h"
49#include "TLegend.h"
50
51//#include "AliLog.h"
52#include "AliHFSystErr.h"
53#include "AliHFPtSpectrum.h"
54
55ClassImp(AliHFPtSpectrum)
56
57//_________________________________________________________________________________________________________
58AliHFPtSpectrum::AliHFPtSpectrum(const char* name, const char* title, Int_t option):
59 TNamed(name,title),
60 fhDirectMCpt(NULL),
61 fhFeedDownMCpt(NULL),
62 fhDirectMCptMax(NULL),
63 fhDirectMCptMin(NULL),
64 fhFeedDownMCptMax(NULL),
65 fhFeedDownMCptMin(NULL),
66 fhDirectEffpt(NULL),
67 fhFeedDownEffpt(NULL),
68 fhRECpt(NULL),
69 fgRECSystematics(NULL),
70 fNevts(1),
71 fLuminosity(),
72 fTrigEfficiency(),
73 fGlobalEfficiencyUncertainties(),
74 fTab(),
75 fhFc(NULL),
76 fhFcMax(NULL),
77 fhFcMin(NULL),
78 fhFcRcb(NULL),
79 fgFcExtreme(NULL),
80 fgFcConservative(NULL),
81 fhYieldCorr(NULL),
82 fhYieldCorrMax(NULL),
83 fhYieldCorrMin(NULL),
84 fhYieldCorrRcb(NULL),
85 fgYieldCorr(NULL),
86 fgYieldCorrExtreme(NULL),
87 fgYieldCorrConservative(NULL),
88 fhSigmaCorr(NULL),
89 fhSigmaCorrMax(NULL),
90 fhSigmaCorrMin(NULL),
91 fhSigmaCorrDataSyst(NULL),
92 fhSigmaCorrRcb(NULL),
93 fgSigmaCorr(NULL),
94 fgSigmaCorrExtreme(NULL),
95 fgSigmaCorrConservative(NULL),
96 fnSigma(NULL),
97 fnHypothesis(NULL),
98 fFeedDownOption(option),
99 fAsymUncertainties(kTRUE),
100 fPbPbElossHypothesis(kFALSE),
101 fIsStatUncEff(kTRUE),
102 fParticleAntiParticle(2),
103 fhStatUncEffcSigma(NULL),
104 fhStatUncEffbSigma(NULL),
105 fhStatUncEffcFD(NULL),
106 fhStatUncEffbFD(NULL)
107{
108 //
109 // Default constructor
110 //
111
112 fLuminosity[0]=1.; fLuminosity[1]=0.;
113 fTrigEfficiency[0]=1.; fTrigEfficiency[1]=0.;
114 fGlobalEfficiencyUncertainties[0]=0.; fGlobalEfficiencyUncertainties[1]=0.;
115 fTab[0]=1.; fTab[1]=0.;
116
117}
118
119//_________________________________________________________________________________________________________
120AliHFPtSpectrum::AliHFPtSpectrum(const AliHFPtSpectrum &rhs):
121 TNamed(rhs),
122 fhDirectMCpt(rhs.fhDirectMCpt),
123 fhFeedDownMCpt(rhs.fhFeedDownMCpt),
124 fhDirectMCptMax(rhs.fhDirectMCptMax),
125 fhDirectMCptMin(rhs.fhDirectMCptMin),
126 fhFeedDownMCptMax(rhs.fhFeedDownMCptMax),
127 fhFeedDownMCptMin(rhs.fhFeedDownMCptMin),
128 fhDirectEffpt(rhs.fhDirectEffpt),
129 fhFeedDownEffpt(rhs.fhFeedDownEffpt),
130 fhRECpt(rhs.fhRECpt),
131 fgRECSystematics(rhs.fgRECSystematics),
132 fNevts(rhs.fNevts),
133 fLuminosity(),
134 fTrigEfficiency(),
135 fGlobalEfficiencyUncertainties(),
136 fTab(),
137 fhFc(rhs.fhFc),
138 fhFcMax(rhs.fhFcMax),
139 fhFcMin(rhs.fhFcMin),
140 fhFcRcb(rhs.fhFcRcb),
141 fgFcExtreme(rhs.fgFcExtreme),
142 fgFcConservative(rhs.fgFcConservative),
143 fhYieldCorr(rhs.fhYieldCorr),
144 fhYieldCorrMax(rhs.fhYieldCorrMax),
145 fhYieldCorrMin(rhs.fhYieldCorrMin),
146 fhYieldCorrRcb(rhs.fhYieldCorrRcb),
147 fgYieldCorr(rhs.fgYieldCorr),
148 fgYieldCorrExtreme(rhs.fgYieldCorrExtreme),
149 fgYieldCorrConservative(rhs.fgYieldCorrConservative),
150 fhSigmaCorr(rhs.fhSigmaCorr),
151 fhSigmaCorrMax(rhs.fhSigmaCorrMax),
152 fhSigmaCorrMin(rhs.fhSigmaCorrMin),
153 fhSigmaCorrDataSyst(rhs.fhSigmaCorrDataSyst),
154 fhSigmaCorrRcb(rhs.fhSigmaCorrRcb),
155 fgSigmaCorr(rhs.fgSigmaCorr),
156 fgSigmaCorrExtreme(rhs.fgSigmaCorrExtreme),
157 fgSigmaCorrConservative(rhs.fgSigmaCorrConservative),
158 fnSigma(rhs.fnSigma),
159 fnHypothesis(rhs.fnHypothesis),
160 fFeedDownOption(rhs.fFeedDownOption),
161 fAsymUncertainties(rhs.fAsymUncertainties),
162 fPbPbElossHypothesis(rhs.fPbPbElossHypothesis),
163 fIsStatUncEff(rhs.fIsStatUncEff),
164 fParticleAntiParticle(rhs.fParticleAntiParticle),
165 fhStatUncEffcSigma(NULL),
166 fhStatUncEffbSigma(NULL),
167 fhStatUncEffcFD(NULL),
168 fhStatUncEffbFD(NULL)
169{
170 //
171 // Copy constructor
172 //
173
174 for(Int_t i=0; i<2; i++){
175 fLuminosity[i] = rhs.fLuminosity[i];
176 fTrigEfficiency[i] = rhs.fTrigEfficiency[i];
177 fGlobalEfficiencyUncertainties[i] = rhs.fGlobalEfficiencyUncertainties[i];
178 fTab[i] = rhs.fTab[i];
179 }
180
181}
182
183//_________________________________________________________________________________________________________
184AliHFPtSpectrum &AliHFPtSpectrum::operator=(const AliHFPtSpectrum &source){
185 //
186 // Assignment operator
187 //
188
189 if (&source == this) return *this;
190
191 fhDirectMCpt = source.fhDirectMCpt;
192 fhFeedDownMCpt = source.fhFeedDownMCpt;
193 fhDirectMCptMax = source.fhDirectMCptMax;
194 fhDirectMCptMin = source.fhDirectMCptMin;
195 fhFeedDownMCptMax = source.fhFeedDownMCptMax;
196 fhFeedDownMCptMin = source.fhFeedDownMCptMin;
197 fhDirectEffpt = source.fhDirectEffpt;
198 fhFeedDownEffpt = source.fhFeedDownEffpt;
199 fhRECpt = source.fhRECpt;
200 fgRECSystematics = source.fgRECSystematics;
201 fNevts = source.fNevts;
202 fhFc = source.fhFc;
203 fhFcMax = source.fhFcMax;
204 fhFcMin = source.fhFcMin;
205 fhFcRcb = source.fhFcRcb;
206 fgFcExtreme = source.fgFcExtreme;
207 fgFcConservative = source.fgFcConservative;
208 fhYieldCorr = source.fhYieldCorr;
209 fhYieldCorrMax = source.fhYieldCorrMax;
210 fhYieldCorrMin = source.fhYieldCorrMin;
211 fhYieldCorrRcb = source.fhYieldCorrRcb;
212 fgYieldCorr = source.fgYieldCorr;
213 fgYieldCorrExtreme = source.fgYieldCorrExtreme;
214 fgYieldCorrConservative = source.fgYieldCorrConservative;
215 fhSigmaCorr = source.fhSigmaCorr;
216 fhSigmaCorrMax = source.fhSigmaCorrMax;
217 fhSigmaCorrMin = source.fhSigmaCorrMin;
218 fhSigmaCorrDataSyst = source.fhSigmaCorrDataSyst;
219 fhSigmaCorrRcb = source.fhSigmaCorrRcb;
220 fgSigmaCorr = source.fgSigmaCorr;
221 fgSigmaCorrExtreme = source.fgSigmaCorrExtreme;
222 fgSigmaCorrConservative = source.fgSigmaCorrConservative;
223 fnSigma = source.fnSigma;
224 fnHypothesis = source.fnHypothesis;
225 fFeedDownOption = source.fFeedDownOption;
226 fAsymUncertainties = source.fAsymUncertainties;
227 fPbPbElossHypothesis = source.fPbPbElossHypothesis;
228 fIsStatUncEff = source.fIsStatUncEff;
229 fParticleAntiParticle = source.fParticleAntiParticle;
230
231 for(Int_t i=0; i<2; i++){
232 fLuminosity[i] = source.fLuminosity[i];
233 fTrigEfficiency[i] = source.fTrigEfficiency[i];
234 fGlobalEfficiencyUncertainties[i] = source.fGlobalEfficiencyUncertainties[i];
235 fTab[i] = source.fTab[i];
236 }
237
238 return *this;
239}
240
241//_________________________________________________________________________________________________________
242AliHFPtSpectrum::~AliHFPtSpectrum(){
243 //
244 // Destructor
245 //
246 if (fhDirectMCpt) delete fhDirectMCpt;
247 if (fhFeedDownMCpt) delete fhFeedDownMCpt;
248 if (fhDirectMCptMax) delete fhDirectMCptMax;
249 if (fhDirectMCptMin) delete fhDirectMCptMin;
250 if (fhFeedDownMCptMax) delete fhFeedDownMCptMax;
251 if (fhFeedDownMCptMin) delete fhFeedDownMCptMin;
252 if (fhDirectEffpt) delete fhDirectEffpt;
253 if (fhFeedDownEffpt) delete fhFeedDownEffpt;
254 if (fhRECpt) delete fhRECpt;
255 if (fgRECSystematics) delete fgRECSystematics;
256 if (fhFc) delete fhFc;
257 if (fhFcMax) delete fhFcMax;
258 if (fhFcMin) delete fhFcMin;
259 if (fhFcRcb) delete fhFcRcb;
260 if (fgFcExtreme) delete fgFcExtreme;
261 if (fgFcConservative) delete fgFcConservative;
262 if (fhYieldCorr) delete fhYieldCorr;
263 if (fhYieldCorrMax) delete fhYieldCorrMax;
264 if (fhYieldCorrMin) delete fhYieldCorrMin;
265 if (fhYieldCorrRcb) delete fhYieldCorrRcb;
266 if (fgYieldCorr) delete fgYieldCorr;
267 if (fgYieldCorrExtreme) delete fgYieldCorrExtreme;
268 if (fgYieldCorrConservative) delete fgYieldCorrConservative;
269 if (fhSigmaCorr) delete fhSigmaCorr;
270 if (fhSigmaCorrMax) delete fhSigmaCorrMax;
271 if (fhSigmaCorrMin) delete fhSigmaCorrMin;
272 if (fhSigmaCorrDataSyst) delete fhSigmaCorrDataSyst;
273 if (fgSigmaCorr) delete fgSigmaCorr;
274 if (fgSigmaCorrExtreme) delete fgSigmaCorrExtreme;
275 if (fgSigmaCorrConservative) delete fgSigmaCorrConservative;
276 if (fnSigma) delete fnSigma;
277 if (fnHypothesis) delete fnHypothesis;
278}
279
280
281//_________________________________________________________________________________________________________
282TH1D * AliHFPtSpectrum::RebinTheoreticalSpectra(TH1D *hTheory, const char *name) {
283 //
284 // Function to rebin the theoretical spectrum
285 // with respect to the real-data reconstructed spectrum binning
286 //
287
288 if (!hTheory || !fhRECpt) {
289 AliError("Feed-down or reconstructed spectra don't exist");
290 return NULL;
291 }
292
293 //
294 // Get the reconstructed spectra bins & limits
295 Int_t nbins = fhRECpt->GetNbinsX();
296 Int_t nbinsMC = hTheory->GetNbinsX();
297 Double_t *limits = new Double_t[nbins+1];
298 Double_t xlow=0., binwidth=0.;
299 for (Int_t i=1; i<=nbins; i++) {
300 binwidth = fhRECpt->GetBinWidth(i);
301 xlow = fhRECpt->GetBinLowEdge(i);
302 limits[i-1] = xlow;
303 }
304 limits[nbins] = xlow + binwidth;
305
306 // Check that the reconstructed spectra binning
307 // is larger than the theoretical one
308 Double_t thbinwidth = hTheory->GetBinWidth(1);
309 for (Int_t i=1; i<=nbins; i++) {
310 binwidth = fhRECpt->GetBinWidth(i);
311 if ( thbinwidth > binwidth ) {
312 AliInfo(" Beware it seems that the reconstructed spectra has a smaller binning than the theoretical predictions !! ");
313 }
314 }
315
316 //
317 // Define a new histogram with the real-data reconstructed spectrum binning
318 TH1D * hTheoryRebin = new TH1D(name," theoretical rebinned prediction",nbins,limits);
319
320 Double_t sum[nbins], items[nbins];
321 for (Int_t ibin=0; ibin<nbins; ibin++) {
322 sum[ibin]=0.; items[ibin]=0.;
323 }
324 for (Int_t ibin=0; ibin<=nbinsMC; ibin++){
325
326 for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++){
327 if (hTheory->GetBinCenter(ibin)>limits[ibinrec] &&
328 hTheory->GetBinCenter(ibin)<limits[ibinrec+1]){
329 sum[ibinrec]+=hTheory->GetBinContent(ibin);
330 items[ibinrec]+=1.;
331 }
332 }
333
334 }
335
336 // set the theoretical rebinned spectra to ( sum-bins / n-bins ) per new bin
337 for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++) {
338 hTheoryRebin->SetBinContent(ibinrec+1,sum[ibinrec]/items[ibinrec]);
339 }
340
341 return (TH1D*)hTheoryRebin;
342}
343
344//_________________________________________________________________________________________________________
345void AliHFPtSpectrum::SetMCptSpectra(TH1D *hDirect, TH1D *hFeedDown){
346 //
347 // Set the MonteCarlo or Theoretical spectra
348 // both for direct and feed-down contributions
349 //
350
351 if (!hDirect || !hFeedDown || !fhRECpt) {
352 AliError("One or both (direct, feed-down) spectra or the reconstructed spectra don't exist");
353 return;
354 }
355
356 Bool_t areconsistent = kTRUE;
357 areconsistent = CheckHistosConsistency(hDirect,hFeedDown);
358 if (!areconsistent) {
359 AliInfo("Histograms are not consistent (bin width, bounds)");
360 return;
361 }
362
363 //
364 // Rebin the theoretical predictions to the reconstructed spectra binning
365 //
366 fhDirectMCpt = RebinTheoreticalSpectra(hDirect,"fhDirectMCpt");
367 fhDirectMCpt->SetNameTitle("fhDirectMCpt"," direct theoretical prediction");
368 fhFeedDownMCpt = RebinTheoreticalSpectra(hFeedDown,"fhFeedDownMCpt");
369 fhFeedDownMCpt->SetNameTitle("fhFeedDownMCpt"," feed-down theoretical prediction");
370
371}
372
373//_________________________________________________________________________________________________________
374void AliHFPtSpectrum::SetFeedDownMCptSpectra(TH1D *hFeedDown){
375 //
376 // Set the MonteCarlo or Theoretical spectra
377 // for feed-down contribution
378 //
379
380 if (!hFeedDown || !fhRECpt) {
381 AliError("Feed-down or reconstructed spectra don't exist");
382 return;
383 }
384
385 //
386 // Rebin the theoretical predictions to the reconstructed spectra binning
387 //
388 fhFeedDownMCpt = RebinTheoreticalSpectra(hFeedDown,"fhFeedDownMCpt");
389 fhFeedDownMCpt->SetNameTitle("fhFeedDownMCpt"," feed-down theoretical prediction");
390
391}
392
393//_________________________________________________________________________________________________________
394void AliHFPtSpectrum::SetMCptDistributionsBounds(TH1D *hDirectMax, TH1D *hDirectMin, TH1D *hFeedDownMax, TH1D *hFeedDownMin){
395 //
396 // Set the maximum and minimum MonteCarlo or Theoretical spectra
397 // both for direct and feed-down contributions
398 // used in case uncertainties are asymmetric and ca not be on the "basic histograms"
399 //
400
401 if (!hDirectMax || !hDirectMin || !hFeedDownMax|| !hFeedDownMin || !fhRECpt) {
402 AliError("One or all of the max/min direct/feed-down or the reconstructed spectra don't exist");
403 return;
404 }
405
406 Bool_t areconsistent = kTRUE;
407 areconsistent &= CheckHistosConsistency(hDirectMax,hDirectMin);
408 areconsistent &= CheckHistosConsistency(hFeedDownMax,hFeedDownMin);
409 areconsistent &= CheckHistosConsistency(hDirectMax,hFeedDownMax);
410 if (!areconsistent) {
411 AliInfo("Histograms are not consistent (bin width, bounds)");
412 return;
413 }
414
415
416 //
417 // Rebin the theoretical predictions to the reconstructed spectra binning
418 //
419 fhDirectMCptMax = RebinTheoreticalSpectra(hDirectMax,"fhDirectMCptMax");
420 fhDirectMCptMax->SetNameTitle("fhDirectMCptMax"," maximum direct theoretical prediction");
421 fhDirectMCptMin = RebinTheoreticalSpectra(hDirectMin,"fhDirectMCptMin");
422 fhDirectMCptMin->SetNameTitle("fhDirectMCptMin"," minimum direct theoretical prediction");
423 fhFeedDownMCptMax = RebinTheoreticalSpectra(hFeedDownMax,"fhFeedDownMCptMax");
424 fhFeedDownMCptMax->SetNameTitle("fhFeedDownMCptMax"," maximum feed-down theoretical prediction");
425 fhFeedDownMCptMin = RebinTheoreticalSpectra(hFeedDownMin,"fhFeedDownMCptMin");
426 fhFeedDownMCptMin->SetNameTitle("fhFeedDownMCptMin"," minimum feed-down theoretical prediction");
427
428}
429
430//_________________________________________________________________________________________________________
431void AliHFPtSpectrum::SetFeedDownMCptDistributionsBounds(TH1D *hFeedDownMax, TH1D *hFeedDownMin){
432 //
433 // Set the maximum and minimum MonteCarlo or Theoretical spectra
434 // for feed-down contributions
435 // used in case uncertainties are asymmetric and can not be on the "basic histogram"
436 //
437
438 if (!hFeedDownMax || !hFeedDownMin || !fhRECpt) {
439 AliError("One or all of the max/min direct/feed-down spectra don't exist");
440 return;
441 }
442
443 Bool_t areconsistent = kTRUE;
444 areconsistent &= CheckHistosConsistency(hFeedDownMax,hFeedDownMin);
445 if (!areconsistent) {
446 AliInfo("Histograms are not consistent (bin width, bounds)");
447 return;
448 }
449
450
451 //
452 // Rebin the theoretical predictions to the reconstructed spectra binning
453 //
454 fhFeedDownMCptMax = RebinTheoreticalSpectra(hFeedDownMax,"fhFeedDownMCptMax");
455 fhFeedDownMCptMax->SetNameTitle("fhFeedDownMCptMax"," maximum feed-down theoretical prediction");
456 fhFeedDownMCptMin = RebinTheoreticalSpectra(hFeedDownMin,"fhFeedDownMCptMin");
457 fhFeedDownMCptMin->SetNameTitle("fhFeedDownMCptMin"," minimum feed-down theoretical prediction");
458
459}
460
461//_________________________________________________________________________________________________________
462void AliHFPtSpectrum::SetDirectAccEffCorrection(TH1D *hDirectEff){
463 //
464 // Set the Acceptance and Efficiency corrections
465 // for the direct contribution
466 //
467
468 if (!hDirectEff) {
469 AliError("The direct acceptance and efficiency corrections doesn't exist");
470 return;
471 }
472
473 fhDirectEffpt = (TH1D*)hDirectEff->Clone();
474 fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance x efficiency correction");
475}
476
477//_________________________________________________________________________________________________________
478void AliHFPtSpectrum::SetAccEffCorrection(TH1D *hDirectEff, TH1D *hFeedDownEff){
479 //
480 // Set the Acceptance and Efficiency corrections
481 // both for direct and feed-down contributions
482 //
483
484 if (!hDirectEff || !hFeedDownEff) {
485 AliError("One or both (direct, feed-down) acceptance and efficiency corrections don't exist");
486 return;
487 }
488
489 Bool_t areconsistent=kTRUE;
490 areconsistent = CheckHistosConsistency(hDirectEff,hFeedDownEff);
491 if (!areconsistent) {
492 AliInfo("Histograms are not consistent (bin width, bounds)");
493 return;
494 }
495
496 fhDirectEffpt = (TH1D*)hDirectEff->Clone();
497 fhFeedDownEffpt = (TH1D*)hFeedDownEff->Clone();
498 fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance x efficiency correction");
499 fhFeedDownEffpt->SetNameTitle("fhFeedDownEffpt"," feed-down acceptance x efficiency correction");
500}
501
502//_________________________________________________________________________________________________________
503void AliHFPtSpectrum::SetReconstructedSpectrum(TH1D *hRec) {
504 //
505 // Set the reconstructed spectrum
506 //
507
508 if (!hRec) {
509 AliError("The reconstructed spectrum doesn't exist");
510 return;
511 }
512
513 fhRECpt = (TH1D*)hRec->Clone();
514 fhRECpt->SetNameTitle("fhRECpt"," reconstructed spectrum");
515}
516
517//_________________________________________________________________________________________________________
518void AliHFPtSpectrum::SetReconstructedSpectrumSystematics(TGraphAsymmErrors *gRec) {
519 //
520 // Set the reconstructed spectrum (uncorrected yield) systematic uncertainties
521 //
522
523 // Check the compatibility with the reconstructed spectrum
524 Double_t gbinwidth = gRec->GetErrorXlow(1) + gRec->GetErrorXhigh(1) ;
525 Double_t hbinwidth = fhRECpt->GetBinWidth(1);
526 Double_t gxbincenter=0., gybincenter=0.;
527 gRec->GetPoint(1,gxbincenter,gybincenter);
528 Double_t hbincenter = fhRECpt->GetBinCenter(1);
529 if ( (gbinwidth != hbinwidth) || (gxbincenter!=hbincenter) ) {
530 AliError(" The reconstructed spectrum and its systematics don't seem compatible");
531 return;
532 }
533
534 fgRECSystematics = gRec;
535}
536
537//_________________________________________________________________________________________________________
538void AliHFPtSpectrum::ComputeHFPtSpectrum(Double_t deltaY, Double_t branchingRatioC, Double_t branchingRatioBintoFinalDecay) {
539 //
540 // Main function to compute the corrected cross-section:
541 // variables : analysed delta_y, BR for the final correction,
542 // BR b --> D --> decay (relative to the input theoretical prediction)
543 //
544 // Sigma = ( 1. / (lumi * delta_y * BR_c * ParticleAntiPartFactor * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
545 //
546 // Uncertainties: (stat) delta_sigma = sigma * sqrt ( (delta_spectra/spectra)^2 )
547 // (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 )
548 // (feed-down syst) delta_sigma = sigma * sqrt ( (delta_spectra_fd/spectra_fd)^2 )
549 //
550 // In HIC the feed-down correction varies with an energy loss hypothesis:
551 // Raa(c-->D) / Raa(b-->D) for the "fc" method, Raa(b-->D) for the "Nb" method (see exact formulas in the functions)
552 //
553
554 //
555 // First: Initialization
556 //
557 Bool_t areHistosOk = Initialize();
558 if (!areHistosOk) {
559 AliInfo(" Histos not properly initialized. Check : inconsistent binning ? missing histos ?");
560 return;
561 }
562 // Reset the statistical uncertainties on the efficiencies if needed
563 if(!fIsStatUncEff) ResetStatUncEff();
564
565 //
566 // Second: Correct for feed-down
567 //
568 if (fFeedDownOption==1) {
569 // Compute the feed-down correction via fc-method
570 CalculateFeedDownCorrectionFc();
571 // Correct the yield for feed-down correction via fc-method
572 CalculateFeedDownCorrectedSpectrumFc();
573 }
574 else if (fFeedDownOption==2) {
575 // Correct the yield for feed-down correction via Nb-method
576 CalculateFeedDownCorrectedSpectrumNb(deltaY,branchingRatioBintoFinalDecay);
577 }
578 else if (fFeedDownOption==0) {
579 // If there is no need for feed-down correction,
580 // the "corrected" yield is equal to the raw yield
581 fhYieldCorr = (TH1D*)fhRECpt->Clone();
582 fhYieldCorr->SetNameTitle("fhYieldCorr","un-corrected yield");
583 fhYieldCorrMax = (TH1D*)fhRECpt->Clone();
584 fhYieldCorrMin = (TH1D*)fhRECpt->Clone();
585 fhYieldCorrMax->SetNameTitle("fhYieldCorrMax","un-corrected yield");
586 fhYieldCorrMin->SetNameTitle("fhYieldCorrMin","un-corrected yield");
587 fAsymUncertainties=kFALSE;
588 }
589 else {
590 AliInfo(" Are you sure the feed-down correction option is right ?");
591 }
592
593
594 // Print out information
595 printf("\n\n Correcting the spectra with : \n luminosity = %2.2e +- %2.2e, trigger efficiency = %2.2e +- %2.2e, \n delta_y = %2.2f, BR_c = %2.2e, BR_b_decay = %2.2e \n %2.2f percent uncertainty on the efficiencies, and %2.2f percent uncertainty on the b/c efficiencies ratio \n\n",fLuminosity[0],fLuminosity[1],fTrigEfficiency[0],fTrigEfficiency[1],deltaY,branchingRatioC,branchingRatioBintoFinalDecay,fGlobalEfficiencyUncertainties[0],fGlobalEfficiencyUncertainties[1]);
596 if (fPbPbElossHypothesis) printf("\n\n The considered Tab is %4.2e +- %2.2e \n\n",fTab[0],fTab[1]);
597
598 //
599 // Finally: Correct from yields to cross-section
600 //
601 Int_t nbins = fhRECpt->GetNbinsX();
602 Double_t binwidth = fhRECpt->GetBinWidth(1);
603 Double_t *limits = new Double_t[nbins+1];
604 Double_t *binwidths = new Double_t[nbins];
605 Double_t xlow=0.;
606 for (Int_t i=1; i<=nbins; i++) {
607 binwidth = fhRECpt->GetBinWidth(i);
608 xlow = fhRECpt->GetBinLowEdge(i);
609 limits[i-1] = xlow;
610 binwidths[i-1] = binwidth;
611 }
612 limits[nbins] = xlow + binwidth;
613
614
615 // declare the output histograms
616 fhSigmaCorr = new TH1D("fhSigmaCorr","corrected sigma",nbins,limits);
617 fhSigmaCorrMax = new TH1D("fhSigmaCorrMax","max corrected sigma",nbins,limits);
618 fhSigmaCorrMin = new TH1D("fhSigmaCorrMin","min corrected sigma",nbins,limits);
619 fhSigmaCorrDataSyst = new TH1D("fhSigmaCorrDataSyst","data syst uncertainties on the corrected sigma",nbins,limits);
620 if (fPbPbElossHypothesis && fFeedDownOption==1) {
621 fhSigmaCorrRcb = new TH2D("fhSigmaCorrRcb","corrected sigma vs Rcb Eloss hypothesis; p_{T} [GeV/c] ; Rcb Eloss hypothesis ; #sigma",nbins,limits,800,0.,4.);
622 fnSigma = new TNtuple("fnSigma"," Sigma ntuple calculation","pt:Signal:Rcb:fc:Yield:Sigma:SigmaStatUnc:SigmaMax:SigmaMin");
623 }
624 if (fPbPbElossHypothesis && fFeedDownOption==2) {
625 fhSigmaCorrRcb = new TH2D("fhSigmaCorrRcb","corrected sigma vs Rb Eloss hypothesis; p_{T} [GeV/c] ; Rb Eloss hypothesis ; #sigma",nbins,limits,800,0.,4.);
626 fnSigma = new TNtuple("fnSigma"," Sigma ntuple calculation","pt:Signal:Rb:fc:Yield:Sigma:SigmaStatUnc:SigmaMax:SigmaMin");
627 }
628 // and the output TGraphAsymmErrors
629 if (fAsymUncertainties){
630 fgSigmaCorr = new TGraphAsymmErrors(nbins+1);
631 fgSigmaCorrExtreme = new TGraphAsymmErrors(nbins+1);
632 fgSigmaCorrConservative = new TGraphAsymmErrors(nbins+1);
633 }
634 fhStatUncEffcSigma = new TH1D("fhStatUncEffcSigma","direct charm stat unc on the cross section",nbins,limits);
635 fhStatUncEffbSigma = new TH1D("fhStatUncEffbSigma","secondary charm stat unc on the cross section",nbins,limits);
636
637
638 // protect against null denominator
639 if (deltaY==0. || fLuminosity[0]==0. || fTrigEfficiency[0]==0. || branchingRatioC==0.) {
640 AliError(" Hey you ! Why luminosity or trigger-efficiency or the c-BR or delta_y are set to zero ?! ");
641 delete [] limits;
642 delete [] binwidths;
643 return ;
644 }
645
646 Double_t value=0, errValue=0, errvalueMax=0., errvalueMin=0.;
647 Double_t errvalueExtremeMax=0., errvalueExtremeMin=0.;
648 Double_t errvalueConservativeMax=0., errvalueConservativeMin=0.;
649 Double_t errvalueStatUncEffc=0., errvalueStatUncEffb=0.;
650 for(Int_t ibin=1; ibin<=nbins; ibin++){
651
652 // Variables initialization
653 value=0.; errValue=0.; errvalueMax=0.; errvalueMin=0.;
654 errvalueExtremeMax=0.; errvalueExtremeMin=0.;
655 errvalueConservativeMax=0.; errvalueConservativeMin=0.;
656 errvalueStatUncEffc=0.; errvalueStatUncEffb=0.;
657
658 // Sigma calculation
659 // Sigma = ( 1. / (lumi * delta_y * BR_c * ParticleAntiPartFactor * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
660 value = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0. && fhRECpt->GetBinContent(ibin)>0.) ?
661 ( fhYieldCorr->GetBinContent(ibin) / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
662 : 0. ;
663
664 // Sigma statistical uncertainty:
665 // delta_sigma = sigma * sqrt ( (delta_spectra/spectra)^2 )
666 errValue = (value!=0.) ? value * (fhYieldCorr->GetBinError(ibin)/fhYieldCorr->GetBinContent(ibin)) : 0. ;
667
668 // cout<< " x "<< fhRECpt->GetBinCenter(ibin) << " sigma " << value << " +- "<< errValue << " (stat)"<<endl;
669
670 //
671 // Sigma systematic uncertainties
672 //
673 if (fAsymUncertainties && value>0.) {
674
675 // (syst but feed-down) delta_sigma = sigma * sqrt ( (delta_spectra_syst/spectra)^2 +
676 // (delta_lumi/lumi)^2 + (delta_eff_trig/eff_trig)^2 + (delta_eff/eff)^2 + (global_eff)^2 )
677 errvalueMax = value * TMath::Sqrt( (fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin))*(fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin)) +
678 (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
679 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) +
680 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
681 fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] );
682 errvalueMin = value * TMath::Sqrt( (fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin))*(fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin)) +
683 (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
684 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) +
685 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
686 fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] );
687
688 // Uncertainties from feed-down
689 // (feed-down syst) delta_sigma = sigma * sqrt ( (delta_spectra_fd/spectra_fd)^2 )
690 // extreme case
691 errvalueExtremeMax = value * (fgYieldCorrExtreme->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin));
692 errvalueExtremeMin = value * (fgYieldCorrExtreme->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin));
693 //
694 // conservative case
695 errvalueConservativeMax = value * (fgYieldCorrConservative->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin));
696 errvalueConservativeMin = value * (fgYieldCorrConservative->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin));
697
698
699 // stat unc of the efficiencies, separately
700 errvalueStatUncEffc = value * (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) ;
701 errvalueStatUncEffb = 0.;
702
703 }
704 else {
705 // protect against null denominator
706 errvalueMax = (value!=0.) ?
707 value * TMath::Sqrt( (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
708 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) +
709 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
710 fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] )
711 : 0. ;
712 errvalueMin = errvalueMax;
713 }
714
715 //
716 // Fill the histograms
717 //
718 fhSigmaCorr->SetBinContent(ibin,value);
719 fhSigmaCorr->SetBinError(ibin,errValue);
720
721 //
722 // Fill the histos and ntuple vs the Eloss hypothesis
723 //
724 if (fPbPbElossHypothesis) {
725
726 // Loop over the Eloss hypothesis
727 if(!fnHypothesis) {
728 AliError("Error reading the fc hypothesis no ntuple, please check !!");
729 delete [] limits;
730 delete [] binwidths;
731 return ;
732 }
733 Int_t entriesHypo = fnHypothesis->GetEntries();
734 Float_t pt=0., Rhy=0., fc=0., fcMin=0., fcMax=0.;
735 fnHypothesis->SetBranchAddress("pt",&pt);
736 if (fFeedDownOption==2) fnHypothesis->SetBranchAddress("Rb",&Rhy);
737 else if (fFeedDownOption==1) fnHypothesis->SetBranchAddress("Rcb",&Rhy);
738 fnHypothesis->SetBranchAddress("fc",&fc);
739 fnHypothesis->SetBranchAddress("fcMax",&fcMax);
740 fnHypothesis->SetBranchAddress("fcMin",&fcMin);
741
742 for (Int_t item=0; item<entriesHypo; item++) {
743
744 fnHypothesis->GetEntry(item);
745 if ( TMath::Abs( pt - fhDirectEffpt->GetBinCenter(ibin) ) > 0.15 ) continue;
746
747 Double_t yieldRcbvalue = (fhRECpt->GetBinContent(ibin) ) ? fhRECpt->GetBinContent(ibin) * fc : 0. ;
748 yieldRcbvalue /= fhRECpt->GetBinWidth(ibin) ;
749 Double_t yieldRcbvalueMax = (fhRECpt->GetBinContent(ibin) ) ? fhRECpt->GetBinContent(ibin) * fcMax : 0. ;
750 yieldRcbvalueMax /= fhRECpt->GetBinWidth(ibin) ;
751 Double_t yieldRcbvalueMin = (fhRECpt->GetBinContent(ibin) ) ? fhRECpt->GetBinContent(ibin) * fcMin : 0. ;
752 yieldRcbvalueMin /= fhRECpt->GetBinWidth(ibin) ;
753
754 // Sigma calculation
755 // Sigma = ( 1. / (lumi * delta_y * BR_c * ParticleAntiPartFactor * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
756 Double_t sigmaRcbvalue = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)>0.) ?
757 ( yieldRcbvalue / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
758 : 0. ;
759 Double_t sigmaRcbvalueMax = (sigmaRcbvalue!=0.) ?
760 ( yieldRcbvalueMax / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
761 : 0. ;
762 Double_t sigmaRcbvalueMin = (sigmaRcbvalue!=0.) ?
763 ( yieldRcbvalueMin / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
764 : 0. ;
765 // Sigma statistical uncertainty:
766 // delta_sigma = sigma * sqrt ( (delta_spectra/spectra)^2 ) = sigma * ( delta_spectra / (spectra-corr * binwidth) )
767 Double_t sigmaRcbvalueStatUnc = (sigmaRcbvalue!=0.) ?
768 sigmaRcbvalue * ( fhRECpt->GetBinError(ibin) / ( yieldRcbvalue * fhRECpt->GetBinWidth(ibin) ) ) : 0. ;
769
770 fhSigmaCorrRcb->Fill( fhSigmaCorr->GetBinCenter(ibin) , Rhy, sigmaRcbvalue );
771 // if(ibin==3)
772 // cout << " pt "<< fhRECpt->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<" fc-value="<< fhFcRcb->GetBinContent(ibin,rbin) <<", yield-fcRbvalue="<<yieldRcbvalue<<", sigma-fcRbvalue="<<sigmaRcbvalue<<endl;
773 fnSigma->Fill(fhRECpt->GetBinCenter(ibin), fhRECpt->GetBinContent(ibin),
774 Rhy, fc, yieldRcbvalue, sigmaRcbvalue, sigmaRcbvalueStatUnc, sigmaRcbvalueMax, sigmaRcbvalueMin );
775 }
776 }
777 //
778 // Fill the TGraphAsymmErrors
779 if (fAsymUncertainties) {
780 Double_t x = fhYieldCorr->GetBinCenter(ibin);
781 fgSigmaCorr->SetPoint(ibin,x,value); // i,x,y
782 fgSigmaCorr->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
783 fhSigmaCorrMax->SetBinContent(ibin,value+errvalueMax);
784 fhSigmaCorrMin->SetBinContent(ibin,value-errvalueMin);
785 fgSigmaCorrExtreme->SetPoint(ibin,x,value); // i,x,y
786 fgSigmaCorrExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
787 fgSigmaCorrConservative->SetPoint(ibin,x,value); // i,x,y
788 fgSigmaCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueConservativeMin,errvalueConservativeMax); // i,xl,xh,yl,yh
789
790 fhStatUncEffcSigma->SetBinContent(ibin,0.);
791 if(value>0.) fhStatUncEffcSigma->SetBinError(ibin,((errvalueStatUncEffc/value)*100.));
792 fhStatUncEffbSigma->SetBinContent(ibin,0.); fhStatUncEffbSigma->SetBinError(ibin,0.);
793 // cout << " pt "<< fhRECpt->GetBinCenter(ibin) <<" bin "<< ibin<<" stat-unc-c-sigma "<< errvalueStatUncEffc/value << endl;
794 }
795
796 }
797 delete [] binwidths;
798 delete [] limits;
799
800}
801
802//_________________________________________________________________________________________________________
803TH1D * AliHFPtSpectrum::EstimateEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco, const char *name) {
804 //
805 // Function that computes the acceptance and efficiency correction
806 // based on the simulated and reconstructed spectra
807 // and using the reconstructed spectra bin width
808 //
809 // eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
810 //
811
812 if(!fhRECpt){
813 AliInfo("Hey, the reconstructed histogram was not set yet !");
814 return NULL;
815 }
816
817 Int_t nbins = fhRECpt->GetNbinsX();
818 Double_t *limits = new Double_t[nbins+1];
819 Double_t xlow=0.,binwidth=0.;
820 for (Int_t i=1; i<=nbins; i++) {
821 binwidth = fhRECpt->GetBinWidth(i);
822 xlow = fhRECpt->GetBinLowEdge(i);
823 limits[i-1] = xlow;
824 }
825 limits[nbins] = xlow + binwidth;
826
827 TH1D * hEfficiency = new TH1D(name," acceptance #times efficiency",nbins,limits);
828
829 Double_t *sumSimu=new Double_t[nbins];
830 Double_t *sumReco=new Double_t[nbins];
831 for (Int_t ibin=0; ibin<nbins; ibin++){
832 sumSimu[ibin]=0.; sumReco[ibin]=0.;
833 }
834 for (Int_t ibin=0; ibin<=hSimu->GetNbinsX(); ibin++){
835
836 for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++){
837 if ( hSimu->GetBinCenter(ibin)>limits[ibinrec] &&
838 hSimu->GetBinCenter(ibin)<limits[ibinrec+1] ) {
839 sumSimu[ibinrec]+=hSimu->GetBinContent(ibin);
840 }
841 if ( hReco->GetBinCenter(ibin)>limits[ibinrec] &&
842 hReco->GetBinCenter(ibin)<limits[ibinrec+1] ) {
843 sumReco[ibinrec]+=hReco->GetBinContent(ibin);
844 }
845 }
846
847 }
848
849
850 // the efficiency is computed as reco/sim (in each bin)
851 // its uncertainty is err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
852 Double_t eff=0., erreff=0.;
853 for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++) {
854 if (sumSimu[ibinrec]!= 0. && sumReco[ibinrec]!=0.) {
855 eff = sumReco[ibinrec] / sumSimu[ibinrec] ;
856 // protection in case eff > 1.0
857 // test calculation (make the argument of the sqrt positive)
858 erreff = TMath::Sqrt( eff * TMath::Abs(1.0 - eff) ) / TMath::Sqrt( sumSimu[ibinrec] );
859 }
860 else { eff=0.0; erreff=0.; }
861 hEfficiency->SetBinContent(ibinrec+1,eff);
862 hEfficiency->SetBinError(ibinrec+1,erreff);
863 }
864
865 delete [] sumSimu;
866 delete [] sumReco;
867
868 return (TH1D*)hEfficiency;
869}
870
871//_________________________________________________________________________________________________________
872void AliHFPtSpectrum::EstimateAndSetDirectEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco) {
873 //
874 // Function that computes the Direct acceptance and efficiency correction
875 // based on the simulated and reconstructed spectra
876 // and using the reconstructed spectra bin width
877 //
878 // eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
879 //
880
881 if(!fhRECpt || !hSimu || !hReco){
882 AliError("Hey, the reconstructed histogram was not set yet !");
883 return;
884 }
885
886 fhDirectEffpt = EstimateEfficiencyRecoBin(hSimu,hReco,"fhDirectEffpt");
887 fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance #times efficiency");
888
889}
890
891//_________________________________________________________________________________________________________
892void AliHFPtSpectrum::EstimateAndSetFeedDownEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco) {
893 //
894 // Function that computes the Feed-Down acceptance and efficiency correction
895 // based on the simulated and reconstructed spectra
896 // and using the reconstructed spectra bin width
897 //
898 // eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
899 //
900
901 if(!fhRECpt || !hSimu || !hReco){
902 AliError("Hey, the reconstructed histogram was not set yet !");
903 return;
904 }
905
906 fhFeedDownEffpt = EstimateEfficiencyRecoBin(hSimu,hReco,"fhFeedDownEffpt");
907 fhFeedDownEffpt->SetNameTitle("fhFeedDownEffpt"," feed-down acceptance #times efficiency");
908
909}
910
911//_________________________________________________________________________________________________________
912Bool_t AliHFPtSpectrum::Initialize(){
913 //
914 // Initialization of the variables (histograms)
915 //
916
917 if (fFeedDownOption==0) {
918 AliInfo("Getting ready for the corrections without feed-down consideration");
919 } else if (fFeedDownOption==1) {
920 AliInfo("Getting ready for the fc feed-down correction calculation");
921 } else if (fFeedDownOption==2) {
922 AliInfo("Getting ready for the Nb feed-down correction calculation");
923 } else { AliError("The calculation option must be <=2"); return kFALSE; }
924
925 // Start checking the input histograms consistency
926 Bool_t areconsistent=kTRUE;
927
928 // General checks
929 if (!fhDirectEffpt || !fhRECpt) {
930 AliError(" Reconstructed spectra and/or the Nc efficiency distributions are not defined");
931 return kFALSE;
932 }
933 areconsistent &= CheckHistosConsistency(fhRECpt,fhDirectEffpt);
934 if (!areconsistent) {
935 AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
936 return kFALSE;
937 }
938 if (fFeedDownOption==0) return kTRUE;
939
940 //
941 // Common checks for options 1 (fc) & 2(Nb)
942 if (!fhFeedDownMCpt || !fhFeedDownEffpt) {
943 AliError(" Theoretical Nb and/or the Nb efficiency distributions are not defined");
944 return kFALSE;
945 }
946 areconsistent &= CheckHistosConsistency(fhRECpt,fhFeedDownMCpt);
947 areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownEffpt);
948 if (fAsymUncertainties) {
949 if (!fhFeedDownMCptMax || !fhFeedDownMCptMin) {
950 AliError(" Max/Min theoretical Nb distributions are not defined");
951 return kFALSE;
952 }
953 areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownMCptMax);
954 }
955 if (!areconsistent) {
956 AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
957 return kFALSE;
958 }
959 if (fFeedDownOption>1) return kTRUE;
960
961 //
962 // Now checks for option 1 (fc correction)
963 if (!fhDirectMCpt) {
964 AliError("Theoretical Nc distributions is not defined");
965 return kFALSE;
966 }
967 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhFeedDownMCpt);
968 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectEffpt);
969 if (fAsymUncertainties) {
970 if (!fhDirectMCptMax || !fhDirectMCptMin) {
971 AliError(" Max/Min theoretical Nc distributions are not defined");
972 return kFALSE;
973 }
974 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectMCptMax);
975 }
976 if (!areconsistent) {
977 AliInfo("Histograms required for fc correction are not consistent (bin width, bounds)");
978 return kFALSE;
979 }
980
981 return kTRUE;
982}
983
984//_________________________________________________________________________________________________________
985Bool_t AliHFPtSpectrum::CheckHistosConsistency(TH1D *h1, TH1D *h2){
986 //
987 // Check the histograms consistency (bins, limits)
988 //
989
990 if (!h1 || !h2) {
991 AliError("One or both histograms don't exist");
992 return kFALSE;
993 }
994
995 Double_t binwidth1 = h1->GetBinWidth(1);
996 Double_t binwidth2 = h2->GetBinWidth(1);
997 Double_t min1 = h1->GetBinCenter(1) - (binwidth1/2.) ;
998// Double_t max1 = h1->GetBinCenter(nbins1) + (binwidth1/2.) ;
999 Double_t min2 = h2->GetBinCenter(1) - (binwidth2/2.) ;
1000// Double_t max2 = h2->GetBinCenter(nbins2) + (binwidth2/2.) ;
1001
1002 if (binwidth1!=binwidth2) {
1003 AliInfo(" histograms with different bin width");
1004 return kFALSE;
1005 }
1006 if (min1!=min2) {
1007 AliInfo(" histograms with different minimum");
1008 return kFALSE;
1009 }
1010// if (max1!=max2) {
1011// AliInfo(" histograms with different maximum");
1012// return kFALSE;
1013// }
1014
1015 return kTRUE;
1016}
1017
1018//_________________________________________________________________________________________________________
1019void AliHFPtSpectrum::CalculateFeedDownCorrectionFc(){
1020 //
1021 // Compute fc factor and its uncertainties bin by bin
1022 // fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
1023 //
1024 // uncertainties: (conservative) combine the upper/lower N_b & N_c predictions together
1025 // (extreme) combine the upper N_b predictions with the lower N_c predictions & viceversa
1026 // systematic uncertainty on the acceptance x efficiency b/c ratio are included
1027 //
1028 // In addition, in HIC the feed-down correction varies with an energy loss hypothesis: Raa(c-->D) / Raa(b-->D) = Rcb
1029 // fc (Rcb) = ( 1. / ( 1 + (eff_b/eff_c)*(N_b/N_c)* (1/Rcb) ) );
1030 //
1031 AliInfo("Calculating the feed-down correction factor (fc method)");
1032
1033 // define the variables
1034 Int_t nbins = fhRECpt->GetNbinsX();
1035 Double_t binwidth = fhRECpt->GetBinWidth(1);
1036 Double_t *limits = new Double_t[nbins+1];
1037 Double_t *binwidths = new Double_t[nbins];
1038 Double_t xlow=0.;
1039 for (Int_t i=1; i<=nbins; i++) {
1040 binwidth = fhRECpt->GetBinWidth(i);
1041 xlow = fhRECpt->GetBinLowEdge(i);
1042 limits[i-1] = xlow;
1043 binwidths[i-1] = binwidth;
1044 }
1045 limits[nbins] = xlow + binwidth;
1046
1047 Double_t correction=1.;
1048 Double_t theoryRatio=1.;
1049 Double_t effRatio=1.;
1050 Double_t correctionExtremeA=1., correctionExtremeB=1.;
1051 Double_t theoryRatioExtremeA=1., theoryRatioExtremeB=1.;
1052 Double_t correctionConservativeA=1., correctionConservativeB=1.;
1053 Double_t theoryRatioConservativeA=1., theoryRatioConservativeB=1.;
1054 Double_t correctionUnc=0.;
1055 Double_t correctionExtremeAUnc=0., correctionExtremeBUnc=0.;
1056 Double_t correctionConservativeAUnc=0., correctionConservativeBUnc=0.;
1057
1058 // declare the output histograms
1059 fhFc = new TH1D("fhFc","fc correction factor",nbins,limits);
1060 fhFcMax = new TH1D("fhFcMax","max fc correction factor",nbins,limits);
1061 fhFcMin = new TH1D("fhFcMin","min fc correction factor",nbins,limits);
1062 if(fPbPbElossHypothesis) {
1063 fhFcRcb = new TH2D("fhFcRcb","fc correction factor vs Rcb Eloss hypothesis; p_{T} [GeV/c] ; Rcb Eloss hypothesis ; fc correction",nbins,limits,800,0.,4.);
1064 fnHypothesis = new TNtuple("fnHypothesis"," Feed-down correction vs hypothesis (fc)","pt:Rcb:fc:fcMax:fcMin");
1065 }
1066 // two local control histograms
1067 TH1D *hTheoryRatio = new TH1D("hTheoryRatio","Theoretical B-->D over c-->D (feed-down/direct) ratio",nbins,limits);
1068 TH1D *hEffRatio = new TH1D("hEffRatio","Efficiency B-->D over c-->D (feed-down/direct) ratio",nbins,limits);
1069 // and the output TGraphAsymmErrors
1070 if (fAsymUncertainties) {
1071 fgFcExtreme = new TGraphAsymmErrors(nbins+1);
1072 fgFcExtreme->SetNameTitle("fgFcExtreme","fgFcExtreme");
1073 fgFcConservative = new TGraphAsymmErrors(nbins+1);
1074 fgFcConservative->SetNameTitle("fgFcConservative","fgFcConservative");
1075 }
1076
1077 fhStatUncEffcFD = new TH1D("fhStatUncEffcFD","direct charm stat unc on the feed-down correction",nbins,limits);
1078 fhStatUncEffbFD = new TH1D("fhStatUncEffbFD","secondary charm stat unc on the feed-down correction",nbins,limits);
1079 Double_t correctionConservativeAUncStatEffc=0., correctionConservativeBUncStatEffc=0.;
1080 Double_t correctionConservativeAUncStatEffb=0., correctionConservativeBUncStatEffb=0.;
1081
1082 //
1083 // Compute fc
1084 //
1085 for (Int_t ibin=1; ibin<=nbins; ibin++) {
1086
1087 // theory_ratio = (N_b/N_c)
1088 theoryRatio = (fhDirectMCpt->GetBinContent(ibin)>0. && fhFeedDownMCpt->GetBinContent(ibin)>0.) ?
1089 fhFeedDownMCpt->GetBinContent(ibin) / fhDirectMCpt->GetBinContent(ibin) : 1.0 ;
1090
1091 //
1092 // Calculate the uncertainty [ considering only the theoretical uncertainties on Nb & Nc for now !!! ]
1093 //
1094 // extreme A = direct-max, feed-down-min
1095 theoryRatioExtremeA = (fhDirectMCptMax->GetBinContent(ibin)>0. && fhFeedDownMCptMin->GetBinContent(ibin)>0.) ?
1096 fhFeedDownMCptMin->GetBinContent(ibin) / fhDirectMCptMax->GetBinContent(ibin) : 1.0 ;
1097 // extreme B = direct-min, feed-down-max
1098 theoryRatioExtremeB = (fhDirectMCptMin->GetBinContent(ibin)>0. && fhDirectMCptMax->GetBinContent(ibin)>0.) ?
1099 fhFeedDownMCptMax->GetBinContent(ibin) / fhDirectMCptMin->GetBinContent(ibin) : 1.0 ;
1100 // conservative A = direct-max, feed-down-max
1101 theoryRatioConservativeA = (fhDirectMCptMax->GetBinContent(ibin)>0. && fhFeedDownMCptMin->GetBinContent(ibin)>0.) ?
1102 fhFeedDownMCptMax->GetBinContent(ibin) / fhDirectMCptMax->GetBinContent(ibin) : 1.0 ;
1103 // conservative B = direct-min, feed-down-min
1104 theoryRatioConservativeB = (fhDirectMCptMin->GetBinContent(ibin)>0. && fhDirectMCptMax->GetBinContent(ibin)>0.) ?
1105 fhFeedDownMCptMin->GetBinContent(ibin) / fhDirectMCptMin->GetBinContent(ibin) : 1.0 ;
1106
1107 // eff_ratio = (eff_b/eff_c)
1108 effRatio = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ?
1109 fhFeedDownEffpt->GetBinContent(ibin) / fhDirectEffpt->GetBinContent(ibin) : 1.0 ;
1110
1111 // fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
1112 if( TMath::Abs(effRatio - 1.0)<0.01 || TMath::Abs(theoryRatio - 1.0)<0.01 ) {
1113 correction = 1.0;
1114 correctionExtremeA = 1.0;
1115 correctionExtremeB = 1.0;
1116 correctionConservativeA = 1.0;
1117 correctionConservativeB = 1.0;
1118 }
1119 else {
1120 correction = ( 1. / ( 1 + ( effRatio * theoryRatio ) ) );
1121 correctionExtremeA = ( 1. / ( 1 + ( effRatio * theoryRatioExtremeA ) ) );
1122 correctionExtremeB = ( 1. / ( 1 + ( effRatio * theoryRatioExtremeB ) ) );
1123 correctionConservativeA = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeA ) ) );
1124 correctionConservativeB = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeB ) ) );
1125 }
1126
1127
1128 // fc uncertainty from (eff_b/eff_c) = fc^2 * (N_b/N_c) * delta(eff_b/eff_c)
1129 // delta(eff_b/eff_c) is a percentage = effRatio * sqrt( fGlobalEfficiencyUncertainties[1]^2 + unc_eff_c ^2 + unc_eff_b ^2 )
1130 Double_t relEffUnc = TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] +
1131 (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
1132 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))
1133 );
1134
1135 correctionUnc = correction*correction * theoryRatio * effRatio * relEffUnc;
1136 correctionExtremeAUnc = correctionExtremeA*correctionExtremeA * theoryRatioExtremeA * effRatio * relEffUnc;
1137 correctionExtremeBUnc = correctionExtremeB*correctionExtremeB * theoryRatioExtremeB * effRatio * relEffUnc;
1138
1139 correctionConservativeAUnc = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA *effRatio * relEffUnc;
1140 //
1141 correctionConservativeAUncStatEffc = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA *effRatio *
1142 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin));
1143 correctionConservativeAUncStatEffb = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA *effRatio *
1144 (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin));
1145
1146 correctionConservativeBUnc = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB *effRatio * relEffUnc;
1147
1148 correctionConservativeBUncStatEffb = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB *effRatio *
1149 (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin));
1150 correctionConservativeBUncStatEffc = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB *effRatio *
1151 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin));
1152
1153
1154 // Fill in the histograms
1155 hTheoryRatio->SetBinContent(ibin,theoryRatio);
1156 hEffRatio->SetBinContent(ibin,effRatio);
1157 fhFc->SetBinContent(ibin,correction);
1158 //
1159 // Estimate how the result varies vs charm/beauty Eloss hypothesis
1160 //
1161 if ( TMath::Abs(correction-1.0)>0.0001 && fPbPbElossHypothesis){
1162 // Loop over the Eloss hypothesis
1163 // Int_t rbin=0;
1164 for (Float_t rval=0.0025; rval<4.0; rval+=0.005){
1165 // Central fc with Eloss hypothesis
1166 Double_t correctionRcb = ( 1. / ( 1 + ( effRatio * theoryRatio * (1/rval) ) ) );
1167 fhFcRcb->Fill( fhFc->GetBinCenter(ibin) , rval, correctionRcb );
1168 // if(ibin==3){
1169 // cout << " pt "<< fhFc->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<", fc-Rcb-value="<<correctionRcb<<endl;
1170 // rbin++;
1171 // }
1172 // Upper / lower fc with up / low FONLL bands and Eloss hypothesis
1173 Double_t correctionConservativeARcb = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeA * (1/rval) ) ) );
1174 Double_t correctionConservativeBRcb = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeB * (1/rval) ) ) );
1175 Double_t correctionConservativeARcbUnc = correctionConservativeARcb*correctionConservativeARcb * theoryRatioConservativeA * (1/rval) *effRatio * relEffUnc;
1176 Double_t correctionConservativeBRcbUnc = correctionConservativeBRcb*correctionConservativeBRcb * theoryRatioConservativeB * (1/rval) *effRatio * relEffUnc;
1177 Double_t consvalRcb[4] = { correctionConservativeARcb - correctionConservativeARcbUnc, correctionConservativeARcb + correctionConservativeARcbUnc,
1178 correctionConservativeBRcb - correctionConservativeBRcbUnc, correctionConservativeBRcb + correctionConservativeBRcbUnc};
1179 Double_t uncConservativeRcbMin = correctionRcb - TMath::MinElement(4,consvalRcb);
1180 Double_t uncConservativeRcbMax = TMath::MaxElement(4,consvalRcb) - correctionRcb;
1181// if(ibin==3)
1182// cout << " pt="<<fhDirectEffpt->GetBinCenter(ibin)<<", hypo="<<rval<<", fc="<<correctionRcb<<" +"<<uncConservativeRcbMax<<" -"<<uncConservativeRcbMin<<endl;
1183 fnHypothesis->Fill( fhDirectEffpt->GetBinCenter(ibin), rval, correctionRcb, correctionRcb+uncConservativeRcbMax, correctionRcb-uncConservativeRcbMin);
1184 }
1185 }
1186 //
1187 // Fill the rest of (asymmetric) histograms
1188 //
1189 if (fAsymUncertainties) {
1190 Double_t x = fhDirectMCpt->GetBinCenter(ibin);
1191 Double_t val[4] = { correctionExtremeA + correctionExtremeAUnc, correctionExtremeA - correctionExtremeAUnc,
1192 correctionExtremeB + correctionExtremeBUnc, correctionExtremeB - correctionExtremeBUnc };
1193 Double_t uncExtremeMin = correction - TMath::MinElement(4,val);
1194 Double_t uncExtremeMax = TMath::MaxElement(4,val) - correction;
1195 fgFcExtreme->SetPoint(ibin,x,correction); // i,x,y
1196 fgFcExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),uncExtremeMin,uncExtremeMax); // i,xl,xh,yl,yh
1197 fhFcMax->SetBinContent(ibin,correction+uncExtremeMax);
1198 fhFcMin->SetBinContent(ibin,correction-uncExtremeMin);
1199 Double_t consval[4] = { correctionConservativeA - correctionConservativeAUnc, correctionConservativeA + correctionConservativeAUnc,
1200 correctionConservativeB - correctionConservativeBUnc, correctionConservativeB + correctionConservativeBUnc};
1201 Double_t uncConservativeMin = correction - TMath::MinElement(4,consval);
1202 Double_t uncConservativeMax = TMath::MaxElement(4,consval) - correction;
1203 fgFcConservative->SetPoint(ibin,x,correction); // i,x,y
1204 fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),uncConservativeMin,uncConservativeMax); // i,xl,xh,yl,yh
1205 if( !(correction>0.) ){
1206 fgFcExtreme->SetPoint(ibin,x,0.); // i,x,y
1207 fgFcExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.); // i,xl,xh,yl,yh
1208 fgFcConservative->SetPoint(ibin,x,0.); // i,x,y
1209 fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.); // i,xl,xh,yl,yh
1210 }
1211
1212 Double_t valStatEffc[2] = { correctionConservativeAUncStatEffc/correctionConservativeA,
1213 correctionConservativeBUncStatEffc/correctionConservativeB };
1214 Double_t valStatEffb[2] = { correctionConservativeAUncStatEffb/correctionConservativeA,
1215 correctionConservativeBUncStatEffb/correctionConservativeB };
1216 Double_t uncConservativeStatEffc = TMath::MaxElement(2,valStatEffc);
1217 Double_t uncConservativeStatEffb = TMath::MaxElement(2,valStatEffb);
1218 fhStatUncEffcFD->SetBinContent(ibin,0.); fhStatUncEffcFD->SetBinError(ibin,uncConservativeStatEffc*100.);
1219 fhStatUncEffbFD->SetBinContent(ibin,0.); fhStatUncEffbFD->SetBinError(ibin,uncConservativeStatEffb*100.);
1220 // cout << " pt "<< fhStatUncEffcFD->GetBinCenter(ibin) <<" bin "<< ibin<<" fc-stat-c ="<<uncConservativeStatEffc<<" fc-stat-b ="<<uncConservativeStatEffb<<endl;
1221 }
1222
1223 }
1224 delete [] binwidths;
1225 delete [] limits;
1226
1227}
1228
1229//_________________________________________________________________________________________________________
1230void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumFc(){
1231 //
1232 // Compute the feed-down corrected spectrum if feed-down correction is done via fc factor (bin by bin)
1233 // physics = reco * fc / bin-width
1234 //
1235 // uncertainty: (stat) delta_physics = physics * sqrt ( (delta_reco/reco)^2 )
1236 // (syst but feed-down) delta_physics = physics * sqrt ( (delta_reco_syst/reco)^2 )
1237 // (feed-down syst) delta_physics = physics * sqrt ( (delta_fc/fc)^2 )
1238 //
1239 // ( Calculation done bin by bin )
1240 //
1241 // In addition, in HIC the feed-down correction varies with an energy loss hypothesis: Raa(c-->D) / Raa(b-->D) = Rcb
1242
1243 AliInfo(" Calculating the feed-down corrected spectrum (fc method)");
1244
1245 if (!fhFc || !fhRECpt) {
1246 AliError(" Reconstructed or fc distributions are not defined");
1247 return;
1248 }
1249
1250 Int_t nbins = fhRECpt->GetNbinsX();
1251 Double_t value = 0., errvalue = 0., errvalueMax= 0., errvalueMin= 0.;
1252 Double_t valueExtremeMax= 0., valueExtremeMin= 0.;
1253 Double_t valueConservativeMax= 0., valueConservativeMin= 0.;
1254 Double_t binwidth = fhRECpt->GetBinWidth(1);
1255 Double_t *limits = new Double_t[nbins+1];
1256 Double_t *binwidths = new Double_t[nbins];
1257 Double_t xlow=0.;
1258 for (Int_t i=1; i<=nbins; i++) {
1259 binwidth = fhRECpt->GetBinWidth(i);
1260 xlow = fhRECpt->GetBinLowEdge(i);
1261 limits[i-1] = xlow;
1262 binwidths[i-1] = binwidth;
1263 }
1264 limits[nbins] = xlow + binwidth;
1265
1266 // declare the output histograms
1267 fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by fc)",nbins,limits);
1268 fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by fc)",nbins,limits);
1269 fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by fc)",nbins,limits);
1270 if(fPbPbElossHypothesis) fhYieldCorrRcb = new TH2D("fhYieldCorrRcb","corrected yield (by fc) vs Rcb Eloss hypothesis; p_{T} [GeV/c] ; Rcb Eloss hypothesis ; corrected yield",nbins,limits,800,0.,4.);
1271 // and the output TGraphAsymmErrors
1272 if (fAsymUncertainties){
1273 fgYieldCorr = new TGraphAsymmErrors(nbins+1);
1274 fgYieldCorrExtreme = new TGraphAsymmErrors(nbins+1);
1275 fgYieldCorrConservative = new TGraphAsymmErrors(nbins+1);
1276 }
1277
1278 //
1279 // Do the calculation
1280 //
1281 for (Int_t ibin=1; ibin<=nbins; ibin++) {
1282
1283 // calculate the value
1284 // physics = reco * fc / bin-width
1285 value = (fhRECpt->GetBinContent(ibin) && fhFc->GetBinContent(ibin)) ?
1286 fhRECpt->GetBinContent(ibin) * fhFc->GetBinContent(ibin) : 0. ;
1287 value /= fhRECpt->GetBinWidth(ibin) ;
1288
1289 // Statistical uncertainty
1290 // (stat) delta_physics = physics * sqrt ( (delta_reco/reco)^2 )
1291 errvalue = (value!=0. && fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) ?
1292 value * (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) : 0. ;
1293
1294 // Calculate the systematic uncertainties
1295 // (syst but feed-down) delta_physics = physics * sqrt ( (delta_reco_syst/reco)^2 )
1296 // (feed-down syst) delta_physics = physics * sqrt ( (delta_fc/fc)^2 )
1297 //
1298 // Protect against null denominator. If so, define uncertainty as null
1299 if (fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) {
1300
1301 if (fAsymUncertainties) {
1302
1303 // Systematics but feed-down
1304 if (fgRECSystematics) {
1305 errvalueMax = value * ( fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinContent(ibin) );
1306 errvalueMin = value * ( fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinContent(ibin) );
1307 }
1308 else { errvalueMax = 0.; errvalueMin = 0.; }
1309
1310 // Extreme feed-down systematics
1311 valueExtremeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcExtreme->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1312 valueExtremeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcExtreme->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1313
1314 // Conservative feed-down systematics
1315 valueConservativeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcConservative->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1316 valueConservativeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcConservative->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1317
1318 }
1319
1320 }
1321 else { errvalueMax = 0.; errvalueMin = 0.; }
1322
1323 //
1324 // Fill in the histograms
1325 //
1326 fhYieldCorr->SetBinContent(ibin,value);
1327 fhYieldCorr->SetBinError(ibin,errvalue);
1328 //
1329 // Fill the histos and ntuple vs the Eloss hypothesis
1330 //
1331 if (fPbPbElossHypothesis) {
1332 // Loop over the Eloss hypothesis
1333 for (Float_t rval=0.0025; rval<4.0; rval+=0.005){
1334 Int_t rbin = FindTH2YBin(fhYieldCorrRcb,rval);
1335 Double_t fcRcbvalue = fhFcRcb->GetBinContent(ibin,rbin);
1336 // physics = reco * fcRcb / bin-width
1337 Double_t Rcbvalue = (fhRECpt->GetBinContent(ibin) && fcRcbvalue) ?
1338 fhRECpt->GetBinContent(ibin) * fcRcbvalue : 0. ;
1339 Rcbvalue /= fhRECpt->GetBinWidth(ibin) ;
1340 fhYieldCorrRcb->Fill( fhYieldCorr->GetBinCenter(ibin) , rval, Rcbvalue );
1341 // cout << " pt "<< fhRECpt->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<" fc-fcRbvalue="<<fcRcbvalue<<", yield="<<Rcbvalue<<endl;
1342 }
1343 }
1344 if (fAsymUncertainties) {
1345 Double_t center = fhYieldCorr->GetBinCenter(ibin);
1346 fgYieldCorr->SetPoint(ibin,center,value); // i,x,y
1347 fgYieldCorr->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
1348 fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax);
1349 fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
1350 fgYieldCorrExtreme->SetPoint(ibin,center,value);
1351 fgYieldCorrExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),value-valueExtremeMin,valueExtremeMax-value);
1352 fgYieldCorrConservative->SetPoint(ibin,center,value);
1353 fgYieldCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),value-valueConservativeMin,valueConservativeMax-value);
1354 }
1355
1356 }
1357 delete [] binwidths;
1358 delete [] limits;
1359
1360}
1361
1362//_________________________________________________________________________________________________________
1363void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumNb(Double_t deltaY, Double_t branchingRatioBintoFinalDecay) {
1364 //
1365 // Compute the feed-down corrected spectrum if feed-down correction is done via Nb (bin by bin)
1366 // physics = [ reco - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) ] / bin-width
1367 //
1368 // uncertainty: (stat) delta_physics = sqrt ( (delta_reco)^2 ) / bin-width
1369 // (syst but feed-down) delta_physics = sqrt ( (delta_reco_syst)^2 ) / bin-width
1370 // (feed-down syst) delta_physics = sqrt ( (k*delta_lumi/lumi)^2 + (k*delta_eff_trig/eff_trig)^2
1371 // + (k*delta_Nb/Nb)^2 + (k*delta_eff/eff)^2 + (k*global_eff_ratio)^2 ) / bin-width
1372 // where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th
1373 //
1374 // In addition, in HIC the feed-down correction varies with an energy loss hypothesis: Raa(b-->D) = Rb
1375 // physics = [ reco - ( Tab * Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th * Rb ) ] / bin-width
1376 //
1377 AliInfo("Calculating the feed-down correction factor and spectrum (Nb method)");
1378
1379 Int_t nbins = fhRECpt->GetNbinsX();
1380 Double_t binwidth = fhRECpt->GetBinWidth(1);
1381 Double_t value = 0., errvalue = 0., errvalueMax = 0., errvalueMin = 0., kfactor = 0.;
1382 Double_t errvalueExtremeMax = 0., errvalueExtremeMin = 0.;
1383 Double_t *limits = new Double_t[nbins+1];
1384 Double_t *binwidths = new Double_t[nbins];
1385 Double_t xlow=0.;
1386 for (Int_t i=1; i<=nbins; i++) {
1387 binwidth = fhRECpt->GetBinWidth(i);
1388 xlow = fhRECpt->GetBinLowEdge(i);
1389 limits[i-1] = xlow;
1390 binwidths[i-1] = binwidth;
1391 }
1392 limits[nbins] = xlow + binwidth;
1393
1394 // declare the output histograms
1395 fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by Nb)",nbins,limits);
1396 fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by Nb)",nbins,limits);
1397 fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by Nb)",nbins,limits);
1398 if(fPbPbElossHypothesis) {
1399 fhFcRcb = new TH2D("fhFcRcb","fc correction factor (Nb method) vs Rb Eloss hypothesis; p_{T} [GeV/c] ; Rb Eloss hypothesis ; fc correction",nbins,limits,800,0.,4.);
1400 fhYieldCorrRcb = new TH2D("fhYieldCorrRcb","corrected yield (by Nb) vs Rb Eloss hypothesis; p_{T} [GeV/c] ; Rb Eloss hypothesis ; corrected yield",nbins,limits,800,0.,4.);
1401 fnHypothesis = new TNtuple("fnHypothesis"," Feed-down correction vs hypothesis (Nb)","pt:Rb:fc:fcMax:fcMin");
1402 }
1403 // and the output TGraphAsymmErrors
1404 if (fAsymUncertainties){
1405 fgYieldCorr = new TGraphAsymmErrors(nbins+1);
1406 fgYieldCorrExtreme = new TGraphAsymmErrors(nbins+1);
1407 fgYieldCorrConservative = new TGraphAsymmErrors(nbins+1);
1408 // Define fc-conservative
1409 fgFcConservative = new TGraphAsymmErrors(nbins+1);
1410 AliInfo(" Beware the conservative & extreme uncertainties are equal by definition !");
1411 }
1412
1413 // variables to define fc-conservative
1414 double correction=0, correctionMax=0., correctionMin=0.;
1415
1416 fhStatUncEffcFD = new TH1D("fhStatUncEffcFD","direct charm stat unc on the feed-down correction",nbins,limits);
1417 fhStatUncEffbFD = new TH1D("fhStatUncEffbFD","secondary charm stat unc on the feed-down correction",nbins,limits);
1418 Double_t correctionUncStatEffc=0.;
1419 Double_t correctionUncStatEffb=0.;
1420
1421
1422 //
1423 // Do the calculation
1424 //
1425 for (Int_t ibin=1; ibin<=nbins; ibin++) {
1426
1427 // Calculate the value
1428 // physics = [ reco - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) ] / bin-width
1429 // In HIC : physics = [ reco - ( Tab * Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th * Rb ) ] / bin-width
1430 //
1431 //
1432 Double_t frac = 1.0, errfrac =0.;
1433 if(fPbPbElossHypothesis) {
1434 frac = fTab[0]*fNevts;
1435 errfrac = frac * TMath::Sqrt( (fTab[1]/fTab[0])*(fTab[1]/fTab[0]) + (1/fNevts) );
1436 } else {
1437 frac = fLuminosity[0];
1438 errfrac = fLuminosity[1];
1439 }
1440
1441 value = ( fhRECpt->GetBinContent(ibin)>0. && fhRECpt->GetBinContent(ibin)!=0. &&
1442 fhFeedDownMCpt->GetBinContent(ibin)>0. && fhFeedDownEffpt->GetBinContent(ibin)>0. ) ?
1443 fhRECpt->GetBinContent(ibin) - frac*(deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) * fhRECpt->GetBinWidth(ibin) )
1444 : 0. ;
1445 value /= fhRECpt->GetBinWidth(ibin);
1446 if (value<0.) value =0.;
1447
1448 // Statistical uncertainty: delta_physics = sqrt ( (delta_reco)^2 ) / bin-width
1449 errvalue = (value!=0. && fhRECpt->GetBinError(ibin) && fhRECpt->GetBinError(ibin)!=0.) ?
1450 fhRECpt->GetBinError(ibin) : 0.;
1451 errvalue /= fhRECpt->GetBinWidth(ibin);
1452
1453 // Correction (fc) : Estimate of the relative amount feed-down subtracted
1454 // correction = [ 1 - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th)/reco ]
1455 // in HIC: correction = [ 1 - ( Tab * Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th)/reco ]
1456 correction = (value>0.) ?
1457 1 - (frac*deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) * fhRECpt->GetBinWidth(ibin) ) / fhRECpt->GetBinContent(ibin) : 0. ;
1458 if (correction<0.) correction = 0.;
1459
1460 // Systematic uncertainties
1461 // (syst but feed-down) delta_physics = sqrt ( (delta_reco_syst)^2 ) / bin-width
1462 // (feed-down syst) delta_physics = sqrt ( (k*delta_lumi/lumi)^2 + (k*delta_eff_trig/eff_trig)^2
1463 // + (k*delta_Nb/Nb)^2 + (k*delta_eff/eff)^2 + (k*global_eff_ratio)^2 ) / bin-width
1464 // where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th * bin-width
1465 kfactor = frac*deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) * fhRECpt->GetBinWidth(ibin) ;
1466 //
1467 if (fAsymUncertainties && value>0.) {
1468 Double_t nb = fhFeedDownMCpt->GetBinContent(ibin);
1469 Double_t nbDmax = fhFeedDownMCptMax->GetBinContent(ibin) - fhFeedDownMCpt->GetBinContent(ibin);
1470 Double_t nbDmin = fhFeedDownMCpt->GetBinContent(ibin) - fhFeedDownMCptMin->GetBinContent(ibin);
1471
1472 // Systematics but feed-down
1473 if (fgRECSystematics){
1474 errvalueMax = fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinWidth(ibin) ;
1475 errvalueMin = fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinWidth(ibin);
1476 }
1477 else { errvalueMax = 0.; errvalueMin = 0.; }
1478
1479 // Feed-down systematics
1480 // min value with the maximum Nb
1481 Double_t errCom = ( (kfactor*errfrac/frac)*(kfactor*errfrac/frac) ) +
1482 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1483 ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
1484 ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) ) ;
1485 errvalueExtremeMin = TMath::Sqrt( errCom + ( (kfactor*nbDmax/nb)*(kfactor*nbDmax/nb) ) ) / fhRECpt->GetBinWidth(ibin);
1486 // max value with the minimum Nb
1487 errvalueExtremeMax = TMath::Sqrt( errCom + ( (kfactor*nbDmin/nb)*(kfactor*nbDmin/nb) ) ) / fhRECpt->GetBinWidth(ibin);
1488
1489 // Correction systematics (fc)
1490 // min value with the maximum Nb
1491 correctionMin = TMath::Sqrt( errCom + ( (kfactor*nbDmax/nb)*(kfactor*nbDmax/nb) ) ) / fhRECpt->GetBinContent(ibin) ;
1492 // max value with the minimum Nb
1493 correctionMax = TMath::Sqrt( errCom + ( (kfactor*nbDmin/nb)*(kfactor*nbDmin/nb) ) ) / fhRECpt->GetBinContent(ibin) ;
1494 //
1495 correctionUncStatEffb = TMath::Sqrt( ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) )
1496 ) / fhRECpt->GetBinContent(ibin) ;
1497 correctionUncStatEffc = 0.;
1498 }
1499 else{ // Don't consider Nb uncertainty in this case [ to be tested!!! ]
1500 errvalueExtremeMax = TMath::Sqrt( ( (kfactor*errfrac/frac)*(kfactor*errfrac/frac) ) +
1501 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1502 ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
1503 ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1504 ) / fhRECpt->GetBinWidth(ibin);
1505 errvalueExtremeMin = errvalueExtremeMax ;
1506 }
1507
1508
1509 // fill in histograms
1510 fhYieldCorr->SetBinContent(ibin,value);
1511 fhYieldCorr->SetBinError(ibin,errvalue);
1512 //
1513 // Estimate how the result varies vs charm/beauty Eloss hypothesis
1514 //
1515 if ( correction>0.0001 && fPbPbElossHypothesis){
1516 // Loop over the Eloss hypothesis
1517 // Int_t rbin=0;
1518 for (Float_t rval=0.0025; rval<4.0; rval+=0.005){
1519 // correction = [ 1 - (Tab *Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th *binwidth )* (rval) /reco ]
1520 Double_t fcRcbvalue = 1 - (fTab[0]*fNevts*deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin)*fhRECpt->GetBinWidth(ibin) * rval ) / fhRECpt->GetBinContent(ibin) ;
1521 if(fcRcbvalue<0.) fcRcbvalue=0.;
1522 fhFcRcb->Fill( fhRECpt->GetBinCenter(ibin) , rval, fcRcbvalue );
1523 // physics = reco * fcRcb / bin-width
1524 Double_t Rcbvalue = (fhRECpt->GetBinContent(ibin) && fcRcbvalue) ?
1525 fhRECpt->GetBinContent(ibin) * fcRcbvalue : 0. ;
1526 Rcbvalue /= fhRECpt->GetBinWidth(ibin) ;
1527 fhYieldCorrRcb->Fill( fhYieldCorr->GetBinCenter(ibin) , rval, Rcbvalue );
1528 // if(ibin==3){
1529 // cout << " pt "<< fhFcRcb->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<", fc-Rb-value="<< fcRcbvalue << ", yield-Rb-value="<< Rcbvalue <<endl;
1530 // cout << " pt "<< fhFcRcb->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", fc-Rb-value="<< fcRcbvalue << ", yield-Rb-value="<< Rcbvalue <<endl;
1531 // rbin++;
1532 // }
1533 Double_t correctionMaxRcb = correctionMax*rval;
1534 Double_t correctionMinRcb = correctionMin*rval;
1535 fnHypothesis->Fill( fhYieldCorr->GetBinCenter(ibin), rval, fcRcbvalue, fcRcbvalue + correctionMaxRcb, fcRcbvalue - correctionMinRcb);
1536// if(ibin==3){
1537// cout << " pt="<< fhFcRcb->GetBinCenter(ibin) <<", rval="<<rval<<", fc="<<fcRcbvalue<<" +"<<correctionMaxRcb<<" -"<<correctionMinRcb<<endl;}
1538 }
1539 }
1540 //
1541 // Fill the rest of (asymmetric) histograms
1542 //
1543 if (fAsymUncertainties) {
1544 Double_t x = fhYieldCorr->GetBinCenter(ibin);
1545 fgYieldCorr->SetPoint(ibin,x,value); // i,x,y
1546 fgYieldCorr->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
1547 fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax);
1548 fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
1549 fgYieldCorrExtreme->SetPoint(ibin,x,value); // i,x,y
1550 fgYieldCorrExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
1551 fgYieldCorrConservative->SetPoint(ibin,x,value); // i,x,y
1552 fgYieldCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
1553 // cout << " bin " << ibin << ", correction " << correction << ", min correction unc " << correctionMin << ", max correction unc " << correctionMax << endl;
1554 if(correction>0.){
1555 fgFcConservative->SetPoint(ibin,x,correction);
1556 fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),correctionMin,correctionMax);
1557
1558 fhStatUncEffbFD->SetBinContent(ibin,0.); fhStatUncEffbFD->SetBinError(ibin,correctionUncStatEffb/correction*100.);
1559 fhStatUncEffcFD->SetBinContent(ibin,0.); fhStatUncEffcFD->SetBinError(ibin,correctionUncStatEffc/correction*100.);
1560 // cout << " pt "<< fhStatUncEffcFD->GetBinCenter(ibin) <<" bin "<< ibin<<" fc-stat-c ="<< correctionUncStatEffc/correction <<" fc-stat-b ="<< correctionUncStatEffb/correction <<endl;
1561 }
1562 else{
1563 fgFcConservative->SetPoint(ibin,x,0.);
1564 fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.);
1565 }
1566 }
1567
1568 }
1569 delete [] binwidths;
1570 delete [] limits;
1571
1572}
1573
1574
1575//_________________________________________________________________________________________________________
1576void AliHFPtSpectrum::ComputeSystUncertainties(AliHFSystErr *systematics, Bool_t combineFeedDown) {
1577 //
1578 // Function that re-calculates the global systematic uncertainties
1579 // by calling the class AliHFSystErr and combining those
1580 // (in quadrature) with the feed-down subtraction uncertainties
1581 //
1582
1583 // Estimate the feed-down uncertainty in percentage
1584 Int_t nentries = fgSigmaCorrConservative->GetN();
1585 TGraphAsymmErrors *grErrFeeddown = new TGraphAsymmErrors(nentries);
1586 Double_t x=0., y=0., errx=0., erryl=0., erryh=0;
1587 for(Int_t i=0; i<nentries; i++) {
1588 x=0.; y=0.; errx=0.; erryl=0.; erryh=0.;
1589 fgSigmaCorrConservative->GetPoint(i,x,y);
1590 if(y>0.){
1591 errx = fgSigmaCorrConservative->GetErrorXlow(i) ;
1592 erryl = fgSigmaCorrConservative->GetErrorYlow(i) / y ;
1593 erryh = fgSigmaCorrConservative->GetErrorYhigh(i) / y ;
1594 }
1595 // cout << " x "<< x << " +- "<<errx<<" , y "<<y<<" + "<<erryh<<" - "<<erryl<<endl;
1596 grErrFeeddown->SetPoint(i,x,0.);
1597 grErrFeeddown->SetPointError(i,errx,errx,erryl,erryh); //i, xl, xh, yl, yh
1598 }
1599
1600 // Draw all the systematics independently
1601 systematics->DrawErrors(grErrFeeddown);
1602
1603 // Set the sigma systematic uncertainties
1604 // possibly combine with the feed-down uncertainties
1605 Double_t errylcomb=0., erryhcomb=0;
1606 for(Int_t i=1; i<nentries; i++) {
1607 fgSigmaCorr->GetPoint(i,x,y);
1608 errx = grErrFeeddown->GetErrorXlow(i) ;
1609 erryl = grErrFeeddown->GetErrorYlow(i);
1610 erryh = grErrFeeddown->GetErrorYhigh(i);
1611 if (combineFeedDown) {
1612 errylcomb = systematics->GetTotalSystErr(x,erryl) * y ;
1613 erryhcomb = systematics->GetTotalSystErr(x,erryh) * y ;
1614 } else {
1615 errylcomb = systematics->GetTotalSystErr(x) * y ;
1616 erryhcomb = systematics->GetTotalSystErr(x) * y ;
1617 }
1618 fgSigmaCorr->SetPointError(i,errx,errx,errylcomb,erryhcomb);
1619 //
1620 fhSigmaCorrDataSyst->SetBinContent(i,y);
1621 erryl = systematics->GetTotalSystErr(x) * y ;
1622 fhSigmaCorrDataSyst->SetBinError(i,erryl);
1623 }
1624
1625}
1626
1627
1628//_________________________________________________________________________________________________________
1629void AliHFPtSpectrum::DrawSpectrum(TGraphAsymmErrors *gPrediction) {
1630 //
1631 // Example method to draw the corrected spectrum & the theoretical prediction
1632 //
1633
1634 TCanvas *csigma = new TCanvas("csigma","Draw the corrected cross-section & the prediction");
1635 csigma->SetFillColor(0);
1636 gPrediction->GetXaxis()->SetTitleSize(0.05);
1637 gPrediction->GetXaxis()->SetTitleOffset(0.95);
1638 gPrediction->GetYaxis()->SetTitleSize(0.05);
1639 gPrediction->GetYaxis()->SetTitleOffset(0.95);
1640 gPrediction->GetXaxis()->SetTitle("p_{T} [GeV]");
1641 gPrediction->GetYaxis()->SetTitle("BR #times #frac{d#sigma}{dp_{T}} |_{|y|<0.5} [pb/GeV]");
1642 gPrediction->SetLineColor(kGreen+2);
1643 gPrediction->SetLineWidth(3);
1644 gPrediction->SetFillColor(kGreen+1);
1645 gPrediction->Draw("3CA");
1646 fgSigmaCorr->SetLineColor(kRed);
1647 fgSigmaCorr->SetLineWidth(1);
1648 fgSigmaCorr->SetFillColor(kRed);
1649 fgSigmaCorr->SetFillStyle(0);
1650 fgSigmaCorr->Draw("2");
1651 fhSigmaCorr->SetMarkerColor(kRed);
1652 fhSigmaCorr->Draw("esame");
1653 csigma->SetLogy();
1654 TLegend * leg = new TLegend(0.7,0.75,0.87,0.5);
1655 leg->SetBorderSize(0);
1656 leg->SetLineColor(0);
1657 leg->SetFillColor(0);
1658 leg->SetTextFont(42);
1659 leg->AddEntry(gPrediction,"FONLL ","fl");
1660 leg->AddEntry(fhSigmaCorr,"data stat. unc.","pl");
1661 leg->AddEntry(fgSigmaCorr,"data syst. unc.","f");
1662 leg->Draw();
1663 csigma->Draw();
1664
1665}
1666
1667//_________________________________________________________________________________________________________
1668TH1D * AliHFPtSpectrum::ReweightHisto(TH1D *hToReweight, TH1D *hReference){
1669 //
1670 // Function to reweight histograms for testing purposes:
1671 // This function takes the histo hToReweight and reweights
1672 // it (its pt shape) with respect to hReference
1673 //
1674
1675 // check histograms consistency
1676 Bool_t areconsistent=kTRUE;
1677 areconsistent &= CheckHistosConsistency(hToReweight,hReference);
1678 if (!areconsistent) {
1679 AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
1680 return NULL;
1681 }
1682
1683 // define a new empty histogram
1684 TH1D *hReweighted = (TH1D*)hToReweight->Clone("hReweighted");
1685 hReweighted->Reset();
1686 Double_t weight=1.0;
1687 Double_t yvalue=1.0;
1688 Double_t integralRef = hReference->Integral();
1689 Double_t integralH = hToReweight->Integral();
1690
1691 // now reweight the spectra
1692 //
1693 // the weight is the relative probability of the given pt bin in the reference histo
1694 // divided by its relative probability (to normalize it) on the histo to re-weight
1695 for (Int_t i=0; i<=hToReweight->GetNbinsX(); i++) {
1696 weight = (hReference->GetBinContent(i)/integralRef) / (hToReweight->GetBinContent(i)/integralH) ;
1697 yvalue = hToReweight->GetBinContent(i);
1698 hReweighted->SetBinContent(i,yvalue*weight);
1699 }
1700
1701 return (TH1D*)hReweighted;
1702}
1703
1704//_________________________________________________________________________________________________________
1705TH1D * AliHFPtSpectrum::ReweightRecHisto(TH1D *hRecToReweight, TH1D *hMCToReweight, TH1D *hMCReference){
1706 //
1707 // Function to reweight histograms for testing purposes:
1708 // This function takes the histo hToReweight and reweights
1709 // it (its pt shape) with respect to hReference /hMCToReweight
1710 //
1711
1712 // check histograms consistency
1713 Bool_t areconsistent=kTRUE;
1714 areconsistent &= CheckHistosConsistency(hMCToReweight,hMCReference);
1715 areconsistent &= CheckHistosConsistency(hRecToReweight,hMCReference);
1716 if (!areconsistent) {
1717 AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
1718 return NULL;
1719 }
1720
1721 // define a new empty histogram
1722 TH1D *hReweighted = (TH1D*)hMCToReweight->Clone("hReweighted");
1723 hReweighted->Reset();
1724 TH1D *hRecReweighted = (TH1D*)hRecToReweight->Clone("hRecReweighted");
1725 hRecReweighted->Reset();
1726 Double_t weight=1.0;
1727 Double_t yvalue=1.0, yrecvalue=1.0;
1728 Double_t integralRef = hMCReference->Integral();
1729 Double_t integralH = hMCToReweight->Integral();
1730
1731 // now reweight the spectra
1732 //
1733 // the weight is the relative probability of the given pt bin
1734 // that should be applied in the MC histo to get the reference histo shape
1735 // Probabilities are properly normalized.
1736 for (Int_t i=0; i<=hMCToReweight->GetNbinsX(); i++) {
1737 weight = (hMCReference->GetBinContent(i)/integralRef) / (hMCToReweight->GetBinContent(i)/integralH) ;
1738 yvalue = hMCToReweight->GetBinContent(i);
1739 hReweighted->SetBinContent(i,yvalue*weight);
1740 yrecvalue = hRecToReweight->GetBinContent(i);
1741 hRecReweighted->SetBinContent(i,yrecvalue*weight);
1742 }
1743
1744 return (TH1D*)hRecReweighted;
1745}
1746
1747
1748
1749//_________________________________________________________________________________________________________
1750Int_t AliHFPtSpectrum::FindTH2YBin(TH2D *histo, Float_t yvalue){
1751 //
1752 // Function to find the y-axis bin of a TH2 for a given y-value
1753 //
1754
1755 Int_t nbins = histo->GetNbinsY();
1756 Int_t ybin=0;
1757 for (int j=0; j<=nbins; j++) {
1758 Float_t value = histo->GetYaxis()->GetBinCenter(j);
1759 Float_t width = histo->GetYaxis()->GetBinWidth(j);
1760 // if( TMath::Abs(yvalue-value)<= width/2. ) {
1761 if( TMath::Abs(yvalue-value)<= width ) {
1762 ybin =j;
1763 // cout <<" value "<<value << ", yval "<< yvalue<<", bin width "<<width/2.<< " y ="<<ybin<<endl;
1764 break;
1765 }
1766 }
1767
1768 return ybin;
1769}
1770
1771//_________________________________________________________________________________________________________
1772void AliHFPtSpectrum::ResetStatUncEff(){
1773
1774 Int_t entries = fhDirectEffpt->GetNbinsX();
1775 for(Int_t i=0; i<=entries; i++){
1776 fhDirectEffpt->SetBinError(i,0.);
1777 fhFeedDownEffpt->SetBinError(i,0.);
1778 }
1779
1780}