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