]> 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
New class to produce corrected dsigma/dpt for charm hadrons, with efficiency and...
[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(),
46 fhDirectMCpt_max(),
47 fhDirectMCpt_min(),
48 fhFeedDownMCpt_max(),
49 fhFeedDownMCpt_min(),
50 fhDirectEffpt(),
51 fhFeedDownEffpt(),
52 fhRECpt(),
53 fLuminosity(),
54 fTrigEfficiency(),
55 fhFc(),
56 fhFc_max(),
57 fhFc_min(),
58 fgFc(),
59 fhYieldCorr(),
60 fhYieldCorr_max(),
61 fhYieldCorr_min(),
62 fgYieldCorr(),
63 fhSigmaCorr(),
64 fhSigmaCorr_max(),
65 fhSigmaCorr_min(),
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),
84 fhDirectMCpt_max(rhs.fhDirectMCpt_max),
85 fhDirectMCpt_min(rhs.fhDirectMCpt_min),
86 fhFeedDownMCpt_max(rhs.fhFeedDownMCpt_max),
87 fhFeedDownMCpt_min(rhs.fhFeedDownMCpt_min),
88 fhDirectEffpt(rhs.fhDirectEffpt),
89 fhFeedDownEffpt(rhs.fhFeedDownEffpt),
90 fhRECpt(rhs.fhRECpt),
91 fLuminosity(),
92 fTrigEfficiency(),
93 fhFc(rhs.fhFc),
94 fhFc_max(rhs.fhFc_max),
95 fhFc_min(rhs.fhFc_min),
96 fgFc(rhs.fgFc),
97 fhYieldCorr(rhs.fhYieldCorr),
98 fhYieldCorr_max(rhs.fhYieldCorr_max),
99 fhYieldCorr_min(rhs.fhYieldCorr_min),
100 fgYieldCorr(rhs.fgYieldCorr),
101 fhSigmaCorr(rhs.fhSigmaCorr),
102 fhSigmaCorr_max(rhs.fhSigmaCorr_max),
103 fhSigmaCorr_min(rhs.fhSigmaCorr_min),
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;
129 fhDirectMCpt_max = source.fhDirectMCpt_max;
130 fhDirectMCpt_min = source.fhDirectMCpt_min;
131 fhFeedDownMCpt_max = source.fhFeedDownMCpt_max;
132 fhFeedDownMCpt_min = source.fhFeedDownMCpt_min;
133 fhDirectEffpt = source.fhDirectEffpt;
134 fhFeedDownEffpt = source.fhFeedDownEffpt;
135 fhRECpt = source.fhRECpt;
136 fhFc = source.fhFc;
137 fhFc_max = source.fhFc_max;
138 fhFc_min = source.fhFc_min;
139 fgFc = source.fgFc;
140 fhYieldCorr = source.fhYieldCorr;
141 fhYieldCorr_max = source.fhYieldCorr_max;
142 fhYieldCorr_min = source.fhYieldCorr_min;
143 fgYieldCorr = source.fgYieldCorr;
144 fhSigmaCorr = source.fhSigmaCorr;
145 fhSigmaCorr_max = source.fhSigmaCorr_max;
146 fhSigmaCorr_min = source.fhSigmaCorr_min;
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
227 fhDirectMCpt_max = hDirectMax;
228 fhDirectMCpt_min = hDirectMin;
229 fhFeedDownMCpt_max = hFeedDownMax;
230 fhFeedDownMCpt_min = hFeedDownMin;
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
253 fhFeedDownMCpt_max = hFeedDownMax;
254 fhFeedDownMCpt_min = hFeedDownMin;
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//_________________________________________________________________________________________________________
310void AliHFPtSpectrum::ComputeHFPtSpectrum(Double_t delta_y, Double_t BR_c, Double_t BR_b){
311 //
312 // Main function to compute the corrected cross-section:
313 //
314 // Sigma = ( 1. / (lumi * delta_y * BR_c * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
315 //
316 // Uncertainties: delta_sigma = sigma * sqrt ( (delta_reco/reco)^2 + (delta_lumi/lumi)^2 + (delta_eff_trig/eff_trig)^2 )
317
318 //
319 // First: Initialization
320 //
321 Bool_t areHistosOk = Initialize();
322 if (!areHistosOk) {
323 AliInfo(" Histos not properly initialized. Check : inconsistent binning ? missing histos ?");
324 return;
325 }
326
327 //
328 // Second: Correct for feed-down
329 //
330 if (fFeedDownOption==1) {
331 // Compute the feed-down correction via fc-method
332 CalculateFeedDownCorrection_fc();
333 // Correct the yield for feed-down correction via fc-method
334 CalculateFeedDownCorrectedSpectrum_fc();
335 }
336 else if (fFeedDownOption==2) {
337 // Correct the yield for feed-down correction via Nb-method
338 CalculateFeedDownCorrectedSpectrum_Nb(delta_y,BR_b);
339 }
340 else if (fFeedDownOption==0) {
341 // If there is no need for feed-down correction,
342 // the "corrected" yield is equal to the raw yield
343 fhYieldCorr = fhRECpt;
344 fhYieldCorr->SetNameTitle("fhYieldCorr","un-corrected yield");
345 fhYieldCorr_max = fhRECpt;
346 fhYieldCorr_min = fhRECpt;
347 fhYieldCorr_max->SetNameTitle("fhYieldCorr_max","un-corrected yield");
348 fhYieldCorr_min->SetNameTitle("fhYieldCorr_min","un-corrected yield");
349 fAsymUncertainties=kFALSE;
350 }
351 else {
352 AliInfo(" Are you sure the feed-down correction option is right ?");
353 }
354
355 // Print out information
356 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],delta_y,BR_c,BR_b);
357
358 //
359 // Finally: Correct from yields to cross-section
360 //
361 Int_t nbins = fhRECpt->GetNbinsX();
362 Double_t binwidth = fhRECpt->GetBinWidth(1);
363 Double_t xmin = fhRECpt->GetBinCenter(1) - (binwidth/2.) ;
364 Double_t xmax = fhRECpt->GetBinCenter(nbins) + (binwidth/2.) ;
365
366 // declare the output histograms
367 TH1D *hSigmaCorr = new TH1D("hSigmaCorr","corrected sigma",nbins,xmin,xmax);
368 TH1D *hSigmaCorr_max = new TH1D("hSigmaCorr_max","max corrected sigma",nbins,xmin,xmax);
369 TH1D *hSigmaCorr_min = new TH1D("hSigmaCorr_min","min corrected sigma",nbins,xmin,xmax);
370 // and the output TGraphAsymmErrors
371 if (fAsymUncertainties & !fgSigmaCorr) fgSigmaCorr = new TGraphAsymmErrors(nbins);
372
373 // protect against null denominator
374 if (delta_y==0. || fLuminosity[0]==0. || fTrigEfficiency[0]==0. || BR_c==0.) {
375 AliError(" Hey you ! Why luminosity or trigger-efficiency or the c-BR or delta_y are set to zero ?! ");
376 return ;
377 }
378
379 Double_t value=0, value_max=0., value_min=0.;
380 for(Int_t ibin=0; ibin<=nbins; ibin++){
381
382 // Sigma = ( 1. / (lumi * delta_y * BR_c * eff_trig * eff_c ) ) * spectra (corrected for feed-down)
383 value = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ?
384 ( fhYieldCorr->GetBinContent(ibin) / ( delta_y * BR_c * fLuminosity[0] * fTrigEfficiency[0] * fhDirectEffpt->GetBinContent(ibin) ) )
385 : 0. ;
386
387 // Uncertainties: delta_sigma = sigma * sqrt ( (delta_reco/reco)^2 + (delta_lumi/lumi)^2 + (delta_eff_trig/eff_trig)^2 )
388 if (fAsymUncertainties) {
389 value_max = value * TMath::Sqrt( (fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin))* (fgYieldCorr->GetErrorYhigh(ibin)/fhYieldCorr->GetBinContent(ibin)) +
390 (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
391 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) );
392 value_min = value * TMath::Sqrt( (fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin))* (fgYieldCorr->GetErrorYlow(ibin)/fhYieldCorr->GetBinContent(ibin)) +
393 (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
394 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) );
395 }
396 else {
397 // protect against null denominator
398 value_max = (value!=0.) ?
399 value * TMath::Sqrt( (fhYieldCorr->GetBinError(ibin)/fhYieldCorr->GetBinContent(ibin))* (fhYieldCorr->GetBinError(ibin)/fhYieldCorr->GetBinContent(ibin)) +
400 (fLuminosity[1]/fLuminosity[0])*(fLuminosity[1]/fLuminosity[0]) +
401 (fTrigEfficiency[1]/fTrigEfficiency[0])*(fTrigEfficiency[1]/fTrigEfficiency[0]) )
402 : 0. ;
403 value_min = value_max;
404 }
405
406 // Fill the histograms
407 hSigmaCorr->SetBinContent(ibin,value);
408 hSigmaCorr_max->SetBinContent(ibin,value_max);
409 hSigmaCorr_min->SetBinContent(ibin,value_min);
410 // Fill the TGraphAsymmErrors
411 if (fAsymUncertainties) {
412 Double_t x = fhYieldCorr->GetBinCenter(ibin);
413 fgSigmaCorr->SetPoint(ibin,x,value); // i,x,y
414 fgSigmaCorr->SetPointError(ibin,(binwidth/2.),(binwidth/2.),value_min,value_max); // i,xl,xh,yl,yh
415 }
416
417 }
418
419 fhSigmaCorr = hSigmaCorr ;
420 fhSigmaCorr_max = hSigmaCorr_max ;
421 fhSigmaCorr_min = hSigmaCorr_min ;
422}
423
424//_________________________________________________________________________________________________________
425Bool_t AliHFPtSpectrum::Initialize(){
426 //
427 // Initialization of the variables (histograms)
428 //
429
430 if (fFeedDownOption==0) {
431 AliInfo("Getting ready for the corrections without feed-down consideration");
432 } else if (fFeedDownOption==1) {
433 AliInfo("Getting ready for the fc feed-down correction calculation");
434 } else if (fFeedDownOption==2) {
435 AliInfo("Getting ready for the Nb feed-down correction calculation");
436 } else { AliError("The calculation option must be <=2"); return kFALSE; }
437
438 // Start checking the input histograms consistency
439 Bool_t areconsistent=kTRUE;
440
441 // General checks
442 if (!fhDirectEffpt || !fhRECpt) {
443 AliError(" Reconstructed spectra and/or the Nc efficiency distributions are not defined");
444 return kFALSE;
445 }
446 areconsistent &= CheckHistosConsistency(fhRECpt,fhDirectEffpt);
447 if (!areconsistent) {
448 AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
449 return kFALSE;
450 }
451 if (fFeedDownOption==0) return kTRUE;
452
453 //
454 // Common checks for options 1 (fc) & 2(Nb)
455 if (!fhFeedDownMCpt || !fhFeedDownEffpt) {
456 AliError(" Theoretical Nb and/or the Nb efficiency distributions are not defined");
457 return kFALSE;
458 }
459 areconsistent &= CheckHistosConsistency(fhRECpt,fhFeedDownMCpt);
460 areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownEffpt);
461 if (fAsymUncertainties) {
462 if (!fhFeedDownMCpt_max || !fhFeedDownMCpt_min) {
463 AliError(" Max/Min theoretical Nb distributions are not defined");
464 return kFALSE;
465 }
466 areconsistent &= CheckHistosConsistency(fhFeedDownMCpt,fhFeedDownMCpt_max);
467 }
468 if (!areconsistent) {
469 AliInfo("Histograms required for Nb correction are not consistent (bin width, bounds)");
470 return kFALSE;
471 }
472 if (fFeedDownOption>1) return kTRUE;
473
474 //
475 // Now checks for option 1 (fc correction)
476 if (!fhDirectMCpt) {
477 AliError("Theoretical Nc distributions is not defined");
478 return kFALSE;
479 }
480 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhFeedDownMCpt);
481 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectEffpt);
482 if (fAsymUncertainties) {
483 if (!fhDirectMCpt_max || !fhDirectMCpt_min) {
484 AliError(" Max/Min theoretical Nc distributions are not defined");
485 return kFALSE;
486 }
487 areconsistent &= CheckHistosConsistency(fhDirectMCpt,fhDirectMCpt_max);
488 }
489 if (!areconsistent) {
490 AliInfo("Histograms required for fc correction are not consistent (bin width, bounds)");
491 return kFALSE;
492 }
493
494 return kTRUE;
495}
496
497//_________________________________________________________________________________________________________
498Bool_t AliHFPtSpectrum::CheckHistosConsistency(TH1 *h1, TH1 *h2){
499 //
500 // Check the histograms consistency (bins, limits)
501 //
502
503 if (!h1 || !h2) {
504 AliError("One or both histograms don't exist");
505 return kFALSE;
506 }
507
508 Double_t binwidth1 = h1->GetBinWidth(1);
509 Double_t binwidth2 = h2->GetBinWidth(1);
510 Double_t min1 = h1->GetBinCenter(1) - (binwidth1/2.) ;
511// Double_t max1 = h1->GetBinCenter(nbins1) + (binwidth1/2.) ;
512 Double_t min2 = h2->GetBinCenter(1) - (binwidth2/2.) ;
513// Double_t max2 = h2->GetBinCenter(nbins2) + (binwidth2/2.) ;
514
515 if (binwidth1!=binwidth2) {
516 AliInfo(" histograms with different bin width");
517 return kFALSE;
518 }
519 if (min1!=min2) {
520 AliInfo(" histograms with different minimum");
521 return kFALSE;
522 }
523// if (max1!=max2) {
524// AliInfo(" histograms with different maximum");
525// return kFALSE;
526// }
527
528 return kTRUE;
529}
530
531//_________________________________________________________________________________________________________
532void AliHFPtSpectrum::CalculateFeedDownCorrection_fc(){
533 //
534 // Compute fc factor and its uncertainties bin by bin
535 // fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
536 //
537
538 // define the variables
539 Int_t nbins = fhRECpt->GetNbinsX();
540 Double_t binwidth = fhRECpt->GetBinWidth(1);
541 Double_t xmin = fhRECpt->GetBinCenter(1) - (binwidth/2.) ;
542 Double_t xmax = fhRECpt->GetBinCenter(nbins) + (binwidth/2.) ;
543 Double_t correction=1.;
544 Double_t correction_max=1., correction_min=1.;
545 Double_t theory_ratio=1.;
546 Double_t eff_ratio=1.;
547
548 // declare the output histograms
549 TH1D *hfc = new TH1D("hfc","fc correction factor",nbins,xmin,xmax);
550 TH1D *hfc_max = new TH1D("hfc_max","max fc correction factor",nbins,xmin,xmax);
551 TH1D *hfc_min = new TH1D("hfc_min","min fc correction factor",nbins,xmin,xmax);
552 // two local control histograms
553 TH1D *hTheoryRatio = new TH1D("hTheoryRatio","Theoretical B-->D over c-->D (feed-down/direct) ratio",nbins,xmin,xmax);
554 TH1D *hEffRatio = new TH1D("hEffRatio","Efficiency B-->D over c-->D (feed-down/direct) ratio",nbins,xmin,xmax);
555 // and the output TGraphAsymmErrors
556 if (fAsymUncertainties & !fgFc) fgFc = new TGraphAsymmErrors(nbins);
557
558 //
559 // Compute fc
560 //
561 for (Int_t ibin=0; ibin<=nbins; ibin++) {
562
563 // theory_ratio = (N_b/N_c)
564 theory_ratio = (fhDirectMCpt->GetBinContent(ibin) && fhDirectMCpt->GetBinContent(ibin)!=0.) ? fhFeedDownMCpt->GetBinContent(ibin) / fhDirectMCpt->GetBinContent(ibin) : 1.0 ;
565 // eff_ratio = (eff_b/eff_c)
566 eff_ratio = (fhDirectEffpt->GetBinContent(ibin) && fhDirectEffpt->GetBinContent(ibin)!=0.) ? fhFeedDownEffpt->GetBinContent(ibin) / fhDirectEffpt->GetBinContent(ibin) : 1.0 ;
567 // fc = 1 / ( 1 + (eff_b/eff_c)*(N_b/N_c) )
568 correction = (eff_ratio && theory_ratio) ? ( 1. / ( 1 + ( eff_ratio * theory_ratio ) ) ) : 1.0 ;
569
570 // Calculate the uncertainty [ considering only the theoretical uncertainties on Nb & Nc for now !!! ]
571 // delta_fc = fc^2 * (Nb/Nc) * sqrt ( (delta_Nb/Nb)^2 + (delta_Nc/Nc)^2 )
572 Double_t delta_Nb_max = fhFeedDownMCpt_max->GetBinContent(ibin) - fhFeedDownMCpt->GetBinContent(ibin) ;
573 Double_t delta_Nb_min = fhFeedDownMCpt->GetBinContent(ibin) - fhFeedDownMCpt_min->GetBinContent(ibin) ;
574 Double_t delta_Nc_max = fhDirectMCpt_max->GetBinContent(ibin) - fhDirectMCpt->GetBinContent(ibin) ;
575 Double_t delta_Nc_min = fhDirectMCpt->GetBinContent(ibin) - fhDirectMCpt_min->GetBinContent(ibin) ;
576
577 // Protect against null denominator. If so, define uncertainty as null
578 if (fhFeedDownMCpt->GetBinContent(ibin) && fhFeedDownMCpt->GetBinContent(ibin)!=0. &&
579 fhDirectMCpt->GetBinContent(ibin) && fhDirectMCpt->GetBinContent(ibin)!=0. ) {
580 correction_max = correction*correction*theory_ratio *
581 TMath::Sqrt(
582 (delta_Nb_max/fhFeedDownMCpt->GetBinContent(ibin))*(delta_Nb_max/fhFeedDownMCpt->GetBinContent(ibin)) +
583 (delta_Nc_max/fhDirectMCpt->GetBinContent(ibin))*(delta_Nc_max/fhDirectMCpt->GetBinContent(ibin))
584 );
585 correction_min = correction*correction*theory_ratio *
586 TMath::Sqrt(
587 (delta_Nb_min/fhFeedDownMCpt->GetBinContent(ibin))*(delta_Nb_min/fhFeedDownMCpt->GetBinContent(ibin)) +
588 (delta_Nc_min/fhDirectMCpt->GetBinContent(ibin))*(delta_Nc_min/fhDirectMCpt->GetBinContent(ibin))
589 );
590 }
591 else { correction_max = 0.; correction_min = 0.; }
592
593
594 // Fill in the histograms
595 hTheoryRatio->SetBinContent(ibin,theory_ratio);
596 hEffRatio->SetBinContent(ibin,eff_ratio);
597 hfc->SetBinContent(ibin,correction);
598 hfc_max->SetBinContent(ibin,correction+correction_max);
599 hfc_min->SetBinContent(ibin,correction-correction_min);
600 if (fAsymUncertainties) {
601 Double_t x = fhDirectMCpt->GetBinCenter(ibin);
602 fgFc->SetPoint(ibin,x,correction); // i,x,y
603 fgFc->SetPointError(ibin,(binwidth/2.),(binwidth/2.),correction_min,correction_max); // i,xl,xh,yl,yh
604 }
605
606 }
607
608 fhFc = hfc;
609 fhFc_max = hfc_max;
610 fhFc_min = hfc_min;
611}
612
613//_________________________________________________________________________________________________________
614void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrum_fc(){
615 //
616 // Compute the feed-down corrected spectrum if feed-down correction is done via fc factor (bin by bin)
617 // physics = reco * fc
618 //
619 // uncertainty: delta_physics = physics * sqrt ( (delta_reco/reco)^2 + (delta_fc/fc)^2 )
620 //
621 // ( Calculation done bin by bin )
622
623 if (!fhFc || !fhRECpt) {
624 AliError(" Reconstructed or fc distributions are not defined");
625 return;
626 }
627
628 Int_t nbins = fhRECpt->GetNbinsX();
629 Double_t value = 0., value_dmax= 0., value_dmin= 0.;
630 Double_t binwidth = fhRECpt->GetBinWidth(1);
631 Double_t xmin = fhRECpt->GetBinCenter(1) - (binwidth/2.) ;
632 Double_t xmax = fhRECpt->GetBinCenter(nbins) + (binwidth/2.) ;
633
634 // declare the output histograms
635 TH1D *hYield = new TH1D("hYield","corrected yield (by fc)",nbins,xmin,xmax);
636 TH1D *hYield_max = new TH1D("hYield_max","max corrected yield (by fc)",nbins,xmin,xmax);
637 TH1D *hYield_min = new TH1D("hYield_min","min corrected yield (by fc)",nbins,xmin,xmax);
638 // and the output TGraphAsymmErrors
639 if (fAsymUncertainties & !fgYieldCorr) fgYieldCorr = new TGraphAsymmErrors(nbins);
640
641 //
642 // Do the calculation
643 //
644 for (Int_t ibin=0; ibin<=nbins; ibin++) {
645
646 // calculate the value
647 value = fhRECpt->GetBinContent(ibin) * fhFc->GetBinContent(ibin) ;
648
649 // calculate the value uncertainty
650 // Protect against null denominator. If so, define uncertainty as null
651 if (fhRECpt->GetBinContent(ibin) && fhRECpt->GetBinContent(ibin)!=0.) {
652
653 if (fAsymUncertainties) {
654
655 if (fhFc->GetBinContent(ibin) && fhFc->GetBinContent(ibin)!=0.) {
656 value_dmax = 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)) ) );
657 value_dmin = 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)) ) );
658 }
659 else { value_dmax = 0.; value_dmin = 0.; }
660
661 }
662 else { // Don't consider fc uncertainty in this case [ to be tested!!! ]
663 value_dmax = value * (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) ;
664 value_dmin = value * (fhRECpt->GetBinError(ibin)/fhRECpt->GetBinContent(ibin)) ;
665 }
666
667 }
668 else { value_dmax = 0.; value_dmin = 0.; }
669
670 // fill in the histograms
671 hYield->SetBinContent(ibin,value);
672 hYield_max->SetBinContent(ibin,value+value_dmax);
673 hYield_min->SetBinContent(ibin,value-value_dmin);
674 if (fAsymUncertainties) {
675 Double_t center = hYield->GetBinCenter(ibin);
676 fgYieldCorr->SetPoint(ibin,center,value); // i,x,y
677 fgYieldCorr->SetPointError(ibin,(binwidth/2.),(binwidth/2.),value_dmin,value_dmax); // i,xl,xh,yl,yh
678 }
679
680 }
681
682 fhYieldCorr = hYield;
683 fhYieldCorr_max = hYield_max;
684 fhYieldCorr_min = hYield_min;
685}
686
687//_________________________________________________________________________________________________________
688void AliHFPtSpectrum::CalculateFeedDownCorrectedSpectrum_Nb(Float_t delta_y, Double_t BR_b){
689 //
690 // Compute the feed-down corrected spectrum if feed-down correction is done via Nb (bin by bin)
691 // physics = reco - (lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th)
692 //
693 // uncertainty: delta_physics = sqrt ( (delta_reco)^2 + (k*delta_lumi/lumi)^2 +
694 // (k*delta_eff_trig/eff_trig)^2 + (k*delta_Nb/Nb)^2 )
695 // where k = lumi * delta_y * BR_b * eff_trig * eff_b * Nb_th
696 //
697
698 Int_t nbins = fhRECpt->GetNbinsX();
699 Double_t binwidth = fhRECpt->GetBinWidth(1);
700 Double_t value = 0., value_dmax= 0., value_dmin= 0., kfactor=0.;
701 Double_t xmin = fhRECpt->GetBinCenter(1) - (binwidth/2.) ;
702 Double_t xmax = fhRECpt->GetBinCenter(nbins) + (binwidth/2.) ;
703
704 // declare the output histograms
705 TH1D *hYield = new TH1D("hYield","corrected yield (by Nb)",nbins,xmin,xmax);
706 TH1D *hYield_max = new TH1D("hYield_max","max corrected yield (by Nb)",nbins,xmin,xmax);
707 TH1D *hYield_min = new TH1D("hYield_min","min corrected yield (by Nb)",nbins,xmin,xmax);
708 // and the output TGraphAsymmErrors
709 if (fAsymUncertainties & !fgYieldCorr) fgYieldCorr = new TGraphAsymmErrors(nbins);
710
711 //
712 // Do the calculation
713 //
714 for (Int_t ibin=0; ibin<=nbins; ibin++) {
715
716 // calculate the value
717 value = fhRECpt->GetBinContent(ibin) - (delta_y*BR_b*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) );
718
719 kfactor = delta_y*BR_b*fLuminosity[0]*fTrigEfficiency[0]*fhFeedDownEffpt->GetBinContent(ibin)*fhFeedDownMCpt->GetBinContent(ibin) ;
720
721 // calculate the value uncertainty
722 if (fAsymUncertainties) {
723 Double_t Nb = fhFeedDownMCpt->GetBinContent(ibin);
724 Double_t Nb_dmax = fhFeedDownMCpt_max->GetBinContent(ibin) - fhFeedDownMCpt->GetBinContent(ibin);
725 Double_t Nb_dmin = fhFeedDownMCpt->GetBinContent(ibin) - fhFeedDownMCpt_min->GetBinContent(ibin);
726 value_dmax = TMath::Sqrt( ( fhRECpt->GetBinError(ibin)*fhRECpt->GetBinError(ibin) ) +
727 ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
728 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) +
729 ( (kfactor*Nb_dmax/Nb)*(kfactor*Nb_dmax/Nb) ) );
730 value_dmin = 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*Nb_dmin/Nb)*(kfactor*Nb_dmin/Nb) ) );
734 }
735 else{ // Don't consider Nb uncertainty in this case [ to be tested!!! ]
736 value_dmax = TMath::Sqrt( ( fhRECpt->GetBinError(ibin)*fhRECpt->GetBinError(ibin) ) +
737 ( (kfactor*fLuminosity[1]/fLuminosity[0])*(kfactor*fLuminosity[1]/fLuminosity[0]) ) +
738 ( (kfactor*fTrigEfficiency[1]/fTrigEfficiency[0])*(kfactor*fTrigEfficiency[1]/fTrigEfficiency[0]) ) );
739 value_dmin = value_dmax ;
740 }
741
742 // fill in histograms
743 hYield->SetBinContent(ibin,value);
744 hYield_max->SetBinContent(ibin,value+value_dmax);
745 hYield_min->SetBinContent(ibin,value-value_dmin);
746 if (fAsymUncertainties) {
747 Double_t x = hYield->GetBinCenter(ibin);
748 fgYieldCorr->SetPoint(ibin,x,value); // i,x,y
749 fgYieldCorr->SetPointError(ibin,(binwidth/2.),(binwidth/2.),value_dmin,value_dmax); // i,xl,xh,yl,yh
750 }
751
752 }
753
754 fhYieldCorr = hYield;
755 fhYieldCorr_max = hYield_max;
756 fhYieldCorr_min = hYield_min;
757}
758
759
760//_________________________________________________________________________________________________________
761TH1 * AliHFPtSpectrum::ReweightHisto(TH1 *hToReweight, TH1 *hReference){
762 //
763 // Function to reweight histograms for testing purposes:
764 // This function takes the histo hToReweight and reweights
765 // it (its pt shape) with respect to hReference
766 //
767
768 // check histograms consistency
769 Bool_t areconsistent=kTRUE;
770 areconsistent &= CheckHistosConsistency(hToReweight,hReference);
771 if (!areconsistent) {
772 AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
773 return NULL;
774 }
775
776 // define a new empty histogram
777 TH1 *hReweighted = (TH1*)hToReweight->Clone("hReweighted");
778 hReweighted->Reset();
779 Double_t weight=1.0;
780 Double_t yvalue=1.0;
781 Double_t integral_ref = hReference->Integral();
782 Double_t integral_h = hToReweight->Integral();
783
784 // now reweight the spectra
785 //
786 // the weight is the relative probability of the given pt bin in the reference histo
787 // divided by its relative probability (to normalize it) on the histo to re-weight
788 for (Int_t i=0; i<=hToReweight->GetNbinsX(); i++) {
789 weight = (hReference->GetBinContent(i)/integral_ref) / (hToReweight->GetBinContent(i)/integral_h) ;
790 yvalue = hToReweight->GetBinContent(i);
791 hReweighted->SetBinContent(i,yvalue*weight);
792 }
793
794 return (TH1*)hReweighted;
795}
796
797//_________________________________________________________________________________________________________
798TH1 * AliHFPtSpectrum::ReweightRecHisto(TH1 *hRecToReweight, TH1 *hMCToReweight, TH1 *hMCReference){
799 //
800 // Function to reweight histograms for testing purposes:
801 // This function takes the histo hToReweight and reweights
802 // it (its pt shape) with respect to hReference /hMCToReweight
803 //
804
805 // check histograms consistency
806 Bool_t areconsistent=kTRUE;
807 areconsistent &= CheckHistosConsistency(hMCToReweight,hMCReference);
808 areconsistent &= CheckHistosConsistency(hRecToReweight,hMCReference);
809 if (!areconsistent) {
810 AliInfo("the histograms to reweight are not consistent (bin width, bounds)");
811 return NULL;
812 }
813
814 // define a new empty histogram
815 TH1 *hReweighted = (TH1*)hMCToReweight->Clone("hReweighted");
816 hReweighted->Reset();
817 TH1 *hRecReweighted = (TH1*)hRecToReweight->Clone("hRecReweighted");
818 hRecReweighted->Reset();
819 Double_t weight=1.0;
820 Double_t yvalue=1.0, yrecvalue=1.0;
821 Double_t integral_ref = hMCReference->Integral();
822 Double_t integral_h = hMCToReweight->Integral();
823
824 // now reweight the spectra
825 //
826 // the weight is the relative probability of the given pt bin
827 // that should be applied in the MC histo to get the reference histo shape
828 // Probabilities are properly normalized.
829 for (Int_t i=0; i<=hMCToReweight->GetNbinsX(); i++) {
830 weight = (hMCReference->GetBinContent(i)/integral_ref) / (hMCToReweight->GetBinContent(i)/integral_h) ;
831 yvalue = hMCToReweight->GetBinContent(i);
832 hReweighted->SetBinContent(i,yvalue*weight);
833 yrecvalue = hRecToReweight->GetBinContent(i);
834 hRecReweighted->SetBinContent(i,yrecvalue*weight);
835 }
836
837 return (TH1*)hRecReweighted;
838}
839