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