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