Last active
August 19, 2022 07:32
-
-
Save JohannesBuchner/b25425ac0bf5ae5fa99b892a3578c04d to your computer and use it in GitHub Desktop.
The bias of rebinning to a minimum of bmin=1 background counts with ftgrouppha, then using wstat to estimate background contribution
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
import matplotlib.pyplot as plt | |
import numpy as np | |
def rebin( | |
Nbins = 40, | |
minimum = 0.1, | |
): | |
bins = np.linspace(0, 1, Nbins) | |
lam = minimum + 0 * bins | |
c = np.random.poisson(lam) | |
edges = np.unique([0] + (np.where(c!=0)[0]+1).tolist()[:-1] + [Nbins]) | |
estimate = np.zeros(Nbins) * np.nan | |
for lo, hi in zip(edges[:-1], edges[1:]): | |
estimate[lo:hi] = c[lo:hi].sum() / (hi - lo) | |
return (estimate - lam) / lam | |
def simulate_many(Nsim, | |
Nbins = 40, | |
minimum = 0.1, | |
): | |
return np.mean([rebin(Nbins=Nbins, minimum=minimum) for i in range(Nsim)], axis=0) | |
for minimum in np.linspace(0.1, 4.1, 11): | |
plt.plot(100 * simulate_many(10000, Nbins=40, minimum=minimum), label='%.2f' % minimum) | |
plt.ylabel('Estimation error (%)') | |
plt.ylim(-25, 25) | |
plt.legend(title='counts per channel') | |
plt.xlabel('Channel') | |
plt.tight_layout() | |
plt.savefig('wstatrebin-counts.pdf') | |
plt.savefig('wstatrebin-counts.png') | |
plt.close() |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment