in case of less centrality bins histograms requested than centrality bins provided...
[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 //
459 // Set the reconstructed spectrum systematic uncertainties
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
65e55bbd 989 }
990
991 }
992
65e55bbd 993}
994
995//_________________________________________________________________________________________________________
a17b17dd 996void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumFc(){
65e55bbd 997 //
998 // Compute the feed-down corrected spectrum if feed-down correction is done via fc factor (bin by bin)
bb427707 999 // physics = reco * fc / bin-width
65e55bbd 1000 //
86bdcd8c 1001 // uncertainty: (stat) delta_physics = physics * sqrt ( (delta_reco/reco)^2 )
1002 // (syst but feed-down) delta_physics = physics * sqrt ( (delta_reco_syst/reco)^2 )
1003 // (feed-down syst) delta_physics = physics * sqrt ( (delta_fc/fc)^2 )
65e55bbd 1004 //
1005 // ( Calculation done bin by bin )
1006
1007 if (!fhFc || !fhRECpt) {
1008 AliError(" Reconstructed or fc distributions are not defined");
1009 return;
1010 }
1011
1012 Int_t nbins = fhRECpt->GetNbinsX();
86bdcd8c 1013 Double_t value = 0., errvalue = 0., errvalueMax= 0., errvalueMin= 0.;
1014 Double_t valueExtremeMax= 0., valueExtremeMin= 0.;
1015 Double_t valueConservativeMax= 0., valueConservativeMin= 0.;
65e55bbd 1016 Double_t binwidth = fhRECpt->GetBinWidth(1);
bb427707 1017 Double_t *limits = new Double_t[nbins+1];
8998180c 1018 Double_t *binwidths = new Double_t[nbins];
bb427707 1019 Double_t xlow=0.;
86bdcd8c 1020 for (Int_t i=1; i<=nbins; i++) {
bb427707 1021 binwidth = fhRECpt->GetBinWidth(i);
1022 xlow = fhRECpt->GetBinLowEdge(i);
1023 limits[i-1] = xlow;
de97013c 1024 binwidths[i] = binwidth;
bb427707 1025 }
86bdcd8c 1026 limits[nbins] = xlow + binwidth;
65e55bbd 1027
1028 // declare the output histograms
86bdcd8c 1029 if (!fhYieldCorr) fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by fc)",nbins,limits);
1030 if (!fhYieldCorrMax) fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by fc)",nbins,limits);
1031 if (!fhYieldCorrMin) fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by fc)",nbins,limits);
65e55bbd 1032 // and the output TGraphAsymmErrors
86bdcd8c 1033 if (fAsymUncertainties & !fgYieldCorr) fgYieldCorr = new TGraphAsymmErrors(nbins+1);
1034 if (fAsymUncertainties & !fgYieldCorrExtreme) fgYieldCorrExtreme = new TGraphAsymmErrors(nbins+1);
1035 if (fAsymUncertainties & !fgYieldCorrConservative) fgYieldCorrConservative = new TGraphAsymmErrors(nbins+1);
65e55bbd 1036
1037 //
1038 // Do the calculation
1039 //
1040 for (Int_t ibin=0; ibin<=nbins; ibin++) {
1041
1042 // calculate the value
86bdcd8c 1043 // physics = reco * fc / bin-width
65e55bbd 1044 value = fhRECpt->GetBinContent(ibin) * fhFc->GetBinContent(ibin) ;
bb427707 1045 value /= fhRECpt->GetBinWidth(ibin) ;
86bdcd8c 1046
1047 // Statistical uncertainty
1048 // (stat) delta_physics = physics * sqrt ( (delta_reco/reco)^2 )
1049 errvalue = value * (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) ;
1050
1051 // Calculate the systematic uncertainties
1052 // (syst but feed-down) delta_physics = physics * sqrt ( (delta_reco_syst/reco)^2 )
1053 // (feed-down syst) delta_physics = physics * sqrt ( (delta_fc/fc)^2 )
1054 //
65e55bbd 1055 // Protect against null denominator. If so, define uncertainty as null
1056 if (fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) {
1057
1058 if (fAsymUncertainties) {
1059
86bdcd8c 1060 // Systematics but feed-down
1061 if (fgRECSystematics) {
1062 errvalueMax = value * ( fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinContent(ibin) );
1063 errvalueMin = value * ( fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinContent(ibin) );
65e55bbd 1064 }
d38da1b0 1065 else { errvalueMax = 0.; errvalueMin = 0.; }
86bdcd8c 1066
1067 // Extreme feed-down systematics
1068 valueExtremeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcExtreme->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1069 valueExtremeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcExtreme->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1070
1071 // Conservative feed-down systematics
1072 valueConservativeMax = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) + fgFcConservative->GetErrorYhigh(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
1073 valueConservativeMin = fhRECpt->GetBinContent(ibin) * ( fhFc->GetBinContent(ibin) - fgFcConservative->GetErrorYlow(ibin) ) / fhRECpt->GetBinWidth(ibin) ;
65e55bbd 1074
1075 }
65e55bbd 1076
1077 }
86bdcd8c 1078 else { errvalueMax = 0.; errvalueMin = 0.; }
65e55bbd 1079
1080 // fill in the histograms
86bdcd8c 1081 fhYieldCorr->SetBinContent(ibin,value);
1082 fhYieldCorr->SetBinError(ibin,errvalue);
65e55bbd 1083 if (fAsymUncertainties) {
86bdcd8c 1084 Double_t center = fhYieldCorr->GetBinCenter(ibin);
65e55bbd 1085 fgYieldCorr->SetPoint(ibin,center,value); // i,x,y
86bdcd8c 1086 fgYieldCorr->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
1087 fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax);
1088 fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
1089 fgYieldCorrExtreme->SetPoint(ibin,center,value);
1090 fgYieldCorrExtreme->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),value-valueExtremeMin,valueExtremeMax-value);
1091 fgYieldCorrConservative->SetPoint(ibin,center,value);
1092 fgYieldCorrConservative->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),value-valueConservativeMin,valueConservativeMax-value);
65e55bbd 1093 }
1094
1095 }
1096
65e55bbd 1097}
1098
1099//_________________________________________________________________________________________________________
a17b17dd 1100void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumNb(Double_t deltaY, Double_t branchingRatioBintoFinalDecay) {
65e55bbd 1101 //
1102 // Compute the feed-down corrected spectrum if feed-down correction is done via Nb (bin by bin)
bb427707 1103 // physics = [ reco - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) ] / bin-width
65e55bbd 1104 //
86bdcd8c 1105 // uncertainty: (stat) delta_physics = sqrt ( (delta_reco)^2 ) / bin-width
1106 // (syst but feed-down) delta_physics = sqrt ( (delta_reco_syst)^2 ) / bin-width
1107 // (feed-down syst) delta_physics = sqrt ( (k*delta_lumi/lumi)^2 + (k*delta_eff_trig/eff_trig)^2
8998180c 1108 // + (k*delta_Nb/Nb)^2 + (k*delta_eff/eff)^2 + (k*global_eff_ratio)^2 ) / bin-width
65e55bbd 1109 // where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th
1110 //
1111
1112 Int_t nbins = fhRECpt->GetNbinsX();
1113 Double_t binwidth = fhRECpt->GetBinWidth(1);
86bdcd8c 1114 Double_t value = 0., errvalue = 0., errvalueMax = 0., errvalueMin = 0., kfactor = 0.;
1115 Double_t errvalueExtremeMax = 0., errvalueExtremeMin = 0.;
bb427707 1116 Double_t *limits = new Double_t[nbins+1];
de97013c 1117 Double_t *binwidths = new Double_t[nbins+1];
bb427707 1118 Double_t xlow=0.;
86bdcd8c 1119 for (Int_t i=1; i<=nbins; i++) {
bb427707 1120 binwidth = fhRECpt->GetBinWidth(i);
1121 xlow = fhRECpt->GetBinLowEdge(i);
1122 limits[i-1] = xlow;
de97013c 1123 binwidths[i] = binwidth;
bb427707 1124 }
86bdcd8c 1125 limits[nbins] = xlow + binwidth;
65e55bbd 1126
1127 // declare the output histograms
86bdcd8c 1128 if (!fhYieldCorr) fhYieldCorr = new TH1D("fhYieldCorr","corrected yield (by Nb)",nbins,limits);
1129 if (!fhYieldCorrMax) fhYieldCorrMax = new TH1D("fhYieldCorrMax","max corrected yield (by Nb)",nbins,limits);
1130 if (!fhYieldCorrMin) fhYieldCorrMin = new TH1D("fhYieldCorrMin","min corrected yield (by Nb)",nbins,limits);
65e55bbd 1131 // and the output TGraphAsymmErrors
86bdcd8c 1132 if (fAsymUncertainties & !fgYieldCorr) fgYieldCorr = new TGraphAsymmErrors(nbins+1);
1133 if (fAsymUncertainties & !fgYieldCorrExtreme) fgYieldCorrExtreme = new TGraphAsymmErrors(nbins+1);
1134 if (fAsymUncertainties & !fgYieldCorrConservative) {
1135 fgYieldCorrConservative = new TGraphAsymmErrors(nbins+1);
1136 AliInfo(" Beware the conservative & extreme uncertainties are equal by definition !");
1137 }
65e55bbd 1138
1139 //
1140 // Do the calculation
1141 //
1142 for (Int_t ibin=0; ibin<=nbins; ibin++) {
1143
86bdcd8c 1144 // Calculate the value
1145 // physics = [ reco - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th) ] / bin-width
a17b17dd 1146 value = fhRECpt->GetBinContent(ibin) - (deltaY*branchingRatioBintoFinalDecay*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) );
bb427707 1147 value /= fhRECpt->GetBinWidth(ibin);
65e55bbd 1148
86bdcd8c 1149 // Statistical uncertainty: delta_physics = sqrt ( (delta_reco)^2 ) / bin-width
1150 errvalue = fhRECpt->GetBinError(ibin);
1151 errvalue /= fhRECpt->GetBinWidth(ibin);
65e55bbd 1152
86bdcd8c 1153 // Systematic uncertainties
1154 // (syst but feed-down) delta_physics = sqrt ( (delta_reco_syst)^2 ) / bin-width
1155 // (feed-down syst) delta_physics = sqrt ( (k*delta_lumi/lumi)^2 + (k*delta_eff_trig/eff_trig)^2
8998180c 1156 // + (k*delta_Nb/Nb)^2 + (k*delta_eff/eff)^2 + (k*global_eff_ratio)^2 ) / bin-width
86bdcd8c 1157 // where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th
1158 kfactor = deltaY*branchingRatioBintoFinalDecay*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) ;
2ee3afe2 1159
86bdcd8c 1160 //
65e55bbd 1161 if (fAsymUncertainties) {
8998180c 1162 Double_t Nb = fhFeedDownMCpt->GetBinContent(ibin);
1163 Double_t NbDmax = fhFeedDownMCptMax->GetBinContent(ibin) - fhFeedDownMCpt->GetBinContent(ibin);
1164 Double_t NbDmin = fhFeedDownMCpt->GetBinContent(ibin) - fhFeedDownMCptMin->GetBinContent(ibin);
86bdcd8c 1165
1166 // Systematics but feed-down
1167 if (fgRECSystematics){
1168 errvalueMax = fgRECSystematics->GetErrorYhigh(ibin) / fhRECpt->GetBinWidth(ibin) ;
1169 errvalueMin = fgRECSystematics->GetErrorYlow(ibin) / fhRECpt->GetBinWidth(ibin);
1170 }
d38da1b0 1171 else { errvalueMax = 0.; errvalueMin = 0.; }
86bdcd8c 1172
1173 // Feed-down systematics
1174 // min value with the maximum Nb
1175 errvalueExtremeMin = TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
1176 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
8998180c 1177 ( (kfactor*NbDmax/Nb)*(kfactor*NbDmax/Nb) ) +
4c7dab0f 1178 ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
8998180c 1179 ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1180 ) / fhRECpt->GetBinWidth(ibin);
86bdcd8c 1181 // max value with the minimum Nb
1182 errvalueExtremeMax = TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
1183 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
8998180c 1184 ( (kfactor*NbDmin/Nb)*(kfactor*NbDmin/Nb) ) +
4c7dab0f 1185 ( (kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin))*(kfactor*fhFeedDownEffpt->GetBinError(ibin)/fhFeedDownEffpt->GetBinContent(ibin)) ) +
8998180c 1186 ( (kfactor*fGlobalEfficiencyUncertainties[1])*(kfactor*fGlobalEfficiencyUncertainties[1]) )
1187 ) / fhRECpt->GetBinWidth(ibin);
65e55bbd 1188 }
1189 else{ // Don't consider Nb uncertainty in this case [ to be tested!!! ]
86bdcd8c 1190 errvalueExtremeMax = TMath::Sqrt( ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
1191 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
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 errvalueExtremeMin = errvalueExtremeMax ;
65e55bbd 1196 }
1197
2ee3afe2 1198
65e55bbd 1199 // fill in histograms
86bdcd8c 1200 fhYieldCorr->SetBinContent(ibin,value);
1201 fhYieldCorr->SetBinError(ibin,errvalue);
65e55bbd 1202 if (fAsymUncertainties) {
86bdcd8c 1203 Double_t x = fhYieldCorr->GetBinCenter(ibin);
65e55bbd 1204 fgYieldCorr->SetPoint(ibin,x,value); // i,x,y
86bdcd8c 1205 fgYieldCorr->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),errvalueMin,errvalueMax); // i,xl,xh,yl,yh
1206 fhYieldCorrMax->SetBinContent(ibin,value+errvalueMax);
1207 fhYieldCorrMin->SetBinContent(ibin,value-errvalueMin);
1208 fgYieldCorrExtreme->SetPoint(ibin,x,value); // i,x,y
1209 fgYieldCorrExtreme->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
1210 fgYieldCorrConservative->SetPoint(ibin,x,value); // i,x,y
1211 fgYieldCorrConservative->SetPointError(ibin,(binwidths[ibin]/2.),(binwidths[ibin]/2.),errvalueExtremeMin,errvalueExtremeMax); // i,xl,xh,yl,yh
65e55bbd 1212 }
1213
1214 }
1215
65e55bbd 1216}
1217
1218
1219//_________________________________________________________________________________________________________
8998180c 1220void AliHFPtSpectrum::ComputeSystUncertainties(Int_t decay, Bool_t combineFeedDown) {
1221
1222 // Call the systematics uncertainty class for a given decay
1223 AliHFSystErr systematics(decay);
1224
1225 // Estimate the feed-down uncertainty in percentage
1226 Int_t nentries = fgSigmaCorrConservative->GetN();
1227 TGraphAsymmErrors *grErrFeeddown = new TGraphAsymmErrors(nentries);
1228 Double_t x=0., y=0., errx=0., erryl=0., erryh=0;
1229 for(Int_t i=0; i<nentries; i++) {
1230 fgSigmaCorrConservative->GetPoint(i,x,y);
1231 errx = fgSigmaCorrConservative->GetErrorXlow(i) ;
1232 erryl = fgSigmaCorrConservative->GetErrorYlow(i) / y ;
1233 erryh = fgSigmaCorrConservative->GetErrorYhigh(i) / y ;
1234 grErrFeeddown->SetPoint(i,x,0.);
1235 grErrFeeddown->SetPointError(i,errx,errx,erryl,erryh); //i, xl, xh, yl, yh
1236 }
1237
1238 // Draw all the systematics independently
1239 systematics.DrawErrors(grErrFeeddown);
1240
1241 // Set the sigma systematic uncertainties
1242 // possibly combine with the feed-down uncertainties
1243 Double_t errylcomb=0., erryhcomb=0;
1244 for(Int_t i=1; i<nentries; i++) {
1245 fgSigmaCorr->GetPoint(i,x,y);
1246 errx = grErrFeeddown->GetErrorXlow(i) ;
1247 erryl = grErrFeeddown->GetErrorYlow(i);
1248 erryh = grErrFeeddown->GetErrorYhigh(i);
1249 if (combineFeedDown) {
1250 errylcomb = systematics.GetTotalSystErr(x,erryl) * y ;
1251 erryhcomb = systematics.GetTotalSystErr(x,erryh) * y ;
1252 } else {
1253 errylcomb = systematics.GetTotalSystErr(x) * y ;
1254 erryhcomb = systematics.GetTotalSystErr(x) * y ;
1255 }
1256 fgSigmaCorr->SetPointError(i,errx,errx,errylcomb,erryhcomb);
1257 }
1258
1259}
1260
1261
1262//_________________________________________________________________________________________________________
1263void AliHFPtSpectrum::DrawSpectrum(TGraphAsymmErrors *gPrediction) {
1264
1265 TCanvas *csigma = new TCanvas("csigma","Draw the corrected cross-section & the prediction");
1266 csigma->SetFillColor(0);
1267 gPrediction->GetXaxis()->SetTitleSize(0.05);
1268 gPrediction->GetXaxis()->SetTitleOffset(0.95);
1269 gPrediction->GetYaxis()->SetTitleSize(0.05);
1270 gPrediction->GetYaxis()->SetTitleOffset(0.95);
1271 gPrediction->GetXaxis()->SetTitle("p_{T} [GeV]");
1272 gPrediction->GetYaxis()->SetTitle("BR #times #frac{d#sigma}{dp_{T}} |_{|y|<0.5} [pb/GeV]");
1273 gPrediction->SetLineColor(kGreen+2);
1274 gPrediction->SetLineWidth(3);
1275 gPrediction->SetFillColor(kGreen+1);
1276 gPrediction->Draw("3CA");
1277 fgSigmaCorr->SetLineColor(kRed);
1278 fgSigmaCorr->SetLineWidth(1);
1279 fgSigmaCorr->SetFillColor(kRed);
1280 fgSigmaCorr->SetFillStyle(0);
1281 fgSigmaCorr->Draw("2");
1282 fhSigmaCorr->SetMarkerColor(kRed);
1283 fhSigmaCorr->Draw("esame");
1284 csigma->SetLogy();
1285 TLegend * leg = new TLegend(0.7,0.75,0.87,0.5);
1286 leg->SetBorderSize(0);
1287 leg->SetLineColor(0);
1288 leg->SetFillColor(0);
1289 leg->SetTextFont(42);
1290 leg->AddEntry(gPrediction,"FONLL ","fl");
1291 leg->AddEntry(fhSigmaCorr,"data stat. unc.","pl");
1292 leg->AddEntry(fgSigmaCorr,"data syst. unc.","f");
1293 leg->Draw();
1294 csigma->Draw();
1295
1296}
1297
1298//_________________________________________________________________________________________________________
86bdcd8c 1299TH1D * AliHFPtSpectrum::ReweightHisto(TH1D *hToReweight, TH1D *hReference){
65e55bbd 1300 //
1301 // Function to reweight histograms for testing purposes:
1302 // This function takes the histo hToReweight and reweights
1303 // it (its pt shape) with respect to hReference
1304 //
1305
1306 // check histograms consistency
1307 Bool_t areconsistent=kTRUE;
1308 areconsistent &= CheckHistosConsistency(hToReweight,hReference);
1309 if (!areconsistent) {
1310 AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
1311 return NULL;
1312 }
1313
1314 // define a new empty histogram
86bdcd8c 1315 TH1D *hReweighted = (TH1D*)hToReweight->Clone("hReweighted");
65e55bbd 1316 hReweighted->Reset();
1317 Double_t weight=1.0;
1318 Double_t yvalue=1.0;
a17b17dd 1319 Double_t integralRef = hReference->Integral();
1320 Double_t integralH = hToReweight->Integral();
65e55bbd 1321
1322 // now reweight the spectra
1323 //
1324 // the weight is the relative probability of the given pt bin in the reference histo
1325 // divided by its relative probability (to normalize it) on the histo to re-weight
1326 for (Int_t i=0; i<=hToReweight->GetNbinsX(); i++) {
a17b17dd 1327 weight = (hReference->GetBinContent(i)/integralRef) / (hToReweight->GetBinContent(i)/integralH) ;
65e55bbd 1328 yvalue = hToReweight->GetBinContent(i);
1329 hReweighted->SetBinContent(i,yvalue*weight);
1330 }
1331
86bdcd8c 1332 return (TH1D*)hReweighted;
65e55bbd 1333}
1334
1335//_________________________________________________________________________________________________________
86bdcd8c 1336TH1D * AliHFPtSpectrum::ReweightRecHisto(TH1D *hRecToReweight, TH1D *hMCToReweight, TH1D *hMCReference){
65e55bbd 1337 //
1338 // Function to reweight histograms for testing purposes:
1339 // This function takes the histo hToReweight and reweights
1340 // it (its pt shape) with respect to hReference /hMCToReweight
1341 //
1342
1343 // check histograms consistency
1344 Bool_t areconsistent=kTRUE;
1345 areconsistent &= CheckHistosConsistency(hMCToReweight,hMCReference);
1346 areconsistent &= CheckHistosConsistency(hRecToReweight,hMCReference);
1347 if (!areconsistent) {
1348 AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
1349 return NULL;
1350 }
1351
1352 // define a new empty histogram
86bdcd8c 1353 TH1D *hReweighted = (TH1D*)hMCToReweight->Clone("hReweighted");
65e55bbd 1354 hReweighted->Reset();
86bdcd8c 1355 TH1D *hRecReweighted = (TH1D*)hRecToReweight->Clone("hRecReweighted");
65e55bbd 1356 hRecReweighted->Reset();
1357 Double_t weight=1.0;
1358 Double_t yvalue=1.0, yrecvalue=1.0;
a17b17dd 1359 Double_t integralRef = hMCReference->Integral();
1360 Double_t integralH = hMCToReweight->Integral();
65e55bbd 1361
1362 // now reweight the spectra
1363 //
1364 // the weight is the relative probability of the given pt bin
1365 // that should be applied in the MC histo to get the reference histo shape
1366 // Probabilities are properly normalized.
1367 for (Int_t i=0; i<=hMCToReweight->GetNbinsX(); i++) {
a17b17dd 1368 weight = (hMCReference->GetBinContent(i)/integralRef) / (hMCToReweight->GetBinContent(i)/integralH) ;
65e55bbd 1369 yvalue = hMCToReweight->GetBinContent(i);
1370 hReweighted->SetBinContent(i,yvalue*weight);
1371 yrecvalue = hRecToReweight->GetBinContent(i);
1372 hRecReweighted->SetBinContent(i,yrecvalue*weight);
1373 }
1374
86bdcd8c 1375 return (TH1D*)hRecReweighted;
65e55bbd 1376}
1377