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