In [1]:
%load_ext autoreload
%autoreload 2
In [2]:
from __future__ import division, print_function
import os
from collections import namedtuple
import numpy as np
import pandas as pd
import pyunfold
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import seaborn as sns
import comptools as comp
sns.set_context(context='paper', font_scale=1.5)
color_dict = comp.color_dict
%matplotlib inline
In [3]:
config = 'IC86.2012'
num_groups = 2
comp_list = comp.get_comp_list(num_groups=num_groups)
energybins = comp.get_energybins(config)
num_ebins = len(energybins.log_energy_midpoints)
unfolding_dir = os.path.join(comp.paths.comp_data_dir,
config,
'unfolding')
In [4]:
feature_list, feature_labels = comp.get_training_features()
Load trained models
In [5]:
energy_pipeline = comp.load_trained_model('linearregression_energy_{}'.format(config))
pipeline_str = 'xgboost_comp_{}_{}-groups'.format(config, num_groups)
comp_pipeline = comp.load_trained_model(pipeline_str)
In [6]:
# Load DataFrame with saved prior distributions
df_file = os.path.join(unfolding_dir,
'unfolding-df_{}-groups.hdf'.format(num_groups))
df_priors = pd.read_hdf(df_file)
df_priors.head()
Out[6]:
In [7]:
# Load fitted efficiencies and calculate effective areas
eff_path = os.path.join(comp.paths.comp_data_dir, config, 'efficiencies',
'efficiency_fit_num_groups_{}.hdf'.format(num_groups))
df_eff = pd.read_hdf(eff_path)
df_eff.head()
Out[7]:
In [8]:
# Get simulation thrown areas for each energy bin
thrown_radii = comp.simfunctions.get_sim_thrown_radius(energybins.log_energy_midpoints)
thrown_area = np.pi * thrown_radii.max()**2
thrown_area
Out[8]:
In [9]:
efficiencies, efficiencies_err = comp.get_detector_efficiencies(config=config,
num_groups=num_groups,
sigmoid='slant',
pyunfold_format=True)
Load and process data
In [10]:
print('Loading data into memory...')
df_data = comp.load_data(config=config,
processed=True,
columns=feature_list,
n_jobs=15,
verbose=True)
# data_dir = os.path.join(comp.paths.comp_data_dir,
# 'IC86.2012',
# 'data')
# filepath_up_shift = os.path.join(data_dir,
# 'data_dataframe_vem_cal_up.hdf')
# df_data = pd.read_hdf(filepath_up_shift,
# columns=feature_list)
In [11]:
X_data = df_data[feature_list].values
In [12]:
print('Running composition classifications...')
df_data['pred_comp_target'] = comp_pipeline.predict(X_data)
print('Running energy regressions...')
df_data['reco_log_energy'] = energy_pipeline.predict(X_data)
In [13]:
counts_observed = {}
counts_observed_err = {}
for comp_idx, composition in enumerate(comp_list):
print(composition)
# Filter out events that don't pass composition & energy mask
pred_comp_mask = df_data['pred_comp_target'] == comp_idx
mask = pred_comp_mask
energies = df_data.loc[mask, 'reco_log_energy'].values
comp_counts, _ = np.histogram(energies, bins=energybins.log_energy_bins)
counts_observed[composition] = comp_counts
counts_observed_err[composition] = np.sqrt(comp_counts)
counts_observed_err['total'] = np.sqrt(np.sum(counts_observed_err[composition]**2 for composition in comp_list))
# Calculate total counts
counts_observed['total'] = np.sum(counts_observed[composition] for composition in comp_list)
# Format observed counts, detection efficiencies, and priors for PyUnfold use
counts_pyunfold = np.empty(num_groups * len(energybins.energy_midpoints))
counts_err_pyunfold = np.empty(num_groups * len(energybins.energy_midpoints))
for idx, composition in enumerate(comp_list):
counts_pyunfold[idx::num_groups] = counts_observed[composition]
counts_err_pyunfold[idx::num_groups] = counts_observed_err[composition]
In [14]:
counts_observed
Out[14]:
In [15]:
# Effective area
eff_area, eff_area_err = {}, {}
for composition in comp_list + ['total']:
eff_area[composition] = df_eff['eff_median_{}'.format(composition)].values * thrown_area
eff_area_err[composition] = df_eff['eff_err_low_{}'.format(composition)].values * thrown_area
# Solid angle
theta_max = 40 if config == 'IC79.2010' else 65
solid_angle = np.pi*np.sin(np.deg2rad(theta_max))**2
# Livetime
livetime, livetime_err = comp.get_detector_livetime(config=config)
In [16]:
def counts_to_flux(counts, counts_err=None, composition=None):
return comp.get_flux(counts, counts_err,
energybins=energybins.energy_bins,
eff_area=thrown_area,
livetime=livetime,
livetime_err=livetime_err,
solid_angle=solid_angle,
scalingindex=2.7)
In [17]:
UnfoldedCounts = namedtuple('UnfoldedCounts', ['counts', 'counts_sys_err', 'counts_stat_err'])
In [18]:
UnfoldedCounts(counts=[1, 2, 3],
counts_sys_err=[4, 5, 6],
counts_stat_err=[7, 8, 9])
Out[18]:
In [22]:
def get_unfolded_counts(df_file):
print('df_file = {}'.format(df_file))
# Load simulation and train composition classifier
df_sim_train, df_sim_test = comp.load_sim(df_file=df_file,
config=config,
energy_reco=False,
log_energy_min=None,
log_energy_max=None,
test_size=0.5,
n_jobs=10,
verbose=True)
print('Running energy reconstruction...')
for df in [df_sim_train, df_sim_test]:
X = df_sim_train[feature_list].values
# Energy reconstruction
df['reco_log_energy'] = energy_pipeline.predict(df[feature_list].values)
df['reco_energy'] = 10**df['reco_log_energy']
print('Running composition classifications...')
pred_target = comp_pipeline.predict(df_sim_test[feature_list].values)
print('Making response matrix...')
log_reco_energy_sim_test = df_sim_test['reco_log_energy']
log_true_energy_sim_test = df_sim_test['MC_log_energy']
true_target = df_sim_test['comp_target_{}'.format(num_groups)].values
response, response_err = comp.response_matrix(true_energy=log_true_energy_sim_test,
reco_energy=log_reco_energy_sim_test,
true_target=true_target,
pred_target=pred_target,
efficiencies=efficiencies,
efficiencies_err=efficiencies_err,
energy_bins=energybins.log_energy_bins)
fig, ax = plt.subplots()
im = ax.imshow(response, origin='lower', cmap='viridis')
ax.plot([0, response.shape[0] - 1], [0, response.shape[1] - 1], marker='None', ls=':', color='white')
comp.plotting.colorbar(im, label='Normalized response matrix')
ax.set_xlabel('True bins')
ax.set_ylabel('Reconstructed bins')
plt.show()
# Run unfolding for each of the priors
prior_name = 'H4a'
prior = None if prior_name == 'uniform' else df_priors['{}_prior'.format(prior_name)]
# priors = 'Jeffreys' if prior_name == 'Jeffreys' else df['{}_prior'.format(prior_name)]
unfolded_result = pyunfold.iterative_unfold(data=counts_pyunfold,
data_err=counts_err_pyunfold,
response=response,
response_err=response_err,
efficiencies=efficiencies,
efficiencies_err=efficiencies_err,
ts='ks',
ts_stopping=0.005,
prior=prior,
return_iterations=True)
# callbacks=[pyunfold.callbacks.SplineRegularizer(degree=1, smooth=256)])
# unfolding_results[prior_name] = df_unfolding_iter
# counts, counts_sys_err, counts_stat_err = comp.unfolded_counts_dist(unfolded_result,
# num_groups=num_groups)
# unfolded_counts = UnfoldedCounts(counts=counts,
# counts_sys_err=counts_sys_err,
# counts_stat_err=counts_stat_err)
return unfolded_result
In [23]:
sim_files = {'up': '/data/user/jbourbeau/composition/IC86.2012/sim/sim_dataframe_vem_cal_down.hdf',
'down': '/data/user/jbourbeau/composition/IC86.2012/sim/sim_dataframe_vem_cal_up.hdf',
'nominal': '/data/user/jbourbeau/composition/IC86.2012/sim/processed_hdf/*.hdf',
}
unfolded_results = {key: get_unfolded_counts(file) for key, file in sim_files.items()}
In [24]:
unfolded_results
Out[24]:
In [27]:
fig, ax = plt.subplots()
# fig, ax = plt.subplots(figsize=(10, 5))
linestyles = iter(['-', ':', '-.', '--'])
markers = iter(['.', '^', 'X'])
for shift_type, df_unfolded in unfolded_results.items():
marker = next(markers)
ls = next(linestyles)
# if prior_name != 'H4a':
# continue
counts, counts_sys_err, counts_stat_err = comp.unfolded_counts_dist(df_unfolded,
num_groups=num_groups)
for composition in comp_list + ['total']:
# print(composition)
# print(counts[composition])
flux, flux_sys_err = counts_to_flux(counts=counts[composition],
counts_err=counts_sys_err[composition],
composition=composition)
flux, flux_stat_err = counts_to_flux(counts=counts[composition],
counts_err=counts_stat_err[composition],
composition=composition)
comp.plot_steps(energybins.log_energy_bins, flux,
yerr=flux_sys_err,
color=color_dict[composition],
ls=ls,
ax=ax)
# label = '{} ({})'.format(composition, prior_name.replace('_', '-'))
if composition == 'total':
label = shift_type.capitalize() if shift_type == 'nominal' else '{} 3\%'.format(shift_type.capitalize())
else:
label = ''
ax.errorbar(energybins.log_energy_midpoints, flux,
yerr=flux_stat_err,
color=color_dict[composition],
ls='None',
marker=marker,
capsize=20,
elinewidth=1,
label=label)
ax.set_yscale('log', nonposy='clip')
ax.set_xlabel('$\mathrm{\log_{10}(E/GeV)}$')
ax.set_ylabel('$\mathrm{ E^{2.7} \ J(E) \ [GeV^{1.7} m^{-2} sr^{-1} s^{-1}]}$')
ax.set_xlim(6.4, 7.8)
ax.set_ylim(1e3, 1e5)
ax.grid(lw=0.5, which='both')
ax.legend(title='$\mathrm{S_{125}}$ shift')
# ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))
# ax.legend(loc='upper center',
# bbox_to_anchor=(0.5, 1.75),
# ncol=2,
# fancybox=True)
outfile = os.path.join(comp.paths.figures_dir,
'systematics',
'vem-calibration-flux-comparison.png')
comp.check_output_dir(outfile)
plt.savefig(outfile)
plt.show()
In [31]:
unfolded_counts = {}
for shift_type, df_unfolded in unfolded_results.items():
counts, counts_sys_err, counts_stat_err = comp.unfolded_counts_dist(df_unfolded,
num_groups=num_groups)
unfolded_counts[shift_type] = {composition: counts[composition] for composition in comp_list + ['total']}
unfolded_counts
Out[31]:
In [77]:
fig, ax = plt.subplots(figsize=(10, 6))
# unfolded_counts = {}
linestyles = iter(['-', ':', '-.', '--'])
markers = iter(['o', '^', 'X'])
vem_systematic = []
for shift_type in ['down', 'up']:
shift_systematic = {}
marker = next(markers)
ls = next(linestyles)
for composition in comp_list + ['total']:
nominal = unfolded_counts['nominal'][composition]
deviation = (unfolded_counts[shift_type][composition] - nominal) / nominal
energy_mask = np.logical_and(energybins.log_energy_midpoints < 7.8,
energybins.log_energy_midpoints > 6.4)
ax.plot(energybins.log_energy_midpoints[energy_mask], deviation[energy_mask],
color=color_dict[composition],
marker=marker,
ls=ls,
label='{} 3\% {}'.format(shift_type, composition))
shift_systematic[composition] = np.mean(np.abs(deviation[energy_mask]))
vem_systematic.append(shift_systematic)
ax.axhline(0.0, marker='None', ls='-.', color='k')
ax.set(xlabel='$\mathrm{\log_{10}(E/GeV)}$',
xlim=(6.4, 7.8),
ylabel='\% Difference w.r.t. Nominal')
ax.grid()
# ax.legend()
ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))
outfile = os.path.join(comp.paths.figures_dir,
'systematics',
'vem-calibration-counts-deviation.png')
comp.check_output_dir(outfile)
plt.savefig(outfile)
plt.show()
In [87]:
np.stack((unfolded_counts['down']['total'], unfolded_counts['up']['total'])).max(axis=0)
Out[87]:
In [98]:
fig, ax = plt.subplots()
# fig, ax = plt.subplots(figsize=(10, 6))
for composition in comp_list + ['total']:
nominal = unfolded_counts['nominal'][composition]
max_counts = np.stack((unfolded_counts['down'][composition], unfolded_counts['up'][composition])).max(axis=0)
min_counts = np.stack((unfolded_counts['down'][composition], unfolded_counts['up'][composition])).min(axis=0)
deviation_upper = (max_counts - nominal) / nominal
deviation_lower = (min_counts - nominal) / nominal
energy_mask = np.logical_and(energybins.log_energy_midpoints < 7.8,
energybins.log_energy_midpoints > 6.4)
ax.fill_between(energybins.log_energy_midpoints[energy_mask],
deviation_upper[energy_mask],
deviation_lower[energy_mask],
color=color_dict[composition],
alpha=0.4,
label=composition.capitalize())
ax.axhline(0.0, marker='None', ls='-.', color='k')
ax.set(xlabel='$\mathrm{\log_{10}(E/GeV)}$',
xlim=(6.4, 7.8),
ylabel='\% Difference w.r.t. Nominal')
ax.grid()
ax.legend(loc='upper left')
# ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))
outfile = os.path.join(comp.paths.figures_dir,
'systematics',
'vem-calibration-counts-deviation-bands.png')
comp.check_output_dir(outfile)
plt.savefig(outfile)
plt.show()
In [126]:
fig, axes = plt.subplots(ncols=1, nrows=3,
sharey=True, sharex=True,
gridspec_kw={'hspace':0.15})
for idx, (composition, ax) in enumerate(zip(comp_list + ['total'], axes)):
nominal = unfolded_counts['nominal'][composition]
max_counts = np.stack((unfolded_counts['down'][composition], unfolded_counts['up'][composition])).max(axis=0)
min_counts = np.stack((unfolded_counts['down'][composition], unfolded_counts['up'][composition])).min(axis=0)
deviation_upper = (max_counts - nominal) / nominal
deviation_lower = (min_counts - nominal) / nominal
energy_mask = np.logical_and(energybins.log_energy_midpoints < 7.8,
energybins.log_energy_midpoints > 6.4)
ax.fill_between(energybins.log_energy_midpoints[energy_mask],
deviation_upper[energy_mask],
deviation_lower[energy_mask],
color=color_dict[composition],
alpha=0.4,
label=composition.capitalize())
ax.axhline(0.0, marker='None', ls='-.', color='k')
ax.set_xlim(6.4, 7.8)
if idx == 1:
ax.set_ylabel('\% Difference w.r.t. Nominal')
ax.grid()
ax.legend(loc='lower left')
ax.set(xlabel='$\mathrm{\log_{10}(E/GeV)}$')
# ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))
outfile = os.path.join(comp.paths.figures_dir,
'systematics',
'vem-calibration-counts-deviation-bands-2.png')
comp.check_output_dir(outfile)
plt.savefig(outfile)
plt.tight_layout()
plt.show()
In [78]:
vem_systematic = pd.DataFrame.from_records(vem_systematic, index=['down', 'up'])
In [115]:
vem_systematic
Out[115]:
In [28]:
# fig, ax = plt.subplots()
# im = ax.imshow(response, origin='lower', cmap='viridis')
# ax.plot([0, response.shape[0] - 1], [0, response.shape[1] - 1], marker='None', ls=':', color='white')
# comp.plotting.colorbar(im, label='Normalized response matrix')
# ax.set_xlabel('True bins')
# ax.set_ylabel('Reconstructed bins')
# plt.show()
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [23]:
# Run unfolding for each of the priors
names = ['uniform', 'H3a', 'H4a', 'simple_power_law']
# names = ['Jeffreys', 'H3a', 'H4a', 'Polygonato']
unfolding_results = {}
for prior_name in names:
prior = None if prior_name == 'uniform' else df_priors['{}_prior'.format(prior_name)]
# priors = 'Jeffreys' if prior_name == 'Jeffreys' else df['{}_prior'.format(prior_name)]
df_unfolding_iter = pyunfold.iterative_unfold(data=counts_pyunfold,
data_err=counts_err_pyunfold,
response=response,
response_err=response_err,
efficiencies=efficiencies,
efficiencies_err=efficiencies_err,
ts='ks',
ts_stopping=0.005,
prior=prior,
return_iterations=True)
# callbacks=[pyunfold.callbacks.SplineRegularizer(degree=1, smooth=256)])
unfolding_results[prior_name] = df_unfolding_iter
In [24]:
names
Out[24]:
In [25]:
fig, ax = plt.subplots()
# fig, ax = plt.subplots(figsize=(10, 5))
linestyles = ['-', ':', '-.', '--']
markers = ['.', '*', '^', 'X']
for prior_name, ls, marker in zip(names, linestyles, markers):
# if prior_name != 'H4a':
# continue
counts, counts_sys_err, counts_stat_err = comp.unfolded_counts_dist(unfolding_results[prior_name],
num_groups=num_groups)
for composition in comp_list + ['total']:
# print(composition)
# print(counts[composition])
flux, flux_sys_err = counts_to_flux(counts=counts[composition],
counts_err=counts_sys_err[composition],
composition=composition)
flux, flux_stat_err = counts_to_flux(counts=counts[composition],
counts_err=counts_stat_err[composition],
composition=composition)
comp.plot_steps(energybins.log_energy_bins, flux,
yerr=flux_sys_err,
color=color_dict[composition],
ls=ls,
ax=ax)
# label = '{} ({})'.format(composition, prior_name.replace('_', '-'))
label = prior_name.replace('_', ' ') if composition == 'total' else ''
ax.errorbar(energybins.log_energy_midpoints, flux,
yerr=flux_stat_err,
color=color_dict[composition],
ls='None',
marker=marker,
capsize=20,
elinewidth=1,
label=label)
ax.set_yscale('log', nonposy='clip')
ax.set_xlabel('$\mathrm{\log_{10}(E/GeV)}$')
ax.set_ylabel('$\mathrm{ E^{2.7} \ J(E) \ [GeV^{1.7} m^{-2} sr^{-1} s^{-1}]}$')
ax.set_xlim(6.4, 7.8)
ax.set_ylim(1e3, 1e5)
ax.grid(lw=1, which='both')
ax.legend()
# ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))
# ax.legend(loc='upper center',
# bbox_to_anchor=(0.5, 1.75),
# ncol=2,
# fancybox=True)
# outfile = os.path.join(comp.paths.figures_dir, 'skymaps', 'allsky.png')
# comp.check_output_dir(outfile)
# plt.savefig(outfile)
plt.show()
In [220]:
unfolding_results[prior_name]
Out[220]:
In [224]:
fig, ax = plt.subplots()
# fig, ax = plt.subplots(figsize=(10, 5))
prior_name = 'H4a'
marker = '.'
ls = '-'
counts, counts_sys_err, counts_stat_err = comp.unfolded_counts_dist(unfolding_results[prior_name],
num_groups=num_groups)
for composition in comp_list + ['total']:
flux, flux_sys_err = counts_to_flux(counts=counts[composition],
counts_err=counts_sys_err[composition],
composition=composition)
flux, flux_stat_err = counts_to_flux(counts=counts[composition],
counts_err=counts_stat_err[composition],
composition=composition)
comp.plot_steps(energybins.log_energy_bins, flux,
yerr=flux_sys_err,
color=color_dict[composition],
fillalpha=0.4,
ls=ls,
label=composition,
ax=ax)
ax.errorbar(energybins.log_energy_midpoints, flux,
yerr=flux_stat_err,
color=color_dict[composition],
ls='None',
marker=marker,
capsize=20,
elinewidth=1)
ax.set_yscale('log', nonposy='clip')
ax.set_xlabel('$\mathrm{\log_{10}(E/GeV)}$')
ax.set_ylabel('$\mathrm{ E^{2.7} \ J(E) \ [GeV^{1.7} m^{-2} sr^{-1} s^{-1}]}$')
ax.set_xlim(6.4, 7.8)
ax.set_ylim(1e3, 1e5)
ax.grid(lw=1, which='both')
ax.legend()
# ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))
# ax.legend(loc='upper center',
# bbox_to_anchor=(0.5, 1.75),
# ncol=2,
# fancybox=True)
outfile = os.path.join(comp.paths.figures_dir, 'skymaps', 'allsky.png')
comp.check_output_dir(outfile)
plt.savefig(outfile)
plt.show()