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