]> 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
Merge branch 'master_patch'
[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.;
6be55058 655 Double_t errvalueStatUncEffc=0.;
b188dc47 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.;
6be55058 662 errvalueStatUncEffc=0.;
b188dc47 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)) ;
b188dc47 707
708 }
709 else {
710 // protect against null denominator
711 errvalueMax = (value!=0.) ?
712 value * TMath::Sqrt( (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
713 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) +
714 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
715 fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] )
716 : 0. ;
717 errvalueMin = errvalueMax;
718 }
719
720 //
721 // Fill the histograms
722 //
723 fhSigmaCorr->SetBinContent(ibin,value);
724 fhSigmaCorr->SetBinError(ibin,errValue);
725
c7d86f5e 726 Double_t x = fhYieldCorr->GetBinCenter(ibin);
727 fgSigmaCorr->SetPoint(ibin,x,value); // i,x,y
728 fgSigmaCorr->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
729 fhSigmaCorrMax->SetBinContent(ibin,value+errvalueMax);
730 fhSigmaCorrMin->SetBinContent(ibin,value-errvalueMin);
731
b188dc47 732 //
733 // Fill the histos and ntuple vs the Eloss hypothesis
734 //
735 if (fPbPbElossHypothesis) {
736
737 // Loop over the Eloss hypothesis
738 if(!fnHypothesis) {
739 AliError("Error reading the fc hypothesis no ntuple, please check !!");
740 delete [] limits;
741 delete [] binwidths;
742 return ;
743 }
744 Int_t entriesHypo = fnHypothesis->GetEntries();
745 Float_t pt=0., Rhy=0., fc=0., fcMin=0., fcMax=0.;
746 fnHypothesis->SetBranchAddress("pt",&pt);
747 if (fFeedDownOption==2) fnHypothesis->SetBranchAddress("Rb",&Rhy);
748 else if (fFeedDownOption==1) fnHypothesis->SetBranchAddress("Rcb",&Rhy);
749 fnHypothesis->SetBranchAddress("fc",&fc);
750 fnHypothesis->SetBranchAddress("fcMax",&fcMax);
751 fnHypothesis->SetBranchAddress("fcMin",&fcMin);
752
753 for (Int_t item=0; item<entriesHypo; item++) {
754
755 fnHypothesis->GetEntry(item);
756 if ( TMath::Abs( pt - fhDirectEffpt->GetBinCenter(ibin) ) > 0.15 ) continue;
757
758 Double_t yieldRcbvalue = (fhRECpt->GetBinContent(ibin) ) ? fhRECpt->GetBinContent(ibin) * fc : 0. ;
759 yieldRcbvalue /= fhRECpt->GetBinWidth(ibin) ;
760 Double_t yieldRcbvalueMax = (fhRECpt->GetBinContent(ibin) ) ? fhRECpt->GetBinContent(ibin) * fcMax : 0. ;
761 yieldRcbvalueMax /= fhRECpt->GetBinWidth(ibin) ;
762 Double_t yieldRcbvalueMin = (fhRECpt->GetBinContent(ibin) ) ? fhRECpt->GetBinContent(ibin) * fcMin : 0. ;
763 yieldRcbvalueMin /= fhRECpt->GetBinWidth(ibin) ;
764
765 // Sigma calculation
766 // Sigma = ( 1. / (lumi * delta_y * BR_c * ParticleAntiPartFactor * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
767 Double_t sigmaRcbvalue = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)>0.) ?
768 ( yieldRcbvalue / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
769 : 0. ;
770 Double_t sigmaRcbvalueMax = (sigmaRcbvalue!=0.) ?
771 ( yieldRcbvalueMax / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
772 : 0. ;
773 Double_t sigmaRcbvalueMin = (sigmaRcbvalue!=0.) ?
774 ( yieldRcbvalueMin / ( deltaY * branchingRatioC * fParticleAntiParticle * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
775 : 0. ;
776 // Sigma statistical uncertainty:
777 // delta_sigma = sigma * sqrt ( (delta_spectra/spectra)^2 ) = sigma * ( delta_spectra / (spectra-corr * binwidth) )
778 Double_t sigmaRcbvalueStatUnc = (sigmaRcbvalue!=0.) ?
779 sigmaRcbvalue * ( fhRECpt->GetBinError(ibin) / ( yieldRcbvalue * fhRECpt->GetBinWidth(ibin) ) ) : 0. ;
780
781 fhSigmaCorrRcb->Fill( fhSigmaCorr->GetBinCenter(ibin) , Rhy, sigmaRcbvalue );
782 // if(ibin==3)
783 // cout << " pt "<< fhRECpt->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<" fc-value="<< fhFcRcb->GetBinContent(ibin,rbin) <<", yield-fcRbvalue="<<yieldRcbvalue<<", sigma-fcRbvalue="<<sigmaRcbvalue<<endl;
784 fnSigma->Fill(fhRECpt->GetBinCenter(ibin), fhRECpt->GetBinContent(ibin),
785 Rhy, fc, yieldRcbvalue, sigmaRcbvalue, sigmaRcbvalueStatUnc, sigmaRcbvalueMax, sigmaRcbvalueMin );
786 }
787 }
788 //
789 // Fill the TGraphAsymmErrors
790 if (fAsymUncertainties) {
b188dc47 791 fgSigmaCorrExtreme->SetPoint(ibin,x,value); // i,x,y
792 fgSigmaCorrExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
793 fgSigmaCorrConservative->SetPoint(ibin,x,value); // i,x,y
794 fgSigmaCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueConservativeMin,errvalueConservativeMax); // i,xl,xh,yl,yh
795
796 fhStatUncEffcSigma->SetBinContent(ibin,0.);
797 if(value>0.) fhStatUncEffcSigma->SetBinError(ibin,((errvalueStatUncEffc/value)*100.));
798 fhStatUncEffbSigma->SetBinContent(ibin,0.); fhStatUncEffbSigma->SetBinError(ibin,0.);
799 // cout << " pt "<< fhRECpt->GetBinCenter(ibin) <<" bin "<< ibin<<" stat-unc-c-sigma "<< errvalueStatUncEffc/value << endl;
800 }
801
802 }
803 delete [] binwidths;
804 delete [] limits;
805
806}
807
808//_________________________________________________________________________________________________________
809TH1D * AliHFPtSpectrum::EstimateEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco, const char *name) {
810 //
811 // Function that computes the acceptance and efficiency correction
812 // based on the simulated and reconstructed spectra
813 // and using the reconstructed spectra bin width
814 //
815 // eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
816 //
817
818 if(!fhRECpt){
819 AliInfo("Hey, the reconstructed histogram was not set yet !");
820 return NULL;
821 }
822
823 Int_t nbins = fhRECpt->GetNbinsX();
824 Double_t *limits = new Double_t[nbins+1];
825 Double_t xlow=0.,binwidth=0.;
826 for (Int_t i=1; i<=nbins; i++) {
827 binwidth = fhRECpt->GetBinWidth(i);
828 xlow = fhRECpt->GetBinLowEdge(i);
829 limits[i-1] = xlow;
830 }
831 limits[nbins] = xlow + binwidth;
832
833 TH1D * hEfficiency = new TH1D(name," acceptance #times efficiency",nbins,limits);
834
835 Double_t *sumSimu=new Double_t[nbins];
836 Double_t *sumReco=new Double_t[nbins];
837 for (Int_t ibin=0; ibin<nbins; ibin++){
838 sumSimu[ibin]=0.; sumReco[ibin]=0.;
839 }
840 for (Int_t ibin=0; ibin<=hSimu->GetNbinsX(); ibin++){
841
842 for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++){
843 if ( hSimu->GetBinCenter(ibin)>limits[ibinrec] &&
844 hSimu->GetBinCenter(ibin)<limits[ibinrec+1] ) {
845 sumSimu[ibinrec]+=hSimu->GetBinContent(ibin);
846 }
847 if ( hReco->GetBinCenter(ibin)>limits[ibinrec] &&
848 hReco->GetBinCenter(ibin)<limits[ibinrec+1] ) {
849 sumReco[ibinrec]+=hReco->GetBinContent(ibin);
850 }
851 }
852
853 }
854
855
856 // the efficiency is computed as reco/sim (in each bin)
857 // its uncertainty is err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
858 Double_t eff=0., erreff=0.;
859 for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++) {
860 if (sumSimu[ibinrec]!= 0. && sumReco[ibinrec]!=0.) {
861 eff = sumReco[ibinrec] / sumSimu[ibinrec] ;
862 // protection in case eff > 1.0
863 // test calculation (make the argument of the sqrt positive)
864 erreff = TMath::Sqrt( eff * TMath::Abs(1.0 - eff) ) / TMath::Sqrt( sumSimu[ibinrec] );
865 }
866 else { eff=0.0; erreff=0.; }
867 hEfficiency->SetBinContent(ibinrec+1,eff);
868 hEfficiency->SetBinError(ibinrec+1,erreff);
869 }
870
871 delete [] sumSimu;
872 delete [] sumReco;
873
874 return (TH1D*)hEfficiency;
875}
876
877//_________________________________________________________________________________________________________
878void AliHFPtSpectrum::EstimateAndSetDirectEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco) {
879 //
880 // Function that computes the Direct acceptance and efficiency correction
881 // based on the simulated and reconstructed spectra
882 // and using the reconstructed spectra bin width
883 //
884 // eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
885 //
886
887 if(!fhRECpt || !hSimu || !hReco){
888 AliError("Hey, the reconstructed histogram was not set yet !");
889 return;
890 }
891
892 fhDirectEffpt = EstimateEfficiencyRecoBin(hSimu,hReco,"fhDirectEffpt");
893 fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance #times efficiency");
894
895}
896
897//_________________________________________________________________________________________________________
898void AliHFPtSpectrum::EstimateAndSetFeedDownEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco) {
899 //
900 // Function that computes the Feed-Down acceptance and efficiency correction
901 // based on the simulated and reconstructed spectra
902 // and using the reconstructed spectra bin width
903 //
904 // eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
905 //
906
907 if(!fhRECpt || !hSimu || !hReco){
908 AliError("Hey, the reconstructed histogram was not set yet !");
909 return;
910 }
911
912 fhFeedDownEffpt = EstimateEfficiencyRecoBin(hSimu,hReco,"fhFeedDownEffpt");
913 fhFeedDownEffpt->SetNameTitle("fhFeedDownEffpt"," feed-down acceptance #times efficiency");
914
915}
916
917//_________________________________________________________________________________________________________
918Bool_t AliHFPtSpectrum::Initialize(){
919 //
920 // Initialization of the variables (histograms)
921 //
922
923 if (fFeedDownOption==0) {
924 AliInfo("Getting ready for the corrections without feed-down consideration");
925 } else if (fFeedDownOption==1) {
926 AliInfo("Getting ready for the fc feed-down correction calculation");
927 } else if (fFeedDownOption==2) {
928 AliInfo("Getting ready for the Nb feed-down correction calculation");
929 } else { AliError("The calculation option must be <=2"); return kFALSE; }
930
931 // Start checking the input histograms consistency
932 Bool_t areconsistent=kTRUE;
933
934 // General checks
935 if (!fhDirectEffpt || !fhRECpt) {
936 AliError(" Reconstructed spectra and/or the Nc efficiency distributions are not defined");
937 return kFALSE;
938 }
939 areconsistent &= CheckHistosConsistency(fhRECpt,fhDirectEffpt);
940 if (!areconsistent) {
941 AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
942 return kFALSE;
943 }
944 if (fFeedDownOption==0) return kTRUE;
945
946 //
947 // Common checks for options 1 (fc) & 2(Nb)
948 if (!fhFeedDownMCpt || !fhFeedDownEffpt) {
949 AliError(" Theoretical Nb and/or the Nb efficiency distributions are not defined");
950 return kFALSE;
951 }
952 areconsistent &= CheckHistosConsistency(fhRECpt,fhFeedDownMCpt);
953 areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownEffpt);
954 if (fAsymUncertainties) {
955 if (!fhFeedDownMCptMax || !fhFeedDownMCptMin) {
956 AliError(" Max/Min theoretical Nb distributions are not defined");
957 return kFALSE;
958 }
959 areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownMCptMax);
960 }
961 if (!areconsistent) {
962 AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
963 return kFALSE;
964 }
965 if (fFeedDownOption>1) return kTRUE;
966
967 //
968 // Now checks for option 1 (fc correction)
969 if (!fhDirectMCpt) {
970 AliError("Theoretical Nc distributions is not defined");
971 return kFALSE;
972 }
973 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhFeedDownMCpt);
974 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectEffpt);
975 if (fAsymUncertainties) {
976 if (!fhDirectMCptMax || !fhDirectMCptMin) {
977 AliError(" Max/Min theoretical Nc distributions are not defined");
978 return kFALSE;
979 }
980 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectMCptMax);
981 }
982 if (!areconsistent) {
983 AliInfo("Histograms required for fc correction are not consistent (bin width, bounds)");
984 return kFALSE;
985 }
986
987 return kTRUE;
988}
989
990//_________________________________________________________________________________________________________
991Bool_t AliHFPtSpectrum::CheckHistosConsistency(TH1D *h1, TH1D *h2){
992 //
993 // Check the histograms consistency (bins, limits)
994 //
995
996 if (!h1 || !h2) {
997 AliError("One or both histograms don't exist");
998 return kFALSE;
999 }
1000
1001 Double_t binwidth1 = h1->GetBinWidth(1);
1002 Double_t binwidth2 = h2->GetBinWidth(1);
1003 Double_t min1 = h1->GetBinCenter(1) - (binwidth1/2.) ;
1004// Double_t max1 = h1->GetBinCenter(nbins1) + (binwidth1/2.) ;
1005 Double_t min2 = h2->GetBinCenter(1) - (binwidth2/2.) ;
1006// Double_t max2 = h2->GetBinCenter(nbins2) + (binwidth2/2.) ;
1007
1008 if (binwidth1!=binwidth2) {
1009 AliInfo(" histograms with different bin width");
1010 return kFALSE;
1011 }
1012 if (min1!=min2) {
1013 AliInfo(" histograms with different minimum");
1014 return kFALSE;
1015 }
1016// if (max1!=max2) {
1017// AliInfo(" histograms with different maximum");
1018// return kFALSE;
1019// }
1020
1021 return kTRUE;
1022}
1023
1024//_________________________________________________________________________________________________________
1025void AliHFPtSpectrum::CalculateFeedDownCorrectionFc(){
1026 //
1027 // Compute fc factor and its uncertainties bin by bin
1028 // fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
1029 //
1030 // uncertainties: (conservative) combine the upper/lower N_b & N_c predictions together
1031 // (extreme) combine the upper N_b predictions with the lower N_c predictions & viceversa
1032 // systematic uncertainty on the acceptance x efficiency b/c ratio are included
1033 //
1034 // In addition, in HIC the feed-down correction varies with an energy loss hypothesis: Raa(c-->D) / Raa(b-->D) = Rcb
1035 // fc (Rcb) = ( 1. / ( 1 + (eff_b/eff_c)*(N_b/N_c)* (1/Rcb) ) );
1036 //
1037 AliInfo("Calculating the feed-down correction factor (fc method)");
1038
1039 // define the variables
1040 Int_t nbins = fhRECpt->GetNbinsX();
1041 Double_t binwidth = fhRECpt->GetBinWidth(1);
1042 Double_t *limits = new Double_t[nbins+1];
1043 Double_t *binwidths = new Double_t[nbins];
1044 Double_t xlow=0.;
1045 for (Int_t i=1; i<=nbins; i++) {
1046 binwidth = fhRECpt->GetBinWidth(i);
1047 xlow = fhRECpt->GetBinLowEdge(i);
1048 limits[i-1] = xlow;
1049 binwidths[i-1] = binwidth;
1050 }
1051 limits[nbins] = xlow + binwidth;
1052
1053 Double_t correction=1.;
1054 Double_t theoryRatio=1.;
1055 Double_t effRatio=1.;
1056 Double_t correctionExtremeA=1., correctionExtremeB=1.;
1057 Double_t theoryRatioExtremeA=1., theoryRatioExtremeB=1.;
1058 Double_t correctionConservativeA=1., correctionConservativeB=1.;
1059 Double_t theoryRatioConservativeA=1., theoryRatioConservativeB=1.;
273fb9a5 1060 // Double_t correctionUnc=0.;
b188dc47 1061 Double_t correctionExtremeAUnc=0., correctionExtremeBUnc=0.;
1062 Double_t correctionConservativeAUnc=0., correctionConservativeBUnc=0.;
1063
1064 // declare the output histograms
1065 fhFc = new TH1D("fhFc","fc correction factor",nbins,limits);
1066 fhFcMax = new TH1D("fhFcMax","max fc correction factor",nbins,limits);
1067 fhFcMin = new TH1D("fhFcMin","min fc correction factor",nbins,limits);
1068 if(fPbPbElossHypothesis) {
1069 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.);
1070 fnHypothesis = new TNtuple("fnHypothesis"," Feed-down correction vs hypothesis (fc)","pt:Rcb:fc:fcMax:fcMin");
1071 }
1072 // two local control histograms
1073 TH1D *hTheoryRatio = new TH1D("hTheoryRatio","Theoretical B-->D over c-->D (feed-down/direct) ratio",nbins,limits);
1074 TH1D *hEffRatio = new TH1D("hEffRatio","Efficiency B-->D over c-->D (feed-down/direct) ratio",nbins,limits);
1075 // and the output TGraphAsymmErrors
1076 if (fAsymUncertainties) {
1077 fgFcExtreme = new TGraphAsymmErrors(nbins+1);
1078 fgFcExtreme->SetNameTitle("fgFcExtreme","fgFcExtreme");
1079 fgFcConservative = new TGraphAsymmErrors(nbins+1);
1080 fgFcConservative->SetNameTitle("fgFcConservative","fgFcConservative");
1081 }
1082
1083 fhStatUncEffcFD = new TH1D("fhStatUncEffcFD","direct charm stat unc on the feed-down correction",nbins,limits);
1084 fhStatUncEffbFD = new TH1D("fhStatUncEffbFD","secondary charm stat unc on the feed-down correction",nbins,limits);
1085 Double_t correctionConservativeAUncStatEffc=0., correctionConservativeBUncStatEffc=0.;
1086 Double_t correctionConservativeAUncStatEffb=0., correctionConservativeBUncStatEffb=0.;
1087
1088 //
1089 // Compute fc
1090 //
1091 for (Int_t ibin=1; ibin<=nbins; ibin++) {
1092
1287e14b 1093 // Variables initialization
1094 correction=1.; theoryRatio=1.; effRatio=1.;
1095 correctionExtremeA=1.; correctionExtremeB=1.;
1096 theoryRatioExtremeA=1.; theoryRatioExtremeB=1.;
1097 correctionConservativeA=1.; correctionConservativeB=1.;
1098 theoryRatioConservativeA=1.; theoryRatioConservativeB=1.;
273fb9a5 1099 // correctionUnc=0.;
1287e14b 1100 correctionExtremeAUnc=0.; correctionExtremeBUnc=0.;
1101 correctionConservativeAUnc=0.; correctionConservativeBUnc=0.;
1102 correctionConservativeAUncStatEffc=0.; correctionConservativeBUncStatEffc=0.;
1103 correctionConservativeAUncStatEffb=0.; correctionConservativeBUncStatEffb=0.;
1104
b188dc47 1105 // theory_ratio = (N_b/N_c)
1106 theoryRatio = (fhDirectMCpt->GetBinContent(ibin)>0. && fhFeedDownMCpt->GetBinContent(ibin)>0.) ?
1107 fhFeedDownMCpt->GetBinContent(ibin) / fhDirectMCpt->GetBinContent(ibin) : 1.0 ;
1108
1109 //
1110 // Calculate the uncertainty [ considering only the theoretical uncertainties on Nb & Nc for now !!! ]
1111 //
1112 // extreme A = direct-max, feed-down-min
1113 theoryRatioExtremeA = (fhDirectMCptMax->GetBinContent(ibin)>0. && fhFeedDownMCptMin->GetBinContent(ibin)>0.) ?
1114 fhFeedDownMCptMin->GetBinContent(ibin) / fhDirectMCptMax->GetBinContent(ibin) : 1.0 ;
1115 // extreme B = direct-min, feed-down-max
1116 theoryRatioExtremeB = (fhDirectMCptMin->GetBinContent(ibin)>0. && fhDirectMCptMax->GetBinContent(ibin)>0.) ?
1117 fhFeedDownMCptMax->GetBinContent(ibin) / fhDirectMCptMin->GetBinContent(ibin) : 1.0 ;
1118 // conservative A = direct-max, feed-down-max
1119 theoryRatioConservativeA = (fhDirectMCptMax->GetBinContent(ibin)>0. && fhFeedDownMCptMin->GetBinContent(ibin)>0.) ?
1120 fhFeedDownMCptMax->GetBinContent(ibin) / fhDirectMCptMax->GetBinContent(ibin) : 1.0 ;
1121 // conservative B = direct-min, feed-down-min
1122 theoryRatioConservativeB = (fhDirectMCptMin->GetBinContent(ibin)>0. && fhDirectMCptMax->GetBinContent(ibin)>0.) ?
1123 fhFeedDownMCptMin->GetBinContent(ibin) / fhDirectMCptMin->GetBinContent(ibin) : 1.0 ;
1124
1125 // eff_ratio = (eff_b/eff_c)
1126 effRatio = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ?
1127 fhFeedDownEffpt->GetBinContent(ibin) / fhDirectEffpt->GetBinContent(ibin) : 1.0 ;
1128
1129 // fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
7a385c9e 1130 if( TMath::Abs(effRatio - 1.0)<0.0001 || TMath::Abs(theoryRatio - 1.0)<0.0001 ) {
b188dc47 1131 correction = 1.0;
1132 correctionExtremeA = 1.0;
1133 correctionExtremeB = 1.0;
1134 correctionConservativeA = 1.0;
1135 correctionConservativeB = 1.0;
1136 }
1137 else {
1138 correction = ( 1. / ( 1 + ( effRatio * theoryRatio ) ) );
1139 correctionExtremeA = ( 1. / ( 1 + ( effRatio * theoryRatioExtremeA ) ) );
1140 correctionExtremeB = ( 1. / ( 1 + ( effRatio * theoryRatioExtremeB ) ) );
1141 correctionConservativeA = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeA ) ) );
1142 correctionConservativeB = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeB ) ) );
1143 }
1144
1145
1146 // fc uncertainty from (eff_b/eff_c) = fc^2 * (N_b/N_c) * delta(eff_b/eff_c)
1147 // delta(eff_b/eff_c) is a percentage = effRatio * sqrt( fGlobalEfficiencyUncertainties[1]^2 + unc_eff_c ^2 + unc_eff_b ^2 )
1148 Double_t relEffUnc = TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] +
1149 (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
1150 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))
1151 );
1152
273fb9a5 1153 // correctionUnc = correction*correction * theoryRatio * effRatio * relEffUnc;
b188dc47 1154 correctionExtremeAUnc = correctionExtremeA*correctionExtremeA * theoryRatioExtremeA * effRatio * relEffUnc;
1155 correctionExtremeBUnc = correctionExtremeB*correctionExtremeB * theoryRatioExtremeB * effRatio * relEffUnc;
1156
1157 correctionConservativeAUnc = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA *effRatio * relEffUnc;
1158 //
1159 correctionConservativeAUncStatEffc = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA *effRatio *
1160 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin));
1161 correctionConservativeAUncStatEffb = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA *effRatio *
1162 (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin));
1163
1164 correctionConservativeBUnc = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB *effRatio * relEffUnc;
1165
1166 correctionConservativeBUncStatEffb = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB *effRatio *
1167 (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin));
1168 correctionConservativeBUncStatEffc = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB *effRatio *
1169 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin));
1170
1171
1172 // Fill in the histograms
1173 hTheoryRatio->SetBinContent(ibin,theoryRatio);
1174 hEffRatio->SetBinContent(ibin,effRatio);
1175 fhFc->SetBinContent(ibin,correction);
1176 //
1177 // Estimate how the result varies vs charm/beauty Eloss hypothesis
1178 //
7a385c9e 1179 if ( correction>1.0e-16 && fPbPbElossHypothesis){
b188dc47 1180 // Loop over the Eloss hypothesis
1181 // Int_t rbin=0;
1182 for (Float_t rval=0.0025; rval<4.0; rval+=0.005){
1183 // Central fc with Eloss hypothesis
1184 Double_t correctionRcb = ( 1. / ( 1 + ( effRatio * theoryRatio * (1/rval) ) ) );
1185 fhFcRcb->Fill( fhFc->GetBinCenter(ibin) , rval, correctionRcb );
1186 // if(ibin==3){
1187 // cout << " pt "<< fhFc->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<", fc-Rcb-value="<<correctionRcb<<endl;
1188 // rbin++;
1189 // }
1190 // Upper / lower fc with up / low FONLL bands and Eloss hypothesis
1191 Double_t correctionConservativeARcb = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeA * (1/rval) ) ) );
1192 Double_t correctionConservativeBRcb = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeB * (1/rval) ) ) );
1193 Double_t correctionConservativeARcbUnc = correctionConservativeARcb*correctionConservativeARcb * theoryRatioConservativeA * (1/rval) *effRatio * relEffUnc;
1194 Double_t correctionConservativeBRcbUnc = correctionConservativeBRcb*correctionConservativeBRcb * theoryRatioConservativeB * (1/rval) *effRatio * relEffUnc;
1195 Double_t consvalRcb[4] = { correctionConservativeARcb - correctionConservativeARcbUnc, correctionConservativeARcb + correctionConservativeARcbUnc,
1196 correctionConservativeBRcb - correctionConservativeBRcbUnc, correctionConservativeBRcb + correctionConservativeBRcbUnc};
1197 Double_t uncConservativeRcbMin = correctionRcb - TMath::MinElement(4,consvalRcb);
1198 Double_t uncConservativeRcbMax = TMath::MaxElement(4,consvalRcb) - correctionRcb;
1199// if(ibin==3)
1200// cout << " pt="<<fhDirectEffpt->GetBinCenter(ibin)<<", hypo="<<rval<<", fc="<<correctionRcb<<" +"<<uncConservativeRcbMax<<" -"<<uncConservativeRcbMin<<endl;
1201 fnHypothesis->Fill( fhDirectEffpt->GetBinCenter(ibin), rval, correctionRcb, correctionRcb+uncConservativeRcbMax, correctionRcb-uncConservativeRcbMin);
1202 }
1203 }
1204 //
1205 // Fill the rest of (asymmetric) histograms
1206 //
1207 if (fAsymUncertainties) {
1208 Double_t x = fhDirectMCpt->GetBinCenter(ibin);
1209 Double_t val[4] = { correctionExtremeA + correctionExtremeAUnc, correctionExtremeA - correctionExtremeAUnc,
1210 correctionExtremeB + correctionExtremeBUnc, correctionExtremeB - correctionExtremeBUnc };
1211 Double_t uncExtremeMin = correction - TMath::MinElement(4,val);
1212 Double_t uncExtremeMax = TMath::MaxElement(4,val) - correction;
1213 fgFcExtreme->SetPoint(ibin,x,correction); // i,x,y
1214 fgFcExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),uncExtremeMin,uncExtremeMax); // i,xl,xh,yl,yh
1215 fhFcMax->SetBinContent(ibin,correction+uncExtremeMax);
1216 fhFcMin->SetBinContent(ibin,correction-uncExtremeMin);
1217 Double_t consval[4] = { correctionConservativeA - correctionConservativeAUnc, correctionConservativeA + correctionConservativeAUnc,
1218 correctionConservativeB - correctionConservativeBUnc, correctionConservativeB + correctionConservativeBUnc};
1219 Double_t uncConservativeMin = correction - TMath::MinElement(4,consval);
1220 Double_t uncConservativeMax = TMath::MaxElement(4,consval) - correction;
1221 fgFcConservative->SetPoint(ibin,x,correction); // i,x,y
1222 fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),uncConservativeMin,uncConservativeMax); // i,xl,xh,yl,yh
1223 if( !(correction>0.) ){
1224 fgFcExtreme->SetPoint(ibin,x,0.); // i,x,y
1225 fgFcExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.); // i,xl,xh,yl,yh
1226 fgFcConservative->SetPoint(ibin,x,0.); // i,x,y
1227 fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.); // i,xl,xh,yl,yh
1228 }
1229
1230 Double_t valStatEffc[2] = { correctionConservativeAUncStatEffc/correctionConservativeA,
1231 correctionConservativeBUncStatEffc/correctionConservativeB };
1232 Double_t valStatEffb[2] = { correctionConservativeAUncStatEffb/correctionConservativeA,
1233 correctionConservativeBUncStatEffb/correctionConservativeB };
1234 Double_t uncConservativeStatEffc = TMath::MaxElement(2,valStatEffc);
1235 Double_t uncConservativeStatEffb = TMath::MaxElement(2,valStatEffb);
1236 fhStatUncEffcFD->SetBinContent(ibin,0.); fhStatUncEffcFD->SetBinError(ibin,uncConservativeStatEffc*100.);
1237 fhStatUncEffbFD->SetBinContent(ibin,0.); fhStatUncEffbFD->SetBinError(ibin,uncConservativeStatEffb*100.);
1238 // cout << " pt "<< fhStatUncEffcFD->GetBinCenter(ibin) <<" bin "<< ibin<<" fc-stat-c ="<<uncConservativeStatEffc<<" fc-stat-b ="<<uncConservativeStatEffb<<endl;
1239 }
1240
1241 }
1242 delete [] binwidths;
1243 delete [] limits;
1244
1245}
1246
1247//_________________________________________________________________________________________________________
1248void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumFc(){
1249 //
1250 // Compute the feed-down corrected spectrum if feed-down correction is done via fc factor (bin by bin)
1251 // physics = reco * fc / bin-width
1252 //
1253 // uncertainty: (stat) delta_physics = physics * sqrt ( (delta_reco/reco)^2 )
1254 // (syst but feed-down) delta_physics = physics * sqrt ( (delta_reco_syst/reco)^2 )
1255 // (feed-down syst) delta_physics = physics * sqrt ( (delta_fc/fc)^2 )
1256 //
1257 // ( Calculation done bin by bin )
1258 //
1259 // In addition, in HIC the feed-down correction varies with an energy loss hypothesis: Raa(c-->D) / Raa(b-->D) = Rcb
1260
1261 AliInfo(" Calculating the feed-down corrected spectrum (fc method)");
1262
1263 if (!fhFc || !fhRECpt) {
1264 AliError(" Reconstructed or fc distributions are not defined");
1265 return;
1266 }
1267
1268 Int_t nbins = fhRECpt->GetNbinsX();
1269 Double_t value = 0., errvalue = 0., errvalueMax= 0., errvalueMin= 0.;
1270 Double_t valueExtremeMax= 0., valueExtremeMin= 0.;
1271 Double_t valueConservativeMax= 0., valueConservativeMin= 0.;
1272 Double_t binwidth = fhRECpt->GetBinWidth(1);
1273 Double_t *limits = new Double_t[nbins+1];
1274 Double_t *binwidths = new Double_t[nbins];
1275 Double_t xlow=0.;
1276 for (Int_t i=1; i<=nbins; i++) {
1277 binwidth = fhRECpt->GetBinWidth(i);
1278 xlow = fhRECpt->GetBinLowEdge(i);
1279 limits[i-1] = xlow;
1280 binwidths[i-1] = binwidth;
1281 }
1282 limits[nbins] = xlow + binwidth;
1283
1284 // declare the output histograms
1285 fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by fc)",nbins,limits);
1286 fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by fc)",nbins,limits);
1287 fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by fc)",nbins,limits);
1288 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.);
1289 // and the output TGraphAsymmErrors
c7d86f5e 1290 fgYieldCorr = new TGraphAsymmErrors(nbins+1);
b188dc47 1291 if (fAsymUncertainties){
b188dc47 1292 fgYieldCorrExtreme = new TGraphAsymmErrors(nbins+1);
1293 fgYieldCorrConservative = new TGraphAsymmErrors(nbins+1);
1294 }
1295
1296 //
1297 // Do the calculation
1298 //
1299 for (Int_t ibin=1; ibin<=nbins; ibin++) {
1300
1287e14b 1301 // Variables initialization
1302 value = 0.; errvalue = 0.; errvalueMax= 0.; errvalueMin= 0.;
1303 valueExtremeMax= 0.; valueExtremeMin= 0.;
1304 valueConservativeMax= 0.; valueConservativeMin= 0.;
1305
1306
b188dc47 1307 // calculate the value
1308 // physics = reco * fc / bin-width
1309 value = (fhRECpt->GetBinContent(ibin) && fhFc->GetBinContent(ibin)) ?
1310 fhRECpt->GetBinContent(ibin) * fhFc->GetBinContent(ibin) : 0. ;
1311 value /= fhRECpt->GetBinWidth(ibin) ;
1312
1313 // Statistical uncertainty
1314 // (stat) delta_physics = physics * sqrt ( (delta_reco/reco)^2 )
1315 errvalue = (value!=0. && fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) ?
1316 value * (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) : 0. ;
1317
1318 // Calculate the systematic uncertainties
1319 // (syst but feed-down) delta_physics = physics * sqrt ( (delta_reco_syst/reco)^2 )
1320 // (feed-down syst) delta_physics = physics * sqrt ( (delta_fc/fc)^2 )
1321 //
1322 // Protect against null denominator. If so, define uncertainty as null
1323 if (fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) {
1324
1325 if (fAsymUncertainties) {
1326
1327 // Systematics but feed-down
1328 if (fgRECSystematics) {
1329 errvalueMax = value * ( fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinContent(ibin) );
1330 errvalueMin = value * ( fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinContent(ibin) );
1331 }
1332 else { errvalueMax = 0.; errvalueMin = 0.; }
1333
1334 // Extreme feed-down systematics
1335 valueExtremeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcExtreme->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1336 valueExtremeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcExtreme->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1337
1338 // Conservative feed-down systematics
1339 valueConservativeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcConservative->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1340 valueConservativeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcConservative->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1341
1342 }
1343
1344 }
1345 else { errvalueMax = 0.; errvalueMin = 0.; }
1346
1347 //
1348 // Fill in the histograms
1349 //
1350 fhYieldCorr->SetBinContent(ibin,value);
1351 fhYieldCorr->SetBinError(ibin,errvalue);
1352 //
1353 // Fill the histos and ntuple vs the Eloss hypothesis
1354 //
1355 if (fPbPbElossHypothesis) {
1356 // Loop over the Eloss hypothesis
1357 for (Float_t rval=0.0025; rval<4.0; rval+=0.005){
1358 Int_t rbin = FindTH2YBin(fhYieldCorrRcb,rval);
1359 Double_t fcRcbvalue = fhFcRcb->GetBinContent(ibin,rbin);
1360 // physics = reco * fcRcb / bin-width
1361 Double_t Rcbvalue = (fhRECpt->GetBinContent(ibin) && fcRcbvalue) ?
1362 fhRECpt->GetBinContent(ibin) * fcRcbvalue : 0. ;
1363 Rcbvalue /= fhRECpt->GetBinWidth(ibin) ;
1364 fhYieldCorrRcb->Fill( fhYieldCorr->GetBinCenter(ibin) , rval, Rcbvalue );
1365 // cout << " pt "<< fhRECpt->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<" fc-fcRbvalue="<<fcRcbvalue<<", yield="<<Rcbvalue<<endl;
1366 }
1367 }
1368 if (fAsymUncertainties) {
1369 Double_t center = fhYieldCorr->GetBinCenter(ibin);
1370 fgYieldCorr->SetPoint(ibin,center,value); // i,x,y
1371 fgYieldCorr->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
1372 fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax);
1373 fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
1374 fgYieldCorrExtreme->SetPoint(ibin,center,value);
1375 fgYieldCorrExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),value-valueExtremeMin,valueExtremeMax-value);
1376 fgYieldCorrConservative->SetPoint(ibin,center,value);
1377 fgYieldCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),value-valueConservativeMin,valueConservativeMax-value);
1378 }
1379
1380 }
1381 delete [] binwidths;
1382 delete [] limits;
1383
1384}
1385
1386//_________________________________________________________________________________________________________
1387void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumNb(Double_t deltaY, Double_t branchingRatioBintoFinalDecay) {
1388 //
1389 // Compute the feed-down corrected spectrum if feed-down correction is done via Nb (bin by bin)
1390 // physics = [ reco - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) ] / bin-width
1391 //
1392 // uncertainty: (stat) delta_physics = sqrt ( (delta_reco)^2 ) / bin-width
1393 // (syst but feed-down) delta_physics = sqrt ( (delta_reco_syst)^2 ) / bin-width
1394 // (feed-down syst) delta_physics = sqrt ( (k*delta_lumi/lumi)^2 + (k*delta_eff_trig/eff_trig)^2
1395 // + (k*delta_Nb/Nb)^2 + (k*delta_eff/eff)^2 + (k*global_eff_ratio)^2 ) / bin-width
1396 // where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th
1397 //
1398 // In addition, in HIC the feed-down correction varies with an energy loss hypothesis: Raa(b-->D) = Rb
1399 // physics = [ reco - ( Tab * Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th * Rb ) ] / bin-width
1400 //
1401 AliInfo("Calculating the feed-down correction factor and spectrum (Nb method)");
1402
1403 Int_t nbins = fhRECpt->GetNbinsX();
1404 Double_t binwidth = fhRECpt->GetBinWidth(1);
1405 Double_t value = 0., errvalue = 0., errvalueMax = 0., errvalueMin = 0., kfactor = 0.;
1406 Double_t errvalueExtremeMax = 0., errvalueExtremeMin = 0.;
1407 Double_t *limits = new Double_t[nbins+1];
1408 Double_t *binwidths = new Double_t[nbins];
1409 Double_t xlow=0.;
1410 for (Int_t i=1; i<=nbins; i++) {
1411 binwidth = fhRECpt->GetBinWidth(i);
1412 xlow = fhRECpt->GetBinLowEdge(i);
1413 limits[i-1] = xlow;
1414 binwidths[i-1] = binwidth;
1415 }
1416 limits[nbins] = xlow + binwidth;
1417
1418 // declare the output histograms
1419 fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by Nb)",nbins,limits);
1420 fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by Nb)",nbins,limits);
1421 fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by Nb)",nbins,limits);
1422 if(fPbPbElossHypothesis) {
1423 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.);
1424 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.);
1425 fnHypothesis = new TNtuple("fnHypothesis"," Feed-down correction vs hypothesis (Nb)","pt:Rb:fc:fcMax:fcMin");
1426 }
1427 // and the output TGraphAsymmErrors
c7d86f5e 1428 fgYieldCorr = new TGraphAsymmErrors(nbins+1);
b188dc47 1429 if (fAsymUncertainties){
b188dc47 1430 fgYieldCorrExtreme = new TGraphAsymmErrors(nbins+1);
1431 fgYieldCorrConservative = new TGraphAsymmErrors(nbins+1);
1432 // Define fc-conservative
1433 fgFcConservative = new TGraphAsymmErrors(nbins+1);
1434 AliInfo(" Beware the conservative & extreme uncertainties are equal by definition !");
1435 }
1436
1437 // variables to define fc-conservative
1438 double correction=0, correctionMax=0., correctionMin=0.;
1439
1440 fhStatUncEffcFD = new TH1D("fhStatUncEffcFD","direct charm stat unc on the feed-down correction",nbins,limits);
1441 fhStatUncEffbFD = new TH1D("fhStatUncEffbFD","secondary charm stat unc on the feed-down correction",nbins,limits);
1442 Double_t correctionUncStatEffc=0.;
1443 Double_t correctionUncStatEffb=0.;
1444
1445
1446 //
1447 // Do the calculation
1448 //
1449 for (Int_t ibin=1; ibin<=nbins; ibin++) {
1450
1451 // Calculate the value
1452 // physics = [ reco - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) ] / bin-width
1453 // In HIC : physics = [ reco - ( Tab * Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th * Rb ) ] / bin-width
1454 //
1455 //
1456 Double_t frac = 1.0, errfrac =0.;
1287e14b 1457
1458 // Variables initialization
1459 value = 0.; errvalue = 0.; errvalueMax = 0.; errvalueMin = 0.; kfactor = 0.;
1460 errvalueExtremeMax = 0.; errvalueExtremeMin = 0.;
1461 correction=0; correctionMax=0.; correctionMin=0.;
1462 correctionUncStatEffc=0.; correctionUncStatEffb=0.;
1463
b188dc47 1464 if(fPbPbElossHypothesis) {
ddd86f95 1465 frac = fTab[0]*fNevts;
1466 if(fIsEventPlane) frac = frac/2.0;
b188dc47 1467 errfrac = frac * TMath::Sqrt( (fTab[1]/fTab[0])*(fTab[1]/fTab[0]) + (1/fNevts) );
1468 } else {
1469 frac = fLuminosity[0];
1470 errfrac = fLuminosity[1];
1471 }
1472
1473 value = ( fhRECpt->GetBinContent(ibin)>0. && fhRECpt->GetBinContent(ibin)!=0. &&
1474 fhFeedDownMCpt->GetBinContent(ibin)>0. && fhFeedDownEffpt->GetBinContent(ibin)>0. ) ?
ddd86f95 1475 fhRECpt->GetBinContent(ibin) - frac*(deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) * fhRECpt->GetBinWidth(ibin) )
b188dc47 1476 : 0. ;
1477 value /= fhRECpt->GetBinWidth(ibin);
1478 if (value<0.) value =0.;
1479
1480 // Statistical uncertainty: delta_physics = sqrt ( (delta_reco)^2 ) / bin-width
1481 errvalue = (value!=0. && fhRECpt->GetBinError(ibin) && fhRECpt->GetBinError(ibin)!=0.) ?
1482 fhRECpt->GetBinError(ibin) : 0.;
1483 errvalue /= fhRECpt->GetBinWidth(ibin);
1484
1485 // Correction (fc) : Estimate of the relative amount feed-down subtracted
1486 // correction = [ 1 - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th)/reco ]
1487 // in HIC: correction = [ 1 - ( Tab * Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th)/reco ]
1488 correction = (value>0.) ?
1489 1 - (frac*deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) * fhRECpt->GetBinWidth(ibin) ) / fhRECpt->GetBinContent(ibin) : 0. ;
1490 if (correction<0.) correction = 0.;
1491
ddd86f95 1492 // 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;
1493
b188dc47 1494 // Systematic uncertainties
1495 // (syst but feed-down) delta_physics = sqrt ( (delta_reco_syst)^2 ) / bin-width
1496 // (feed-down syst) delta_physics = sqrt ( (k*delta_lumi/lumi)^2 + (k*delta_eff_trig/eff_trig)^2
1497 // + (k*delta_Nb/Nb)^2 + (k*delta_eff/eff)^2 + (k*global_eff_ratio)^2 ) / bin-width
1498 // where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th * bin-width
1499 kfactor = frac*deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) * fhRECpt->GetBinWidth(ibin) ;
1500 //
1501 if (fAsymUncertainties && value>0.) {
1502 Double_t nb = fhFeedDownMCpt->GetBinContent(ibin);
1503 Double_t nbDmax = fhFeedDownMCptMax->GetBinContent(ibin) - fhFeedDownMCpt->GetBinContent(ibin);
1504 Double_t nbDmin = fhFeedDownMCpt->GetBinContent(ibin) - fhFeedDownMCptMin->GetBinContent(ibin);
1505
1506 // Systematics but feed-down
1507 if (fgRECSystematics){
1508 errvalueMax = fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinWidth(ibin) ;
1509 errvalueMin = fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinWidth(ibin);
1510 }
1511 else { errvalueMax = 0.; errvalueMin = 0.; }
1512
1513 // Feed-down systematics
1514 // min value with the maximum Nb
1515 Double_t errCom = ( (kfactor*errfrac/frac)*(kfactor*errfrac/frac) ) +
1516 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1517 ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
1518 ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) ) ;
1519 errvalueExtremeMin = TMath::Sqrt( errCom + ( (kfactor*nbDmax/nb)*(kfactor*nbDmax/nb) ) ) / fhRECpt->GetBinWidth(ibin);
1520 // max value with the minimum Nb
1521 errvalueExtremeMax = TMath::Sqrt( errCom + ( (kfactor*nbDmin/nb)*(kfactor*nbDmin/nb) ) ) / fhRECpt->GetBinWidth(ibin);
1522
1523 // Correction systematics (fc)
1524 // min value with the maximum Nb
1525 correctionMin = TMath::Sqrt( errCom + ( (kfactor*nbDmax/nb)*(kfactor*nbDmax/nb) ) ) / fhRECpt->GetBinContent(ibin) ;
1526 // max value with the minimum Nb
1527 correctionMax = TMath::Sqrt( errCom + ( (kfactor*nbDmin/nb)*(kfactor*nbDmin/nb) ) ) / fhRECpt->GetBinContent(ibin) ;
1528 //
1529 correctionUncStatEffb = TMath::Sqrt( ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) )
1530 ) / fhRECpt->GetBinContent(ibin) ;
1531 correctionUncStatEffc = 0.;
1532 }
1533 else{ // Don't consider Nb uncertainty in this case [ to be tested!!! ]
1534 errvalueExtremeMax = TMath::Sqrt( ( (kfactor*errfrac/frac)*(kfactor*errfrac/frac) ) +
1535 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1536 ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
1537 ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1538 ) / fhRECpt->GetBinWidth(ibin);
1539 errvalueExtremeMin = errvalueExtremeMax ;
1540 }
1541
1542
1543 // fill in histograms
1544 fhYieldCorr->SetBinContent(ibin,value);
1545 fhYieldCorr->SetBinError(ibin,errvalue);
1546 //
1547 // Estimate how the result varies vs charm/beauty Eloss hypothesis
1548 //
ddd86f95 1549 if ( correction>1.0e-16 && fPbPbElossHypothesis){
b188dc47 1550 // Loop over the Eloss hypothesis
1551 // Int_t rbin=0;
1552 for (Float_t rval=0.0025; rval<4.0; rval+=0.005){
1553 // correction = [ 1 - (Tab *Nevt * delta_y * BR_b * eff_trig * eff_b * Nb_th *binwidth )* (rval) /reco ]
38a106d5 1554 Double_t fcRcbvalue = 1 - (frac*deltaY*branchingRatioBintoFinalDecay*fParticleAntiParticle*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin)*fhRECpt->GetBinWidth(ibin) * rval ) / fhRECpt->GetBinContent(ibin) ;
ddd86f95 1555 // cout << " rval="<<rval <<" , fc="<<fcRcbvalue<<" ";
1556 if(fcRcbvalue<1.e-16) fcRcbvalue=0.;
b188dc47 1557 fhFcRcb->Fill( fhRECpt->GetBinCenter(ibin) , rval, fcRcbvalue );
1558 // physics = reco * fcRcb / bin-width
1559 Double_t Rcbvalue = (fhRECpt->GetBinContent(ibin) && fcRcbvalue) ?
1560 fhRECpt->GetBinContent(ibin) * fcRcbvalue : 0. ;
1561 Rcbvalue /= fhRECpt->GetBinWidth(ibin) ;
1562 fhYieldCorrRcb->Fill( fhYieldCorr->GetBinCenter(ibin) , rval, Rcbvalue );
1563 // if(ibin==3){
1564 // cout << " pt "<< fhFcRcb->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", rbin="<<rbin<<", fc-Rb-value="<< fcRcbvalue << ", yield-Rb-value="<< Rcbvalue <<endl;
1565 // cout << " pt "<< fhFcRcb->GetBinCenter(ibin) <<" bin "<< ibin<<" rval="<<rval<<", fc-Rb-value="<< fcRcbvalue << ", yield-Rb-value="<< Rcbvalue <<endl;
1566 // rbin++;
1567 // }
1568 Double_t correctionMaxRcb = correctionMax*rval;
1569 Double_t correctionMinRcb = correctionMin*rval;
1570 fnHypothesis->Fill( fhYieldCorr->GetBinCenter(ibin), rval, fcRcbvalue, fcRcbvalue + correctionMaxRcb, fcRcbvalue - correctionMinRcb);
1571// if(ibin==3){
1572// cout << " pt="<< fhFcRcb->GetBinCenter(ibin) <<", rval="<<rval<<", fc="<<fcRcbvalue<<" +"<<correctionMaxRcb<<" -"<<correctionMinRcb<<endl;}
1573 }
1574 }
1575 //
1576 // Fill the rest of (asymmetric) histograms
1577 //
1578 if (fAsymUncertainties) {
1579 Double_t x = fhYieldCorr->GetBinCenter(ibin);
1580 fgYieldCorr->SetPoint(ibin,x,value); // i,x,y
1581 fgYieldCorr->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
1582 fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax);
1583 fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
1584 fgYieldCorrExtreme->SetPoint(ibin,x,value); // i,x,y
1585 fgYieldCorrExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
1586 fgYieldCorrConservative->SetPoint(ibin,x,value); // i,x,y
1587 fgYieldCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
1588 // cout << " bin " << ibin << ", correction " << correction << ", min correction unc " << correctionMin << ", max correction unc " << correctionMax << endl;
1589 if(correction>0.){
1590 fgFcConservative->SetPoint(ibin,x,correction);
1591 fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),correctionMin,correctionMax);
1592
1593 fhStatUncEffbFD->SetBinContent(ibin,0.); fhStatUncEffbFD->SetBinError(ibin,correctionUncStatEffb/correction*100.);
1594 fhStatUncEffcFD->SetBinContent(ibin,0.); fhStatUncEffcFD->SetBinError(ibin,correctionUncStatEffc/correction*100.);
1595 // cout << " pt "<< fhStatUncEffcFD->GetBinCenter(ibin) <<" bin "<< ibin<<" fc-stat-c ="<< correctionUncStatEffc/correction <<" fc-stat-b ="<< correctionUncStatEffb/correction <<endl;
1596 }
1597 else{
1598 fgFcConservative->SetPoint(ibin,x,0.);
1599 fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.);
1600 }
1601 }
1602
1603 }
1604 delete [] binwidths;
1605 delete [] limits;
1606
1607}
1608
1609
1610//_________________________________________________________________________________________________________
1611void AliHFPtSpectrum::ComputeSystUncertainties(AliHFSystErr *systematics, Bool_t combineFeedDown) {
1612 //
1613 // Function that re-calculates the global systematic uncertainties
1614 // by calling the class AliHFSystErr and combining those
1615 // (in quadrature) with the feed-down subtraction uncertainties
1616 //
1617
1618 // Estimate the feed-down uncertainty in percentage
c7d86f5e 1619 Int_t nentries = 0;
59fca4fa 1620 TGraphAsymmErrors *grErrFeeddown = 0;
b188dc47 1621 Double_t x=0., y=0., errx=0., erryl=0., erryh=0;
c7d86f5e 1622 if(fFeedDownOption!=0) {
1623 nentries = fgSigmaCorrConservative->GetN();
1624 grErrFeeddown = new TGraphAsymmErrors(nentries);
1625 for(Int_t i=0; i<nentries; i++) {
1626 x=0.; y=0.; errx=0.; erryl=0.; erryh=0.;
1627 fgSigmaCorrConservative->GetPoint(i,x,y);
1628 if(y>0.){
1629 errx = fgSigmaCorrConservative->GetErrorXlow(i) ;
1630 erryl = fgSigmaCorrConservative->GetErrorYlow(i) / y ;
1631 erryh = fgSigmaCorrConservative->GetErrorYhigh(i) / y ;
1632 }
1633 // cout << " x "<< x << " +- "<<errx<<" , y "<<y<<" + "<<erryh<<" - "<<erryl<<endl;
1634 grErrFeeddown->SetPoint(i,x,0.);
1635 grErrFeeddown->SetPointError(i,errx,errx,erryl,erryh); //i, xl, xh, yl, yh
b188dc47 1636 }
b188dc47 1637 }
1638
1639 // Draw all the systematics independently
1640 systematics->DrawErrors(grErrFeeddown);
1641
1642 // Set the sigma systematic uncertainties
1643 // possibly combine with the feed-down uncertainties
1644 Double_t errylcomb=0., erryhcomb=0;
1645 for(Int_t i=1; i<nentries; i++) {
1646 fgSigmaCorr->GetPoint(i,x,y);
1647 errx = grErrFeeddown->GetErrorXlow(i) ;
1648 erryl = grErrFeeddown->GetErrorYlow(i);
1649 erryh = grErrFeeddown->GetErrorYhigh(i);
1650 if (combineFeedDown) {
1651 errylcomb = systematics->GetTotalSystErr(x,erryl) * y ;
1652 erryhcomb = systematics->GetTotalSystErr(x,erryh) * y ;
1653 } else {
1654 errylcomb = systematics->GetTotalSystErr(x) * y ;
1655 erryhcomb = systematics->GetTotalSystErr(x) * y ;
1656 }
1657 fgSigmaCorr->SetPointError(i,errx,errx,errylcomb,erryhcomb);
1658 //
1659 fhSigmaCorrDataSyst->SetBinContent(i,y);
1660 erryl = systematics->GetTotalSystErr(x) * y ;
1661 fhSigmaCorrDataSyst->SetBinError(i,erryl);
1662 }
1663
1664}
1665
1666
1667//_________________________________________________________________________________________________________
1668void AliHFPtSpectrum::DrawSpectrum(TGraphAsymmErrors *gPrediction) {
1669 //
1670 // Example method to draw the corrected spectrum & the theoretical prediction
1671 //
1672
1673 TCanvas *csigma = new TCanvas("csigma","Draw the corrected cross-section & the prediction");
1674 csigma->SetFillColor(0);
1675 gPrediction->GetXaxis()->SetTitleSize(0.05);
1676 gPrediction->GetXaxis()->SetTitleOffset(0.95);
1677 gPrediction->GetYaxis()->SetTitleSize(0.05);
1678 gPrediction->GetYaxis()->SetTitleOffset(0.95);
1679 gPrediction->GetXaxis()->SetTitle("p_{T} [GeV]");
1680 gPrediction->GetYaxis()->SetTitle("BR #times #frac{d#sigma}{dp_{T}} |_{|y|<0.5} [pb/GeV]");
1681 gPrediction->SetLineColor(kGreen+2);
1682 gPrediction->SetLineWidth(3);
1683 gPrediction->SetFillColor(kGreen+1);
1684 gPrediction->Draw("3CA");
1685 fgSigmaCorr->SetLineColor(kRed);
1686 fgSigmaCorr->SetLineWidth(1);
1687 fgSigmaCorr->SetFillColor(kRed);
1688 fgSigmaCorr->SetFillStyle(0);
1689 fgSigmaCorr->Draw("2");
1690 fhSigmaCorr->SetMarkerColor(kRed);
1691 fhSigmaCorr->Draw("esame");
1692 csigma->SetLogy();
1693 TLegend * leg = new TLegend(0.7,0.75,0.87,0.5);
1694 leg->SetBorderSize(0);
1695 leg->SetLineColor(0);
1696 leg->SetFillColor(0);
1697 leg->SetTextFont(42);
1698 leg->AddEntry(gPrediction,"FONLL ","fl");
1699 leg->AddEntry(fhSigmaCorr,"data stat. unc.","pl");
1700 leg->AddEntry(fgSigmaCorr,"data syst. unc.","f");
1701 leg->Draw();
1702 csigma->Draw();
1703
1704}
1705
1706//_________________________________________________________________________________________________________
1707TH1D * AliHFPtSpectrum::ReweightHisto(TH1D *hToReweight, TH1D *hReference){
1708 //
1709 // Function to reweight histograms for testing purposes:
1710 // This function takes the histo hToReweight and reweights
1711 // it (its pt shape) with respect to hReference
1712 //
1713
1714 // check histograms consistency
1715 Bool_t areconsistent=kTRUE;
1716 areconsistent &= CheckHistosConsistency(hToReweight,hReference);
1717 if (!areconsistent) {
1718 AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
1719 return NULL;
1720 }
1721
1722 // define a new empty histogram
1723 TH1D *hReweighted = (TH1D*)hToReweight->Clone("hReweighted");
1724 hReweighted->Reset();
1725 Double_t weight=1.0;
1726 Double_t yvalue=1.0;
1727 Double_t integralRef = hReference->Integral();
1728 Double_t integralH = hToReweight->Integral();
1729
1730 // now reweight the spectra
1731 //
1732 // the weight is the relative probability of the given pt bin in the reference histo
1733 // divided by its relative probability (to normalize it) on the histo to re-weight
1734 for (Int_t i=0; i<=hToReweight->GetNbinsX(); i++) {
1735 weight = (hReference->GetBinContent(i)/integralRef) / (hToReweight->GetBinContent(i)/integralH) ;
1736 yvalue = hToReweight->GetBinContent(i);
1737 hReweighted->SetBinContent(i,yvalue*weight);
1738 }
1739
1740 return (TH1D*)hReweighted;
1741}
1742
1743//_________________________________________________________________________________________________________
1744TH1D * AliHFPtSpectrum::ReweightRecHisto(TH1D *hRecToReweight, TH1D *hMCToReweight, TH1D *hMCReference){
1745 //
1746 // Function to reweight histograms for testing purposes:
1747 // This function takes the histo hToReweight and reweights
1748 // it (its pt shape) with respect to hReference /hMCToReweight
1749 //
1750
1751 // check histograms consistency
1752 Bool_t areconsistent=kTRUE;
1753 areconsistent &= CheckHistosConsistency(hMCToReweight,hMCReference);
1754 areconsistent &= CheckHistosConsistency(hRecToReweight,hMCReference);
1755 if (!areconsistent) {
1756 AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
1757 return NULL;
1758 }
1759
1760 // define a new empty histogram
1761 TH1D *hReweighted = (TH1D*)hMCToReweight->Clone("hReweighted");
1762 hReweighted->Reset();
1763 TH1D *hRecReweighted = (TH1D*)hRecToReweight->Clone("hRecReweighted");
1764 hRecReweighted->Reset();
1765 Double_t weight=1.0;
1766 Double_t yvalue=1.0, yrecvalue=1.0;
1767 Double_t integralRef = hMCReference->Integral();
1768 Double_t integralH = hMCToReweight->Integral();
1769
1770 // now reweight the spectra
1771 //
1772 // the weight is the relative probability of the given pt bin
1773 // that should be applied in the MC histo to get the reference histo shape
1774 // Probabilities are properly normalized.
1775 for (Int_t i=0; i<=hMCToReweight->GetNbinsX(); i++) {
1776 weight = (hMCReference->GetBinContent(i)/integralRef) / (hMCToReweight->GetBinContent(i)/integralH) ;
1777 yvalue = hMCToReweight->GetBinContent(i);
1778 hReweighted->SetBinContent(i,yvalue*weight);
1779 yrecvalue = hRecToReweight->GetBinContent(i);
1780 hRecReweighted->SetBinContent(i,yrecvalue*weight);
1781 }
1782
1783 return (TH1D*)hRecReweighted;
1784}
1785
1786
1787
1788//_________________________________________________________________________________________________________
1789Int_t AliHFPtSpectrum::FindTH2YBin(TH2D *histo, Float_t yvalue){
1790 //
1791 // Function to find the y-axis bin of a TH2 for a given y-value
1792 //
1793
1794 Int_t nbins = histo->GetNbinsY();
1795 Int_t ybin=0;
1796 for (int j=0; j<=nbins; j++) {
1797 Float_t value = histo->GetYaxis()->GetBinCenter(j);
1798 Float_t width = histo->GetYaxis()->GetBinWidth(j);
1799 // if( TMath::Abs(yvalue-value)<= width/2. ) {
1800 if( TMath::Abs(yvalue-value)<= width ) {
1801 ybin =j;
1802 // cout <<" value "<<value << ", yval "<< yvalue<<", bin width "<<width/2.<< " y ="<<ybin<<endl;
1803 break;
1804 }
1805 }
1806
1807 return ybin;
1808}
1809
1810//_________________________________________________________________________________________________________
1811void AliHFPtSpectrum::ResetStatUncEff(){
1812
1813 Int_t entries = fhDirectEffpt->GetNbinsX();
1814 for(Int_t i=0; i<=entries; i++){
1815 fhDirectEffpt->SetBinError(i,0.);
1816 fhFeedDownEffpt->SetBinError(i,0.);
1817 }
1818
1819}