]>
Commit | Line | Data |
---|---|---|
1f1d6be7 | 1 | /* definition of the fields in the histogram returned */ |
2 | enum EValue_t { | |
3 | kYield = 1, | |
4 | kYieldStat, | |
5 | kYieldSysHi, | |
6 | kYieldSysLo, | |
7 | kMean, | |
8 | kMeanStat, | |
9 | kMeanSysHi, | |
10 | kMeanSysLo | |
11 | }; | |
12 | ||
13 | TH1 * | |
14 | YieldMean(TH1 *hstat, TH1 *hsys, TF1 *f = NULL, Double_t min = 0., Double_t max = 10., Double_t loprecision = 0.01, Double_t hiprecision = 0.1, Option_t *opt = "0q") | |
15 | { | |
16 | /* set many iterations when fitting the data so we don't | |
17 | stop minimization with MAX_CALLS */ | |
18 | TVirtualFitter::SetMaxIterations(1000000); | |
19 | ||
20 | /* create output histo */ | |
21 | Double_t integral, mean; | |
22 | TH1 *hout = new TH1D("hout", "", 8, 0, 8); | |
23 | TH1 *hlo, *hhi; | |
24 | ||
25 | /* create histo with stat+sys errors */ | |
26 | TH1 *htot = (TH1 *)hstat->Clone("htot"); | |
27 | for (Int_t ibin = 0; ibin < htot->GetNbinsX(); ibin++) { | |
28 | htot->SetBinError(ibin + 1, TMath::Sqrt(hsys->GetBinError(ibin + 1) * hsys->GetBinError(ibin + 1) + hstat->GetBinError(ibin + 1) * hstat->GetBinError(ibin + 1))); | |
29 | } | |
30 | ||
31 | /* | |
32 | * measure the central value | |
33 | */ | |
34 | ||
35 | do Int_t fitres = htot->Fit(f, opt); | |
36 | while (fitres != 0); | |
37 | hlo = YieldMean_LowExtrapolationHisto(htot, f, min, loprecision); | |
38 | hhi = YieldMean_HighExtrapolationHisto(htot, f, max, hiprecision); | |
39 | YieldMean_IntegralMean(htot, hlo, hhi, integral, mean); | |
40 | hout->SetBinContent(kYield, integral); | |
41 | hout->SetBinContent(kMean, mean); | |
42 | ||
43 | /* | |
44 | * STATISTICS | |
45 | */ | |
46 | ||
47 | TCanvas *cCanvasStat = new TCanvas("cCanvasStat"); | |
48 | cCanvasStat->Divide(2, 1); | |
49 | ||
50 | /* | |
51 | * measure statistical error | |
52 | */ | |
53 | ||
54 | /* fit with stat error */ | |
55 | do Int_t fitres = hstat->Fit(f, opt); | |
56 | while (fitres != 0); | |
57 | hlo = YieldMean_LowExtrapolationHisto(hstat, f, min, loprecision); | |
58 | hhi = YieldMean_HighExtrapolationHisto(hstat, f, max, hiprecision); | |
59 | ||
60 | /* random generation with integration (coarse) */ | |
61 | TH1 *hIntegral_tmp = new TH1F("hIntegral_tmp", "", 1000, 0.75 * integral, 1.25 * integral); | |
62 | TH1 *hMean_tmp = new TH1F("hMean_tmp", "", 1000, 0.75 * mean, 1.25 * mean); | |
63 | for (Int_t irnd = 0; irnd < 100; irnd++) { | |
64 | /* get random histogram */ | |
65 | TH1 *hrnd = YieldMean_ReturnRandom(hstat); | |
66 | /* fit */ | |
67 | TH1 *hrndlo = YieldMean_ReturnCoherentRandom(hlo); | |
68 | TH1 *hrndhi = YieldMean_ReturnCoherentRandom(hhi); | |
69 | /* integrate */ | |
70 | YieldMean_IntegralMean(hrnd, hrndlo, hrndhi, integral, mean); | |
71 | hIntegral_tmp->Fill(integral); | |
72 | hMean_tmp->Fill(mean); | |
73 | delete hrnd; | |
74 | delete hrndlo; | |
75 | delete hrndhi; | |
76 | } | |
77 | /* random generation with integration (fine) */ | |
78 | TH1 *hIntegral = new TH1F("hIntegral", "", 100, | |
79 | hIntegral_tmp->GetMean() - 10. * hIntegral_tmp->GetRMS(), | |
80 | hIntegral_tmp->GetMean() + 10. * hIntegral_tmp->GetRMS()); | |
81 | TH1 *hMean = new TH1F("hMean", "", 100, | |
82 | hMean_tmp->GetMean() - 10. * hMean_tmp->GetRMS(), | |
83 | hMean_tmp->GetMean() + 10. * hMean_tmp->GetRMS()); | |
84 | for (Int_t irnd = 0; irnd < 1000; irnd++) { | |
85 | /* get random histogram */ | |
86 | TH1 *hrnd = YieldMean_ReturnRandom(hstat); | |
87 | /* fit */ | |
88 | TH1 *hrndlo = YieldMean_ReturnCoherentRandom(hlo); | |
89 | TH1 *hrndhi = YieldMean_ReturnCoherentRandom(hhi); | |
90 | /* integrate */ | |
91 | YieldMean_IntegralMean(hrnd, hrndlo, hrndhi, integral, mean); | |
92 | hIntegral->Fill(integral); | |
93 | hMean->Fill(mean); | |
94 | delete hrnd; | |
95 | delete hrndlo; | |
96 | delete hrndhi; | |
97 | } | |
98 | TF1 *gaus = (TF1 *)gROOT->GetFunction("gaus"); | |
99 | ||
100 | cCanvasStat->cd(1); | |
101 | hIntegral->Fit(gaus, "q"); | |
102 | integral = hout->GetBinContent(kYield) * gaus->GetParameter(2) / gaus->GetParameter(1); | |
103 | hout->SetBinContent(kYieldStat, integral); | |
104 | ||
105 | cCanvasStat->cd(2); | |
106 | hMean->Fit(gaus, "q"); | |
107 | mean = hout->GetBinContent(kMean) * gaus->GetParameter(2) / gaus->GetParameter(1); | |
108 | hout->SetBinContent(kMeanStat, mean); | |
109 | ||
110 | /* | |
111 | * SYSTEMATICS | |
112 | */ | |
113 | ||
114 | TCanvas *cCanvasSys = new TCanvas("cCanvasYieldSys"); | |
115 | cCanvasSys->Divide(2, 1); | |
116 | cCanvasSys->cd(1)->DrawFrame(min, 1.e-3, max, 1.e3); | |
117 | hsys->SetMarkerStyle(20); | |
118 | hsys->SetMarkerColor(1); | |
119 | hsys->SetMarkerSize(1); | |
120 | hsys->Draw("same"); | |
121 | cCanvasSys->cd(2)->DrawFrame(min, 1.e-3, max, 1.e3); | |
122 | hsys->Draw("same"); | |
123 | ||
124 | /* | |
125 | * systematic error high | |
126 | */ | |
127 | ||
128 | TH1 *hhigh = YieldMean_ReturnExtremeHighHisto(hsys); | |
129 | do Int_t fitres = hhigh->Fit(f, opt); | |
130 | while (fitres != 0); | |
131 | hlo = YieldMean_LowExtrapolationHisto(hhigh, f, min, loprecision); | |
132 | hhi = YieldMean_HighExtrapolationHisto(hhigh, f, max, hiprecision); | |
133 | YieldMean_IntegralMean(hhigh, hlo, hhi, integral, mean); | |
134 | integral = TMath::Abs(integral - hout->GetBinContent(kYield)); | |
135 | hout->SetBinContent(kYieldSysHi, integral); | |
136 | ||
137 | cCanvasSys->cd(1); | |
138 | f->SetLineColor(2); | |
139 | f->DrawCopy("same"); | |
140 | ||
141 | /* | |
142 | * systematic error hard | |
143 | */ | |
144 | ||
145 | TH1 *hhard = YieldMean_ReturnExtremeHardHisto(hsys); | |
146 | do Int_t fitres = hhard->Fit(f, opt); | |
147 | while (fitres != 0); | |
148 | hlo = YieldMean_LowExtrapolationHisto(hhard, f, min, loprecision); | |
149 | hhi = YieldMean_HighExtrapolationHisto(hhard, f, max, hiprecision); | |
150 | YieldMean_IntegralMean(hhard, hlo, hhi, integral, mean); | |
151 | mean = TMath::Abs(mean - hout->GetBinContent(kMean)); | |
152 | hout->SetBinContent(kMeanSysHi, mean); | |
153 | ||
154 | cCanvasSys->cd(2); | |
155 | f->SetLineColor(2); | |
156 | f->DrawCopy("same"); | |
157 | ||
158 | /* | |
159 | * systematic error low | |
160 | */ | |
161 | ||
162 | TH1 *hlow = YieldMean_ReturnExtremeLowHisto(hsys); | |
163 | do Int_t fitres = hlow->Fit(f, opt); | |
164 | while (fitres != 0); | |
165 | hlo = YieldMean_LowExtrapolationHisto(hlow, f, min, loprecision); | |
166 | hhi = YieldMean_HighExtrapolationHisto(hlow, f, max, hiprecision); | |
167 | YieldMean_IntegralMean(hlow, hlo, hhi, integral, mean); | |
168 | integral = TMath::Abs(integral - hout->GetBinContent(kYield)); | |
169 | hout->SetBinContent(kYieldSysLo, integral); | |
170 | ||
171 | cCanvasSys->cd(1); | |
172 | f->SetLineColor(4); | |
173 | f->DrawCopy("same"); | |
174 | ||
175 | /* | |
176 | * systematic error soft | |
177 | */ | |
178 | ||
179 | TH1 *hsoft = YieldMean_ReturnExtremeSoftHisto(hsys); | |
180 | do Int_t fitres = hsoft->Fit(f, opt); | |
181 | while (fitres != 0); | |
182 | hlo = YieldMean_LowExtrapolationHisto(hsoft, f, min, loprecision); | |
183 | hhi = YieldMean_HighExtrapolationHisto(hsoft, f, max, hiprecision); | |
184 | YieldMean_IntegralMean(hsoft, hlo, hhi, integral, mean); | |
185 | mean = TMath::Abs(mean - hout->GetBinContent(kMean)); | |
186 | hout->SetBinContent(kMeanSysLo, mean); | |
187 | ||
188 | cCanvasSys->cd(2); | |
189 | f->SetLineColor(4); | |
190 | f->DrawCopy("same"); | |
191 | ||
192 | return hout; | |
193 | } | |
194 | ||
195 | TH1 * | |
196 | YieldMean_LowExtrapolationHisto(TH1 *h, TF1 *f, Double_t min, Double_t binwidth = 0.01) | |
197 | { | |
198 | /* find lowest edge in histo */ | |
199 | Int_t binlo; | |
200 | Double_t lo; | |
201 | for (Int_t ibin = 1; ibin < h->GetNbinsX() + 1; ibin++) { | |
202 | if (h->GetBinContent(ibin) != 0.) { | |
203 | binlo = ibin; | |
204 | lo = h->GetBinLowEdge(ibin); | |
205 | break; | |
206 | } | |
207 | } | |
208 | ||
209 | Int_t nbins = (lo - min) / binwidth; | |
210 | TH1 *hlo = new TH1F("hlo", "", nbins, min, lo); | |
211 | ||
212 | /* integrate function in histogram bins */ | |
213 | Double_t cont, err, width; | |
214 | for (Int_t ibin = 0; ibin < hlo->GetNbinsX(); ibin++) { | |
215 | width = hlo->GetBinWidth(ibin + 1); | |
216 | cont = f->Integral(hlo->GetBinLowEdge(ibin + 1), hlo->GetBinLowEdge(ibin + 2), (Double_t *)0, 1.e-6); | |
217 | err = f->IntegralError(hlo->GetBinLowEdge(ibin + 1), hlo->GetBinLowEdge(ibin + 2), (Double_t *)0, (Double_t *)0, 1.e-6); | |
218 | hlo->SetBinContent(ibin + 1, cont / width); | |
219 | hlo->SetBinError(ibin + 1, err / width); | |
220 | } | |
221 | ||
222 | return hlo; | |
223 | } | |
224 | ||
225 | TH1 * | |
226 | YieldMean_HighExtrapolationHisto(TH1 *h, TF1 *f, Double_t max, Double_t binwidth = 0.1) | |
227 | { | |
228 | /* find highest edge in histo */ | |
229 | Int_t binhi; | |
230 | Double_t hi; | |
231 | for (Int_t ibin = h->GetNbinsX(); ibin > 0; ibin--) { | |
232 | if (h->GetBinContent(ibin) != 0.) { | |
233 | binhi = ibin + 1; | |
234 | hi = h->GetBinLowEdge(ibin + 1); | |
235 | break; | |
236 | } | |
237 | } | |
238 | ||
239 | Int_t nbins = (max - hi) / binwidth; | |
240 | TH1 *hhi = new TH1F("hhi", "", nbins, hi, max); | |
241 | ||
242 | /* integrate function in histogram bins */ | |
243 | Double_t cont, err, width; | |
244 | for (Int_t ibin = 0; ibin < hhi->GetNbinsX(); ibin++) { | |
245 | width = hhi->GetBinWidth(ibin + 1); | |
246 | cont = f->Integral(hhi->GetBinLowEdge(ibin + 1), hhi->GetBinLowEdge(ibin + 2), (Double_t *)0, 1.e-6); | |
247 | err = f->IntegralError(hhi->GetBinLowEdge(ibin + 1), hhi->GetBinLowEdge(ibin + 2), (Double_t *)0, (Double_t *)0, 1.e-6); | |
248 | hhi->SetBinContent(ibin + 1, cont / width); | |
249 | hhi->SetBinError(ibin + 1, err / width); | |
250 | } | |
251 | ||
252 | return hhi; | |
253 | } | |
254 | ||
255 | TH1 * | |
256 | YieldMean_ReturnRandom(TH1 *hin) | |
257 | { | |
258 | TH1 *hout = (TH1 *)hin->Clone("hout"); | |
259 | hout->Reset(); | |
260 | Double_t cont, err; | |
261 | for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) { | |
262 | if (hin->GetBinError(ibin + 1) <= 0.) continue; | |
263 | cont = hin->GetBinContent(ibin + 1); | |
264 | err = hin->GetBinError(ibin + 1); | |
265 | hout->SetBinContent(ibin + 1, gRandom->Gaus(cont, err)); | |
266 | hout->SetBinError(ibin + 1, err); | |
267 | } | |
268 | return hout; | |
269 | } | |
270 | ||
271 | TH1 * | |
272 | YieldMean_ReturnCoherentRandom(TH1 *hin) | |
273 | { | |
274 | TH1 *hout = (TH1 *)hin->Clone("hout"); | |
275 | hout->Reset(); | |
276 | Double_t cont, err, cohe; | |
277 | cohe = gRandom->Gaus(0., 1.); | |
278 | for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) { | |
279 | if (hin->GetBinError(ibin + 1) <= 0.) continue; | |
280 | cont = hin->GetBinContent(ibin + 1); | |
281 | err = hin->GetBinError(ibin + 1); | |
282 | hout->SetBinContent(ibin + 1, cont + cohe * err); | |
283 | hout->SetBinError(ibin + 1, err); | |
284 | } | |
285 | return hout; | |
286 | } | |
287 | ||
288 | TH1 * | |
289 | YieldMean_ReturnExtremeHighHisto(TH1 *hin) | |
290 | { | |
291 | TH1 *hout = (TH1 *)hin->Clone(Form("%s_extremehigh", hin->GetName())); | |
292 | for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) { | |
293 | if (hin->GetBinError(ibin + 1) <= 0.) continue; | |
294 | Double_t val = hin->GetBinContent(ibin + 1); | |
295 | Double_t err = hin->GetBinError(ibin + 1); | |
296 | hout->SetBinContent(ibin + 1, val + err); | |
297 | } | |
298 | return hout; | |
299 | } | |
300 | ||
301 | TH1 * | |
302 | YieldMean_ReturnExtremeLowHisto(TH1 *hin) | |
303 | { | |
304 | TH1 *hout = (TH1 *)hin->Clone(Form("%s_extremelow", hin->GetName())); | |
305 | for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) { | |
306 | if (hin->GetBinError(ibin + 1) <= 0.) continue; | |
307 | Double_t val = hin->GetBinContent(ibin + 1); | |
308 | Double_t err = hin->GetBinError(ibin + 1); | |
309 | hout->SetBinContent(ibin + 1, val - err); | |
310 | } | |
311 | return hout; | |
312 | } | |
313 | ||
314 | TH1 * | |
315 | YieldMean_ReturnExtremeSoftHisto(TH1 *hin) | |
316 | { | |
317 | return YieldMean_ReturnExtremeHisto(hin, -1.); | |
318 | } | |
319 | ||
320 | TH1 * | |
321 | YieldMean_ReturnExtremeHardHisto(TH1 *hin) | |
322 | { | |
323 | return YieldMean_ReturnExtremeHisto(hin, 1.); | |
324 | } | |
325 | ||
326 | TH1 * | |
327 | YieldMean_ReturnExtremeHisto(TH1 *hin, Float_t sign = 1.) | |
328 | { | |
329 | Double_t ptlow, pthigh; | |
330 | for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) { | |
331 | if (hin->GetBinError(ibin + 1) <= 0.) continue; | |
332 | ptlow = hin->GetBinLowEdge(ibin + 1); | |
333 | break; | |
334 | } | |
335 | for (Int_t ibin = hin->GetNbinsX(); ibin >= 0; ibin--) { | |
336 | if (hin->GetBinError(ibin + 1) <= 0.) continue; | |
337 | pthigh = hin->GetBinLowEdge(ibin + 2); | |
338 | break; | |
339 | } | |
340 | ||
341 | Double_t mean = hin->GetMean(); | |
342 | Double_t maxdiff = 0.; | |
343 | TH1 *hmax = NULL; | |
344 | for (Int_t inode = 0; inode < hin->GetNbinsX(); inode++) { | |
345 | ||
346 | Double_t ptnode = hin->GetBinCenter(inode + 1); | |
347 | TH1 *hout = (TH1 *)hin->Clone(Form("%s_extremehard", hin->GetName())); | |
348 | ||
349 | for (Int_t ibin = 0; ibin < hin->GetNbinsX(); ibin++) { | |
350 | if (hin->GetBinError(ibin + 1) <= 0.) continue; | |
351 | Double_t val = hin->GetBinContent(ibin + 1); | |
352 | Double_t err = hin->GetBinError(ibin + 1); | |
353 | Double_t cen = hin->GetBinCenter(ibin + 1); | |
354 | if (cen < ptnode) | |
355 | err *= -1. + (cen - ptlow) / (ptnode - ptlow); | |
356 | else | |
357 | err *= (cen - ptnode) / (pthigh - ptnode); | |
358 | ||
359 | hout->SetBinContent(ibin + 1, val + sign * err); | |
360 | } | |
361 | ||
362 | Double_t diff = TMath::Abs(mean - hout->GetMean()); | |
363 | if (diff > maxdiff) { | |
364 | // printf("found max at %f\n", ptnode); | |
365 | if (hmax) delete hmax; | |
366 | hmax = (TH1 *)hout->Clone("hmax"); | |
367 | maxdiff = diff; | |
368 | } | |
369 | delete hout; | |
370 | } | |
371 | return hmax; | |
372 | } | |
373 | ||
374 | YieldMean_IntegralMean(TH1 *hdata, TH1 *hlo, TH1 *hhi, Double_t &integral, Double_t &mean) | |
375 | { | |
376 | ||
377 | /* | |
378 | * compute integrals | |
379 | */ | |
380 | ||
381 | Double_t cont, err, width, cent; | |
382 | Double_t I = 0., IX = 0., Ierr = 0., IXerr = 0., Ilerr = 0., IXlerr = 0.; | |
383 | Double_t M = 0., Merr = 0., Mlerr = 0., C; | |
384 | ||
385 | /* integrate the data */ | |
386 | for (Int_t ibin = 0; ibin < hdata->GetNbinsX(); ibin++) { | |
387 | cent = hdata->GetBinCenter(ibin + 1); | |
388 | width = hdata->GetBinWidth(ibin + 1); | |
389 | cont = width * hdata->GetBinContent(ibin + 1); | |
390 | err = width * hdata->GetBinError(ibin + 1); | |
391 | if (err <= 0.) continue; | |
392 | I += cont; | |
393 | IX += cont * cent; | |
394 | } | |
395 | /* integrate low */ | |
396 | for (Int_t ibin = 0; ibin < hlo->GetNbinsX(); ibin++) { | |
397 | cent = hlo->GetBinCenter(ibin + 1); | |
398 | width = hlo->GetBinWidth(ibin + 1); | |
399 | cont = width * hlo->GetBinContent(ibin + 1); | |
400 | err = width * hlo->GetBinError(ibin + 1); | |
401 | if (err <= 0.) continue; | |
402 | I += cont; | |
403 | IX += cont * cent; | |
404 | } | |
405 | /* integrate high */ | |
406 | for (Int_t ibin = 0; ibin < hhi->GetNbinsX(); ibin++) { | |
407 | cent = hhi->GetBinCenter(ibin + 1); | |
408 | width = hhi->GetBinWidth(ibin + 1); | |
409 | cont = width * hhi->GetBinContent(ibin + 1); | |
410 | err = width * hhi->GetBinError(ibin + 1); | |
411 | if (err <= 0.) continue; | |
412 | I += cont; | |
413 | IX += cont * cent; | |
414 | } | |
415 | ||
416 | /* set values */ | |
417 | integral = I; | |
418 | mean = IX / I; | |
419 | } |