'''
This is the hard working code in order to calculate ULs, sensitivities,
and time to detections.
'''
import matplotlib.pyplot as plt
from scipy.interpolate import interpolate
from scipy.optimize import brentq
from scipy import integrate
import numpy as np
import corner
[docs]def upper_limit(t_obs, lambda_lim, a_eff, e_0, plot_resolution=30):
'''
This function generates all plots for the command 'ul' from
input data. It takes:
t_obs in seconds
lambda_lim from Rolke, Knoetig, ...
a_eff A path to the file with effective area over true energy after cuts
plot_resolution a parameter for running tests faster
It returns a dictionary with results.
'''
a_eff_interpol = get_effective_area(a_eff)
# make the figures
phasespace_figure = get_ul_phasespace_figure(
t_obs,
lambda_lim,
a_eff_interpol,
e_0,
pixels_per_line=plot_resolution)
spectrum_figure, energy_x, dn_de_y = get_ul_spectrum_figure(
t_obs,
lambda_lim,
a_eff_interpol,
e_0,
n_points_to_plot=plot_resolution)
sensitive_energy_figure, gamma_s, e_sens_s = get_sensitive_energy_figure(
a_eff_interpol
)
a_eff_figure = get_effective_area_figure(a_eff_interpol)
figures = {
'ul_phasespace': phasespace_figure,
'ul_integral_spectral_exclusion_zone': spectrum_figure,
'ul_sensitive_energy': sensitive_energy_figure,
'ul_effective_area': a_eff_figure
}
dictionary = {
'plots': figures,
'data': {
'ul_integral_spectral_exclusion_zone':
np.transpose((energy_x, dn_de_y)),
'ul_sensitive_energy':
np.transpose((gamma_s, e_sens_s))
}
}
return dictionary
[docs]def sensitivity(sigma_bg, alpha, t_obs, a_eff, e_0, plot_resolution=30):
'''
This function generates all plots for the command 'sens' from
input data. It takes:
sigma_bg in 1/seconds
alpha on/off exposure ratio
t_obs observation time in seconds
a_eff A path to the file with effective area over true energy after cuts
It returns a dictionary with results.
'''
a_eff_interpol = get_effective_area(a_eff)
# make the figures
phasespace_figure = get_sens_phasespace_figure(
sigma_bg,
alpha,
t_obs,
a_eff_interpol,
pixels_per_line=plot_resolution)
spectrum_figure, energy_x, dn_de_y = get_sens_spectrum_figure(
sigma_bg,
alpha,
t_obs,
a_eff_interpol,
e_0,
n_points_to_plot=plot_resolution
)
sensitive_energy_figure, gamma_s, e_sens_s = get_sensitive_energy_figure(
a_eff_interpol
)
a_eff_figure = get_effective_area_figure(a_eff_interpol)
figures = {
'sens_phasespace': phasespace_figure,
'sens_integral_spectral_exclusion_zone': spectrum_figure,
'sens_sensitive_energy': sensitive_energy_figure,
'sens_effective_area': a_eff_figure
}
dictionary = {
'plots': figures,
'data': {
'sens_integral_spectral_exclusion_zone':
np.transpose((energy_x, dn_de_y)),
'sens_sensitive_energy':
np.transpose((gamma_s, e_sens_s))
}
}
return dictionary
[docs]def predict(
sigma_bg,
alpha,
f_0,
df_0,
gamma,
dgamma,
e_0,
a_eff,
plot_resolution=30
):
'''
This function generates all plots for the command 'predict' from
input data. It takes:
sigma_bg in 1/seconds
alpha on/off exposure ratio
f_0 flux normalization in 1/(cm^2 s TeV)
df_0 flux normalization error 1 sigma 1/(cm^2 s TeV)
gamma power law index (<0)
dgamma power law index error 1 sigma
e_0 reference energy in TeV
a_eff A path to the file with effective area over true energy after cuts
It returns a dictionary with results.
'''
a_eff_interpol = get_effective_area(a_eff)
# make the figures
phasespace_figure = get_predict_phasespace_figure(
sigma_bg,
alpha,
f_0,
df_0,
gamma,
dgamma,
e_0,
a_eff_interpol,
pixels_per_line=plot_resolution)
t_obs_samples = get_t_obs_samples(
sigma_bg,
alpha,
f_0,
df_0,
gamma,
dgamma,
e_0,
a_eff_interpol,
n_samples=plot_resolution*3000
)
t_obs_est = get_t_obs_est_from_samples(t_obs_samples)
spectrum_figure, energy_x, dn_de_ys = get_predict_spectrum_figure(
sigma_bg,
alpha,
t_obs_est,
f_0,
df_0,
gamma,
dgamma,
e_0,
a_eff_interpol,
n_points_to_plot=plot_resolution
)
sensitive_energy_figure, gamma_s, e_sens_s = get_sensitive_energy_figure(
a_eff_interpol
)
a_eff_figure = get_effective_area_figure(a_eff_interpol)
figures = {
'predict_phasespace': phasespace_figure,
'predict_integral_spectral_exclusion_zone': spectrum_figure,
'predict_sensitive_energy': sensitive_energy_figure,
'predict_effective_area': a_eff_figure
}
dictionary = {
'plots': figures,
'data': {
'predict_integral_spectral_exclusion_zone':
np.vstack((energy_x, dn_de_ys)).T,
'predict_sensitive_energy':
np.transpose((gamma_s, e_sens_s)),
'predict_t_obs_est':
np.transpose(t_obs_est)
}
}
return dictionary
[docs]def get_effective_area_test_relative_paths():
'''
Helper function to get the paths of stored effective areas
'''
a_eff_test_relative_paths = [
'/resources/A_eff/MAGIC_lowZd_Ecut_300GeV.dat',
'/resources/A_eff/MAGIC_medZd_Ecut_300GeV.dat',
'/resources/A_eff/VERITAS_V5_lowZd_McCutcheon.dat',
'/resources/A_eff/uhe_test_aperture.dat',
'/resources/A_eff/AeffEnergy_P8R2_OnAxis_Total.dat'
]
return a_eff_test_relative_paths
[docs]def get_effective_area(a_eff_path):
'''
Function to get the interpolated effective area
from a file path
'''
a_eff_data = np.loadtxt(a_eff_path, delimiter=',')
# interpolate the data points, every energy outside definition range
# from the data file is assumed to have 0 effective area
a_eff_interpol = interpolate.interp1d(
a_eff_data[:, 0],
a_eff_data[:, 1],
bounds_error=False,
fill_value=0.
)
return a_eff_interpol
[docs]def prepare_phasespace_meshes(
t_obs, lambda_lim, a_eff_interpol, e_0, pixels_per_line, gamma=-2.6):
'''
determine parameter plot ranges
'''
f_0 = get_ul_f_0(t_obs, lambda_lim, a_eff_interpol, e_0, gamma)
f_0_limits, gamma_limits = get_f_0_gamma_limits(f_0, gamma)
f_0_mesh, gamma_mesh = get_f_0_gamma_mesh(
f_0_limits,
gamma_limits,
pixels_per_line)
return f_0_mesh, gamma_mesh
[docs]def plot_predict_contours_from_phasespace_parameters(
f_0,
df_0,
gamma,
dgamma,
n_samples=100000
):
'''
This function draws contours from phase space parameters
'''
random_data = get_phasespace_samples(
f_0,
df_0,
gamma,
dgamma,
n_samples=n_samples
)
plot_contours_from_sample(random_data)
[docs]def get_phasespace_samples(
f_0,
df_0,
gamma,
dgamma,
n_samples,
):
'''
Function to return a sample from the phase space according to
the given parameters and parameter uncertainties
assuming there can only be::
f_0 > 0
gamma < 0
'''
mean = [f_0*1e12, gamma] # scale f_0, later downscale again
# diagonal covariance
cov = [[df_0*1e12*df_0*1e12, 0.], [0., dgamma*dgamma]]
random_data = np.random.multivariate_normal(mean, cov, n_samples*2)
random_data[:, 0] = random_data[:, 0]/1e12
# truncate to physical values
random_data = random_data[random_data[:, 0] > 0] # f_0
random_data = random_data[random_data[:, 1] < 0] # gamma
if len(random_data) < n_samples:
raise IndexError(
'less than 50 percent of the sample from the highest density' +
'confidence interval in the physical region f_0>0 && Gamma<0' +
'-> time to detection is uncertain/ source may be undetectable')
return random_data[:n_samples]
[docs]def get_t_obs_samples(
sigma_bg,
alpha,
f_0,
df_0,
gamma,
dgamma,
e_0,
a_eff_interpol,
n_samples=100000,
thinning=1./400.
):
'''
Function to calculate samples of the time to detection
'''
phasespace_samples = get_phasespace_samples(
f_0,
df_0,
gamma,
dgamma,
n_samples
)
t_obs_samples = np.array([
t_obs_li_ma_criterion(
point[0] * effective_area_averaged_flux(
point[1],
e_0,
a_eff_interpol
),
sigma_bg,
alpha)
for point
in phasespace_samples[:(int(n_samples*thinning)+1)]
])
return t_obs_samples
[docs]def get_t_obs_est_from_samples(
t_obs_samples,
lower_percentile=16.,
upper_percentile=84.
):
'''
Function to calculate the asymmetrical 68% CI
from a sample of times to detection
'''
t_obs_est = np.array([
np.percentile(t_obs_samples, lower_percentile),
np.percentile(t_obs_samples, 50.),
np.percentile(t_obs_samples, upper_percentile),
])
return t_obs_est
[docs]def get_energy_range(a_eff_interpol):
'''
Get the definition energy range of the effective area
for integration and plotting purposes.
'''
return np.power(10, np.array([
a_eff_interpol.x.min()*0.999,
a_eff_interpol.x.max()*1.001
]))
[docs]def get_f_0_gamma_limits(f_0, gamma):
'''
Get a nice box power law phase space box for plotting
'''
f_0_limits = [f_0*0.1, f_0*1.9]
gamma_limits = [gamma-1., gamma+1.]
if gamma_limits[1] > 0:
gamma_limits[1] = 0.
return f_0_limits, gamma_limits
[docs]def get_f_0_gamma_mesh(f_0_limits, gamma_limits, pixels_per_line):
'''
Generate two numpy.meshgrids for 2D plotting
'''
f_0_stepsize = (f_0_limits[1]-f_0_limits[0])/pixels_per_line
gamma_stepsize = (gamma_limits[1]-gamma_limits[0])/pixels_per_line
f_0_stepsize = f_0_stepsize+f_0_stepsize*1e-9
gamma_stepsize = gamma_stepsize+gamma_stepsize*1e-9
f_0_buf = np.arange(f_0_limits[0], f_0_limits[1], f_0_stepsize)
gamma_buf = np.arange(gamma_limits[1], gamma_limits[0], -gamma_stepsize)
f_0_mesh, gamma_mesh = np.meshgrid(f_0_buf, gamma_buf)
return f_0_mesh, gamma_mesh
[docs]def get_ul_f_0(t_obs, lambda_lim, a_eff_interpol, e_0, gamma):
'''
Calculate f_0 on the exclusion line from solving the boundary condition
lambda_lim = lambda_s
'''
return lambda_lim / t_obs / effective_area_averaged_flux(
gamma, e_0, a_eff_interpol)
[docs]def sensitive_energy(gamma, a_eff_interpol):
'''
Function returning the sensitive energy, given gamma
and the effective area
'''
mu = ln_energy_weighted(
gamma,
a_eff_interpol
)/effective_area_averaged_flux(
gamma,
1.,
a_eff_interpol
)
return np.exp(mu)
[docs]def get_useful_e_0(a_eff_interpol):
'''
This function will try to yield a useful e_0
value for anchoring the power law on.
This is calculated from the effective area definition range
'''
# round to the next order of magnitude
log10_e_0_suggestion = round(mean_log10_energy(a_eff_interpol), 0)
return 10**(log10_e_0_suggestion) / 10
[docs]def get_gamma_from_sensitive_energy(E_sens, a_eff_interpol):
'''
numerical inverse of the sensitive energy
'''
gamma_min = -30.
gamma_max = -0.05
gamma_num = brentq(lambda x: (sensitive_energy(
gamma=x,
a_eff_interpol=a_eff_interpol
) - E_sens
), gamma_min, gamma_max
)
return gamma_num
[docs]def effective_area_averaged_flux(gamma, e_0, a_eff_interpol):
'''
I define this in the paper as c(gamma)
'''
energy_range = get_energy_range(a_eff_interpol)
integrand = lambda x: power_law(
x,
f_0=1.,
gamma=gamma,
e_0=e_0
)*a_eff_interpol(np.log10(x))
return integrate.quad(
integrand,
energy_range[0],
energy_range[1],
limit=10000,
full_output=1,
points=[energy_range[0], energy_range[0]*10]
)[0]
[docs]def ln_energy_weighted(gamma, a_eff_interpol):
'''
For calculating the sensitive energy.
This function gets the unnormalized mean ln(energy)
integrated over the power law flux and sensitive area
'''
f_0 = 1.
energy_range = get_energy_range(a_eff_interpol)
integrand = lambda x: power_law(
x,
f_0=f_0,
gamma=gamma,
e_0=1.
)*a_eff_interpol(np.log10(x))*np.log(x)
return integrate.quad(
integrand,
energy_range[0],
energy_range[1],
limit=10000,
full_output=1,
points=[energy_range[0], energy_range[0]*10]
)[0]
[docs]def mean_log10_energy(a_eff_interpol):
'''
For calculating mean of the log10 of the energy
over the effective area, to estimate a useful e_0
'''
energy_range = get_energy_range(a_eff_interpol)
integrand_numerator = lambda x: a_eff_interpol(np.log10(x))*np.log10(x)
integrand_denominator = lambda x: a_eff_interpol(np.log10(x))
numerator = integrate.quad(
integrand_numerator,
energy_range[0],
energy_range[1],
limit=10000,
full_output=1,
points=[energy_range[0], energy_range[0]*10]
)[0]
denominator = integrate.quad(
integrand_denominator,
energy_range[0],
energy_range[1],
limit=10000,
full_output=1,
points=[energy_range[0], energy_range[0]*10]
)[0]
return numerator/denominator
[docs]def power_law(energy, f_0, gamma, e_0=1.):
'''
A power law function as defined in the paper
'''
return f_0*(energy/e_0)**(gamma)
[docs]def plot_effective_area(a_eff_interpol, style='k', label=''):
'''
fill a plot with the effective energy from the supplied
interpolated data
'''
start = a_eff_interpol.x.min()
stop = a_eff_interpol.x.max()
samples = 1000
energy_samples = np.linspace(start, stop, samples)
area_samples = np.array([
a_eff_interpol(energy)
for energy
in energy_samples
])
plt.plot(np.power(10, energy_samples), area_samples/10000.,
style,
label=label)
plt.loglog()
plt.title('Effective Area')
plt.xlabel('Energy / TeV')
plt.ylabel('A$_{eff}$ / m$^2$')
return
[docs]def plot_sensitive_energy(a_eff_interpol):
'''
fill a sensitive energy plot figure
'''
gamma_range = [-5., -0.5]
stepsize = 0.1
gammas = np.arange(gamma_range[0], gamma_range[1]+stepsize, stepsize)
e_sens = np.array([sensitive_energy(i, a_eff_interpol) for i in gammas])
plt.plot(gammas, e_sens, 'k')
plt.title('sensitive energy E$_{sens}$($\\Gamma$)')
plt.semilogy()
plt.ylabel('E$_{sens}$ / TeV')
plt.xlabel('$\\Gamma$')
return gammas, e_sens
[docs]def plot_lambda_s_mesh(
t_obs,
lambda_lim,
a_eff_interpol,
e_0,
f_0_mesh,
gamma_mesh,
n_levels=9,
linestyles='dashed',
linewidths=1,
colors='k'
):
'''
Function to get the lambda_s plot in the phase space of the power law
'''
pixels_per_line = np.shape(f_0_mesh)[0]
lambda_s = np.array([[t_obs*f_0_mesh[i, j]*effective_area_averaged_flux(
gamma_mesh[i, j],
e_0=e_0,
a_eff_interpol=a_eff_interpol
) for j in range(pixels_per_line)] for i in range(pixels_per_line)])
levels = np.array([lambda_lim/((1.5)**(int(n_levels/2)-i))
for i in range(n_levels)])
limit_index = np.where(levels == lambda_lim)[0][0]
linestyles = [linestyles for i in range(n_levels)]
linestyles[limit_index] = 'solid'
linewidths = [linewidths for i in range(n_levels)]
linewidths[limit_index] = 2
cset = plt.contour(
f_0_mesh,
gamma_mesh,
lambda_s,
levels=levels,
linestyles=linestyles,
linewidths=linewidths,
colors=colors
)
plt.clabel(cset, inline=True, fmt='%1.1f', fontsize=10)
plt.title(
'signal counts per {0:1.1e} h, E$_0$={1:1.1e} TeV assuming power law'.
format(t_obs/3600., e_0)
)
plt.xlabel('$f_0$ / [(cm$^2$ s TeV)$^{-1}$]')
plt.ylabel('$\\Gamma$')
return lambda_s
[docs]def plot_t_obs_mesh(
sigma_bg,
alpha,
a_eff_interpol,
e_0,
f_0_mesh,
gamma_mesh,
n_levels=9,
linestyles='dashed',
linewidths=1,
colors='k',
t_obs=None,
):
'''
This function puts the times until detection into
a plot in phase space (f_0, Gamma) for given telescope
analysis parameters
'''
pixels_per_line = np.shape(f_0_mesh)[0]
t_obs_s = np.array([[
t_obs_li_ma_criterion(
# calculate the sigma_s from c(Gamma)*f_0
f_0_mesh[i, j] *
effective_area_averaged_flux(
gamma_mesh[i, j],
e_0,
a_eff_interpol),
sigma_bg,
alpha)/3600.
for j in range(pixels_per_line)] for i in range(pixels_per_line)])
# if there is no predefined info about t_obs
# in prediction mode -> get it from the t_obs mesh
print_solid = True
if t_obs is None:
print_solid = False
t_obs = np.median(t_obs_s.flatten())
else:
t_obs = t_obs/3600.
levels = np.array([t_obs/((1.5)**(int(n_levels/2)-i))
for i in range(n_levels)])
limit_index = np.where(
np.isclose(levels, t_obs)
)[0][0]
linestyles = [linestyles for i in range(n_levels)]
if print_solid:
linestyles[limit_index] = 'solid'
linewidths = [linewidths for i in range(n_levels)]
if print_solid:
linewidths[limit_index] = 2
cset = plt.contour(
f_0_mesh,
gamma_mesh,
t_obs_s,
levels=levels,
linestyles=linestyles,
linewidths=linewidths,
colors=colors
)
plt.clabel(cset, inline=True, fmt='%1.1e', fontsize=10)
plt.title(
'time to detection in h, E$_0$={0:1.1e} TeV assuming power law'.
format(e_0)
)
plt.xlabel('$f_0$ / [(cm$^2$ s TeV)$^{-1}$]')
plt.ylabel('$\\Gamma$')
return t_obs_s
[docs]def plot_contours_from_sample(
samples,
levels=None,
smooth=False
):
'''
This function draws a contour plot into the current
figure showing selected highest density 2D confidence intervals
The standard is: 1 sigma, 2 sigma, 3 sigma in 2D
'''
# one, two, and three sigma in 2D
if levels is None:
levels = 1 - np.exp(-0.5*(np.arange(1, 3.1, 1)**2))
axis = plt.gca()
x_lim = axis.get_xlim()
y_lim = axis.get_ylim()
corner.hist2d(
samples[:, 0],
samples[:, 1],
plot_datapoints=False,
plot_density=False,
no_fill_contours=True,
ax=axis,
smooth=smooth,
levels=levels)
axis.set_xlim(x_lim)
axis.set_ylim(y_lim)
[docs]def plot_power_law(
f_0,
gamma,
e_0,
energy_range,
fmt='k:',
label='',
alpha_plot=0.7
):
'''
This function generates a power law plot in
the current figure
'''
e_x = 10**np.arange(
np.log10(energy_range[0]),
np.log10(energy_range[1])+0.05,
0.05)
e_y = power_law(e_x, f_0, gamma, e_0=e_0)
plt.plot(e_x, e_y, fmt, label=label, alpha=alpha_plot)
plt.loglog()
plt.xlabel("E / TeV")
plt.ylabel("dN/dE / [(cm$^2$ s TeV)$^{-1}$]")
[docs]def plot_source_emission_spectrum_with_uncertainties(
phasespace_samples,
e_0,
energy_range,
label=''
):
'''
This function draws 100 power laws as an illustration for
the source emission uncertainty
'''
f0_mc, gamma_mc = zip(*np.percentile(
phasespace_samples,
[16, 50, 84],
axis=0)
)
for f_0, gamma in phasespace_samples[
np.random.randint(len(phasespace_samples), size=100)
]:
plot_power_law(
f_0,
gamma,
e_0=e_0,
energy_range=energy_range,
fmt='k',
alpha_plot=0.03
)
plot_power_law(
f0_mc[1],
gamma_mc[1],
e_0=e_0,
energy_range=energy_range,
fmt='r',
label=label,
alpha_plot=0.8
)
[docs]def integral_spectral_exclusion_zone(
energy, lambda_lim, a_eff_interpol, t_obs, e_0
):
'''
This function returns the integral spectral exclusion zone value
at one point in energy for given lambda_lim, a_eff_interpol, and t_obs
'''
f_0, gamma = integral_spectral_exclusion_zone_parameters(
energy,
lambda_lim,
a_eff_interpol,
t_obs,
e_0
)
return power_law(energy, f_0, gamma, e_0)
[docs]def integral_spectral_exclusion_zone_parameters(
energy,
lambda_lim,
a_eff_interpol,
t_obs,
e_0=1.
):
'''
This function calculates the integral spectral exclusion zone parameters
f_0 and gamma at at a given energy in order to draw it into spectral plots.
It is done by utilizing the Lagrangian results from the paper.
'''
gamma_calc = get_gamma_from_sensitive_energy(energy, a_eff_interpol)
f_0_calc = get_ul_f_0(t_obs, lambda_lim, a_eff_interpol, e_0, gamma_calc)
return f_0_calc, gamma_calc
[docs]def t_obs_li_ma_criterion(sigma_s, sigma_bg, alpha, threshold=5.):
'''
This function calculates the limit average time to detection,
given a signal rate, background rate, and alpha
'''
estimated_rate_in_off = sigma_bg / alpha
estimated_rate_in_on = sigma_s + sigma_bg
t_obs_min = 0.
t_obs_max = 36e16
try:
t_obs = brentq(lambda x: (
li_ma_significance(
estimated_rate_in_on*x,
estimated_rate_in_off*x,
alpha
) - threshold
), t_obs_min, t_obs_max)
except ValueError:
raise ValueError('The time to detection could not be calculated. '
'This can be the case when sigma_bg is overestimated.'
' Please check that your stated background rate'
' is actually given as: sigma_bg / s')
return t_obs
[docs]def sigma_lim_li_ma_criterion(sigma_bg, alpha, t_obs, threshold=5.):
'''
This function returns the limit signal count rate
using the LiMa criterion
'''
estimated_bg_counts_in_off = sigma_bg * t_obs / alpha
estimated_bg_counts_in_on = sigma_bg * t_obs
sigma_lim_min = 0.
sigma_lim_max = 1e5 # more than 100k gamma / s is likely unrealistic
sigma_lim = brentq(lambda x: (
li_ma_significance(
x*t_obs + estimated_bg_counts_in_on,
estimated_bg_counts_in_off,
alpha
) - threshold
), sigma_lim_min, sigma_lim_max)
# implement the low statistics limit which most authors use:
# "at least 10 excess counts"
if sigma_lim*t_obs < 10.:
sigma_lim = 10./t_obs
return sigma_lim
[docs]def li_ma_significance(n_on, n_off, alpha):
'''
A function to calculate the significance, according to :
Li, T-P., and Y-Q. Ma.
"Analysis methods for results in gamma-ray astronomy."
The Astrophysical Journal 272 (1983): 317-324.
'''
n_on = np.array(n_on, copy=False, ndmin=1)
n_off = np.array(n_off, copy=False, ndmin=1)
with np.errstate(divide='ignore', invalid='ignore'):
p_on = n_on / (n_on + n_off)
p_off = n_off / (n_on + n_off)
t1 = n_on * np.log(((1 + alpha) / alpha) * p_on)
t2 = n_off * np.log((1 + alpha) * p_off)
ts = (t1 + t2)
significance = np.sqrt(ts * 2)
significance[np.isnan(significance)] = 0
significance[n_on < alpha * n_off] = 0
return significance