In [1]:
import sys
sys.path.append('/home/jbourbeau/cr-composition')
sys.path
Out[1]:
In [2]:
import os
import sys
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
import argparse
import seaborn.apionly as sns
from icecube import ShowerLLH
import composition as comp
# from composition.analysis.load_sim import load_sim
# import composition.analysis.data_functions as datafunc
# import composition.analysis.plotting_functions as plotting
In [3]:
sns.set_palette('muted')
sns.set_color_codes()
In [5]:
df, cut_dict = comp.load_sim(return_cut_dict=True)
df_old, cut_dict_old = comp.load_sim(old=True, return_cut_dict=True)
selection_mask = np.array([True] * len(df))
selection_mask_old = np.array([True] * len(df_old))
standard_cut_keys = ['reco_exists', 'reco_zenith', 'num_hits_1_60', 'IT_signal',
'StationDensity', 'max_qfrac_1_60', 'reco_containment',
'energy_range']
for key in standard_cut_keys:
selection_mask *= cut_dict[key]
selection_mask_old *= cut_dict_old[key]
df = df[selection_mask]
df_old = df_old[selection_mask_old]
In [4]:
MC_log_energy = df.MC_log_energy
reco_log_energy = df.reco_log_energy
MC_log_energy_old = df_old.MC_log_energy
reco_log_energy_old = df_old.reco_log_energy
In [12]:
energy_bins = np.arange(6.2, 8, 0.05)
fig, ax = plt.subplots()
plotting.histogram_2D(MC_log_energy, reco_log_energy, energy_bins, log_counts=True)
plt.plot([0,10], [0,10], marker='None', linestyle='-.')
plt.xlim([6.2, 8])
plt.ylim([6.2, 8])
plt.xlabel('$\log_{10}(E_{\mathrm{MC}}/\mathrm{GeV})$')
plt.ylabel('$\log_{10}(E_{\mathrm{reco}}/\mathrm{GeV})$')
plt.title('True vs. ShowerLLH Energy')
plt.show()
In [5]:
energy_res = reco_log_energy - MC_log_energy
energy_res_old = reco_log_energy_old - MC_log_energy_old
MC_proton_mask = (df['MC_comp'] == 'P')
MC_iron_mask = (df['MC_comp'] == 'Fe')
MC_proton_mask_old = (df_old['MC_comp'] == 'P')
MC_iron_mask_old = (df_old['MC_comp'] == 'Fe')
energy_bins = np.linspace(6.2, 8, 30)
bin_centers, bin_medians_proton, error_proton = datafunc.get_medians(MC_log_energy[MC_proton_mask],
energy_res[MC_proton_mask],
energy_bins)
bin_centers, bin_medians_iron, error_iron = datafunc.get_medians(MC_log_energy[MC_iron_mask],
energy_res[MC_iron_mask],
energy_bins)
bin_centers, bin_medians_proton_old, error_proton_old = datafunc.get_medians(MC_log_energy_old[MC_proton_mask_old],
energy_res_old[MC_proton_mask_old],
energy_bins)
bin_centers, bin_medians_iron_old, error_iron_old = datafunc.get_medians(MC_log_energy_old[MC_iron_mask_old],
energy_res_old[MC_iron_mask_old],
energy_bins)
In [6]:
fig, ax = plt.subplots()
ax.errorbar(bin_centers, bin_medians_proton, yerr=error_proton, marker='.', label='MC Proton', alpha=0.75)
ax.errorbar(bin_centers, bin_medians_iron, yerr=error_iron, marker='.', label='MC Iron', alpha=0.75)
ax.errorbar(bin_centers, bin_medians_proton_old, yerr=error_proton_old, marker='.', label='MC Proton (old)', alpha=0.75)
ax.errorbar(bin_centers, bin_medians_iron_old, yerr=error_iron_old, marker='.', label='MC Iron (old)', alpha=0.75)
ax.axhline(0, marker='None', linestyle='-.')
ax.set_xlabel('$\log_{10}(E_{\mathrm{MC}}/\mathrm{GeV})$')
ax.set_ylabel('$\log_{10}(E_{\mathrm{reco}}/E_{\mathrm{MC}})$')
plt.legend()
plt.show()
In [13]:
fig, ax = plt.subplots()
energy_mask = (df['MC_log_energy'] >= 6.4)
weights = df['MC_energy']**-2
print(weights[MC_proton_mask & energy_mask].values)
n, bins, patches = ax.hist(energy_res[MC_proton_mask & energy_mask].values, bins=np.arange(-0.4, 0.4, 0.01),
weights=weights[MC_proton_mask & energy_mask].values, histtype='step', label='Proton', normed=False)
n, bins, patches = ax.hist(energy_res[MC_iron_mask & energy_mask].values, bins=np.arange(-0.4, 0.4, 0.01),
weights=weights[MC_iron_mask & energy_mask].values, histtype='step', label='Iron', normed=False)
n, bins, patches = ax.hist(e_h4a_res[energy_mask].values, bins=np.arange(-0.4, 0.4, 0.01),
weights=weights[energy_mask].values, histtype='step', label='H4a', normed=False)
ax.set_xlim([-0.4, 0.4])
ax.set_xlabel('$\log_{10}(E_{\mathrm{reco}}/E_{\mathrm{MC}})$')
ax.set_ylabel('Counts')
ax.set_title('$\log_{10}(E_{\mathrm{MC}}/\mathrm{GeV}) \geq 6.4$')
plt.legend()
plt.show()
In [8]:
MC_x = df.MC_x
MC_y = df.MC_y
reco_x = df.reco_x
reco_y = df.reco_y
MC_x_old = df_old.MC_x
MC_y_old = df_old.MC_y
reco_x_old = df_old.reco_x
reco_y_old = df_old.reco_y
In [9]:
core_res = np.sqrt((reco_x - MC_x)**2+(reco_y - MC_y)**2)
bin_centers, reco_bin_medians, reco_error = datafunc.get_medians(MC_log_energy, core_res, energy_bins)
core_res_old = np.sqrt((reco_x_old - MC_x_old)**2+(reco_y_old - MC_y_old)**2)
bin_centers, reco_bin_medians_old, reco_error_old = datafunc.get_medians(MC_log_energy_old, core_res_old, energy_bins)
In [10]:
fig, ax = plt.subplots()
ax.errorbar(bin_centers, reco_bin_medians, yerr=reco_error, marker='.', label='ShowerLLH')
ax.errorbar(bin_centers, reco_bin_medians_old, yerr=reco_error_old, marker='.', label='ShowerLLH old')
ax.set_xlabel('$\log_{10}(E_{\mathrm{MC}}/\mathrm{GeV})$')
ax.set_ylabel('$\\vec{x}_{\mathrm{reco}}-\\vec{x}_{\mathrm{MC}} \ [m]$')
plt.legend()
plt.show()
In [23]:
core_pos = np.sqrt(MC_x**2+MC_y**2).values
core_pos_old = np.sqrt(MC_x_old**2+MC_y_old**2).values
pos_bins = np.linspace(0, 800, 30)
bin_centers, bin_medians, error = datafunc.get_medians(core_pos, energy_res, pos_bins)
bin_centers, bin_medians_old, error_old = datafunc.get_medians(core_pos_old, energy_res_old, pos_bins)
In [24]:
fig, ax = plt.subplots()
ax.errorbar(bin_centers, bin_medians, yerr=error, marker='.', label='ShowerLLH', alpha=0.75)
ax.errorbar(bin_centers, bin_medians_old, yerr=error_old, marker='.', label='ShowerLLH old', alpha=0.75)
ax.axhline(0, marker='None', linestyle='-.')
ax.set_xlabel('Core position [m]')
ax.set_ylabel('$\log_{10}(E_{\mathrm{reco}}/E_{\mathrm{MC}})$')
plt.legend()
plt.show()
In [31]:
LLH_grid = ShowerLLH.LLHGrid('7006')
tank_xy = LLH_grid.tank_xy
tank_x, tank_y = np.transpose(tank_xy)
In [32]:
dist_bins = np.arange(-1700, 1700, 25)
fig, ax = plt.subplots()
plotting.histogram_2D(reco_x, reco_y, dist_bins)
ax.plot(tank_x, tank_y, marker='o', markersize=5,
color='white', alpha=0.5)
ax.set_xlabel('X [m]')
ax.set_ylabel('Y [m]')
ax.set_title('ShowerLLH reconstructed positions')
ax.set_xlim([-800,800])
ax.set_ylim([-800,800])
plt.show()
In [33]:
fig, ax = plt.subplots()
plotting.histogram_2D(lap_x, lap_y, dist_bins)
ax.plot(tank_x, tank_y, marker='o', markersize=5,
color='white', alpha=0.5)
ax.set_xlabel('X [m]')
ax.set_ylabel('Y [m]')
ax.set_title('Laputop reconstructed positions')
ax.set_xlim([-800,800])
ax.set_ylim([-800,800])
plt.show()
In [ ]: