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