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