]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/vertexingHF/AliHFPtSpectrum.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fix compilation error (Chiara, Giacomo)
[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
16//***********************************************************************
17// Class AliHFPtSpectrum
18// Base class for feed-down corrections on heavy-flavour decays
19// computes the cross-section via one of the three implemented methods:
20// 0) Consider no feed-down prediction
21// 1) Subtract the feed-down with the "fc" method
22// Yield = Reco * fc; where fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) ) ;
23// 2) Subtract the feed-down with the "Nb" method
24// Yield = Reco - Feed-down (exact formula on the function implementation)
25//
bb427707 26// (the corrected yields per bin are divided by the bin-width)
27//
65e55bbd 28// Author: Z.Conesa, zconesa@in2p3.fr
29//***********************************************************************
30
31#include <Riostream.h>
8998180c 32
33#include "TMath.h"
34#include "TH1.h"
35#include "TH1D.h"
36#include "TGraphAsymmErrors.h"
bb427707 37#include "TNamed.h"
8998180c 38#include "TCanvas.h"
39#include "TLegend.h"
65e55bbd 40
41#include "AliLog.h"
8998180c 42#include "AliHFSystErr.h"
65e55bbd 43#include "AliHFPtSpectrum.h"
44
45ClassImp(AliHFPtSpectrum)
46
47//_________________________________________________________________________________________________________
48AliHFPtSpectrum::AliHFPtSpectrum(const char* name, const char* title, Int_t option):
49 TNamed(name,title),
86bdcd8c 50 fhDirectMCpt(NULL),
51 fhFeedDownMCpt(NULL),
52 fhDirectMCptMax(NULL),
53 fhDirectMCptMin(NULL),
54 fhFeedDownMCptMax(NULL),
55 fhFeedDownMCptMin(NULL),
56 fhDirectEffpt(NULL),
57 fhFeedDownEffpt(NULL),
58 fhRECpt(NULL),
59 fgRECSystematics(NULL),
65e55bbd 60 fLuminosity(),
61 fTrigEfficiency(),
8998180c 62 fGlobalEfficiencyUncertainties(),
86bdcd8c 63 fhFc(NULL),
64 fhFcMax(NULL),
65 fhFcMin(NULL),
66 fgFcExtreme(NULL),
67 fgFcConservative(NULL),
68 fhYieldCorr(NULL),
69 fhYieldCorrMax(NULL),
70 fhYieldCorrMin(NULL),
71 fgYieldCorr(NULL),
72 fgYieldCorrExtreme(NULL),
73 fgYieldCorrConservative(NULL),
74 fhSigmaCorr(NULL),
75 fhSigmaCorrMax(NULL),
76 fhSigmaCorrMin(NULL),
77 fgSigmaCorr(NULL),
78 fgSigmaCorrExtreme(NULL),
79 fgSigmaCorrConservative(NULL),
65e55bbd 80 fFeedDownOption(option),
81 fAsymUncertainties(kTRUE)
82{
83 //
84 // Default constructor
85 //
86
87 fLuminosity[0]=1.; fLuminosity[1]=0.;
88 fTrigEfficiency[0]=1.; fTrigEfficiency[1]=0.;
8998180c 89 fGlobalEfficiencyUncertainties[0]=0.; fGlobalEfficiencyUncertainties[1]=0.;
65e55bbd 90
91}
92
93//_________________________________________________________________________________________________________
94AliHFPtSpectrum::AliHFPtSpectrum(const AliHFPtSpectrum &rhs):
95 TNamed(rhs),
96 fhDirectMCpt(rhs.fhDirectMCpt),
97 fhFeedDownMCpt(rhs.fhFeedDownMCpt),
a17b17dd 98 fhDirectMCptMax(rhs.fhDirectMCptMax),
99 fhDirectMCptMin(rhs.fhDirectMCptMin),
100 fhFeedDownMCptMax(rhs.fhFeedDownMCptMax),
101 fhFeedDownMCptMin(rhs.fhFeedDownMCptMin),
65e55bbd 102 fhDirectEffpt(rhs.fhDirectEffpt),
103 fhFeedDownEffpt(rhs.fhFeedDownEffpt),
104 fhRECpt(rhs.fhRECpt),
86bdcd8c 105 fgRECSystematics(rhs.fgRECSystematics),
65e55bbd 106 fLuminosity(),
107 fTrigEfficiency(),
8998180c 108 fGlobalEfficiencyUncertainties(),
65e55bbd 109 fhFc(rhs.fhFc),
a17b17dd 110 fhFcMax(rhs.fhFcMax),
111 fhFcMin(rhs.fhFcMin),
86bdcd8c 112 fgFcExtreme(rhs.fgFcExtreme),
113 fgFcConservative(rhs.fgFcConservative),
65e55bbd 114 fhYieldCorr(rhs.fhYieldCorr),
a17b17dd 115 fhYieldCorrMax(rhs.fhYieldCorrMax),
116 fhYieldCorrMin(rhs.fhYieldCorrMin),
65e55bbd 117 fgYieldCorr(rhs.fgYieldCorr),
86bdcd8c 118 fgYieldCorrExtreme(rhs.fgYieldCorrExtreme),
119 fgYieldCorrConservative(rhs.fgYieldCorrConservative),
65e55bbd 120 fhSigmaCorr(rhs.fhSigmaCorr),
a17b17dd 121 fhSigmaCorrMax(rhs.fhSigmaCorrMax),
122 fhSigmaCorrMin(rhs.fhSigmaCorrMin),
65e55bbd 123 fgSigmaCorr(rhs.fgSigmaCorr),
86bdcd8c 124 fgSigmaCorrExtreme(rhs.fgSigmaCorrExtreme),
125 fgSigmaCorrConservative(rhs.fgSigmaCorrConservative),
65e55bbd 126 fFeedDownOption(rhs.fFeedDownOption),
127 fAsymUncertainties(rhs.fAsymUncertainties)
128{
129 //
130 // Copy constructor
131 //
132
133 for(Int_t i=0; i<2; i++){
134 fLuminosity[i] = rhs.fLuminosity[i];
135 fTrigEfficiency[i] = rhs.fTrigEfficiency[i];
8998180c 136 fGlobalEfficiencyUncertainties[i] = rhs.fGlobalEfficiencyUncertainties[i];
65e55bbd 137 }
138
139}
140
141//_________________________________________________________________________________________________________
142AliHFPtSpectrum &AliHFPtSpectrum::operator=(const AliHFPtSpectrum &source){
143 //
144 // Assignment operator
145 //
146
147 if (&source == this) return *this;
148
149 fhDirectMCpt = source.fhDirectMCpt;
150 fhFeedDownMCpt = source.fhFeedDownMCpt;
a17b17dd 151 fhDirectMCptMax = source.fhDirectMCptMax;
152 fhDirectMCptMin = source.fhDirectMCptMin;
153 fhFeedDownMCptMax = source.fhFeedDownMCptMax;
154 fhFeedDownMCptMin = source.fhFeedDownMCptMin;
65e55bbd 155 fhDirectEffpt = source.fhDirectEffpt;
156 fhFeedDownEffpt = source.fhFeedDownEffpt;
157 fhRECpt = source.fhRECpt;
86bdcd8c 158 fgRECSystematics = source.fgRECSystematics;
65e55bbd 159 fhFc = source.fhFc;
a17b17dd 160 fhFcMax = source.fhFcMax;
161 fhFcMin = source.fhFcMin;
86bdcd8c 162 fgFcExtreme = source.fgFcExtreme;
163 fgFcConservative = source.fgFcConservative;
65e55bbd 164 fhYieldCorr = source.fhYieldCorr;
a17b17dd 165 fhYieldCorrMax = source.fhYieldCorrMax;
166 fhYieldCorrMin = source.fhYieldCorrMin;
65e55bbd 167 fgYieldCorr = source.fgYieldCorr;
86bdcd8c 168 fgYieldCorrExtreme = source.fgYieldCorrExtreme;
169 fgYieldCorrConservative = source.fgYieldCorrConservative;
65e55bbd 170 fhSigmaCorr = source.fhSigmaCorr;
a17b17dd 171 fhSigmaCorrMax = source.fhSigmaCorrMax;
172 fhSigmaCorrMin = source.fhSigmaCorrMin;
65e55bbd 173 fgSigmaCorr = source.fgSigmaCorr;
86bdcd8c 174 fgSigmaCorrExtreme = source.fgSigmaCorrExtreme;
175 fgSigmaCorrConservative = source.fgSigmaCorrConservative;
65e55bbd 176 fFeedDownOption = source.fFeedDownOption;
177 fAsymUncertainties = source.fAsymUncertainties;
178
179 for(Int_t i=0; i<2; i++){
180 fLuminosity[i] = source.fLuminosity[i];
181 fTrigEfficiency[i] = source.fTrigEfficiency[i];
8998180c 182 fGlobalEfficiencyUncertainties[i] = source.fGlobalEfficiencyUncertainties[i];
65e55bbd 183 }
184
185 return *this;
186}
187
188//_________________________________________________________________________________________________________
189AliHFPtSpectrum::~AliHFPtSpectrum(){
190 //
191 // Destructor
192 //
86bdcd8c 193 if (fhDirectMCpt) delete fhDirectMCpt;
194 if (fhFeedDownMCpt) delete fhFeedDownMCpt;
195 if (fhDirectMCptMax) delete fhDirectMCptMax;
196 if (fhDirectMCptMin) delete fhDirectMCptMin;
197 if (fhFeedDownMCptMax) delete fhFeedDownMCptMax;
198 if (fhFeedDownMCptMin) delete fhFeedDownMCptMin;
199 if (fhDirectEffpt) delete fhDirectEffpt;
200 if (fhFeedDownEffpt) delete fhFeedDownEffpt;
201 if (fhRECpt) delete fhRECpt;
202 if (fgRECSystematics) delete fgRECSystematics;
203 if (fhFc) delete fhFc;
204 if (fhFcMax) delete fhFcMax;
205 if (fhFcMin) delete fhFcMin;
206 if (fgFcExtreme) delete fgFcExtreme;
207 if (fgFcConservative) delete fgFcConservative;
208 if (fhYieldCorr) delete fhYieldCorr;
209 if (fhYieldCorrMax) delete fhYieldCorrMax;
210 if (fhYieldCorrMin) delete fhYieldCorrMin;
211 if (fgYieldCorr) delete fgYieldCorr;
212 if (fgYieldCorrExtreme) delete fgYieldCorrExtreme;
213 if (fgYieldCorrConservative) delete fgYieldCorrConservative;
214 if (fhSigmaCorr) delete fhSigmaCorr;
215 if (fhSigmaCorrMax) delete fhSigmaCorrMax;
216 if (fhSigmaCorrMin) delete fhSigmaCorrMin;
217 if (fgSigmaCorr) delete fgSigmaCorr;
218 if (fgSigmaCorrExtreme) delete fgSigmaCorrExtreme;
219 if (fgSigmaCorrConservative) delete fgSigmaCorrConservative;
65e55bbd 220}
221
222
5f3c1b97 223//_________________________________________________________________________________________________________
86bdcd8c 224TH1D * AliHFPtSpectrum::RebinTheoreticalSpectra(TH1D *hTheory, const char *name) {
5f3c1b97 225 //
226 // Function to rebin the theoretical spectrum
227 // with respect to the real-data reconstructed spectrum binning
228 //
229
230 if (!hTheory || !fhRECpt) {
231 AliError("Feed-down or reconstructed spectra don't exist");
232 return NULL;
233 }
234
235 //
236 // Get the reconstructed spectra bins & limits
237 Int_t nbins = fhRECpt->GetNbinsX();
238 Int_t nbinsMC = hTheory->GetNbinsX();
239 Double_t *limits = new Double_t[nbins+1];
240 Double_t xlow=0., binwidth=0.;
86bdcd8c 241 for (Int_t i=1; i<=nbins; i++) {
5f3c1b97 242 binwidth = fhRECpt->GetBinWidth(i);
243 xlow = fhRECpt->GetBinLowEdge(i);
244 limits[i-1] = xlow;
245 }
246 limits[nbins] = xlow + binwidth;
247
248 //
249 // Define a new histogram with the real-data reconstructed spectrum binning
86bdcd8c 250 TH1D * hTheoryRebin = new TH1D(name," theoretical rebinned prediction",nbins,limits);
5f3c1b97 251
252 Double_t sum[nbins], items[nbins];
253 for (Int_t ibin=0; ibin<=nbinsMC; ibin++){
254
255 for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++){
86bdcd8c 256 if (hTheory->GetBinCenter(ibin)>limits[ibinrec] &&
5f3c1b97 257 hTheory->GetBinCenter(ibin)<limits[ibinrec+1]){
258 sum[ibinrec]+=hTheory->GetBinContent(ibin);
259 items[ibinrec]+=1.;
260 }
261 }
262
263 }
264
265 // set the theoretical rebinned spectra to ( sum-bins / n-bins ) per new bin
266 for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++) {
267 hTheoryRebin->SetBinContent(ibinrec+1,sum[ibinrec]/items[ibinrec]);
268 }
269
86bdcd8c 270 return (TH1D*)hTheoryRebin;
5f3c1b97 271}
272
65e55bbd 273//_________________________________________________________________________________________________________
86bdcd8c 274void AliHFPtSpectrum::SetMCptSpectra(TH1D *hDirect, TH1D *hFeedDown){
65e55bbd 275 //
276 // Set the MonteCarlo or Theoretical spectra
277 // both for direct and feed-down contributions
278 //
279
bb427707 280 if (!hDirect || !hFeedDown || !fhRECpt) {
281 AliError("One or both (direct, feed-down) spectra or the reconstructed spectra don't exist");
65e55bbd 282 return;
283 }
284
285 Bool_t areconsistent = kTRUE;
286 areconsistent = CheckHistosConsistency(hDirect,hFeedDown);
287 if (!areconsistent) {
288 AliInfo("Histograms are not consistent (bin width, bounds)");
289 return;
290 }
291
bb427707 292 //
293 // If the predictions shape do not have the same
294 // bin width (check via number of bins) as the reconstructed spectra, change them
8998180c 295 if (hDirect->GetBinWidth(1) == fhRECpt->GetBinWidth(1)) {
bb427707 296
86bdcd8c 297 fhDirectMCpt = (TH1D*)hDirect->Clone();
298 fhDirectMCpt->SetNameTitle("fhDirectMCpt"," direct theoretical prediction");
299 fhFeedDownMCpt = (TH1D*)hFeedDown->Clone();
300 fhFeedDownMCpt->SetNameTitle("fhFeedDownMCpt"," feed-down theoretical prediction");
bb427707 301
86bdcd8c 302 }
bb427707 303 else {
304
86bdcd8c 305 fhDirectMCpt = RebinTheoreticalSpectra(hDirect,"fhDirectMCpt");
5f3c1b97 306 fhDirectMCpt->SetNameTitle("fhDirectMCpt"," direct theoretical prediction");
86bdcd8c 307 fhFeedDownMCpt = RebinTheoreticalSpectra(hFeedDown,"fhFeedDownMCpt");
5f3c1b97 308 fhFeedDownMCpt->SetNameTitle("fhFeedDownMCpt"," feed-down theoretical prediction");
bb427707 309
bb427707 310 }
311
65e55bbd 312}
313
314//_________________________________________________________________________________________________________
86bdcd8c 315void AliHFPtSpectrum::SetFeedDownMCptSpectra(TH1D *hFeedDown){
65e55bbd 316 //
317 // Set the MonteCarlo or Theoretical spectra
318 // for feed-down contribution
319 //
320
bb427707 321 if (!hFeedDown || !fhRECpt) {
322 AliError("Feed-down or reconstructed spectra don't exist");
65e55bbd 323 return;
324 }
bb427707 325
bb427707 326 //
327 // If the predictions shape do not have the same
328 // bin width (check via number of bins) as the reconstructed spectra, change them
8998180c 329 if (hFeedDown->GetBinWidth(1) == fhRECpt->GetBinWidth(1)) {
bb427707 330
86bdcd8c 331 fhFeedDownMCpt = (TH1D*)hFeedDown->Clone();
332 fhFeedDownMCpt->SetNameTitle("fhFeedDownMCpt"," feed-down theoretical prediction");
bb427707 333
334 }
335 else {
336
86bdcd8c 337 fhFeedDownMCpt = RebinTheoreticalSpectra(hFeedDown,"fhFeedDownMCpt");
5f3c1b97 338 fhFeedDownMCpt->SetNameTitle("fhFeedDownMCpt"," feed-down theoretical prediction");
bb427707 339
bb427707 340 }
341
65e55bbd 342}
343
344//_________________________________________________________________________________________________________
86bdcd8c 345void AliHFPtSpectrum::SetMCptDistributionsBounds(TH1D *hDirectMax, TH1D *hDirectMin, TH1D *hFeedDownMax, TH1D *hFeedDownMin){
65e55bbd 346 //
347 // Set the maximum and minimum MonteCarlo or Theoretical spectra
348 // both for direct and feed-down contributions
349 // used in case uncertainties are asymmetric and ca not be on the "basic histograms"
350 //
351
bb427707 352 if (!hDirectMax || !hDirectMin || !hFeedDownMax|| !hFeedDownMin || !fhRECpt) {
353 AliError("One or all of the max/min direct/feed-down or the reconstructed spectra don't exist");
65e55bbd 354 return;
355 }
356
357 Bool_t areconsistent = kTRUE;
358 areconsistent &= CheckHistosConsistency(hDirectMax,hDirectMin);
359 areconsistent &= CheckHistosConsistency(hFeedDownMax,hFeedDownMin);
360 areconsistent &= CheckHistosConsistency(hDirectMax,hFeedDownMax);
361 if (!areconsistent) {
362 AliInfo("Histograms are not consistent (bin width, bounds)");
363 return;
364 }
365
bb427707 366
367 //
368 // If the predictions shape do not have the same
369 // bin width (check via number of bins) as the reconstructed spectra, change them
8998180c 370 if (hFeedDownMax->GetBinWidth(1) == fhRECpt->GetBinWidth(1)) {
bb427707 371
86bdcd8c 372 fhDirectMCptMax = (TH1D*)hDirectMax->Clone();
373 fhDirectMCptMin = (TH1D*)hDirectMin->Clone();
374 fhFeedDownMCptMax = (TH1D*)hFeedDownMax->Clone();
375 fhFeedDownMCptMin = (TH1D*)hFeedDownMin->Clone();
376 fhDirectMCptMax->SetNameTitle("fhDirectMCptMax"," maximum direct theoretical prediction");
377 fhDirectMCptMin->SetNameTitle("fhDirectMCptMin"," minimum direct theoretical prediction");
378 fhFeedDownMCptMax->SetNameTitle("fhFeedDownMCptMax"," maximum feed-down theoretical prediction");
379 fhFeedDownMCptMin->SetNameTitle("fhFeedDownMCptMin"," minimum feed-down theoretical prediction");
bb427707 380
381 }
382 else {
383
86bdcd8c 384 fhDirectMCptMax = RebinTheoreticalSpectra(hDirectMax,"fhDirectMCptMax");
5f3c1b97 385 fhDirectMCptMax->SetNameTitle("fhDirectMCptMax"," maximum direct theoretical prediction");
86bdcd8c 386 fhDirectMCptMin = RebinTheoreticalSpectra(hDirectMin,"fhDirectMCptMin");
5f3c1b97 387 fhDirectMCptMin->SetNameTitle("fhDirectMCptMin"," minimum direct theoretical prediction");
86bdcd8c 388 fhFeedDownMCptMax = RebinTheoreticalSpectra(hFeedDownMax,"fhFeedDownMCptMax");
5f3c1b97 389 fhFeedDownMCptMax->SetNameTitle("fhFeedDownMCptMax"," maximum feed-down theoretical prediction");
86bdcd8c 390 fhFeedDownMCptMin = RebinTheoreticalSpectra(hFeedDownMin,"fhFeedDownMCptMin");
5f3c1b97 391 fhFeedDownMCptMin->SetNameTitle("fhFeedDownMCptMin"," minimum feed-down theoretical prediction");
bb427707 392
bb427707 393 }
394
65e55bbd 395}
396
397//_________________________________________________________________________________________________________
86bdcd8c 398void AliHFPtSpectrum::SetFeedDownMCptDistributionsBounds(TH1D *hFeedDownMax, TH1D *hFeedDownMin){
65e55bbd 399 //
400 // Set the maximum and minimum MonteCarlo or Theoretical spectra
401 // for feed-down contributions
402 // used in case uncertainties are asymmetric and can not be on the "basic histogram"
403 //
404
bb427707 405 if (!hFeedDownMax || !hFeedDownMin || !fhRECpt) {
65e55bbd 406 AliError("One or all of the max/min direct/feed-down spectra don't exist");
407 return;
408 }
409
410 Bool_t areconsistent = kTRUE;
411 areconsistent &= CheckHistosConsistency(hFeedDownMax,hFeedDownMin);
412 if (!areconsistent) {
413 AliInfo("Histograms are not consistent (bin width, bounds)");
414 return;
415 }
416
bb427707 417
418 //
419 // If the predictions shape do not have the same
420 // bin width (check via number of bins) as the reconstructed spectra, change them
8998180c 421 if (hFeedDownMax->GetBinWidth(1) == fhRECpt->GetBinWidth(1)) {
bb427707 422
86bdcd8c 423 fhFeedDownMCptMax = (TH1D*)hFeedDownMax->Clone();
424 fhFeedDownMCptMin = (TH1D*)hFeedDownMin->Clone();
425 fhFeedDownMCptMax->SetNameTitle("fhFeedDownMCptMax"," maximum feed-down theoretical prediction");
426 fhFeedDownMCptMin->SetNameTitle("fhFeedDownMCptMin"," minimum feed-down theoretical prediction");
bb427707 427
428 }
429 else {
430
86bdcd8c 431 fhFeedDownMCptMax = RebinTheoreticalSpectra(hFeedDownMax,"fhFeedDownMCptMax");
5f3c1b97 432 fhFeedDownMCptMax->SetNameTitle("fhFeedDownMCptMax"," maximum feed-down theoretical prediction");
86bdcd8c 433 fhFeedDownMCptMin = RebinTheoreticalSpectra(hFeedDownMin,"fhFeedDownMCptMin");
5f3c1b97 434 fhFeedDownMCptMin->SetNameTitle("fhFeedDownMCptMin"," minimum feed-down theoretical prediction");
bb427707 435
bb427707 436 }
437
65e55bbd 438}
439
440//_________________________________________________________________________________________________________
86bdcd8c 441void AliHFPtSpectrum::SetDirectAccEffCorrection(TH1D *hDirectEff){
65e55bbd 442 //
443 // Set the Acceptance and Efficiency corrections
444 // for the direct contribution
445 //
446
447 if (!hDirectEff) {
448 AliError("The direct acceptance and efficiency corrections doesn't exist");
449 return;
450 }
451
86bdcd8c 452 fhDirectEffpt = (TH1D*)hDirectEff->Clone();
453 fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance x efficiency correction");
65e55bbd 454}
455
456//_________________________________________________________________________________________________________
86bdcd8c 457void AliHFPtSpectrum::SetAccEffCorrection(TH1D *hDirectEff, TH1D *hFeedDownEff){
65e55bbd 458 //
459 // Set the Acceptance and Efficiency corrections
460 // both for direct and feed-down contributions
461 //
462
463 if (!hDirectEff || !hFeedDownEff) {
464 AliError("One or both (direct, feed-down) acceptance and efficiency corrections don't exist");
465 return;
466 }
467
468 Bool_t areconsistent=kTRUE;
469 areconsistent = CheckHistosConsistency(hDirectEff,hFeedDownEff);
470 if (!areconsistent) {
471 AliInfo("Histograms are not consistent (bin width, bounds)");
472 return;
473 }
474
86bdcd8c 475 fhDirectEffpt = (TH1D*)hDirectEff->Clone();
476 fhFeedDownEffpt = (TH1D*)hFeedDownEff->Clone();
477 fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance x efficiency correction");
478 fhFeedDownEffpt->SetNameTitle("fhFeedDownEffpt"," feed-down acceptance x efficiency correction");
65e55bbd 479}
480
481//_________________________________________________________________________________________________________
86bdcd8c 482void AliHFPtSpectrum::SetReconstructedSpectrum(TH1D *hRec) {
65e55bbd 483 //
484 // Set the reconstructed spectrum
485 //
486
487 if (!hRec) {
488 AliError("The reconstructed spectrum doesn't exist");
489 return;
490 }
491
86bdcd8c 492 fhRECpt = (TH1D*)hRec->Clone();
493 fhRECpt->SetNameTitle("fhRECpt"," reconstructed spectrum");
494}
495
496//_________________________________________________________________________________________________________
497void AliHFPtSpectrum::SetReconstructedSpectrumSystematics(TGraphAsymmErrors *gRec) {
498 //
499 // Set the reconstructed spectrum systematic uncertainties
500 //
501
502 // Check the compatibility with the reconstructed spectrum
d38da1b0 503 Double_t gbinwidth = gRec->GetErrorXlow(1) + gRec->GetErrorXhigh(1) ;
86bdcd8c 504 Double_t hbinwidth = fhRECpt->GetBinWidth(1);
505 Double_t gxbincenter=0., gybincenter=0.;
d38da1b0 506 gRec->GetPoint(1,gxbincenter,gybincenter);
86bdcd8c 507 Double_t hbincenter = fhRECpt->GetBinCenter(1);
8998180c 508 if ( (gbinwidth != hbinwidth) || (gxbincenter!=hbincenter) ) {
86bdcd8c 509 AliError(" The reconstructed spectrum and its systematics don't seem compatible");
510 return;
511 }
512
513 fgRECSystematics = gRec;
65e55bbd 514}
515
516//_________________________________________________________________________________________________________
a17b17dd 517void AliHFPtSpectrum::ComputeHFPtSpectrum(Double_t deltaY, Double_t branchingRatioC, Double_t branchingRatioBintoFinalDecay) {
65e55bbd 518 //
519 // Main function to compute the corrected cross-section:
a17b17dd 520 // variables : analysed delta_y, BR for the final correction,
521 // BR b --> D --> decay (relative to the input theoretical prediction)
65e55bbd 522 //
523 // Sigma = ( 1. / (lumi * delta_y * BR_c * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
524 //
86bdcd8c 525 // Uncertainties: (stat) delta_sigma = sigma * sqrt ( (delta_spectra/spectra)^2 )
526 // (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 )
527 // (feed-down syst) delta_sigma = sigma * sqrt ( (delta_spectra_fd/spectra_fd)^2 )
65e55bbd 528
529 //
530 // First: Initialization
531 //
532 Bool_t areHistosOk = Initialize();
533 if (!areHistosOk) {
534 AliInfo(" Histos not properly initialized. Check : inconsistent binning ? missing histos ?");
535 return;
536 }
537
538 //
539 // Second: Correct for feed-down
540 //
541 if (fFeedDownOption==1) {
542 // Compute the feed-down correction via fc-method
a17b17dd 543 CalculateFeedDownCorrectionFc();
65e55bbd 544 // Correct the yield for feed-down correction via fc-method
a17b17dd 545 CalculateFeedDownCorrectedSpectrumFc();
65e55bbd 546 }
547 else if (fFeedDownOption==2) {
548 // Correct the yield for feed-down correction via Nb-method
a17b17dd 549 CalculateFeedDownCorrectedSpectrumNb(deltaY,branchingRatioBintoFinalDecay);
65e55bbd 550 }
551 else if (fFeedDownOption==0) {
552 // If there is no need for feed-down correction,
553 // the "corrected" yield is equal to the raw yield
86bdcd8c 554 fhYieldCorr = (TH1D*)fhRECpt->Clone();
65e55bbd 555 fhYieldCorr->SetNameTitle("fhYieldCorr","un-corrected yield");
86bdcd8c 556 fhYieldCorrMax = (TH1D*)fhRECpt->Clone();
557 fhYieldCorrMin = (TH1D*)fhRECpt->Clone();
a17b17dd 558 fhYieldCorrMax->SetNameTitle("fhYieldCorrMax","un-corrected yield");
559 fhYieldCorrMin->SetNameTitle("fhYieldCorrMin","un-corrected yield");
65e55bbd 560 fAsymUncertainties=kFALSE;
561 }
562 else {
563 AliInfo(" Are you sure the feed-down correction option is right ?");
564 }
565
566 // Print out information
8998180c 567 // 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\n",fLuminosity[0],fLuminosity[1],fTrigEfficiency[0],fTrigEfficiency[1],deltaY,branchingRatioC,branchingRatioBintoFinalDecay);
568 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]);
65e55bbd 569
570 //
571 // Finally: Correct from yields to cross-section
572 //
573 Int_t nbins = fhRECpt->GetNbinsX();
574 Double_t binwidth = fhRECpt->GetBinWidth(1);
bb427707 575 Double_t *limits = new Double_t[nbins+1];
de97013c 576 Double_t *binwidths = new Double_t[nbins+1];
bb427707 577 Double_t xlow=0.;
86bdcd8c 578 for (Int_t i=1; i<=nbins; i++) {
bb427707 579 binwidth = fhRECpt->GetBinWidth(i);
580 xlow = fhRECpt->GetBinLowEdge(i);
581 limits[i-1] = xlow;
de97013c 582 binwidths[i] = binwidth;
bb427707 583 }
584 limits[nbins] = xlow + binwidth;
585
65e55bbd 586
587 // declare the output histograms
86bdcd8c 588 if (!fhSigmaCorr) fhSigmaCorr = new TH1D("fhSigmaCorr","corrected sigma",nbins,limits);
589 if (!fhSigmaCorrMax) fhSigmaCorrMax = new TH1D("fhSigmaCorrMax","max corrected sigma",nbins,limits);
590 if (!fhSigmaCorrMin) fhSigmaCorrMin = new TH1D("fhSigmaCorrMin","min corrected sigma",nbins,limits);
65e55bbd 591 // and the output TGraphAsymmErrors
86bdcd8c 592 if (fAsymUncertainties & !fgSigmaCorr) fgSigmaCorr = new TGraphAsymmErrors(nbins+1);
593 if (fAsymUncertainties & !fgSigmaCorrExtreme) fgSigmaCorrExtreme = new TGraphAsymmErrors(nbins+1);
594 if (fAsymUncertainties & !fgSigmaCorrConservative) fgSigmaCorrConservative = new TGraphAsymmErrors(nbins+1);
65e55bbd 595
596 // protect against null denominator
8998180c 597 if (deltaY==0. || fLuminosity[0]==0. || fTrigEfficiency[0]==0. || branchingRatioC==0.) {
65e55bbd 598 AliError(" Hey you ! Why luminosity or trigger-efficiency or the c-BR or delta_y are set to zero ?! ");
599 return ;
600 }
601
86bdcd8c 602 Double_t value=0, errValue=0, errvalueMax=0., errvalueMin=0.;
603 Double_t errvalueExtremeMax=0., errvalueExtremeMin=0.;
604 Double_t errvalueConservativeMax=0., errvalueConservativeMin=0.;
65e55bbd 605 for(Int_t ibin=0; ibin<=nbins; ibin++){
606
86bdcd8c 607 // Sigma calculation
65e55bbd 608 // Sigma = ( 1. / (lumi * delta_y * BR_c * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
609 value = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ?
a17b17dd 610 ( fhYieldCorr->GetBinContent(ibin) / ( deltaY * branchingRatioC * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
65e55bbd 611 : 0. ;
86bdcd8c 612
613 // Sigma statistical uncertainty:
614 // delta_sigma = sigma * sqrt ( (delta_spectra/spectra)^2 )
615 errValue = (value!=0.) ? value * (fhYieldCorr->GetBinError(ibin)/fhYieldCorr->GetBinContent(ibin)) : 0. ;
616
617 //
618 // Sigma systematic uncertainties
619 //
65e55bbd 620 if (fAsymUncertainties) {
86bdcd8c 621
622 // (syst but feed-down) delta_sigma = sigma * sqrt ( (delta_spectra_syst/spectra)^2 +
8998180c 623 // (delta_lumi/lumi)^2 + (delta_eff_trig/eff_trig)^2 + (delta_eff/eff)^2 + (global_eff)^2 )
86bdcd8c 624 errvalueMax = value * TMath::Sqrt( (fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin))*(fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin)) +
625 (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
626 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) +
8998180c 627 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
628 fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] );
86bdcd8c 629 errvalueMin = value * TMath::Sqrt( (fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin))*(fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin)) +
630 (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
631 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) +
8998180c 632 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
633 fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] );
86bdcd8c 634
635 // Uncertainties from feed-down
636 // (feed-down syst) delta_sigma = sigma * sqrt ( (delta_spectra_fd/spectra_fd)^2 )
637 // extreme case
638 errvalueExtremeMax = value * (fgYieldCorrExtreme->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin));
639 errvalueExtremeMin = value * (fgYieldCorrExtreme->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin));
640 //
641 // conservative case
642 errvalueConservativeMax = value * (fgYieldCorrConservative->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin));
643 errvalueConservativeMin = value * (fgYieldCorrConservative->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin));
644
65e55bbd 645 }
646 else {
647 // protect against null denominator
86bdcd8c 648 errvalueMax = (value!=0.) ?
649 value * TMath::Sqrt( (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
650 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) +
8998180c 651 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
652 fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] )
65e55bbd 653 : 0. ;
86bdcd8c 654 errvalueMin = errvalueMax;
65e55bbd 655 }
656
657 // Fill the histograms
86bdcd8c 658 fhSigmaCorr->SetBinContent(ibin,value);
659 fhSigmaCorr->SetBinError(ibin,errValue);
65e55bbd 660 // Fill the TGraphAsymmErrors
661 if (fAsymUncertainties) {
662 Double_t x = fhYieldCorr->GetBinCenter(ibin);
663 fgSigmaCorr->SetPoint(ibin,x,value); // i,x,y
86bdcd8c 664 fgSigmaCorr->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
665 fhSigmaCorrMax->SetBinContent(ibin,value+errvalueMax);
666 fhSigmaCorrMin->SetBinContent(ibin,value-errvalueMin);
667 fgSigmaCorrExtreme->SetPoint(ibin,x,value); // i,x,y
668 fgSigmaCorrExtreme->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
669 fgSigmaCorrConservative->SetPoint(ibin,x,value); // i,x,y
670 fgSigmaCorrConservative->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),errvalueConservativeMin,errvalueConservativeMax); // i,xl,xh,yl,yh
671
65e55bbd 672 }
673
674 }
675
65e55bbd 676}
677
bb427707 678//_________________________________________________________________________________________________________
86bdcd8c 679TH1D * AliHFPtSpectrum::EstimateEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco, const char *name) {
bb427707 680 //
5f3c1b97 681 // Function that computes the acceptance and efficiency correction
bb427707 682 // based on the simulated and reconstructed spectra
683 // and using the reconstructed spectra bin width
684 //
685 // eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
686 //
687
688 if(!fhRECpt){
689 AliInfo("Hey, the reconstructed histogram was not set yet !");
690 return NULL;
691 }
692
693 Int_t nbins = fhRECpt->GetNbinsX();
694 Double_t *limits = new Double_t[nbins+1];
695 Double_t xlow=0.,binwidth=0.;
86bdcd8c 696 for (Int_t i=1; i<=nbins; i++) {
bb427707 697 binwidth = fhRECpt->GetBinWidth(i);
698 xlow = fhRECpt->GetBinLowEdge(i);
699 limits[i-1] = xlow;
700 }
701 limits[nbins] = xlow + binwidth;
702
86bdcd8c 703 TH1D * hEfficiency = new TH1D(name," acceptance #times efficiency",nbins,limits);
bb427707 704
705 Double_t sumSimu[nbins], sumReco[nbins];
706 for (Int_t ibin=0; ibin<=hSimu->GetNbinsX(); ibin++){
707
708 for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++){
709 if ( hSimu->GetBinCenter(ibin)>limits[ibinrec] &&
710 hSimu->GetBinCenter(ibin)<limits[ibinrec+1] ) {
711 sumSimu[ibinrec]+=hSimu->GetBinContent(ibin);
712 }
713 if ( hReco->GetBinCenter(ibin)>limits[ibinrec] &&
714 hReco->GetBinCenter(ibin)<limits[ibinrec+1] ) {
715 sumReco[ibinrec]+=hReco->GetBinContent(ibin);
716 }
717 }
718
719 }
720
86bdcd8c 721
bb427707 722 // the efficiency is computed as reco/sim (in each bin)
723 // its uncertainty is err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
724 Double_t eff=0., erreff=0.;
86bdcd8c 725 for (Int_t ibinrec=0; ibinrec<=nbins; ibinrec++) {
726 if (sumSimu[ibinrec]!= 0. && sumReco[ibinrec]!=0.) {
727 eff = sumReco[ibinrec] / sumSimu[ibinrec] ;
728 erreff = TMath::Sqrt( eff * (1.0 - eff) ) / TMath::Sqrt( sumSimu[ibinrec] );
729 }
730 else { eff=0.0; erreff=0.; }
5f3c1b97 731 hEfficiency->SetBinContent(ibinrec+1,eff);
732 hEfficiency->SetBinError(ibinrec+1,erreff);
bb427707 733 }
734
86bdcd8c 735 return (TH1D*)hEfficiency;
bb427707 736}
737
738//_________________________________________________________________________________________________________
86bdcd8c 739void AliHFPtSpectrum::EstimateAndSetDirectEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco) {
bb427707 740 //
5f3c1b97 741 // Function that computes the Direct acceptance and efficiency correction
bb427707 742 // based on the simulated and reconstructed spectra
743 // and using the reconstructed spectra bin width
744 //
745 // eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
746 //
5f3c1b97 747
86bdcd8c 748 if(!fhRECpt || !hSimu || !hReco){
749 AliError("Hey, the reconstructed histogram was not set yet !");
750 return;
bb427707 751 }
752
86bdcd8c 753 fhDirectEffpt = EstimateEfficiencyRecoBin(hSimu,hReco,"fhDirectEffpt");
5f3c1b97 754 fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance #times efficiency");
bb427707 755
5f3c1b97 756}
bb427707 757
5f3c1b97 758//_________________________________________________________________________________________________________
86bdcd8c 759void AliHFPtSpectrum::EstimateAndSetFeedDownEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco) {
5f3c1b97 760 //
761 // Function that computes the Feed-Down acceptance and efficiency correction
762 // based on the simulated and reconstructed spectra
763 // and using the reconstructed spectra bin width
764 //
765 // eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
766 //
bb427707 767
86bdcd8c 768 if(!fhRECpt || !hSimu || !hReco){
769 AliError("Hey, the reconstructed histogram was not set yet !");
770 return;
bb427707 771 }
772
86bdcd8c 773 fhFeedDownEffpt = EstimateEfficiencyRecoBin(hSimu,hReco,"fhFeedDownEffpt");
5f3c1b97 774 fhFeedDownEffpt->SetNameTitle("fhFeedDownEffpt"," feed-down acceptance #times efficiency");
775
bb427707 776}
777
65e55bbd 778//_________________________________________________________________________________________________________
779Bool_t AliHFPtSpectrum::Initialize(){
780 //
781 // Initialization of the variables (histograms)
782 //
783
784 if (fFeedDownOption==0) {
785 AliInfo("Getting ready for the corrections without feed-down consideration");
786 } else if (fFeedDownOption==1) {
787 AliInfo("Getting ready for the fc feed-down correction calculation");
788 } else if (fFeedDownOption==2) {
789 AliInfo("Getting ready for the Nb feed-down correction calculation");
790 } else { AliError("The calculation option must be <=2"); return kFALSE; }
791
792 // Start checking the input histograms consistency
793 Bool_t areconsistent=kTRUE;
794
795 // General checks
796 if (!fhDirectEffpt || !fhRECpt) {
797 AliError(" Reconstructed spectra and/or the Nc efficiency distributions are not defined");
798 return kFALSE;
799 }
800 areconsistent &= CheckHistosConsistency(fhRECpt,fhDirectEffpt);
801 if (!areconsistent) {
802 AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
803 return kFALSE;
804 }
805 if (fFeedDownOption==0) return kTRUE;
806
807 //
808 // Common checks for options 1 (fc) & 2(Nb)
809 if (!fhFeedDownMCpt || !fhFeedDownEffpt) {
810 AliError(" Theoretical Nb and/or the Nb efficiency distributions are not defined");
811 return kFALSE;
812 }
813 areconsistent &= CheckHistosConsistency(fhRECpt,fhFeedDownMCpt);
814 areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownEffpt);
815 if (fAsymUncertainties) {
a17b17dd 816 if (!fhFeedDownMCptMax || !fhFeedDownMCptMin) {
65e55bbd 817 AliError(" Max/Min theoretical Nb distributions are not defined");
818 return kFALSE;
819 }
a17b17dd 820 areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownMCptMax);
65e55bbd 821 }
822 if (!areconsistent) {
823 AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
824 return kFALSE;
825 }
826 if (fFeedDownOption>1) return kTRUE;
827
828 //
829 // Now checks for option 1 (fc correction)
830 if (!fhDirectMCpt) {
831 AliError("Theoretical Nc distributions is not defined");
832 return kFALSE;
833 }
834 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhFeedDownMCpt);
835 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectEffpt);
836 if (fAsymUncertainties) {
a17b17dd 837 if (!fhDirectMCptMax || !fhDirectMCptMin) {
65e55bbd 838 AliError(" Max/Min theoretical Nc distributions are not defined");
839 return kFALSE;
840 }
a17b17dd 841 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectMCptMax);
65e55bbd 842 }
843 if (!areconsistent) {
844 AliInfo("Histograms required for fc correction are not consistent (bin width, bounds)");
845 return kFALSE;
846 }
847
848 return kTRUE;
849}
850
851//_________________________________________________________________________________________________________
86bdcd8c 852Bool_t AliHFPtSpectrum::CheckHistosConsistency(TH1D *h1, TH1D *h2){
65e55bbd 853 //
854 // Check the histograms consistency (bins, limits)
855 //
856
857 if (!h1 || !h2) {
858 AliError("One or both histograms don't exist");
859 return kFALSE;
860 }
861
862 Double_t binwidth1 = h1->GetBinWidth(1);
863 Double_t binwidth2 = h2->GetBinWidth(1);
864 Double_t min1 = h1->GetBinCenter(1) - (binwidth1/2.) ;
865// Double_t max1 = h1->GetBinCenter(nbins1) + (binwidth1/2.) ;
866 Double_t min2 = h2->GetBinCenter(1) - (binwidth2/2.) ;
867// Double_t max2 = h2->GetBinCenter(nbins2) + (binwidth2/2.) ;
868
869 if (binwidth1!=binwidth2) {
870 AliInfo(" histograms with different bin width");
871 return kFALSE;
872 }
873 if (min1!=min2) {
874 AliInfo(" histograms with different minimum");
875 return kFALSE;
876 }
877// if (max1!=max2) {
878// AliInfo(" histograms with different maximum");
879// return kFALSE;
880// }
881
882 return kTRUE;
883}
884
885//_________________________________________________________________________________________________________
a17b17dd 886void AliHFPtSpectrum::CalculateFeedDownCorrectionFc(){
65e55bbd 887 //
888 // Compute fc factor and its uncertainties bin by bin
889 // fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
890 //
86bdcd8c 891 // uncertainties: (conservative) combine the upper/lower N_b & N_c predictions together
892 // (extreme) combine the upper N_b predictions with the lower N_c predictions & viceversa
8998180c 893 // systematic uncertainty on the acceptance x efficiency b/c ratio are included
86bdcd8c 894 //
65e55bbd 895
896 // define the variables
897 Int_t nbins = fhRECpt->GetNbinsX();
898 Double_t binwidth = fhRECpt->GetBinWidth(1);
bb427707 899 Double_t *limits = new Double_t[nbins+1];
8998180c 900 Double_t *binwidths = new Double_t[nbins];
bb427707 901 Double_t xlow=0.;
86bdcd8c 902 for (Int_t i=1; i<=nbins; i++) {
bb427707 903 binwidth = fhRECpt->GetBinWidth(i);
904 xlow = fhRECpt->GetBinLowEdge(i);
905 limits[i-1] = xlow;
de97013c 906 binwidths[i] = binwidth;
bb427707 907 }
908 limits[nbins] = xlow + binwidth;
909
65e55bbd 910 Double_t correction=1.;
a17b17dd 911 Double_t theoryRatio=1.;
912 Double_t effRatio=1.;
86bdcd8c 913 Double_t correctionExtremeA=1., correctionExtremeB=1.;
914 Double_t theoryRatioExtremeA=1., theoryRatioExtremeB=1.;
915 Double_t correctionConservativeA=1., correctionConservativeB=1.;
916 Double_t theoryRatioConservativeA=1., theoryRatioConservativeB=1.;
8998180c 917 Double_t correctionUnc=0.;
918 Double_t correctionExtremeAUnc=0., correctionExtremeBUnc=0.;
919 Double_t correctionConservativeAUnc=0., correctionConservativeBUnc=0.;
86bdcd8c 920
65e55bbd 921 // declare the output histograms
86bdcd8c 922 if (!fhFc) fhFc = new TH1D("fhFc","fc correction factor",nbins,limits);
923 if (!fhFcMax) fhFcMax = new TH1D("fhFcMax","max fc correction factor",nbins,limits);
924 if (!fhFcMin) fhFcMin = new TH1D("fhFcMin","min fc correction factor",nbins,limits);
65e55bbd 925 // two local control histograms
bb427707 926 TH1D *hTheoryRatio = new TH1D("hTheoryRatio","Theoretical B-->D over c-->D (feed-down/direct) ratio",nbins,limits);
927 TH1D *hEffRatio = new TH1D("hEffRatio","Efficiency B-->D over c-->D (feed-down/direct) ratio",nbins,limits);
65e55bbd 928 // and the output TGraphAsymmErrors
86bdcd8c 929 if (fAsymUncertainties & !fgFcExtreme) {
930 fgFcExtreme = new TGraphAsymmErrors(nbins+1);
931 fgFcExtreme->SetNameTitle("fgFcExtreme","fgFcExtreme");
932 }
933 if (fAsymUncertainties & !fgFcConservative) {
934 fgFcConservative = new TGraphAsymmErrors(nbins+1);
935 fgFcConservative->SetNameTitle("fgFcConservative","fgFcConservative");
936 }
937
65e55bbd 938
939 //
940 // Compute fc
941 //
942 for (Int_t ibin=0; ibin<=nbins; ibin++) {
943
944 // theory_ratio = (N_b/N_c)
86bdcd8c 945 theoryRatio = (fhDirectMCpt->GetBinContent(ibin) && fhDirectMCpt->GetBinContent(ibin)!=0.) ?
946 fhFeedDownMCpt->GetBinContent(ibin) / fhDirectMCpt->GetBinContent(ibin) : 1.0 ;
947 //
948 // Calculate the uncertainty [ considering only the theoretical uncertainties on Nb & Nc for now !!! ]
949 //
950 // extreme A = direct-max, feed-down-min
951 theoryRatioExtremeA = (fhDirectMCptMax->GetBinContent(ibin) && fhDirectMCptMax->GetBinContent(ibin)!=0.) ?
952 fhFeedDownMCptMin->GetBinContent(ibin) / fhDirectMCptMax->GetBinContent(ibin) : 1.0 ;
953 // extreme B = direct-min, feed-down-max
954 theoryRatioExtremeB = (fhDirectMCptMin->GetBinContent(ibin) && fhDirectMCptMin->GetBinContent(ibin)!=0.) ?
955 fhFeedDownMCptMax->GetBinContent(ibin) / fhDirectMCptMin->GetBinContent(ibin) : 1.0 ;
956 // conservative A = direct-max, feed-down-max
957 theoryRatioConservativeA = (fhDirectMCptMax->GetBinContent(ibin) && fhDirectMCptMax->GetBinContent(ibin)!=0.) ?
958 fhFeedDownMCptMax->GetBinContent(ibin) / fhDirectMCptMax->GetBinContent(ibin) : 1.0 ;
959 // conservative B = direct-min, feed-down-min
960 theoryRatioConservativeB = (fhDirectMCptMin->GetBinContent(ibin) && fhDirectMCptMin->GetBinContent(ibin)!=0.) ?
961 fhFeedDownMCptMin->GetBinContent(ibin) / fhDirectMCptMin->GetBinContent(ibin) : 1.0 ;
bb427707 962
65e55bbd 963 // eff_ratio = (eff_b/eff_c)
86bdcd8c 964 effRatio = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ?
965 fhFeedDownEffpt->GetBinContent(ibin) / fhDirectEffpt->GetBinContent(ibin) : 1.0 ;
bb427707 966
65e55bbd 967 // fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
86bdcd8c 968 correction = ( 1. / ( 1 + ( effRatio * theoryRatio ) ) );
969 correctionExtremeA = ( 1. / ( 1 + ( effRatio * theoryRatioExtremeA ) ) );
970 correctionExtremeB = ( 1. / ( 1 + ( effRatio * theoryRatioExtremeB ) ) );
971 correctionConservativeA = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeA ) ) );
972 correctionConservativeB = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeB ) ) );
65e55bbd 973
8998180c 974 // fc uncertainty from (eff_b/eff_c) = fc^2 * (N_b/N_c) * delta(eff_b/eff_c)
975 // delta(eff_b/eff_c) is a percentage = fGlobalEfficiencyUncertainties[1]*effRatio
976 correctionUnc = correction*correction * theoryRatio * fGlobalEfficiencyUncertainties[1]*effRatio;
977 correctionExtremeAUnc = correctionExtremeA*correctionExtremeA * theoryRatioExtremeA * fGlobalEfficiencyUncertainties[1]*effRatio;
978 correctionExtremeBUnc = correctionExtremeB*correctionExtremeB * theoryRatioExtremeB * fGlobalEfficiencyUncertainties[1]*effRatio;
979 correctionConservativeAUnc = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA * fGlobalEfficiencyUncertainties[1]*effRatio;
980 correctionConservativeBUnc = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB * fGlobalEfficiencyUncertainties[1]*effRatio;
981
982
65e55bbd 983 // Fill in the histograms
a17b17dd 984 hTheoryRatio->SetBinContent(ibin,theoryRatio);
985 hEffRatio->SetBinContent(ibin,effRatio);
86bdcd8c 986 fhFc->SetBinContent(ibin,correction);
65e55bbd 987 if (fAsymUncertainties) {
988 Double_t x = fhDirectMCpt->GetBinCenter(ibin);
8998180c 989 Double_t val[4] = { correctionExtremeA + correctionExtremeAUnc, correctionExtremeA - correctionExtremeAUnc,
990 correctionExtremeB + correctionExtremeBUnc, correctionExtremeB - correctionExtremeBUnc };
991 Double_t uncExtremeMin = correction - TMath::MinElement(4,val);
992 Double_t uncExtremeMax = TMath::MaxElement(4,val) - correction;
86bdcd8c 993 fgFcExtreme->SetPoint(ibin,x,correction); // i,x,y
994 fgFcExtreme->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),uncExtremeMin,uncExtremeMax); // i,xl,xh,yl,yh
995 fhFcMax->SetBinContent(ibin,correction+uncExtremeMax);
996 fhFcMin->SetBinContent(ibin,correction-uncExtremeMin);
8998180c 997 Double_t consval[4] = { correctionConservativeA - correctionConservativeAUnc, correctionConservativeA + correctionConservativeAUnc,
998 correctionConservativeB - correctionConservativeBUnc, correctionConservativeB + correctionConservativeBUnc};
999// Double_t uncConservativeMin = correction - (correctionConservativeA - correctionConservativeAUnc);
1000// Double_t uncConservativeMax = (correctionConservativeB + correctionConservativeBUnc) - correction;
1001 Double_t uncConservativeMin = correction - TMath::MinElement(4,consval);
1002 Double_t uncConservativeMax = TMath::MaxElement(4,consval) - correction;
86bdcd8c 1003 fgFcConservative->SetPoint(ibin,x,correction); // i,x,y
1004 fgFcConservative->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),uncConservativeMin,uncConservativeMax); // i,xl,xh,yl,yh
65e55bbd 1005 }
1006
1007 }
1008
65e55bbd 1009}
1010
1011//_________________________________________________________________________________________________________
a17b17dd 1012void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumFc(){
65e55bbd 1013 //
1014 // Compute the feed-down corrected spectrum if feed-down correction is done via fc factor (bin by bin)
bb427707 1015 // physics = reco * fc / bin-width
65e55bbd 1016 //
86bdcd8c 1017 // uncertainty: (stat) delta_physics = physics * sqrt ( (delta_reco/reco)^2 )
1018 // (syst but feed-down) delta_physics = physics * sqrt ( (delta_reco_syst/reco)^2 )
1019 // (feed-down syst) delta_physics = physics * sqrt ( (delta_fc/fc)^2 )
65e55bbd 1020 //
1021 // ( Calculation done bin by bin )
1022
1023 if (!fhFc || !fhRECpt) {
1024 AliError(" Reconstructed or fc distributions are not defined");
1025 return;
1026 }
1027
1028 Int_t nbins = fhRECpt->GetNbinsX();
86bdcd8c 1029 Double_t value = 0., errvalue = 0., errvalueMax= 0., errvalueMin= 0.;
1030 Double_t valueExtremeMax= 0., valueExtremeMin= 0.;
1031 Double_t valueConservativeMax= 0., valueConservativeMin= 0.;
65e55bbd 1032 Double_t binwidth = fhRECpt->GetBinWidth(1);
bb427707 1033 Double_t *limits = new Double_t[nbins+1];
8998180c 1034 Double_t *binwidths = new Double_t[nbins];
bb427707 1035 Double_t xlow=0.;
86bdcd8c 1036 for (Int_t i=1; i<=nbins; i++) {
bb427707 1037 binwidth = fhRECpt->GetBinWidth(i);
1038 xlow = fhRECpt->GetBinLowEdge(i);
1039 limits[i-1] = xlow;
de97013c 1040 binwidths[i] = binwidth;
bb427707 1041 }
86bdcd8c 1042 limits[nbins] = xlow + binwidth;
65e55bbd 1043
1044 // declare the output histograms
86bdcd8c 1045 if (!fhYieldCorr) fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by fc)",nbins,limits);
1046 if (!fhYieldCorrMax) fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by fc)",nbins,limits);
1047 if (!fhYieldCorrMin) fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by fc)",nbins,limits);
65e55bbd 1048 // and the output TGraphAsymmErrors
86bdcd8c 1049 if (fAsymUncertainties & !fgYieldCorr) fgYieldCorr = new TGraphAsymmErrors(nbins+1);
1050 if (fAsymUncertainties & !fgYieldCorrExtreme) fgYieldCorrExtreme = new TGraphAsymmErrors(nbins+1);
1051 if (fAsymUncertainties & !fgYieldCorrConservative) fgYieldCorrConservative = new TGraphAsymmErrors(nbins+1);
65e55bbd 1052
1053 //
1054 // Do the calculation
1055 //
1056 for (Int_t ibin=0; ibin<=nbins; ibin++) {
1057
1058 // calculate the value
86bdcd8c 1059 // physics = reco * fc / bin-width
65e55bbd 1060 value = fhRECpt->GetBinContent(ibin) * fhFc->GetBinContent(ibin) ;
bb427707 1061 value /= fhRECpt->GetBinWidth(ibin) ;
86bdcd8c 1062
1063 // Statistical uncertainty
1064 // (stat) delta_physics = physics * sqrt ( (delta_reco/reco)^2 )
1065 errvalue = value * (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) ;
1066
1067 // Calculate the systematic uncertainties
1068 // (syst but feed-down) delta_physics = physics * sqrt ( (delta_reco_syst/reco)^2 )
1069 // (feed-down syst) delta_physics = physics * sqrt ( (delta_fc/fc)^2 )
1070 //
65e55bbd 1071 // Protect against null denominator. If so, define uncertainty as null
1072 if (fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) {
1073
1074 if (fAsymUncertainties) {
1075
86bdcd8c 1076 // Systematics but feed-down
1077 if (fgRECSystematics) {
1078 errvalueMax = value * ( fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinContent(ibin) );
1079 errvalueMin = value * ( fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinContent(ibin) );
65e55bbd 1080 }
d38da1b0 1081 else { errvalueMax = 0.; errvalueMin = 0.; }
86bdcd8c 1082
1083 // Extreme feed-down systematics
1084 valueExtremeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcExtreme->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1085 valueExtremeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcExtreme->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1086
1087 // Conservative feed-down systematics
1088 valueConservativeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcConservative->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1089 valueConservativeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcConservative->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
65e55bbd 1090
1091 }
65e55bbd 1092
1093 }
86bdcd8c 1094 else { errvalueMax = 0.; errvalueMin = 0.; }
65e55bbd 1095
1096 // fill in the histograms
86bdcd8c 1097 fhYieldCorr->SetBinContent(ibin,value);
1098 fhYieldCorr->SetBinError(ibin,errvalue);
65e55bbd 1099 if (fAsymUncertainties) {
86bdcd8c 1100 Double_t center = fhYieldCorr->GetBinCenter(ibin);
65e55bbd 1101 fgYieldCorr->SetPoint(ibin,center,value); // i,x,y
86bdcd8c 1102 fgYieldCorr->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
1103 fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax);
1104 fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
1105 fgYieldCorrExtreme->SetPoint(ibin,center,value);
1106 fgYieldCorrExtreme->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),value-valueExtremeMin,valueExtremeMax-value);
1107 fgYieldCorrConservative->SetPoint(ibin,center,value);
1108 fgYieldCorrConservative->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),value-valueConservativeMin,valueConservativeMax-value);
65e55bbd 1109 }
1110
1111 }
1112
65e55bbd 1113}
1114
1115//_________________________________________________________________________________________________________
a17b17dd 1116void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumNb(Double_t deltaY, Double_t branchingRatioBintoFinalDecay) {
65e55bbd 1117 //
1118 // Compute the feed-down corrected spectrum if feed-down correction is done via Nb (bin by bin)
bb427707 1119 // physics = [ reco - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) ] / bin-width
65e55bbd 1120 //
86bdcd8c 1121 // uncertainty: (stat) delta_physics = sqrt ( (delta_reco)^2 ) / bin-width
1122 // (syst but feed-down) delta_physics = sqrt ( (delta_reco_syst)^2 ) / bin-width
1123 // (feed-down syst) delta_physics = sqrt ( (k*delta_lumi/lumi)^2 + (k*delta_eff_trig/eff_trig)^2
8998180c 1124 // + (k*delta_Nb/Nb)^2 + (k*delta_eff/eff)^2 + (k*global_eff_ratio)^2 ) / bin-width
65e55bbd 1125 // where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th
1126 //
1127
1128 Int_t nbins = fhRECpt->GetNbinsX();
1129 Double_t binwidth = fhRECpt->GetBinWidth(1);
86bdcd8c 1130 Double_t value = 0., errvalue = 0., errvalueMax = 0., errvalueMin = 0., kfactor = 0.;
1131 Double_t errvalueExtremeMax = 0., errvalueExtremeMin = 0.;
bb427707 1132 Double_t *limits = new Double_t[nbins+1];
de97013c 1133 Double_t *binwidths = new Double_t[nbins+1];
bb427707 1134 Double_t xlow=0.;
86bdcd8c 1135 for (Int_t i=1; i<=nbins; i++) {
bb427707 1136 binwidth = fhRECpt->GetBinWidth(i);
1137 xlow = fhRECpt->GetBinLowEdge(i);
1138 limits[i-1] = xlow;
de97013c 1139 binwidths[i] = binwidth;
bb427707 1140 }
86bdcd8c 1141 limits[nbins] = xlow + binwidth;
65e55bbd 1142
1143 // declare the output histograms
86bdcd8c 1144 if (!fhYieldCorr) fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by Nb)",nbins,limits);
1145 if (!fhYieldCorrMax) fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by Nb)",nbins,limits);
1146 if (!fhYieldCorrMin) fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by Nb)",nbins,limits);
65e55bbd 1147 // and the output TGraphAsymmErrors
86bdcd8c 1148 if (fAsymUncertainties & !fgYieldCorr) fgYieldCorr = new TGraphAsymmErrors(nbins+1);
1149 if (fAsymUncertainties & !fgYieldCorrExtreme) fgYieldCorrExtreme = new TGraphAsymmErrors(nbins+1);
1150 if (fAsymUncertainties & !fgYieldCorrConservative) {
1151 fgYieldCorrConservative = new TGraphAsymmErrors(nbins+1);
1152 AliInfo(" Beware the conservative & extreme uncertainties are equal by definition !");
1153 }
65e55bbd 1154
1155 //
1156 // Do the calculation
1157 //
1158 for (Int_t ibin=0; ibin<=nbins; ibin++) {
1159
86bdcd8c 1160 // Calculate the value
1161 // physics = [ reco - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) ] / bin-width
a17b17dd 1162 value = fhRECpt->GetBinContent(ibin) - (deltaY*branchingRatioBintoFinalDecay*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) );
bb427707 1163 value /= fhRECpt->GetBinWidth(ibin);
65e55bbd 1164
86bdcd8c 1165 // Statistical uncertainty: delta_physics = sqrt ( (delta_reco)^2 ) / bin-width
1166 errvalue = fhRECpt->GetBinError(ibin);
1167 errvalue /= fhRECpt->GetBinWidth(ibin);
65e55bbd 1168
86bdcd8c 1169 // Systematic uncertainties
1170 // (syst but feed-down) delta_physics = sqrt ( (delta_reco_syst)^2 ) / bin-width
1171 // (feed-down syst) delta_physics = sqrt ( (k*delta_lumi/lumi)^2 + (k*delta_eff_trig/eff_trig)^2
8998180c 1172 // + (k*delta_Nb/Nb)^2 + (k*delta_eff/eff)^2 + (k*global_eff_ratio)^2 ) / bin-width
86bdcd8c 1173 // where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th
1174 kfactor = deltaY*branchingRatioBintoFinalDecay*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) ;
1175 //
65e55bbd 1176 if (fAsymUncertainties) {
8998180c 1177 Double_t Nb = fhFeedDownMCpt->GetBinContent(ibin);
1178 Double_t NbDmax = fhFeedDownMCptMax->GetBinContent(ibin) - fhFeedDownMCpt->GetBinContent(ibin);
1179 Double_t NbDmin = fhFeedDownMCpt->GetBinContent(ibin) - fhFeedDownMCptMin->GetBinContent(ibin);
86bdcd8c 1180
1181 // Systematics but feed-down
1182 if (fgRECSystematics){
1183 errvalueMax = fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinWidth(ibin) ;
1184 errvalueMin = fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinWidth(ibin);
1185 }
d38da1b0 1186 else { errvalueMax = 0.; errvalueMin = 0.; }
86bdcd8c 1187
1188 // Feed-down systematics
1189 // min value with the maximum Nb
1190 errvalueExtremeMin = TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
1191 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
8998180c 1192 ( (kfactor*NbDmax/Nb)*(kfactor*NbDmax/Nb) ) +
1193 ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
1194 ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1195 ) / fhRECpt->GetBinWidth(ibin);
86bdcd8c 1196 // max value with the minimum Nb
1197 errvalueExtremeMax = TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
1198 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
8998180c 1199 ( (kfactor*NbDmin/Nb)*(kfactor*NbDmin/Nb) ) +
1200 ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
1201 ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1202 ) / fhRECpt->GetBinWidth(ibin);
65e55bbd 1203 }
1204 else{ // Don't consider Nb uncertainty in this case [ to be tested!!! ]
86bdcd8c 1205 errvalueExtremeMax = TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
1206 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
8998180c 1207 ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
1208 ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1209 ) / fhRECpt->GetBinWidth(ibin);
86bdcd8c 1210 errvalueExtremeMin = errvalueExtremeMax ;
65e55bbd 1211 }
1212
1213 // fill in histograms
86bdcd8c 1214 fhYieldCorr->SetBinContent(ibin,value);
1215 fhYieldCorr->SetBinError(ibin,errvalue);
65e55bbd 1216 if (fAsymUncertainties) {
86bdcd8c 1217 Double_t x = fhYieldCorr->GetBinCenter(ibin);
65e55bbd 1218 fgYieldCorr->SetPoint(ibin,x,value); // i,x,y
86bdcd8c 1219 fgYieldCorr->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
1220 fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax);
1221 fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
1222 fgYieldCorrExtreme->SetPoint(ibin,x,value); // i,x,y
1223 fgYieldCorrExtreme->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
1224 fgYieldCorrConservative->SetPoint(ibin,x,value); // i,x,y
1225 fgYieldCorrConservative->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
65e55bbd 1226 }
1227
1228 }
1229
65e55bbd 1230}
1231
1232
8998180c 1233//_________________________________________________________________________________________________________
1234void AliHFPtSpectrum::ComputeSystUncertainties(Int_t decay, Bool_t combineFeedDown) {
1235
1236 // Call the systematics uncertainty class for a given decay
1237 AliHFSystErr systematics(decay);
1238
1239 // Estimate the feed-down uncertainty in percentage
1240 Int_t nentries = fgSigmaCorrConservative->GetN();
1241 TGraphAsymmErrors *grErrFeeddown = new TGraphAsymmErrors(nentries);
1242 Double_t x=0., y=0., errx=0., erryl=0., erryh=0;
1243 for(Int_t i=0; i<nentries; i++) {
1244 fgSigmaCorrConservative->GetPoint(i,x,y);
1245 errx = fgSigmaCorrConservative->GetErrorXlow(i) ;
1246 erryl = fgSigmaCorrConservative->GetErrorYlow(i) / y ;
1247 erryh = fgSigmaCorrConservative->GetErrorYhigh(i) / y ;
1248 grErrFeeddown->SetPoint(i,x,0.);
1249 grErrFeeddown->SetPointError(i,errx,errx,erryl,erryh); //i, xl, xh, yl, yh
1250 }
1251
1252 // Draw all the systematics independently
1253 systematics.DrawErrors(grErrFeeddown);
1254
1255 // Set the sigma systematic uncertainties
1256 // possibly combine with the feed-down uncertainties
1257 Double_t errylcomb=0., erryhcomb=0;
1258 for(Int_t i=1; i<nentries; i++) {
1259 fgSigmaCorr->GetPoint(i,x,y);
1260 errx = grErrFeeddown->GetErrorXlow(i) ;
1261 erryl = grErrFeeddown->GetErrorYlow(i);
1262 erryh = grErrFeeddown->GetErrorYhigh(i);
1263 if (combineFeedDown) {
1264 errylcomb = systematics.GetTotalSystErr(x,erryl) * y ;
1265 erryhcomb = systematics.GetTotalSystErr(x,erryh) * y ;
1266 } else {
1267 errylcomb = systematics.GetTotalSystErr(x) * y ;
1268 erryhcomb = systematics.GetTotalSystErr(x) * y ;
1269 }
1270 fgSigmaCorr->SetPointError(i,errx,errx,errylcomb,erryhcomb);
1271 }
1272
1273}
1274
1275
1276//_________________________________________________________________________________________________________
1277void AliHFPtSpectrum::DrawSpectrum(TGraphAsymmErrors *gPrediction) {
1278
1279 TCanvas *csigma = new TCanvas("csigma","Draw the corrected cross-section & the prediction");
1280 csigma->SetFillColor(0);
1281 gPrediction->GetXaxis()->SetTitleSize(0.05);
1282 gPrediction->GetXaxis()->SetTitleOffset(0.95);
1283 gPrediction->GetYaxis()->SetTitleSize(0.05);
1284 gPrediction->GetYaxis()->SetTitleOffset(0.95);
1285 gPrediction->GetXaxis()->SetTitle("p_{T} [GeV]");
1286 gPrediction->GetYaxis()->SetTitle("BR #times #frac{d#sigma}{dp_{T}} |_{|y|<0.5} [pb/GeV]");
1287 gPrediction->SetLineColor(kGreen+2);
1288 gPrediction->SetLineWidth(3);
1289 gPrediction->SetFillColor(kGreen+1);
1290 gPrediction->Draw("3CA");
1291 fgSigmaCorr->SetLineColor(kRed);
1292 fgSigmaCorr->SetLineWidth(1);
1293 fgSigmaCorr->SetFillColor(kRed);
1294 fgSigmaCorr->SetFillStyle(0);
1295 fgSigmaCorr->Draw("2");
1296 fhSigmaCorr->SetMarkerColor(kRed);
1297 fhSigmaCorr->Draw("esame");
1298 csigma->SetLogy();
1299 TLegend * leg = new TLegend(0.7,0.75,0.87,0.5);
1300 leg->SetBorderSize(0);
1301 leg->SetLineColor(0);
1302 leg->SetFillColor(0);
1303 leg->SetTextFont(42);
1304 leg->AddEntry(gPrediction,"FONLL ","fl");
1305 leg->AddEntry(fhSigmaCorr,"data stat. unc.","pl");
1306 leg->AddEntry(fgSigmaCorr,"data syst. unc.","f");
1307 leg->Draw();
1308 csigma->Draw();
1309
1310}
1311
65e55bbd 1312//_________________________________________________________________________________________________________
86bdcd8c 1313TH1D * AliHFPtSpectrum::ReweightHisto(TH1D *hToReweight, TH1D *hReference){
65e55bbd 1314 //
1315 // Function to reweight histograms for testing purposes:
1316 // This function takes the histo hToReweight and reweights
1317 // it (its pt shape) with respect to hReference
1318 //
1319
1320 // check histograms consistency
1321 Bool_t areconsistent=kTRUE;
1322 areconsistent &= CheckHistosConsistency(hToReweight,hReference);
1323 if (!areconsistent) {
1324 AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
1325 return NULL;
1326 }
1327
1328 // define a new empty histogram
86bdcd8c 1329 TH1D *hReweighted = (TH1D*)hToReweight->Clone("hReweighted");
65e55bbd 1330 hReweighted->Reset();
1331 Double_t weight=1.0;
1332 Double_t yvalue=1.0;
a17b17dd 1333 Double_t integralRef = hReference->Integral();
1334 Double_t integralH = hToReweight->Integral();
65e55bbd 1335
1336 // now reweight the spectra
1337 //
1338 // the weight is the relative probability of the given pt bin in the reference histo
1339 // divided by its relative probability (to normalize it) on the histo to re-weight
1340 for (Int_t i=0; i<=hToReweight->GetNbinsX(); i++) {
a17b17dd 1341 weight = (hReference->GetBinContent(i)/integralRef) / (hToReweight->GetBinContent(i)/integralH) ;
65e55bbd 1342 yvalue = hToReweight->GetBinContent(i);
1343 hReweighted->SetBinContent(i,yvalue*weight);
1344 }
1345
86bdcd8c 1346 return (TH1D*)hReweighted;
65e55bbd 1347}
1348
1349//_________________________________________________________________________________________________________
86bdcd8c 1350TH1D * AliHFPtSpectrum::ReweightRecHisto(TH1D *hRecToReweight, TH1D *hMCToReweight, TH1D *hMCReference){
65e55bbd 1351 //
1352 // Function to reweight histograms for testing purposes:
1353 // This function takes the histo hToReweight and reweights
1354 // it (its pt shape) with respect to hReference /hMCToReweight
1355 //
1356
1357 // check histograms consistency
1358 Bool_t areconsistent=kTRUE;
1359 areconsistent &= CheckHistosConsistency(hMCToReweight,hMCReference);
1360 areconsistent &= CheckHistosConsistency(hRecToReweight,hMCReference);
1361 if (!areconsistent) {
1362 AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
1363 return NULL;
1364 }
1365
1366 // define a new empty histogram
86bdcd8c 1367 TH1D *hReweighted = (TH1D*)hMCToReweight->Clone("hReweighted");
65e55bbd 1368 hReweighted->Reset();
86bdcd8c 1369 TH1D *hRecReweighted = (TH1D*)hRecToReweight->Clone("hRecReweighted");
65e55bbd 1370 hRecReweighted->Reset();
1371 Double_t weight=1.0;
1372 Double_t yvalue=1.0, yrecvalue=1.0;
a17b17dd 1373 Double_t integralRef = hMCReference->Integral();
1374 Double_t integralH = hMCToReweight->Integral();
65e55bbd 1375
1376 // now reweight the spectra
1377 //
1378 // the weight is the relative probability of the given pt bin
1379 // that should be applied in the MC histo to get the reference histo shape
1380 // Probabilities are properly normalized.
1381 for (Int_t i=0; i<=hMCToReweight->GetNbinsX(); i++) {
a17b17dd 1382 weight = (hMCReference->GetBinContent(i)/integralRef) / (hMCToReweight->GetBinContent(i)/integralH) ;
65e55bbd 1383 yvalue = hMCToReweight->GetBinContent(i);
1384 hReweighted->SetBinContent(i,yvalue*weight);
1385 yrecvalue = hRecToReweight->GetBinContent(i);
1386 hRecReweighted->SetBinContent(i,yrecvalue*weight);
1387 }
1388
86bdcd8c 1389 return (TH1D*)hRecReweighted;
65e55bbd 1390}
1391