]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG3/vertexingHF/AliHFPtSpectrum.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fixed memory leak (Renu)
[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//
26// Author: Z.Conesa, zconesa@in2p3.fr
27//***********************************************************************
28
29#include <Riostream.h>
30
31#include "TMath.h"
32#include "TH1.h"
33#include "TH1D.h"
34#include "TGraphAsymmErrors.h"
35
36#include "AliLog.h"
37#include "AliHFPtSpectrum.h"
38
39ClassImp(AliHFPtSpectrum)
40
41//_________________________________________________________________________________________________________
42AliHFPtSpectrum::AliHFPtSpectrum(const char* name, const char* title, Int_t option):
43 TNamed(name,title),
44 fhDirectMCpt(),
45 fhFeedDownMCpt(),
a17b17dd 46 fhDirectMCptMax(),
47 fhDirectMCptMin(),
48 fhFeedDownMCptMax(),
49 fhFeedDownMCptMin(),
65e55bbd 50 fhDirectEffpt(),
51 fhFeedDownEffpt(),
52 fhRECpt(),
53 fLuminosity(),
54 fTrigEfficiency(),
55 fhFc(),
a17b17dd 56 fhFcMax(),
57 fhFcMin(),
65e55bbd 58 fgFc(),
59 fhYieldCorr(),
a17b17dd 60 fhYieldCorrMax(),
61 fhYieldCorrMin(),
65e55bbd 62 fgYieldCorr(),
63 fhSigmaCorr(),
a17b17dd 64 fhSigmaCorrMax(),
65 fhSigmaCorrMin(),
65e55bbd 66 fgSigmaCorr(),
67 fFeedDownOption(option),
68 fAsymUncertainties(kTRUE)
69{
70 //
71 // Default constructor
72 //
73
74 fLuminosity[0]=1.; fLuminosity[1]=0.;
75 fTrigEfficiency[0]=1.; fTrigEfficiency[1]=0.;
76
77}
78
79//_________________________________________________________________________________________________________
80AliHFPtSpectrum::AliHFPtSpectrum(const AliHFPtSpectrum &rhs):
81 TNamed(rhs),
82 fhDirectMCpt(rhs.fhDirectMCpt),
83 fhFeedDownMCpt(rhs.fhFeedDownMCpt),
a17b17dd 84 fhDirectMCptMax(rhs.fhDirectMCptMax),
85 fhDirectMCptMin(rhs.fhDirectMCptMin),
86 fhFeedDownMCptMax(rhs.fhFeedDownMCptMax),
87 fhFeedDownMCptMin(rhs.fhFeedDownMCptMin),
65e55bbd 88 fhDirectEffpt(rhs.fhDirectEffpt),
89 fhFeedDownEffpt(rhs.fhFeedDownEffpt),
90 fhRECpt(rhs.fhRECpt),
91 fLuminosity(),
92 fTrigEfficiency(),
93 fhFc(rhs.fhFc),
a17b17dd 94 fhFcMax(rhs.fhFcMax),
95 fhFcMin(rhs.fhFcMin),
65e55bbd 96 fgFc(rhs.fgFc),
97 fhYieldCorr(rhs.fhYieldCorr),
a17b17dd 98 fhYieldCorrMax(rhs.fhYieldCorrMax),
99 fhYieldCorrMin(rhs.fhYieldCorrMin),
65e55bbd 100 fgYieldCorr(rhs.fgYieldCorr),
101 fhSigmaCorr(rhs.fhSigmaCorr),
a17b17dd 102 fhSigmaCorrMax(rhs.fhSigmaCorrMax),
103 fhSigmaCorrMin(rhs.fhSigmaCorrMin),
65e55bbd 104 fgSigmaCorr(rhs.fgSigmaCorr),
105 fFeedDownOption(rhs.fFeedDownOption),
106 fAsymUncertainties(rhs.fAsymUncertainties)
107{
108 //
109 // Copy constructor
110 //
111
112 for(Int_t i=0; i<2; i++){
113 fLuminosity[i] = rhs.fLuminosity[i];
114 fTrigEfficiency[i] = rhs.fTrigEfficiency[i];
115 }
116
117}
118
119//_________________________________________________________________________________________________________
120AliHFPtSpectrum &AliHFPtSpectrum::operator=(const AliHFPtSpectrum &source){
121 //
122 // Assignment operator
123 //
124
125 if (&source == this) return *this;
126
127 fhDirectMCpt = source.fhDirectMCpt;
128 fhFeedDownMCpt = source.fhFeedDownMCpt;
a17b17dd 129 fhDirectMCptMax = source.fhDirectMCptMax;
130 fhDirectMCptMin = source.fhDirectMCptMin;
131 fhFeedDownMCptMax = source.fhFeedDownMCptMax;
132 fhFeedDownMCptMin = source.fhFeedDownMCptMin;
65e55bbd 133 fhDirectEffpt = source.fhDirectEffpt;
134 fhFeedDownEffpt = source.fhFeedDownEffpt;
135 fhRECpt = source.fhRECpt;
136 fhFc = source.fhFc;
a17b17dd 137 fhFcMax = source.fhFcMax;
138 fhFcMin = source.fhFcMin;
65e55bbd 139 fgFc = source.fgFc;
140 fhYieldCorr = source.fhYieldCorr;
a17b17dd 141 fhYieldCorrMax = source.fhYieldCorrMax;
142 fhYieldCorrMin = source.fhYieldCorrMin;
65e55bbd 143 fgYieldCorr = source.fgYieldCorr;
144 fhSigmaCorr = source.fhSigmaCorr;
a17b17dd 145 fhSigmaCorrMax = source.fhSigmaCorrMax;
146 fhSigmaCorrMin = source.fhSigmaCorrMin;
65e55bbd 147 fgSigmaCorr = source.fgSigmaCorr;
148 fFeedDownOption = source.fFeedDownOption;
149 fAsymUncertainties = source.fAsymUncertainties;
150
151 for(Int_t i=0; i<2; i++){
152 fLuminosity[i] = source.fLuminosity[i];
153 fTrigEfficiency[i] = source.fTrigEfficiency[i];
154 }
155
156 return *this;
157}
158
159//_________________________________________________________________________________________________________
160AliHFPtSpectrum::~AliHFPtSpectrum(){
161 //
162 // Destructor
163 //
164 ;
165}
166
167
168//_________________________________________________________________________________________________________
169void AliHFPtSpectrum::SetMCptSpectra(TH1 *hDirect, TH1 *hFeedDown){
170 //
171 // Set the MonteCarlo or Theoretical spectra
172 // both for direct and feed-down contributions
173 //
174
175 if (!hDirect || !hFeedDown) {
176 AliError("One or both (direct, feed-down) spectra don't exist");
177 return;
178 }
179
180 Bool_t areconsistent = kTRUE;
181 areconsistent = CheckHistosConsistency(hDirect,hFeedDown);
182 if (!areconsistent) {
183 AliInfo("Histograms are not consistent (bin width, bounds)");
184 return;
185 }
186
187 fhDirectMCpt = hDirect;
188 fhFeedDownMCpt = hFeedDown;
189}
190
191//_________________________________________________________________________________________________________
192void AliHFPtSpectrum::SetFeedDownMCptSpectra(TH1 *hFeedDown){
193 //
194 // Set the MonteCarlo or Theoretical spectra
195 // for feed-down contribution
196 //
197
198 if (!hFeedDown) {
199 AliError("Feed-down spectra don't exist");
200 return;
201 }
202 fhFeedDownMCpt = hFeedDown;
203}
204
205//_________________________________________________________________________________________________________
206void AliHFPtSpectrum::SetMCptDistributionsBounds(TH1 *hDirectMax, TH1 *hDirectMin, TH1 *hFeedDownMax, TH1 *hFeedDownMin){
207 //
208 // Set the maximum and minimum MonteCarlo or Theoretical spectra
209 // both for direct and feed-down contributions
210 // used in case uncertainties are asymmetric and ca not be on the "basic histograms"
211 //
212
213 if (!hDirectMax || !hDirectMin || !hFeedDownMax|| !hFeedDownMin) {
214 AliError("One or all of the max/min direct/feed-down spectra don't exist");
215 return;
216 }
217
218 Bool_t areconsistent = kTRUE;
219 areconsistent &= CheckHistosConsistency(hDirectMax,hDirectMin);
220 areconsistent &= CheckHistosConsistency(hFeedDownMax,hFeedDownMin);
221 areconsistent &= CheckHistosConsistency(hDirectMax,hFeedDownMax);
222 if (!areconsistent) {
223 AliInfo("Histograms are not consistent (bin width, bounds)");
224 return;
225 }
226
a17b17dd 227 fhDirectMCptMax = hDirectMax;
228 fhDirectMCptMin = hDirectMin;
229 fhFeedDownMCptMax = hFeedDownMax;
230 fhFeedDownMCptMin = hFeedDownMin;
65e55bbd 231}
232
233//_________________________________________________________________________________________________________
234void AliHFPtSpectrum::SetFeedDownMCptDistributionsBounds(TH1 *hFeedDownMax, TH1 *hFeedDownMin){
235 //
236 // Set the maximum and minimum MonteCarlo or Theoretical spectra
237 // for feed-down contributions
238 // used in case uncertainties are asymmetric and can not be on the "basic histogram"
239 //
240
241 if (!hFeedDownMax || !hFeedDownMin) {
242 AliError("One or all of the max/min direct/feed-down spectra don't exist");
243 return;
244 }
245
246 Bool_t areconsistent = kTRUE;
247 areconsistent &= CheckHistosConsistency(hFeedDownMax,hFeedDownMin);
248 if (!areconsistent) {
249 AliInfo("Histograms are not consistent (bin width, bounds)");
250 return;
251 }
252
a17b17dd 253 fhFeedDownMCptMax = hFeedDownMax;
254 fhFeedDownMCptMin = hFeedDownMin;
65e55bbd 255}
256
257//_________________________________________________________________________________________________________
258void AliHFPtSpectrum::SetDirectAccEffCorrection(TH1 *hDirectEff){
259 //
260 // Set the Acceptance and Efficiency corrections
261 // for the direct contribution
262 //
263
264 if (!hDirectEff) {
265 AliError("The direct acceptance and efficiency corrections doesn't exist");
266 return;
267 }
268
269 fhDirectEffpt = hDirectEff;
270}
271
272//_________________________________________________________________________________________________________
273void AliHFPtSpectrum::SetAccEffCorrection(TH1 *hDirectEff, TH1 *hFeedDownEff){
274 //
275 // Set the Acceptance and Efficiency corrections
276 // both for direct and feed-down contributions
277 //
278
279 if (!hDirectEff || !hFeedDownEff) {
280 AliError("One or both (direct, feed-down) acceptance and efficiency corrections don't exist");
281 return;
282 }
283
284 Bool_t areconsistent=kTRUE;
285 areconsistent = CheckHistosConsistency(hDirectEff,hFeedDownEff);
286 if (!areconsistent) {
287 AliInfo("Histograms are not consistent (bin width, bounds)");
288 return;
289 }
290
291 fhDirectEffpt = hDirectEff;
292 fhFeedDownEffpt = hFeedDownEff;
293}
294
295//_________________________________________________________________________________________________________
296void AliHFPtSpectrum::SetReconstructedSpectrum(TH1 *hRec) {
297 //
298 // Set the reconstructed spectrum
299 //
300
301 if (!hRec) {
302 AliError("The reconstructed spectrum doesn't exist");
303 return;
304 }
305
306 fhRECpt = hRec;
307}
308
309//_________________________________________________________________________________________________________
a17b17dd 310void AliHFPtSpectrum::ComputeHFPtSpectrum(Double_t deltaY, Double_t branchingRatioC, Double_t branchingRatioBintoFinalDecay) {
65e55bbd 311 //
312 // Main function to compute the corrected cross-section:
a17b17dd 313 // variables : analysed delta_y, BR for the final correction,
314 // BR b --> D --> decay (relative to the input theoretical prediction)
65e55bbd 315 //
316 // Sigma = ( 1. / (lumi * delta_y * BR_c * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
317 //
318 // Uncertainties: delta_sigma = sigma * sqrt ( (delta_reco/reco)^2 + (delta_lumi/lumi)^2 + (delta_eff_trig/eff_trig)^2 )
319
320 //
321 // First: Initialization
322 //
323 Bool_t areHistosOk = Initialize();
324 if (!areHistosOk) {
325 AliInfo(" Histos not properly initialized. Check : inconsistent binning ? missing histos ?");
326 return;
327 }
328
329 //
330 // Second: Correct for feed-down
331 //
332 if (fFeedDownOption==1) {
333 // Compute the feed-down correction via fc-method
a17b17dd 334 CalculateFeedDownCorrectionFc();
65e55bbd 335 // Correct the yield for feed-down correction via fc-method
a17b17dd 336 CalculateFeedDownCorrectedSpectrumFc();
65e55bbd 337 }
338 else if (fFeedDownOption==2) {
339 // Correct the yield for feed-down correction via Nb-method
a17b17dd 340 CalculateFeedDownCorrectedSpectrumNb(deltaY,branchingRatioBintoFinalDecay);
65e55bbd 341 }
342 else if (fFeedDownOption==0) {
343 // If there is no need for feed-down correction,
344 // the "corrected" yield is equal to the raw yield
345 fhYieldCorr = fhRECpt;
346 fhYieldCorr->SetNameTitle("fhYieldCorr","un-corrected yield");
a17b17dd 347 fhYieldCorrMax = fhRECpt;
348 fhYieldCorrMin = fhRECpt;
349 fhYieldCorrMax->SetNameTitle("fhYieldCorrMax","un-corrected yield");
350 fhYieldCorrMin->SetNameTitle("fhYieldCorrMin","un-corrected yield");
65e55bbd 351 fAsymUncertainties=kFALSE;
352 }
353 else {
354 AliInfo(" Are you sure the feed-down correction option is right ?");
355 }
356
357 // Print out information
a17b17dd 358 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);
65e55bbd 359
360 //
361 // Finally: Correct from yields to cross-section
362 //
363 Int_t nbins = fhRECpt->GetNbinsX();
364 Double_t binwidth = fhRECpt->GetBinWidth(1);
365 Double_t xmin = fhRECpt->GetBinCenter(1) - (binwidth/2.) ;
366 Double_t xmax = fhRECpt->GetBinCenter(nbins) + (binwidth/2.) ;
367
368 // declare the output histograms
369 TH1D *hSigmaCorr = new TH1D("hSigmaCorr","corrected sigma",nbins,xmin,xmax);
a17b17dd 370 TH1D *hSigmaCorrMax = new TH1D("hSigmaCorrMax","max corrected sigma",nbins,xmin,xmax);
371 TH1D *hSigmaCorrMin = new TH1D("hSigmaCorrMin","min corrected sigma",nbins,xmin,xmax);
65e55bbd 372 // and the output TGraphAsymmErrors
373 if (fAsymUncertainties & !fgSigmaCorr) fgSigmaCorr = new TGraphAsymmErrors(nbins);
374
375 // protect against null denominator
a17b17dd 376 if (deltaY==0. || fLuminosity[0]==0. || fTrigEfficiency[0]==0. || branchingRatioC==0.) {
65e55bbd 377 AliError(" Hey you ! Why luminosity or trigger-efficiency or the c-BR or delta_y are set to zero ?! ");
378 return ;
379 }
380
a17b17dd 381 Double_t value=0, valueMax=0., valueMin=0.;
65e55bbd 382 for(Int_t ibin=0; ibin<=nbins; ibin++){
383
384 // Sigma = ( 1. / (lumi * delta_y * BR_c * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
385 value = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ?
a17b17dd 386 ( fhYieldCorr->GetBinContent(ibin) / ( deltaY * branchingRatioC * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
65e55bbd 387 : 0. ;
388
389 // Uncertainties: delta_sigma = sigma * sqrt ( (delta_reco/reco)^2 + (delta_lumi/lumi)^2 + (delta_eff_trig/eff_trig)^2 )
390 if (fAsymUncertainties) {
a17b17dd 391 valueMax = value * TMath::Sqrt( (fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin))* (fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin)) +
65e55bbd 392 (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
393 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) );
a17b17dd 394 valueMin = value * TMath::Sqrt( (fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin))* (fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin)) +
65e55bbd 395 (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
396 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) );
397 }
398 else {
399 // protect against null denominator
a17b17dd 400 valueMax = (value!=0.) ?
65e55bbd 401 value * TMath::Sqrt( (fhYieldCorr->GetBinError(ibin)/fhYieldCorr->GetBinContent(ibin))* (fhYieldCorr->GetBinError(ibin)/fhYieldCorr->GetBinContent(ibin)) +
402 (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
403 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) )
404 : 0. ;
a17b17dd 405 valueMin = valueMax;
65e55bbd 406 }
407
408 // Fill the histograms
409 hSigmaCorr->SetBinContent(ibin,value);
a17b17dd 410 hSigmaCorr->SetBinError(ibin,valueMax);
411 hSigmaCorrMax->SetBinContent(ibin,valueMax);
412 hSigmaCorrMin->SetBinContent(ibin,valueMin);
65e55bbd 413 // Fill the TGraphAsymmErrors
414 if (fAsymUncertainties) {
415 Double_t x = fhYieldCorr->GetBinCenter(ibin);
416 fgSigmaCorr->SetPoint(ibin,x,value); // i,x,y
a17b17dd 417 fgSigmaCorr->SetPointError(ibin,(binwidth/2.),(binwidth/2.),valueMin,valueMax); // i,xl,xh,yl,yh
65e55bbd 418 }
419
420 }
421
422 fhSigmaCorr = hSigmaCorr ;
a17b17dd 423 fhSigmaCorrMax = hSigmaCorrMax ;
424 fhSigmaCorrMin = hSigmaCorrMin ;
65e55bbd 425}
426
427//_________________________________________________________________________________________________________
428Bool_t AliHFPtSpectrum::Initialize(){
429 //
430 // Initialization of the variables (histograms)
431 //
432
433 if (fFeedDownOption==0) {
434 AliInfo("Getting ready for the corrections without feed-down consideration");
435 } else if (fFeedDownOption==1) {
436 AliInfo("Getting ready for the fc feed-down correction calculation");
437 } else if (fFeedDownOption==2) {
438 AliInfo("Getting ready for the Nb feed-down correction calculation");
439 } else { AliError("The calculation option must be <=2"); return kFALSE; }
440
441 // Start checking the input histograms consistency
442 Bool_t areconsistent=kTRUE;
443
444 // General checks
445 if (!fhDirectEffpt || !fhRECpt) {
446 AliError(" Reconstructed spectra and/or the Nc efficiency distributions are not defined");
447 return kFALSE;
448 }
449 areconsistent &= CheckHistosConsistency(fhRECpt,fhDirectEffpt);
450 if (!areconsistent) {
451 AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
452 return kFALSE;
453 }
454 if (fFeedDownOption==0) return kTRUE;
455
456 //
457 // Common checks for options 1 (fc) & 2(Nb)
458 if (!fhFeedDownMCpt || !fhFeedDownEffpt) {
459 AliError(" Theoretical Nb and/or the Nb efficiency distributions are not defined");
460 return kFALSE;
461 }
462 areconsistent &= CheckHistosConsistency(fhRECpt,fhFeedDownMCpt);
463 areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownEffpt);
464 if (fAsymUncertainties) {
a17b17dd 465 if (!fhFeedDownMCptMax || !fhFeedDownMCptMin) {
65e55bbd 466 AliError(" Max/Min theoretical Nb distributions are not defined");
467 return kFALSE;
468 }
a17b17dd 469 areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownMCptMax);
65e55bbd 470 }
471 if (!areconsistent) {
472 AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
473 return kFALSE;
474 }
475 if (fFeedDownOption>1) return kTRUE;
476
477 //
478 // Now checks for option 1 (fc correction)
479 if (!fhDirectMCpt) {
480 AliError("Theoretical Nc distributions is not defined");
481 return kFALSE;
482 }
483 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhFeedDownMCpt);
484 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectEffpt);
485 if (fAsymUncertainties) {
a17b17dd 486 if (!fhDirectMCptMax || !fhDirectMCptMin) {
65e55bbd 487 AliError(" Max/Min theoretical Nc distributions are not defined");
488 return kFALSE;
489 }
a17b17dd 490 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectMCptMax);
65e55bbd 491 }
492 if (!areconsistent) {
493 AliInfo("Histograms required for fc correction are not consistent (bin width, bounds)");
494 return kFALSE;
495 }
496
497 return kTRUE;
498}
499
500//_________________________________________________________________________________________________________
501Bool_t AliHFPtSpectrum::CheckHistosConsistency(TH1 *h1, TH1 *h2){
502 //
503 // Check the histograms consistency (bins, limits)
504 //
505
506 if (!h1 || !h2) {
507 AliError("One or both histograms don't exist");
508 return kFALSE;
509 }
510
511 Double_t binwidth1 = h1->GetBinWidth(1);
512 Double_t binwidth2 = h2->GetBinWidth(1);
513 Double_t min1 = h1->GetBinCenter(1) - (binwidth1/2.) ;
514// Double_t max1 = h1->GetBinCenter(nbins1) + (binwidth1/2.) ;
515 Double_t min2 = h2->GetBinCenter(1) - (binwidth2/2.) ;
516// Double_t max2 = h2->GetBinCenter(nbins2) + (binwidth2/2.) ;
517
518 if (binwidth1!=binwidth2) {
519 AliInfo(" histograms with different bin width");
520 return kFALSE;
521 }
522 if (min1!=min2) {
523 AliInfo(" histograms with different minimum");
524 return kFALSE;
525 }
526// if (max1!=max2) {
527// AliInfo(" histograms with different maximum");
528// return kFALSE;
529// }
530
531 return kTRUE;
532}
533
534//_________________________________________________________________________________________________________
a17b17dd 535void AliHFPtSpectrum::CalculateFeedDownCorrectionFc(){
65e55bbd 536 //
537 // Compute fc factor and its uncertainties bin by bin
538 // fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
539 //
540
541 // define the variables
542 Int_t nbins = fhRECpt->GetNbinsX();
543 Double_t binwidth = fhRECpt->GetBinWidth(1);
544 Double_t xmin = fhRECpt->GetBinCenter(1) - (binwidth/2.) ;
545 Double_t xmax = fhRECpt->GetBinCenter(nbins) + (binwidth/2.) ;
546 Double_t correction=1.;
a17b17dd 547 Double_t correctionMax=1., correctionMin=1.;
548 Double_t theoryRatio=1.;
549 Double_t effRatio=1.;
65e55bbd 550
551 // declare the output histograms
552 TH1D *hfc = new TH1D("hfc","fc correction factor",nbins,xmin,xmax);
a17b17dd 553 TH1D *hfcMax = new TH1D("hfcMax","max fc correction factor",nbins,xmin,xmax);
554 TH1D *hfcMin = new TH1D("hfcMin","min fc correction factor",nbins,xmin,xmax);
65e55bbd 555 // two local control histograms
556 TH1D *hTheoryRatio = new TH1D("hTheoryRatio","Theoretical B-->D over c-->D (feed-down/direct) ratio",nbins,xmin,xmax);
557 TH1D *hEffRatio = new TH1D("hEffRatio","Efficiency B-->D over c-->D (feed-down/direct) ratio",nbins,xmin,xmax);
558 // and the output TGraphAsymmErrors
559 if (fAsymUncertainties & !fgFc) fgFc = new TGraphAsymmErrors(nbins);
560
561 //
562 // Compute fc
563 //
564 for (Int_t ibin=0; ibin<=nbins; ibin++) {
565
566 // theory_ratio = (N_b/N_c)
a17b17dd 567 theoryRatio = (fhDirectMCpt->GetBinContent(ibin) && fhDirectMCpt->GetBinContent(ibin)!=0.) ? fhFeedDownMCpt->GetBinContent(ibin) / fhDirectMCpt->GetBinContent(ibin) : 1.0 ;
65e55bbd 568 // eff_ratio = (eff_b/eff_c)
a17b17dd 569 effRatio = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ? fhFeedDownEffpt->GetBinContent(ibin) / fhDirectEffpt->GetBinContent(ibin) : 1.0 ;
65e55bbd 570 // fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
a17b17dd 571 correction = (effRatio && theoryRatio) ? ( 1. / ( 1 + ( effRatio * theoryRatio ) ) ) : 1.0 ;
65e55bbd 572
573 // Calculate the uncertainty [ considering only the theoretical uncertainties on Nb & Nc for now !!! ]
574 // delta_fc = fc^2 * (Nb/Nc) * sqrt ( (delta_Nb/Nb)^2 + (delta_Nc/Nc)^2 )
a17b17dd 575 Double_t deltaNbMax = fhFeedDownMCptMax->GetBinContent(ibin) - fhFeedDownMCpt->GetBinContent(ibin) ;
576 Double_t deltaNbMin = fhFeedDownMCpt->GetBinContent(ibin) - fhFeedDownMCptMin->GetBinContent(ibin) ;
577 Double_t deltaNcMax = fhDirectMCptMax->GetBinContent(ibin) - fhDirectMCpt->GetBinContent(ibin) ;
578 Double_t deltaNcMin = fhDirectMCpt->GetBinContent(ibin) - fhDirectMCptMin->GetBinContent(ibin) ;
65e55bbd 579
580 // Protect against null denominator. If so, define uncertainty as null
581 if (fhFeedDownMCpt->GetBinContent(ibin) && fhFeedDownMCpt->GetBinContent(ibin)!=0. &&
582 fhDirectMCpt->GetBinContent(ibin) && fhDirectMCpt->GetBinContent(ibin)!=0. ) {
a17b17dd 583 correctionMax = correction*correction*theoryRatio *
65e55bbd 584 TMath::Sqrt(
a17b17dd 585 (deltaNbMax/fhFeedDownMCpt->GetBinContent(ibin))*(deltaNbMax/fhFeedDownMCpt->GetBinContent(ibin)) +
586 (deltaNcMax/fhDirectMCpt->GetBinContent(ibin))*(deltaNcMax/fhDirectMCpt->GetBinContent(ibin))
65e55bbd 587 );
a17b17dd 588 correctionMin = correction*correction*theoryRatio *
65e55bbd 589 TMath::Sqrt(
a17b17dd 590 (deltaNbMin/fhFeedDownMCpt->GetBinContent(ibin))*(deltaNbMin/fhFeedDownMCpt->GetBinContent(ibin)) +
591 (deltaNcMin/fhDirectMCpt->GetBinContent(ibin))*(deltaNcMin/fhDirectMCpt->GetBinContent(ibin))
65e55bbd 592 );
593 }
a17b17dd 594 else { correctionMax = 0.; correctionMin = 0.; }
65e55bbd 595
596
597 // Fill in the histograms
a17b17dd 598 hTheoryRatio->SetBinContent(ibin,theoryRatio);
599 hEffRatio->SetBinContent(ibin,effRatio);
65e55bbd 600 hfc->SetBinContent(ibin,correction);
a17b17dd 601 hfcMax->SetBinContent(ibin,correction+correctionMax);
602 hfcMin->SetBinContent(ibin,correction-correctionMin);
65e55bbd 603 if (fAsymUncertainties) {
604 Double_t x = fhDirectMCpt->GetBinCenter(ibin);
605 fgFc->SetPoint(ibin,x,correction); // i,x,y
a17b17dd 606 fgFc->SetPointError(ibin,(binwidth/2.),(binwidth/2.),correctionMin,correctionMax); // i,xl,xh,yl,yh
65e55bbd 607 }
608
609 }
610
611 fhFc = hfc;
a17b17dd 612 fhFcMax = hfcMax;
613 fhFcMin = hfcMin;
65e55bbd 614}
615
616//_________________________________________________________________________________________________________
a17b17dd 617void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumFc(){
65e55bbd 618 //
619 // Compute the feed-down corrected spectrum if feed-down correction is done via fc factor (bin by bin)
620 // physics = reco * fc
621 //
622 // uncertainty: delta_physics = physics * sqrt ( (delta_reco/reco)^2 + (delta_fc/fc)^2 )
623 //
624 // ( Calculation done bin by bin )
625
626 if (!fhFc || !fhRECpt) {
627 AliError(" Reconstructed or fc distributions are not defined");
628 return;
629 }
630
631 Int_t nbins = fhRECpt->GetNbinsX();
a17b17dd 632 Double_t value = 0., valueDmax= 0., valueDmin= 0.;
65e55bbd 633 Double_t binwidth = fhRECpt->GetBinWidth(1);
634 Double_t xmin = fhRECpt->GetBinCenter(1) - (binwidth/2.) ;
635 Double_t xmax = fhRECpt->GetBinCenter(nbins) + (binwidth/2.) ;
636
637 // declare the output histograms
638 TH1D *hYield = new TH1D("hYield","corrected yield (by fc)",nbins,xmin,xmax);
a17b17dd 639 TH1D *hYieldMax = new TH1D("hYieldMax","max corrected yield (by fc)",nbins,xmin,xmax);
640 TH1D *hYieldMin = new TH1D("hYieldMin","min corrected yield (by fc)",nbins,xmin,xmax);
65e55bbd 641 // and the output TGraphAsymmErrors
642 if (fAsymUncertainties & !fgYieldCorr) fgYieldCorr = new TGraphAsymmErrors(nbins);
643
644 //
645 // Do the calculation
646 //
647 for (Int_t ibin=0; ibin<=nbins; ibin++) {
648
649 // calculate the value
650 value = fhRECpt->GetBinContent(ibin) * fhFc->GetBinContent(ibin) ;
651
652 // calculate the value uncertainty
653 // Protect against null denominator. If so, define uncertainty as null
654 if (fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) {
655
656 if (fAsymUncertainties) {
657
658 if (fhFc->GetBinContent(ibin) && fhFc->GetBinContent(ibin)!=0.) {
a17b17dd 659 valueDmax = value * TMath::Sqrt( ( (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin))*(fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) ) + ( (fgFc->GetErrorYhigh(ibin)/fhFc->GetBinContent(ibin))*(fgFc->GetErrorYhigh(ibin)/fhFc->GetBinContent(ibin)) ) );
660 valueDmin = value * TMath::Sqrt( ( (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin))*(fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) ) + ( (fgFc->GetErrorYlow(ibin)/fhFc->GetBinContent(ibin))*(fgFc->GetErrorYlow(ibin)/fhFc->GetBinContent(ibin)) ) );
65e55bbd 661 }
a17b17dd 662 else { valueDmax = 0.; valueDmin = 0.; }
65e55bbd 663
664 }
665 else { // Don't consider fc uncertainty in this case [ to be tested!!! ]
a17b17dd 666 valueDmax = value * (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) ;
667 valueDmin = value * (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) ;
65e55bbd 668 }
669
670 }
a17b17dd 671 else { valueDmax = 0.; valueDmin = 0.; }
65e55bbd 672
673 // fill in the histograms
674 hYield->SetBinContent(ibin,value);
a17b17dd 675 hYield->SetBinError(ibin,valueDmax);
676 hYieldMax->SetBinContent(ibin,value+valueDmax);
677 hYieldMin->SetBinContent(ibin,value-valueDmin);
65e55bbd 678 if (fAsymUncertainties) {
679 Double_t center = hYield->GetBinCenter(ibin);
680 fgYieldCorr->SetPoint(ibin,center,value); // i,x,y
a17b17dd 681 fgYieldCorr->SetPointError(ibin,(binwidth/2.),(binwidth/2.),valueDmin,valueDmax); // i,xl,xh,yl,yh
65e55bbd 682 }
683
684 }
685
686 fhYieldCorr = hYield;
a17b17dd 687 fhYieldCorrMax = hYieldMax;
688 fhYieldCorrMin = hYieldMin;
65e55bbd 689}
690
691//_________________________________________________________________________________________________________
a17b17dd 692void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrumNb(Double_t deltaY, Double_t branchingRatioBintoFinalDecay) {
65e55bbd 693 //
694 // Compute the feed-down corrected spectrum if feed-down correction is done via Nb (bin by bin)
695 // physics = reco - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th)
696 //
697 // uncertainty: delta_physics = sqrt ( (delta_reco)^2 + (k*delta_lumi/lumi)^2 +
698 // (k*delta_eff_trig/eff_trig)^2 + (k*delta_Nb/Nb)^2 )
699 // where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th
700 //
701
702 Int_t nbins = fhRECpt->GetNbinsX();
703 Double_t binwidth = fhRECpt->GetBinWidth(1);
a17b17dd 704 Double_t value = 0., valueDmax = 0., valueDmin = 0., kfactor = 0.;
65e55bbd 705 Double_t xmin = fhRECpt->GetBinCenter(1) - (binwidth/2.) ;
706 Double_t xmax = fhRECpt->GetBinCenter(nbins) + (binwidth/2.) ;
707
708 // declare the output histograms
709 TH1D *hYield = new TH1D("hYield","corrected yield (by Nb)",nbins,xmin,xmax);
a17b17dd 710 TH1D *hYieldMax = new TH1D("hYieldMax","max corrected yield (by Nb)",nbins,xmin,xmax);
711 TH1D *hYieldMin = new TH1D("hYieldMin","min corrected yield (by Nb)",nbins,xmin,xmax);
65e55bbd 712 // and the output TGraphAsymmErrors
713 if (fAsymUncertainties & !fgYieldCorr) fgYieldCorr = new TGraphAsymmErrors(nbins);
714
715 //
716 // Do the calculation
717 //
718 for (Int_t ibin=0; ibin<=nbins; ibin++) {
719
720 // calculate the value
a17b17dd 721 value = fhRECpt->GetBinContent(ibin) - (deltaY*branchingRatioBintoFinalDecay*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) );
65e55bbd 722
a17b17dd 723 kfactor = deltaY*branchingRatioBintoFinalDecay*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) ;
65e55bbd 724
725 // calculate the value uncertainty
726 if (fAsymUncertainties) {
727 Double_t Nb = fhFeedDownMCpt->GetBinContent(ibin);
a17b17dd 728 Double_t NbDmax = fhFeedDownMCptMax->GetBinContent(ibin) - fhFeedDownMCpt->GetBinContent(ibin);
729 Double_t NbDmin = fhFeedDownMCpt->GetBinContent(ibin) - fhFeedDownMCptMin->GetBinContent(ibin);
730 valueDmax = TMath::Sqrt( ( fhRECpt->GetBinError(ibin)*fhRECpt->GetBinError(ibin) ) +
731 ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
732 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
733 ( (kfactor*NbDmax/Nb)*(kfactor*NbDmax/Nb) ) );
734 valueDmin = TMath::Sqrt( ( fhRECpt->GetBinError(ibin)*fhRECpt->GetBinError(ibin) ) +
65e55bbd 735 ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
736 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
a17b17dd 737 ( (kfactor*NbDmin/Nb)*(kfactor*NbDmin/Nb) ) );
65e55bbd 738 }
739 else{ // Don't consider Nb uncertainty in this case [ to be tested!!! ]
a17b17dd 740 valueDmax = TMath::Sqrt( ( fhRECpt->GetBinError(ibin)*fhRECpt->GetBinError(ibin) ) +
65e55bbd 741 ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
742 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) );
a17b17dd 743 valueDmin = valueDmax ;
65e55bbd 744 }
745
746 // fill in histograms
747 hYield->SetBinContent(ibin,value);
a17b17dd 748 hYield->SetBinError(ibin,valueDmax);
749 hYieldMax->SetBinContent(ibin,value+valueDmax);
750 hYieldMin->SetBinContent(ibin,value-valueDmin);
65e55bbd 751 if (fAsymUncertainties) {
752 Double_t x = hYield->GetBinCenter(ibin);
753 fgYieldCorr->SetPoint(ibin,x,value); // i,x,y
a17b17dd 754 fgYieldCorr->SetPointError(ibin,(binwidth/2.),(binwidth/2.),valueDmin,valueDmax); // i,xl,xh,yl,yh
65e55bbd 755 }
756
757 }
758
759 fhYieldCorr = hYield;
a17b17dd 760 fhYieldCorrMax = hYieldMax;
761 fhYieldCorrMin = hYieldMin;
65e55bbd 762}
763
764
765//_________________________________________________________________________________________________________
766TH1 * AliHFPtSpectrum::ReweightHisto(TH1 *hToReweight, TH1 *hReference){
767 //
768 // Function to reweight histograms for testing purposes:
769 // This function takes the histo hToReweight and reweights
770 // it (its pt shape) with respect to hReference
771 //
772
773 // check histograms consistency
774 Bool_t areconsistent=kTRUE;
775 areconsistent &= CheckHistosConsistency(hToReweight,hReference);
776 if (!areconsistent) {
777 AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
778 return NULL;
779 }
780
781 // define a new empty histogram
782 TH1 *hReweighted = (TH1*)hToReweight->Clone("hReweighted");
783 hReweighted->Reset();
784 Double_t weight=1.0;
785 Double_t yvalue=1.0;
a17b17dd 786 Double_t integralRef = hReference->Integral();
787 Double_t integralH = hToReweight->Integral();
65e55bbd 788
789 // now reweight the spectra
790 //
791 // the weight is the relative probability of the given pt bin in the reference histo
792 // divided by its relative probability (to normalize it) on the histo to re-weight
793 for (Int_t i=0; i<=hToReweight->GetNbinsX(); i++) {
a17b17dd 794 weight = (hReference->GetBinContent(i)/integralRef) / (hToReweight->GetBinContent(i)/integralH) ;
65e55bbd 795 yvalue = hToReweight->GetBinContent(i);
796 hReweighted->SetBinContent(i,yvalue*weight);
797 }
798
799 return (TH1*)hReweighted;
800}
801
802//_________________________________________________________________________________________________________
803TH1 * AliHFPtSpectrum::ReweightRecHisto(TH1 *hRecToReweight, TH1 *hMCToReweight, TH1 *hMCReference){
804 //
805 // Function to reweight histograms for testing purposes:
806 // This function takes the histo hToReweight and reweights
807 // it (its pt shape) with respect to hReference /hMCToReweight
808 //
809
810 // check histograms consistency
811 Bool_t areconsistent=kTRUE;
812 areconsistent &= CheckHistosConsistency(hMCToReweight,hMCReference);
813 areconsistent &= CheckHistosConsistency(hRecToReweight,hMCReference);
814 if (!areconsistent) {
815 AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
816 return NULL;
817 }
818
819 // define a new empty histogram
820 TH1 *hReweighted = (TH1*)hMCToReweight->Clone("hReweighted");
821 hReweighted->Reset();
822 TH1 *hRecReweighted = (TH1*)hRecToReweight->Clone("hRecReweighted");
823 hRecReweighted->Reset();
824 Double_t weight=1.0;
825 Double_t yvalue=1.0, yrecvalue=1.0;
a17b17dd 826 Double_t integralRef = hMCReference->Integral();
827 Double_t integralH = hMCToReweight->Integral();
65e55bbd 828
829 // now reweight the spectra
830 //
831 // the weight is the relative probability of the given pt bin
832 // that should be applied in the MC histo to get the reference histo shape
833 // Probabilities are properly normalized.
834 for (Int_t i=0; i<=hMCToReweight->GetNbinsX(); i++) {
a17b17dd 835 weight = (hMCReference->GetBinContent(i)/integralRef) / (hMCToReweight->GetBinContent(i)/integralH) ;
65e55bbd 836 yvalue = hMCToReweight->GetBinContent(i);
837 hReweighted->SetBinContent(i,yvalue*weight);
838 yrecvalue = hRecToReweight->GetBinContent(i);
839 hRecReweighted->SetBinContent(i,yrecvalue*weight);
840 }
841
842 return (TH1*)hRecReweighted;
843}
844