Coverity fixes
[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 ?! ");
86df816f 564 delete [] limits;
565 delete [] binwidths;
65e55bbd 566 return ;
567 }
568
86bdcd8c 569 Double_t value=0, errValue=0, errvalueMax=0., errvalueMin=0.;
570 Double_t errvalueExtremeMax=0., errvalueExtremeMin=0.;
571 Double_t errvalueConservativeMax=0., errvalueConservativeMin=0.;
b890d92c 572 for(Int_t ibin=1; ibin<=nbins; ibin++){
65e55bbd 573
86bdcd8c 574 // Sigma calculation
65e55bbd 575 // Sigma = ( 1. / (lumi * delta_y * BR_c * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
576 value = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ?
a17b17dd 577 ( fhYieldCorr->GetBinContent(ibin) / ( deltaY * branchingRatioC * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
65e55bbd 578 : 0. ;
86bdcd8c 579
580 // Sigma statistical uncertainty:
581 // delta_sigma = sigma * sqrt ( (delta_spectra/spectra)^2 )
582 errValue = (value!=0.) ? value * (fhYieldCorr->GetBinError(ibin)/fhYieldCorr->GetBinContent(ibin)) : 0. ;
583
b890d92c 584
86bdcd8c 585 //
586 // Sigma systematic uncertainties
587 //
65e55bbd 588 if (fAsymUncertainties) {
86bdcd8c 589
590 // (syst but feed-down) delta_sigma = sigma * sqrt ( (delta_spectra_syst/spectra)^2 +
8998180c 591 // (delta_lumi/lumi)^2 + (delta_eff_trig/eff_trig)^2 + (delta_eff/eff)^2 + (global_eff)^2 )
86bdcd8c 592 errvalueMax = value * TMath::Sqrt( (fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin))*(fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin)) +
593 (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
594 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) +
4c7dab0f 595 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
8998180c 596 fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] );
86bdcd8c 597 errvalueMin = value * TMath::Sqrt( (fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin))*(fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin)) +
598 (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
599 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) +
4c7dab0f 600 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
8998180c 601 fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] );
86bdcd8c 602
603 // Uncertainties from feed-down
604 // (feed-down syst) delta_sigma = sigma * sqrt ( (delta_spectra_fd/spectra_fd)^2 )
605 // extreme case
606 errvalueExtremeMax = value * (fgYieldCorrExtreme->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin));
607 errvalueExtremeMin = value * (fgYieldCorrExtreme->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin));
608 //
609 // conservative case
610 errvalueConservativeMax = value * (fgYieldCorrConservative->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin));
611 errvalueConservativeMin = value * (fgYieldCorrConservative->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin));
612
65e55bbd 613 }
614 else {
615 // protect against null denominator
86bdcd8c 616 errvalueMax = (value!=0.) ?
617 value * TMath::Sqrt( (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
618 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) +
4c7dab0f 619 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin)) +
8998180c 620 fGlobalEfficiencyUncertainties[0]*fGlobalEfficiencyUncertainties[0] )
65e55bbd 621 : 0. ;
86bdcd8c 622 errvalueMin = errvalueMax;
65e55bbd 623 }
624
625 // Fill the histograms
86bdcd8c 626 fhSigmaCorr->SetBinContent(ibin,value);
627 fhSigmaCorr->SetBinError(ibin,errValue);
65e55bbd 628 // Fill the TGraphAsymmErrors
629 if (fAsymUncertainties) {
630 Double_t x = fhYieldCorr->GetBinCenter(ibin);
631 fgSigmaCorr->SetPoint(ibin,x,value); // i,x,y
b890d92c 632 fgSigmaCorr->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
86bdcd8c 633 fhSigmaCorrMax->SetBinContent(ibin,value+errvalueMax);
634 fhSigmaCorrMin->SetBinContent(ibin,value-errvalueMin);
635 fgSigmaCorrExtreme->SetPoint(ibin,x,value); // i,x,y
b890d92c 636 fgSigmaCorrExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
86bdcd8c 637 fgSigmaCorrConservative->SetPoint(ibin,x,value); // i,x,y
b890d92c 638 fgSigmaCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueConservativeMin,errvalueConservativeMax); // i,xl,xh,yl,yh
86bdcd8c 639
65e55bbd 640 }
641
642 }
86df816f 643 delete [] binwidths;
644 delete [] limits;
65e55bbd 645
65e55bbd 646}
647
648//_________________________________________________________________________________________________________
86bdcd8c 649TH1D * AliHFPtSpectrum::EstimateEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco, const char *name) {
bb427707 650 //
5f3c1b97 651 // Function that computes the acceptance and efficiency correction
bb427707 652 // based on the simulated and reconstructed spectra
653 // and using the reconstructed spectra bin width
654 //
655 // eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
656 //
657
658 if(!fhRECpt){
659 AliInfo("Hey, the reconstructed histogram was not set yet !");
660 return NULL;
661 }
662
663 Int_t nbins = fhRECpt->GetNbinsX();
664 Double_t *limits = new Double_t[nbins+1];
665 Double_t xlow=0.,binwidth=0.;
86bdcd8c 666 for (Int_t i=1; i<=nbins; i++) {
bb427707 667 binwidth = fhRECpt->GetBinWidth(i);
668 xlow = fhRECpt->GetBinLowEdge(i);
669 limits[i-1] = xlow;
670 }
671 limits[nbins] = xlow + binwidth;
672
86bdcd8c 673 TH1D * hEfficiency = new TH1D(name," acceptance #times efficiency",nbins,limits);
bb427707 674
675 Double_t sumSimu[nbins], sumReco[nbins];
b890d92c 676 for (Int_t ibin=0; ibin<nbins; ibin++){
677 sumSimu[ibin]=0.; sumReco[ibin]=0.;
678 }
bb427707 679 for (Int_t ibin=0; ibin<=hSimu->GetNbinsX(); ibin++){
680
681 for (Int_t ibinrec=0; ibinrec<nbins; ibinrec++){
682 if ( hSimu->GetBinCenter(ibin)>limits[ibinrec] &&
683 hSimu->GetBinCenter(ibin)<limits[ibinrec+1] ) {
684 sumSimu[ibinrec]+=hSimu->GetBinContent(ibin);
685 }
686 if ( hReco->GetBinCenter(ibin)>limits[ibinrec] &&
687 hReco->GetBinCenter(ibin)<limits[ibinrec+1] ) {
688 sumReco[ibinrec]+=hReco->GetBinContent(ibin);
689 }
690 }
691
692 }
693
86bdcd8c 694
bb427707 695 // the efficiency is computed as reco/sim (in each bin)
696 // its uncertainty is err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
697 Double_t eff=0., erreff=0.;
86bdcd8c 698 for (Int_t ibinrec=0; ibinrec<=nbins; ibinrec++) {
699 if (sumSimu[ibinrec]!= 0. && sumReco[ibinrec]!=0.) {
700 eff = sumReco[ibinrec] / sumSimu[ibinrec] ;
b890d92c 701 // protection in case eff > 1.0
702 // test calculation (make the argument of the sqrt positive)
703 erreff = TMath::Sqrt( eff * TMath::Abs(1.0 - eff) ) / TMath::Sqrt( sumSimu[ibinrec] );
86bdcd8c 704 }
705 else { eff=0.0; erreff=0.; }
5f3c1b97 706 hEfficiency->SetBinContent(ibinrec+1,eff);
707 hEfficiency->SetBinError(ibinrec+1,erreff);
bb427707 708 }
709
86bdcd8c 710 return (TH1D*)hEfficiency;
bb427707 711}
712
713//_________________________________________________________________________________________________________
86bdcd8c 714void AliHFPtSpectrum::EstimateAndSetDirectEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco) {
bb427707 715 //
5f3c1b97 716 // Function that computes the Direct acceptance and efficiency correction
bb427707 717 // based on the simulated and reconstructed spectra
718 // and using the reconstructed spectra bin width
719 //
720 // eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
721 //
5f3c1b97 722
86bdcd8c 723 if(!fhRECpt || !hSimu || !hReco){
724 AliError("Hey, the reconstructed histogram was not set yet !");
725 return;
bb427707 726 }
727
86bdcd8c 728 fhDirectEffpt = EstimateEfficiencyRecoBin(hSimu,hReco,"fhDirectEffpt");
5f3c1b97 729 fhDirectEffpt->SetNameTitle("fhDirectEffpt"," direct acceptance #times efficiency");
bb427707 730
5f3c1b97 731}
bb427707 732
5f3c1b97 733//_________________________________________________________________________________________________________
86bdcd8c 734void AliHFPtSpectrum::EstimateAndSetFeedDownEfficiencyRecoBin(TH1D *hSimu, TH1D *hReco) {
5f3c1b97 735 //
736 // Function that computes the Feed-Down acceptance and efficiency correction
737 // based on the simulated and reconstructed spectra
738 // and using the reconstructed spectra bin width
739 //
740 // eff = reco/sim ; err_eff = sqrt( eff*(1-eff) )/ sqrt( sim )
741 //
bb427707 742
86bdcd8c 743 if(!fhRECpt || !hSimu || !hReco){
744 AliError("Hey, the reconstructed histogram was not set yet !");
745 return;
bb427707 746 }
747
86bdcd8c 748 fhFeedDownEffpt = EstimateEfficiencyRecoBin(hSimu,hReco,"fhFeedDownEffpt");
5f3c1b97 749 fhFeedDownEffpt->SetNameTitle("fhFeedDownEffpt"," feed-down acceptance #times efficiency");
750
bb427707 751}
752
753//_________________________________________________________________________________________________________
65e55bbd 754Bool_t AliHFPtSpectrum::Initialize(){
755 //
756 // Initialization of the variables (histograms)
757 //
758
759 if (fFeedDownOption==0) {
760 AliInfo("Getting ready for the corrections without feed-down consideration");
761 } else if (fFeedDownOption==1) {
762 AliInfo("Getting ready for the fc feed-down correction calculation");
763 } else if (fFeedDownOption==2) {
764 AliInfo("Getting ready for the Nb feed-down correction calculation");
765 } else { AliError("The calculation option must be <=2"); return kFALSE; }
766
767 // Start checking the input histograms consistency
768 Bool_t areconsistent=kTRUE;
769
770 // General checks
771 if (!fhDirectEffpt || !fhRECpt) {
772 AliError(" Reconstructed spectra and/or the Nc efficiency distributions are not defined");
773 return kFALSE;
774 }
775 areconsistent &= CheckHistosConsistency(fhRECpt,fhDirectEffpt);
776 if (!areconsistent) {
777 AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
778 return kFALSE;
779 }
780 if (fFeedDownOption==0) return kTRUE;
781
782 //
783 // Common checks for options 1 (fc) & 2(Nb)
784 if (!fhFeedDownMCpt || !fhFeedDownEffpt) {
785 AliError(" Theoretical Nb and/or the Nb efficiency distributions are not defined");
786 return kFALSE;
787 }
788 areconsistent &= CheckHistosConsistency(fhRECpt,fhFeedDownMCpt);
789 areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownEffpt);
790 if (fAsymUncertainties) {
a17b17dd 791 if (!fhFeedDownMCptMax || !fhFeedDownMCptMin) {
65e55bbd 792 AliError(" Max/Min theoretical Nb distributions are not defined");
793 return kFALSE;
794 }
a17b17dd 795 areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownMCptMax);
65e55bbd 796 }
797 if (!areconsistent) {
798 AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
799 return kFALSE;
800 }
801 if (fFeedDownOption>1) return kTRUE;
802
803 //
804 // Now checks for option 1 (fc correction)
805 if (!fhDirectMCpt) {
806 AliError("Theoretical Nc distributions is not defined");
807 return kFALSE;
808 }
809 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhFeedDownMCpt);
810 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectEffpt);
811 if (fAsymUncertainties) {
a17b17dd 812 if (!fhDirectMCptMax || !fhDirectMCptMin) {
65e55bbd 813 AliError(" Max/Min theoretical Nc distributions are not defined");
814 return kFALSE;
815 }
a17b17dd 816 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectMCptMax);
65e55bbd 817 }
818 if (!areconsistent) {
819 AliInfo("Histograms required for fc correction are not consistent (bin width, bounds)");
820 return kFALSE;
821 }
822
823 return kTRUE;
824}
825
826//_________________________________________________________________________________________________________
86bdcd8c 827Bool_t AliHFPtSpectrum::CheckHistosConsistency(TH1D *h1, TH1D *h2){
65e55bbd 828 //
829 // Check the histograms consistency (bins, limits)
830 //
831
832 if (!h1 || !h2) {
833 AliError("One or both histograms don't exist");
834 return kFALSE;
835 }
836
837 Double_t binwidth1 = h1->GetBinWidth(1);
838 Double_t binwidth2 = h2->GetBinWidth(1);
839 Double_t min1 = h1->GetBinCenter(1) - (binwidth1/2.) ;
840// Double_t max1 = h1->GetBinCenter(nbins1) + (binwidth1/2.) ;
841 Double_t min2 = h2->GetBinCenter(1) - (binwidth2/2.) ;
842// Double_t max2 = h2->GetBinCenter(nbins2) + (binwidth2/2.) ;
843
844 if (binwidth1!=binwidth2) {
845 AliInfo(" histograms with different bin width");
846 return kFALSE;
847 }
848 if (min1!=min2) {
849 AliInfo(" histograms with different minimum");
850 return kFALSE;
851 }
852// if (max1!=max2) {
853// AliInfo(" histograms with different maximum");
854// return kFALSE;
855// }
856
857 return kTRUE;
858}
859
860//_________________________________________________________________________________________________________
a17b17dd 861void AliHFPtSpectrum::CalculateFeedDownCorrectionFc(){
65e55bbd 862 //
863 // Compute fc factor and its uncertainties bin by bin
864 // fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
865 //
86bdcd8c 866 // uncertainties: (conservative) combine the upper/lower N_b & N_c predictions together
867 // (extreme) combine the upper N_b predictions with the lower N_c predictions & viceversa
8998180c 868 // systematic uncertainty on the acceptance x efficiency b/c ratio are included
86bdcd8c 869 //
65e55bbd 870
871 // define the variables
872 Int_t nbins = fhRECpt->GetNbinsX();
873 Double_t binwidth = fhRECpt->GetBinWidth(1);
bb427707 874 Double_t *limits = new Double_t[nbins+1];
8998180c 875 Double_t *binwidths = new Double_t[nbins];
bb427707 876 Double_t xlow=0.;
86bdcd8c 877 for (Int_t i=1; i<=nbins; i++) {
bb427707 878 binwidth = fhRECpt->GetBinWidth(i);
879 xlow = fhRECpt->GetBinLowEdge(i);
880 limits[i-1] = xlow;
b890d92c 881 binwidths[i-1] = binwidth;
bb427707 882 }
883 limits[nbins] = xlow + binwidth;
884
65e55bbd 885 Double_t correction=1.;
a17b17dd 886 Double_t theoryRatio=1.;
887 Double_t effRatio=1.;
86bdcd8c 888 Double_t correctionExtremeA=1., correctionExtremeB=1.;
889 Double_t theoryRatioExtremeA=1., theoryRatioExtremeB=1.;
890 Double_t correctionConservativeA=1., correctionConservativeB=1.;
891 Double_t theoryRatioConservativeA=1., theoryRatioConservativeB=1.;
8998180c 892 Double_t correctionUnc=0.;
893 Double_t correctionExtremeAUnc=0., correctionExtremeBUnc=0.;
894 Double_t correctionConservativeAUnc=0., correctionConservativeBUnc=0.;
86bdcd8c 895
65e55bbd 896 // declare the output histograms
b890d92c 897 fhFc = new TH1D("fhFc","fc correction factor",nbins,limits);
898 fhFcMax = new TH1D("fhFcMax","max fc correction factor",nbins,limits);
899 fhFcMin = new TH1D("fhFcMin","min fc correction factor",nbins,limits);
65e55bbd 900 // two local control histograms
bb427707 901 TH1D *hTheoryRatio = new TH1D("hTheoryRatio","Theoretical B-->D over c-->D (feed-down/direct) ratio",nbins,limits);
902 TH1D *hEffRatio = new TH1D("hEffRatio","Efficiency B-->D over c-->D (feed-down/direct) ratio",nbins,limits);
65e55bbd 903 // and the output TGraphAsymmErrors
b890d92c 904 if (fAsymUncertainties) {
86bdcd8c 905 fgFcExtreme = new TGraphAsymmErrors(nbins+1);
906 fgFcExtreme->SetNameTitle("fgFcExtreme","fgFcExtreme");
86bdcd8c 907 fgFcConservative = new TGraphAsymmErrors(nbins+1);
908 fgFcConservative->SetNameTitle("fgFcConservative","fgFcConservative");
909 }
910
65e55bbd 911
912 //
913 // Compute fc
914 //
b890d92c 915 for (Int_t ibin=1; ibin<=nbins; ibin++) {
65e55bbd 916
917 // theory_ratio = (N_b/N_c)
b890d92c 918 theoryRatio = (fhDirectMCpt->GetBinContent(ibin)>0. && fhFeedDownMCpt->GetBinContent(ibin)>0.) ?
86bdcd8c 919 fhFeedDownMCpt->GetBinContent(ibin) / fhDirectMCpt->GetBinContent(ibin) : 1.0 ;
b890d92c 920
86bdcd8c 921 //
922 // Calculate the uncertainty [ considering only the theoretical uncertainties on Nb & Nc for now !!! ]
923 //
924 // extreme A = direct-max, feed-down-min
b890d92c 925 theoryRatioExtremeA = (fhDirectMCptMax->GetBinContent(ibin)>0. && fhFeedDownMCptMin->GetBinContent(ibin)>0.) ?
86bdcd8c 926 fhFeedDownMCptMin->GetBinContent(ibin) / fhDirectMCptMax->GetBinContent(ibin) : 1.0 ;
927 // extreme B = direct-min, feed-down-max
b890d92c 928 theoryRatioExtremeB = (fhDirectMCptMin->GetBinContent(ibin)>0. && fhDirectMCptMax->GetBinContent(ibin)>0.) ?
86bdcd8c 929 fhFeedDownMCptMax->GetBinContent(ibin) / fhDirectMCptMin->GetBinContent(ibin) : 1.0 ;
930 // conservative A = direct-max, feed-down-max
b890d92c 931 theoryRatioConservativeA = (fhDirectMCptMax->GetBinContent(ibin)>0. && fhFeedDownMCptMin->GetBinContent(ibin)>0.) ?
86bdcd8c 932 fhFeedDownMCptMax->GetBinContent(ibin) / fhDirectMCptMax->GetBinContent(ibin) : 1.0 ;
933 // conservative B = direct-min, feed-down-min
b890d92c 934 theoryRatioConservativeB = (fhDirectMCptMin->GetBinContent(ibin)>0. && fhDirectMCptMax->GetBinContent(ibin)>0.) ?
86bdcd8c 935 fhFeedDownMCptMin->GetBinContent(ibin) / fhDirectMCptMin->GetBinContent(ibin) : 1.0 ;
bb427707 936
65e55bbd 937 // eff_ratio = (eff_b/eff_c)
86bdcd8c 938 effRatio = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ?
939 fhFeedDownEffpt->GetBinContent(ibin) / fhDirectEffpt->GetBinContent(ibin) : 1.0 ;
bb427707 940
65e55bbd 941 // fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
b890d92c 942 if( TMath::Abs(effRatio - 1.0)<0.01 || TMath::Abs(theoryRatio - 1.0)<0.01 ) {
943 correction = 1.0;
944 correctionExtremeA = 1.0;
945 correctionExtremeB = 1.0;
946 correctionConservativeA = 1.0;
947 correctionConservativeB = 1.0;
948 }
949 else {
950 correction = ( 1. / ( 1 + ( effRatio * theoryRatio ) ) );
951 correctionExtremeA = ( 1. / ( 1 + ( effRatio * theoryRatioExtremeA ) ) );
952 correctionExtremeB = ( 1. / ( 1 + ( effRatio * theoryRatioExtremeB ) ) );
953 correctionConservativeA = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeA ) ) );
954 correctionConservativeB = ( 1. / ( 1 + ( effRatio * theoryRatioConservativeB ) ) );
955 }
65e55bbd 956
2ee3afe2 957
8998180c 958 // fc uncertainty from (eff_b/eff_c) = fc^2 * (N_b/N_c) * delta(eff_b/eff_c)
959 // delta(eff_b/eff_c) is a percentage = fGlobalEfficiencyUncertainties[1]*effRatio
4c7dab0f 960 // correctionUnc = correction*correction * theoryRatio * fGlobalEfficiencyUncertainties[1]*effRatio;
961 // correctionExtremeAUnc = correctionExtremeA*correctionExtremeA * theoryRatioExtremeA * fGlobalEfficiencyUncertainties[1]*effRatio;
962 // correctionExtremeBUnc = correctionExtremeB*correctionExtremeB * theoryRatioExtremeB * fGlobalEfficiencyUncertainties[1]*effRatio;
963 // correctionConservativeAUnc = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA * fGlobalEfficiencyUncertainties[1]*effRatio;
964 // correctionConservativeBUnc = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB * fGlobalEfficiencyUncertainties[1]*effRatio;
965 correctionUnc = correction*correction * theoryRatio * effRatio *
966 TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] +
967 (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
968 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))
969 );
970 correctionExtremeAUnc = correctionExtremeA*correctionExtremeA * theoryRatioExtremeA * effRatio *
971 TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] +
972 (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
973 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))
974 );
975 correctionExtremeBUnc = correctionExtremeB*correctionExtremeB * theoryRatioExtremeB * effRatio *
976 TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] +
977 (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
978 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))
979 );
980 correctionConservativeAUnc = correctionConservativeA*correctionConservativeA * theoryRatioConservativeA *effRatio *
981 TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] +
982 (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
983 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))
984 );
985 correctionConservativeBUnc = correctionConservativeB*correctionConservativeB * theoryRatioConservativeB *effRatio *
986 TMath::Sqrt( fGlobalEfficiencyUncertainties[1]*fGlobalEfficiencyUncertainties[1] +
987 (fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) +
988 (fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))*(fhDirectEffpt->GetBinError(ibin)/fhDirectEffpt->GetBinContent(ibin))
989 );
8998180c 990
991
65e55bbd 992 // Fill in the histograms
a17b17dd 993 hTheoryRatio->SetBinContent(ibin,theoryRatio);
994 hEffRatio->SetBinContent(ibin,effRatio);
86bdcd8c 995 fhFc->SetBinContent(ibin,correction);
65e55bbd 996 if (fAsymUncertainties) {
997 Double_t x = fhDirectMCpt->GetBinCenter(ibin);
8998180c 998 Double_t val[4] = { correctionExtremeA + correctionExtremeAUnc, correctionExtremeA - correctionExtremeAUnc,
999 correctionExtremeB + correctionExtremeBUnc, correctionExtremeB - correctionExtremeBUnc };
1000 Double_t uncExtremeMin = correction - TMath::MinElement(4,val);
1001 Double_t uncExtremeMax = TMath::MaxElement(4,val) - correction;
86bdcd8c 1002 fgFcExtreme->SetPoint(ibin,x,correction); // i,x,y
b890d92c 1003 fgFcExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),uncExtremeMin,uncExtremeMax); // i,xl,xh,yl,yh
86bdcd8c 1004 fhFcMax->SetBinContent(ibin,correction+uncExtremeMax);
1005 fhFcMin->SetBinContent(ibin,correction-uncExtremeMin);
8998180c 1006 Double_t consval[4] = { correctionConservativeA - correctionConservativeAUnc, correctionConservativeA + correctionConservativeAUnc,
1007 correctionConservativeB - correctionConservativeBUnc, correctionConservativeB + correctionConservativeBUnc};
8998180c 1008 Double_t uncConservativeMin = correction - TMath::MinElement(4,consval);
1009 Double_t uncConservativeMax = TMath::MaxElement(4,consval) - correction;
86bdcd8c 1010 fgFcConservative->SetPoint(ibin,x,correction); // i,x,y
b890d92c 1011 fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),uncConservativeMin,uncConservativeMax); // i,xl,xh,yl,yh
066998f0 1012 if( !(correction>0.) ){
1013 fgFcExtreme->SetPoint(ibin,x,0.); // i,x,y
b890d92c 1014 fgFcExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.); // i,xl,xh,yl,yh
066998f0 1015 fgFcConservative->SetPoint(ibin,x,0.); // i,x,y
b890d92c 1016 fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.); // i,xl,xh,yl,yh
066998f0 1017 }
65e55bbd 1018 }
1019
1020 }
86df816f 1021 delete [] binwidths;
1022 delete [] limits;
65e55bbd 1023
65e55bbd 1024}
1025
1026//_________________________________________________________________________________________________________
a17b17dd 1027void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumFc(){
65e55bbd 1028 //
1029 // Compute the feed-down corrected spectrum if feed-down correction is done via fc factor (bin by bin)
bb427707 1030 // physics = reco * fc / bin-width
65e55bbd 1031 //
86bdcd8c 1032 // uncertainty: (stat) delta_physics = physics * sqrt ( (delta_reco/reco)^2 )
1033 // (syst but feed-down) delta_physics = physics * sqrt ( (delta_reco_syst/reco)^2 )
1034 // (feed-down syst) delta_physics = physics * sqrt ( (delta_fc/fc)^2 )
65e55bbd 1035 //
1036 // ( Calculation done bin by bin )
1037
1038 if (!fhFc || !fhRECpt) {
1039 AliError(" Reconstructed or fc distributions are not defined");
1040 return;
1041 }
1042
1043 Int_t nbins = fhRECpt->GetNbinsX();
86bdcd8c 1044 Double_t value = 0., errvalue = 0., errvalueMax= 0., errvalueMin= 0.;
1045 Double_t valueExtremeMax= 0., valueExtremeMin= 0.;
1046 Double_t valueConservativeMax= 0., valueConservativeMin= 0.;
65e55bbd 1047 Double_t binwidth = fhRECpt->GetBinWidth(1);
bb427707 1048 Double_t *limits = new Double_t[nbins+1];
8998180c 1049 Double_t *binwidths = new Double_t[nbins];
bb427707 1050 Double_t xlow=0.;
86bdcd8c 1051 for (Int_t i=1; i<=nbins; i++) {
bb427707 1052 binwidth = fhRECpt->GetBinWidth(i);
1053 xlow = fhRECpt->GetBinLowEdge(i);
1054 limits[i-1] = xlow;
b890d92c 1055 binwidths[i-1] = binwidth;
bb427707 1056 }
86bdcd8c 1057 limits[nbins] = xlow + binwidth;
65e55bbd 1058
1059 // declare the output histograms
b890d92c 1060 fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by fc)",nbins,limits);
1061 fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by fc)",nbins,limits);
1062 fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by fc)",nbins,limits);
65e55bbd 1063 // and the output TGraphAsymmErrors
b890d92c 1064 if (fAsymUncertainties){
1065 fgYieldCorr = new TGraphAsymmErrors(nbins+1);
1066 fgYieldCorrExtreme = new TGraphAsymmErrors(nbins+1);
1067 fgYieldCorrConservative = new TGraphAsymmErrors(nbins+1);
1068 }
65e55bbd 1069
1070 //
1071 // Do the calculation
1072 //
b890d92c 1073 for (Int_t ibin=1; ibin<=nbins; ibin++) {
65e55bbd 1074
1075 // calculate the value
86bdcd8c 1076 // physics = reco * fc / bin-width
b890d92c 1077 value = (fhRECpt->GetBinContent(ibin) && fhFc->GetBinContent(ibin)) ?
1078 fhRECpt->GetBinContent(ibin) * fhFc->GetBinContent(ibin) : 0. ;
bb427707 1079 value /= fhRECpt->GetBinWidth(ibin) ;
86bdcd8c 1080
1081 // Statistical uncertainty
1082 // (stat) delta_physics = physics * sqrt ( (delta_reco/reco)^2 )
b890d92c 1083 errvalue = (value!=0. && fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) ?
1084 value * (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) : 0. ;
86bdcd8c 1085
1086 // Calculate the systematic uncertainties
1087 // (syst but feed-down) delta_physics = physics * sqrt ( (delta_reco_syst/reco)^2 )
1088 // (feed-down syst) delta_physics = physics * sqrt ( (delta_fc/fc)^2 )
1089 //
65e55bbd 1090 // Protect against null denominator. If so, define uncertainty as null
1091 if (fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) {
1092
1093 if (fAsymUncertainties) {
1094
86bdcd8c 1095 // Systematics but feed-down
1096 if (fgRECSystematics) {
1097 errvalueMax = value * ( fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinContent(ibin) );
1098 errvalueMin = value * ( fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinContent(ibin) );
65e55bbd 1099 }
d38da1b0 1100 else { errvalueMax = 0.; errvalueMin = 0.; }
86bdcd8c 1101
1102 // Extreme feed-down systematics
1103 valueExtremeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcExtreme->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1104 valueExtremeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcExtreme->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1105
1106 // Conservative feed-down systematics
1107 valueConservativeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcConservative->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1108 valueConservativeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcConservative->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
65e55bbd 1109
1110 }
65e55bbd 1111
1112 }
86bdcd8c 1113 else { errvalueMax = 0.; errvalueMin = 0.; }
65e55bbd 1114
1115 // fill in the histograms
86bdcd8c 1116 fhYieldCorr->SetBinContent(ibin,value);
1117 fhYieldCorr->SetBinError(ibin,errvalue);
65e55bbd 1118 if (fAsymUncertainties) {
86bdcd8c 1119 Double_t center = fhYieldCorr->GetBinCenter(ibin);
65e55bbd 1120 fgYieldCorr->SetPoint(ibin,center,value); // i,x,y
b890d92c 1121 fgYieldCorr->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
86bdcd8c 1122 fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax);
1123 fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
1124 fgYieldCorrExtreme->SetPoint(ibin,center,value);
b890d92c 1125 fgYieldCorrExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),value-valueExtremeMin,valueExtremeMax-value);
86bdcd8c 1126 fgYieldCorrConservative->SetPoint(ibin,center,value);
b890d92c 1127 fgYieldCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),value-valueConservativeMin,valueConservativeMax-value);
65e55bbd 1128 }
1129
1130 }
86df816f 1131 delete [] binwidths;
1132 delete [] limits;
65e55bbd 1133
65e55bbd 1134}
1135
1136//_________________________________________________________________________________________________________
a17b17dd 1137void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumNb(Double_t deltaY, Double_t branchingRatioBintoFinalDecay) {
65e55bbd 1138 //
1139 // Compute the feed-down corrected spectrum if feed-down correction is done via Nb (bin by bin)
bb427707 1140 // physics = [ reco - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) ] / bin-width
65e55bbd 1141 //
86bdcd8c 1142 // uncertainty: (stat) delta_physics = sqrt ( (delta_reco)^2 ) / bin-width
1143 // (syst but feed-down) delta_physics = sqrt ( (delta_reco_syst)^2 ) / bin-width
1144 // (feed-down syst) delta_physics = sqrt ( (k*delta_lumi/lumi)^2 + (k*delta_eff_trig/eff_trig)^2
8998180c 1145 // + (k*delta_Nb/Nb)^2 + (k*delta_eff/eff)^2 + (k*global_eff_ratio)^2 ) / bin-width
65e55bbd 1146 // where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th
1147 //
1148
1149 Int_t nbins = fhRECpt->GetNbinsX();
1150 Double_t binwidth = fhRECpt->GetBinWidth(1);
86bdcd8c 1151 Double_t value = 0., errvalue = 0., errvalueMax = 0., errvalueMin = 0., kfactor = 0.;
1152 Double_t errvalueExtremeMax = 0., errvalueExtremeMin = 0.;
bb427707 1153 Double_t *limits = new Double_t[nbins+1];
b890d92c 1154 Double_t *binwidths = new Double_t[nbins];
bb427707 1155 Double_t xlow=0.;
86bdcd8c 1156 for (Int_t i=1; i<=nbins; i++) {
bb427707 1157 binwidth = fhRECpt->GetBinWidth(i);
1158 xlow = fhRECpt->GetBinLowEdge(i);
1159 limits[i-1] = xlow;
b890d92c 1160 binwidths[i-1] = binwidth;
bb427707 1161 }
86bdcd8c 1162 limits[nbins] = xlow + binwidth;
65e55bbd 1163
1164 // declare the output histograms
b890d92c 1165 fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by Nb)",nbins,limits);
1166 fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by Nb)",nbins,limits);
1167 fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by Nb)",nbins,limits);
65e55bbd 1168 // and the output TGraphAsymmErrors
b890d92c 1169 if (fAsymUncertainties){
1170 fgYieldCorr = new TGraphAsymmErrors(nbins+1);
1171 fgYieldCorrExtreme = new TGraphAsymmErrors(nbins+1);
86bdcd8c 1172 fgYieldCorrConservative = new TGraphAsymmErrors(nbins+1);
b890d92c 1173 // Define fc-conservative
1174 fgFcConservative = new TGraphAsymmErrors(nbins+1);
86bdcd8c 1175 AliInfo(" Beware the conservative & extreme uncertainties are equal by definition !");
1176 }
65e55bbd 1177
b890d92c 1178 // variables to define fc-conservative
066998f0 1179 double correction=0, correctionMax=0., correctionMin=0.;
1180
65e55bbd 1181 //
1182 // Do the calculation
1183 //
b890d92c 1184 for (Int_t ibin=1; ibin<=nbins; ibin++) {
65e55bbd 1185
86bdcd8c 1186 // Calculate the value
1187 // physics = [ reco - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) ] / bin-width
e1015a30 1188 value = ( fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0. &&
1189 fhFeedDownMCpt->GetBinContent(ibin)>0. && fhFeedDownEffpt->GetBinContent(ibin)>0. ) ?
b890d92c 1190 fhRECpt->GetBinContent(ibin) - (deltaY*branchingRatioBintoFinalDecay*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) )
1191 : 0. ;
bb427707 1192 value /= fhRECpt->GetBinWidth(ibin);
65e55bbd 1193
86bdcd8c 1194 // Statistical uncertainty: delta_physics = sqrt ( (delta_reco)^2 ) / bin-width
e1015a30 1195 errvalue = (value!=0. && fhRECpt->GetBinError(ibin) && fhRECpt->GetBinError(ibin)!=0.) ?
b890d92c 1196 fhRECpt->GetBinError(ibin) : 0.;
86bdcd8c 1197 errvalue /= fhRECpt->GetBinWidth(ibin);
65e55bbd 1198
066998f0 1199 // Correction (fc) : Estimate of the relative amount feed-down subtracted
1200 // correction = [ 1 - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th)/reco ]
1201 correction = 1 - (deltaY*branchingRatioBintoFinalDecay*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) ) / fhRECpt->GetBinContent(ibin) ;
1202
86bdcd8c 1203 // Systematic uncertainties
1204 // (syst but feed-down) delta_physics = sqrt ( (delta_reco_syst)^2 ) / bin-width
1205 // (feed-down syst) delta_physics = sqrt ( (k*delta_lumi/lumi)^2 + (k*delta_eff_trig/eff_trig)^2
8998180c 1206 // + (k*delta_Nb/Nb)^2 + (k*delta_eff/eff)^2 + (k*global_eff_ratio)^2 ) / bin-width
86bdcd8c 1207 // where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th
1208 kfactor = deltaY*branchingRatioBintoFinalDecay*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) ;
2ee3afe2 1209
86bdcd8c 1210 //
65e55bbd 1211 if (fAsymUncertainties) {
e52da743 1212 Double_t nb = fhFeedDownMCpt->GetBinContent(ibin);
1213 Double_t nbDmax = fhFeedDownMCptMax->GetBinContent(ibin) - fhFeedDownMCpt->GetBinContent(ibin);
1214 Double_t nbDmin = fhFeedDownMCpt->GetBinContent(ibin) - fhFeedDownMCptMin->GetBinContent(ibin);
86bdcd8c 1215
1216 // Systematics but feed-down
1217 if (fgRECSystematics){
1218 errvalueMax = fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinWidth(ibin) ;
1219 errvalueMin = fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinWidth(ibin);
1220 }
d38da1b0 1221 else { errvalueMax = 0.; errvalueMin = 0.; }
86bdcd8c 1222
1223 // Feed-down systematics
1224 // min value with the maximum Nb
1225 errvalueExtremeMin = TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
1226 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
e52da743 1227 ( (kfactor*nbDmax/nb)*(kfactor*nbDmax/nb) ) +
4c7dab0f 1228 ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
8998180c 1229 ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1230 ) / fhRECpt->GetBinWidth(ibin);
86bdcd8c 1231 // max value with the minimum Nb
1232 errvalueExtremeMax = TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
1233 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
e52da743 1234 ( (kfactor*nbDmin/nb)*(kfactor*nbDmin/nb) ) +
4c7dab0f 1235 ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
8998180c 1236 ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1237 ) / fhRECpt->GetBinWidth(ibin);
066998f0 1238
1239 // Correction systematics (fc)
1240 // min value with the maximum Nb
1241 correctionMin = TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
1242 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1243 ( (kfactor*nbDmax/nb)*(kfactor*nbDmax/nb) ) +
1244 ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
1245 ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1246 ) / fhRECpt->GetBinContent(ibin) ;
1247 // max value with the minimum Nb
1248 correctionMax = TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
1249 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
1250 ( (kfactor*nbDmin/nb)*(kfactor*nbDmin/nb) ) +
1251 ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
1252 ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1253 ) / fhRECpt->GetBinContent(ibin) ;
65e55bbd 1254 }
1255 else{ // Don't consider Nb uncertainty in this case [ to be tested!!! ]
86bdcd8c 1256 errvalueExtremeMax = TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
1257 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
4c7dab0f 1258 ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
8998180c 1259 ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1260 ) / fhRECpt->GetBinWidth(ibin);
86bdcd8c 1261 errvalueExtremeMin = errvalueExtremeMax ;
65e55bbd 1262 }
1263
2ee3afe2 1264
65e55bbd 1265 // fill in histograms
86bdcd8c 1266 fhYieldCorr->SetBinContent(ibin,value);
1267 fhYieldCorr->SetBinError(ibin,errvalue);
65e55bbd 1268 if (fAsymUncertainties) {
86bdcd8c 1269 Double_t x = fhYieldCorr->GetBinCenter(ibin);
65e55bbd 1270 fgYieldCorr->SetPoint(ibin,x,value); // i,x,y
b890d92c 1271 fgYieldCorr->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
86bdcd8c 1272 fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax);
1273 fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
1274 fgYieldCorrExtreme->SetPoint(ibin,x,value); // i,x,y
b890d92c 1275 fgYieldCorrExtreme->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
86bdcd8c 1276 fgYieldCorrConservative->SetPoint(ibin,x,value); // i,x,y
b890d92c 1277 fgYieldCorrConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
066998f0 1278 // cout << " bin " << ibin << ", correction " << correction << ", min correction unc " << correctionMin << ", max correction unc " << correctionMax << endl;
1279 if(correction>0.){
1280 fgFcConservative->SetPoint(ibin,x,correction);
b890d92c 1281 fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),correctionMin,correctionMax);
066998f0 1282 }
1283 else{
1284 fgFcConservative->SetPoint(ibin,x,0.);
b890d92c 1285 fgFcConservative->SetPointError(ibin,(binwidths[ibin-1]/2.),(binwidths[ibin-1]/2.),0.,0.);
066998f0 1286 }
65e55bbd 1287 }
1288
1289 }
86df816f 1290 delete [] binwidths;
1291 delete [] limits;
65e55bbd 1292
65e55bbd 1293}
1294
1295
1296//_________________________________________________________________________________________________________
8998180c 1297void AliHFPtSpectrum::ComputeSystUncertainties(Int_t decay, Bool_t combineFeedDown) {
e52da743 1298 //
1299 // Function that re-calculates the global systematic uncertainties
1300 // by calling the class AliHFSystErr and combining those
1301 // (in quadrature) with the feed-down subtraction uncertainties
1302 //
8998180c 1303
1304 // Call the systematics uncertainty class for a given decay
1305 AliHFSystErr systematics(decay);
1306
1307 // Estimate the feed-down uncertainty in percentage
1308 Int_t nentries = fgSigmaCorrConservative->GetN();
1309 TGraphAsymmErrors *grErrFeeddown = new TGraphAsymmErrors(nentries);
1310 Double_t x=0., y=0., errx=0., erryl=0., erryh=0;
1311 for(Int_t i=0; i<nentries; i++) {
1312 fgSigmaCorrConservative->GetPoint(i,x,y);
1313 errx = fgSigmaCorrConservative->GetErrorXlow(i) ;
1314 erryl = fgSigmaCorrConservative->GetErrorYlow(i) / y ;
1315 erryh = fgSigmaCorrConservative->GetErrorYhigh(i) / y ;
1316 grErrFeeddown->SetPoint(i,x,0.);
1317 grErrFeeddown->SetPointError(i,errx,errx,erryl,erryh); //i, xl, xh, yl, yh
1318 }
1319
1320 // Draw all the systematics independently
1321 systematics.DrawErrors(grErrFeeddown);
1322
1323 // Set the sigma systematic uncertainties
1324 // possibly combine with the feed-down uncertainties
1325 Double_t errylcomb=0., erryhcomb=0;
1326 for(Int_t i=1; i<nentries; i++) {
1327 fgSigmaCorr->GetPoint(i,x,y);
1328 errx = grErrFeeddown->GetErrorXlow(i) ;
1329 erryl = grErrFeeddown->GetErrorYlow(i);
1330 erryh = grErrFeeddown->GetErrorYhigh(i);
1331 if (combineFeedDown) {
1332 errylcomb = systematics.GetTotalSystErr(x,erryl) * y ;
1333 erryhcomb = systematics.GetTotalSystErr(x,erryh) * y ;
1334 } else {
1335 errylcomb = systematics.GetTotalSystErr(x) * y ;
1336 erryhcomb = systematics.GetTotalSystErr(x) * y ;
1337 }
1338 fgSigmaCorr->SetPointError(i,errx,errx,errylcomb,erryhcomb);
1339 }
1340
1341}
1342
1343
1344//_________________________________________________________________________________________________________
1345void AliHFPtSpectrum::DrawSpectrum(TGraphAsymmErrors *gPrediction) {
e52da743 1346 //
1347 // Example method to draw the corrected spectrum & the theoretical prediction
1348 //
8998180c 1349
1350 TCanvas *csigma = new TCanvas("csigma","Draw the corrected cross-section & the prediction");
1351 csigma->SetFillColor(0);
1352 gPrediction->GetXaxis()->SetTitleSize(0.05);
1353 gPrediction->GetXaxis()->SetTitleOffset(0.95);
1354 gPrediction->GetYaxis()->SetTitleSize(0.05);
1355 gPrediction->GetYaxis()->SetTitleOffset(0.95);
1356 gPrediction->GetXaxis()->SetTitle("p_{T} [GeV]");
1357 gPrediction->GetYaxis()->SetTitle("BR #times #frac{d#sigma}{dp_{T}} |_{|y|<0.5} [pb/GeV]");
1358 gPrediction->SetLineColor(kGreen+2);
1359 gPrediction->SetLineWidth(3);
1360 gPrediction->SetFillColor(kGreen+1);
1361 gPrediction->Draw("3CA");
1362 fgSigmaCorr->SetLineColor(kRed);
1363 fgSigmaCorr->SetLineWidth(1);
1364 fgSigmaCorr->SetFillColor(kRed);
1365 fgSigmaCorr->SetFillStyle(0);
1366 fgSigmaCorr->Draw("2");
1367 fhSigmaCorr->SetMarkerColor(kRed);
1368 fhSigmaCorr->Draw("esame");
1369 csigma->SetLogy();
1370 TLegend * leg = new TLegend(0.7,0.75,0.87,0.5);
1371 leg->SetBorderSize(0);
1372 leg->SetLineColor(0);
1373 leg->SetFillColor(0);
1374 leg->SetTextFont(42);
1375 leg->AddEntry(gPrediction,"FONLL ","fl");
1376 leg->AddEntry(fhSigmaCorr,"data stat. unc.","pl");
1377 leg->AddEntry(fgSigmaCorr,"data syst. unc.","f");
1378 leg->Draw();
1379 csigma->Draw();
1380
1381}
1382
1383//_________________________________________________________________________________________________________
86bdcd8c 1384TH1D * AliHFPtSpectrum::ReweightHisto(TH1D *hToReweight, TH1D *hReference){
65e55bbd 1385 //
1386 // Function to reweight histograms for testing purposes:
1387 // This function takes the histo hToReweight and reweights
1388 // it (its pt shape) with respect to hReference
1389 //
1390
1391 // check histograms consistency
1392 Bool_t areconsistent=kTRUE;
1393 areconsistent &= CheckHistosConsistency(hToReweight,hReference);
1394 if (!areconsistent) {
1395 AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
1396 return NULL;
1397 }
1398
1399 // define a new empty histogram
86bdcd8c 1400 TH1D *hReweighted = (TH1D*)hToReweight->Clone("hReweighted");
65e55bbd 1401 hReweighted->Reset();
1402 Double_t weight=1.0;
1403 Double_t yvalue=1.0;
a17b17dd 1404 Double_t integralRef = hReference->Integral();
1405 Double_t integralH = hToReweight->Integral();
65e55bbd 1406
1407 // now reweight the spectra
1408 //
1409 // the weight is the relative probability of the given pt bin in the reference histo
1410 // divided by its relative probability (to normalize it) on the histo to re-weight
1411 for (Int_t i=0; i<=hToReweight->GetNbinsX(); i++) {
a17b17dd 1412 weight = (hReference->GetBinContent(i)/integralRef) / (hToReweight->GetBinContent(i)/integralH) ;
65e55bbd 1413 yvalue = hToReweight->GetBinContent(i);
1414 hReweighted->SetBinContent(i,yvalue*weight);
1415 }
1416
86bdcd8c 1417 return (TH1D*)hReweighted;
65e55bbd 1418}
1419
1420//_________________________________________________________________________________________________________
86bdcd8c 1421TH1D * AliHFPtSpectrum::ReweightRecHisto(TH1D *hRecToReweight, TH1D *hMCToReweight, TH1D *hMCReference){
65e55bbd 1422 //
1423 // Function to reweight histograms for testing purposes:
1424 // This function takes the histo hToReweight and reweights
1425 // it (its pt shape) with respect to hReference /hMCToReweight
1426 //
1427
1428 // check histograms consistency
1429 Bool_t areconsistent=kTRUE;
1430 areconsistent &= CheckHistosConsistency(hMCToReweight,hMCReference);
1431 areconsistent &= CheckHistosConsistency(hRecToReweight,hMCReference);
1432 if (!areconsistent) {
1433 AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
1434 return NULL;
1435 }
1436
1437 // define a new empty histogram
86bdcd8c 1438 TH1D *hReweighted = (TH1D*)hMCToReweight->Clone("hReweighted");
65e55bbd 1439 hReweighted->Reset();
86bdcd8c 1440 TH1D *hRecReweighted = (TH1D*)hRecToReweight->Clone("hRecReweighted");
65e55bbd 1441 hRecReweighted->Reset();
1442 Double_t weight=1.0;
1443 Double_t yvalue=1.0, yrecvalue=1.0;
a17b17dd 1444 Double_t integralRef = hMCReference->Integral();
1445 Double_t integralH = hMCToReweight->Integral();
65e55bbd 1446
1447 // now reweight the spectra
1448 //
1449 // the weight is the relative probability of the given pt bin
1450 // that should be applied in the MC histo to get the reference histo shape
1451 // Probabilities are properly normalized.
1452 for (Int_t i=0; i<=hMCToReweight->GetNbinsX(); i++) {
a17b17dd 1453 weight = (hMCReference->GetBinContent(i)/integralRef) / (hMCToReweight->GetBinContent(i)/integralH) ;
65e55bbd 1454 yvalue = hMCToReweight->GetBinContent(i);
1455 hReweighted->SetBinContent(i,yvalue*weight);
1456 yrecvalue = hRecToReweight->GetBinContent(i);
1457 hRecReweighted->SetBinContent(i,yrecvalue*weight);
1458 }
1459
86bdcd8c 1460 return (TH1D*)hRecReweighted;
65e55bbd 1461}
1462