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