Error propagation failure for TH1 objects with the TH1::kPoisson bin error option

[stempl]

Dear experts,

I have two TH1F histograms with Poisson errors (obtained via hist->SetBinErrorOption(TH1::kPoisson)). The Sumw2() is not explicitly called for them on purpose.

Now, I would like to perform a histogram division that correctly propagates the errors. While this works on a technical level (insofar as no errors or warnings occur when calling hist_num->Divide(hist_den)), the resulting bin errors turn out to be incorrect. Most strikingly, the lower bin errors are zero if the resulting histogram bins are smaller than one.

Please see an attached minimal example that reproduces this behavior and its pasted output below.

Could you point me towards a possible solution and tell me how to obtain the result that I need (i.e., a TH1 object with correct bin errors that is the ratio of two histograms with Poisson errors).

I haven’t been able to find a workaround that produces a suited TH1 object as the result, in particular because needed functions like TH1::SetBinErrorLow(Int_t, Double_t) and TH1::SetBinErrorUp(Int_t, Double_t) do not exist (in contrast to their getter functions) and there seem to be no obvious way how to set lower and upper errors manually. Further, any initial call of the Sumw2() function does not give me the desired result as this leads to asymmetric errors being ignored altogether further downstream.

I am aware of the existing functionality of TGraphAsymmErrors, but due to my current setup I would prefer to find a solution that allows me to work with TH1 objects, if at all possible.

Many thanks in advance, best regards,
Sebastian

Minimal example:

#include <iostream>
void printHistogram(TH1F*);

void poissonerrors() {
	TH1F *num = new TH1F("num", "", 2, 0, 2);
	TH1F *den = new TH1F("den", "", 2, 0, 2);

	// setting the following bin error option will break error propagation later on
	num->SetBinErrorOption(TH1::kPoisson);
	den->SetBinErrorOption(TH1::kPoisson);

	// fill some dummy data
	num->SetBinContent(1, 2.0);
	num->SetBinContent(2, 1.0);
	den->SetBinContent(1, 1.0);
	den->SetBinContent(2, 2.0);

	// bin errors are reasonable (see the output)
	printHistogram(num);
	printHistogram(den);

	// the following silently spoils the bin errors
	TH1F *ratio = (TH1F*)num->Clone("ratio");
	ratio->Divide(den);

	// inspect the computed bin errors (see the output)
	printHistogram(ratio);
}

void printHistogram(TH1F* h) {
	for(Int_t i=1; i<=h->GetNbinsX(); i++) {
		std::cout << "[" <<
			h->GetName() << "] bin "
			<< i << ": "
			<< h->GetBinContent(i)  << " (-"
			<< h->GetBinErrorLow(i) << ", +"
			<< h->GetBinErrorUp(i) << ")"
			<< std::endl;
	}
}

Produced output

[num] bin 1: 2 (-1.29181, +2.63786)
[num] bin 2: 1 (-0.827246, +2.29953)
[den] bin 1: 1 (-0.827246, +2.29953)
[den] bin 2: 2 (-1.29181, +2.63786)
[ratio] bin 1: 2 (-1.29181, +2.63786)
[ratio] bin 2: 0.5 (-0, +1.34102)

[couet]

I think @moneta can help you.

[moneta]

Hi,

For computing error of Poisson ratios you should use TGraphAsymmErrors::Divide. If one histogram is a sub-sample of t he other (binomial case), you can use also Efficiency, but you cannot use TH1::Divide. This works only in the normal approximation, so not when you have Poisson errors

Lorenzo

[stempl]

Hi Lorenzo,

Thanks a lot for the prompt answer!

Okay, I see. So in other words - given that there are no functions for setting lower/upper errors in TH1 objects - there is simply no way of assembling a histogram object that contains the desired data, not even by playing some tricks? The only way is to fully switch to TGraphAsymmErrors objects for all related analysis routines?

In any case, this is good to know. Thank you for the help!

Cheers,
Sebastian

PS: Would there be any objections to adding the functions TH1::SetBinErrorLow(Int_t bin; Double_t error) and TH1::SetBinErrorUp(Int_t bin; Double_t error) in a future ROOT release? This way, one could at least do the error computation via the TGraphAsymmErrors methods but then, if really needed, set up a TH1 object containing the relevant data.

[moneta]

Hi,
The problem is that histograms do not store lower/upper errors. The Poisson errors are just computed when requested. For this reason one has to use the TGraphAsymErrors .Adding the storage of both errors will be a large change to the histogram classes and we prefer doing that for the new redesigned histograms planned for ROOT version 7

Cheers
Lorenzo

[stempl]

Hi Lorenzo,

Thanks for this background information! Now I see that this proposed change would be much more work than just adding a handful of simple functions and that it makes a lot of sense to keep this for the next major ROOT version.

Thanks again, cheers,
Sebastian