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