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