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