diff --git a/coffea/hist/plot.py b/coffea/hist/plot.py index ae8a428c1..3f326ad00 100644 --- a/coffea/hist/plot.py +++ b/coffea/hist/plot.py @@ -78,6 +78,24 @@ def clopper_pearson_interval(num, denom, coverage=_coverage1sd): return interval +def normal_interval(pw, tw, pw2, tw2, coverage=_coverage1sd): + # see https://root.cern.ch/doc/master/TEfficiency_8cxx_source.html#l02515 + + eff = pw / tw + + variance = (pw2 * (1 - 2 * eff) + tw2 * eff**2) / (tw**2) + sigma = np.sqrt(variance) + + prob = 0.5 * (1 - coverage) + delta = np.zeros_like(sigma) + delta[sigma != 0] = scipy.stats.norm.ppf(prob, scale=sigma[sigma != 0]) + + lo = eff - np.minimum(eff + delta, np.ones_like(eff)) + hi = np.maximum(eff - delta, np.zeros_like(eff)) - eff + + return np.array([lo, hi]) + + def plot1d(hist, ax=None, clear=True, overlay=None, stack=False, overflow='none', line_opts=None, fill_opts=None, error_opts=None, legend_opts={}, overlay_overflow='none', density=False, binwnorm=None): """Create a 1D plot from a 1D or 2D `Hist` object @@ -352,6 +370,8 @@ def plotratio(num, denom, ax=None, clear=True, overflow='none', error_opts=None, rsumw_err = np.abs(clopper_pearson_interval(sumw_num, sumw_num + sumw_denom) - rsumw) elif unc == 'num': rsumw_err = np.abs(poisson_interval(rsumw, sumw2_num / sumw_denom**2) - rsumw) + elif unc == "normal": + rsumw_err = np.abs(normal_interval(sumw_num, sumw_denom, sumw2_num, sumw2_denom)) else: raise ValueError("Unrecognized uncertainty option: %r" % unc) diff --git a/tests/test_hist_plot.py b/tests/test_hist_plot.py index a2cae2c6b..ce98916a4 100644 --- a/tests/test_hist_plot.py +++ b/tests/test_hist_plot.py @@ -241,3 +241,47 @@ def test_clopper_pearson_interval(): threshold = 1e-6 assert(all((interval[1, :] / ref_hi) - 1 < threshold)) assert(all((interval[0, :] / ref_lo) - 1 < threshold)) + +def test_normal_interval(): + from coffea.hist.plot import normal_interval + + # Reference weighted efficiency and error from ROOTs TEfficiency + + denom = np.array([ 89.01457591590004, 2177.066076428943 , 6122.5256890981855 , + 0. , 100.27757990710668]) + num = np.array([ 75.14287743709515, 2177.066076428943 , 5193.454723043864 , + 0. , 84.97723540536361]) + denom_sumw2 = np.array([ 94.37919737476827, 10000. , 6463.46795877633 , + 0. , 105.90898005417333]) + num_sumw2 = np.array([ 67.2202147680005 , 10000. , 4647.983931785646 , + 0. , 76.01275761253757]) + ref_hi = np.array([0.0514643476600107, 0. , 0.0061403263960343, + np.nan, 0.0480731185500146]) + ref_lo = np.array([0.0514643476600107, 0. , 0.0061403263960343, + np.nan, 0.0480731185500146]) + + interval = normal_interval(num, denom, num_sumw2, denom_sumw2) + threshold = 1e-6 + + lo, hi = interval + + assert len(ref_hi) == len(hi) + assert len(ref_lo) == len(lo) + + for i in range(len(ref_hi)): + if np.isnan(ref_hi[i]): + assert np.isnan(ref_hi[i]) + elif ref_hi[i] == 0.0: + assert hi[i] == 0.0 + else: + assert np.abs(hi[i] / ref_hi[i] - 1) < threshold + + if np.isnan(ref_lo[i]): + assert np.isnan(ref_lo[i]) + elif ref_lo[i] == 0.0: + assert lo[i] == 0.0 + else: + assert np.abs(lo[i] / ref_lo[i] - 1) < threshold + + +