Author: James Bourbeau
In [1]:
%load_ext watermark
%watermark -u -d -v -p numpy,matplotlib,scipy,pandas,sklearn,mlxtend
In [2]:
%matplotlib inline
from __future__ import division, print_function
from collections import defaultdict
import itertools
import numpy as np
from scipy import interp
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap
import seaborn.apionly as sns
import matplotlib as mpl
from sklearn.metrics import accuracy_score, confusion_matrix, roc_curve, auc, classification_report
from sklearn.model_selection import cross_val_score, StratifiedShuffleSplit, KFold, StratifiedKFold
from mlxtend.feature_selection import SequentialFeatureSelector as SFS
import composition as comp
import composition.analysis.plotting as plotting
color_dict = comp.analysis.get_color_dict()
[ back to top ]
In [5]:
bin_midpoints, _, counts, counts_err = comp.get1d('/home/jbourbeau/PyUnfold/unfolded_output_h3a.root', 'NC', 'Unf_ks_ACM/bin0')
Whether or not to train on 'light' and 'heavy' composition classes, or the individual compositions
In [6]:
comp_class = True
comp_list = ['light', 'heavy'] if comp_class else ['P', 'He', 'O', 'Fe']
Get composition classifier pipeline
In [7]:
pipeline_str = 'GBDT'
pipeline = comp.get_pipeline(pipeline_str)
Define energy binning for this analysis
In [8]:
energybins = comp.analysis.get_energybins()
[ back to top ]
In [9]:
sim_train, sim_test = comp.preprocess_sim(comp_class=comp_class, return_energy=True)
In [10]:
# Compute the correlation matrix
df_sim = comp.load_dataframe(datatype='sim', config='IC79')
feature_list, feature_labels = comp.analysis.get_training_features()
In [11]:
fig, ax = plt.subplots()
df_sim[df_sim.MC_comp_class == 'light'].avg_inice_radius.plot(kind='hist', bins=50, ax=ax, alpha=0.75)
df_sim[df_sim.MC_comp_class == 'heavy'].avg_inice_radius.plot(kind='hist', bins=50, ax=ax, alpha=0.75)
ax.grid()
plt.show()
In [12]:
fig, ax = plt.subplots()
df_sim[df_sim.MC_comp_class == 'light'].invcharge_inice_radius.plot(kind='hist', bins=50, ax=ax, alpha=0.75)
df_sim[df_sim.MC_comp_class == 'heavy'].invcharge_inice_radius.plot(kind='hist', bins=50, ax=ax, alpha=0.75)
ax.grid()
plt.show()
In [13]:
fig, ax = plt.subplots()
df_sim[df_sim.MC_comp_class == 'light'].max_inice_radius.plot(kind='hist', bins=50, ax=ax, alpha=0.75)
df_sim[df_sim.MC_comp_class == 'heavy'].max_inice_radius.plot(kind='hist', bins=50, ax=ax, alpha=0.75)
ax.grid()
plt.show()
In [14]:
corr = df_sim[feature_list].corr()
# Generate a mask for the upper triangle
mask = np.zeros_like(corr, dtype=np.bool)
mask[np.triu_indices_from(mask)] = True
fig, ax = plt.subplots()
sns.heatmap(corr, mask=mask, cmap='RdBu_r', center=0,
square=True, xticklabels=feature_labels, yticklabels=feature_labels,
linewidths=.5, cbar_kws={'label': 'Covariance'}, annot=True, ax=ax)
# outfile = args.outdir + '/feature_covariance.png'
# plt.savefig(outfile)
plt.show()
In [15]:
label_replacement = {feature: labels for feature, labels in zip(feature_list, feature_labels)}
with plt.rc_context({'text.usetex': False}):
g = sns.pairplot(df_sim.sample(frac=1)[:1000], vars=feature_list, hue='MC_comp_class',
plot_kws={'alpha': 0.5, 'linewidth': 0},
diag_kws={'histtype': 'step', 'linewidth': 2, 'fill': True, 'alpha': 0.75, 'bins': 15})
for i in range(len(feature_list)):
for j in range(len(feature_list)):
xlabel = g.axes[i][j].get_xlabel()
ylabel = g.axes[i][j].get_ylabel()
if xlabel in label_replacement.keys():
g.axes[i][j].set_xlabel(label_replacement[xlabel])
if ylabel in label_replacement.keys():
g.axes[i][j].set_ylabel(label_replacement[ylabel])
g.fig.get_children()[-1].set_title('Comp class')
# g.fig.get_children()[-1].set_bbox_to_anchor((1.1, 0.5, 0, 0))
In [16]:
data = comp.preprocess_data(comp_class=comp_class, return_energy=True)
In [17]:
is_finite_mask = np.isfinite(data.X)
not_finite_mask = np.logical_not(is_finite_mask)
finite_data_mask = np.logical_not(np.any(not_finite_mask, axis=1))
data = data[finite_data_mask]
Run classifier over training and testing sets to get an idea of the degree of overfitting
In [12]:
clf_name = pipeline.named_steps['classifier'].__class__.__name__
print('=' * 30)
print(clf_name)
weights = sim_train.energy**-1.7
pipeline.fit(sim_train.X, sim_train.y)
# pipeline.fit(sim_train.X, sim_train.y, classifier__sample_weight=weights)
train_pred = pipeline.predict(sim_train.X)
train_acc = accuracy_score(sim_train.y, train_pred)
print('Training accuracy = {:.2%}'.format(train_acc))
test_pred = pipeline.predict(sim_test.X)
test_acc = accuracy_score(sim_test.y, test_pred)
print('Testing accuracy = {:.2%}'.format(test_acc))
print('=' * 30)
In [13]:
num_features = len(feature_list)
importances = pipeline.named_steps['classifier'].feature_importances_
indices = np.argsort(importances)[::-1]
fig, ax = plt.subplots()
for f in range(num_features):
print('{}) {}'.format(f + 1, importances[indices[f]]))
plt.ylabel('Feature Importances')
plt.bar(range(num_features),
importances[indices],
align='center')
plt.xticks(range(num_features),
feature_labels[indices], rotation=90)
plt.xlim([-1, len(feature_list)])
plt.show()
[ back to top ]
In [13]:
def get_frac_correct(train, test, pipeline, comp_list):
assert isinstance(train, comp.analysis.DataSet), 'train dataset must be a DataSet'
assert isinstance(test, comp.analysis.DataSet), 'test dataset must be a DataSet'
assert train.y is not None, 'train must have true y values'
assert test.y is not None, 'test must have true y values'
pipeline.fit(train.X, train.y)
test_predictions = pipeline.predict(test.X)
correctly_identified_mask = (test_predictions == test.y)
# Construct MC composition masks
MC_comp_mask = {}
for composition in comp_list:
MC_comp_mask[composition] = (test.le.inverse_transform(test.y) == composition)
MC_comp_mask['total'] = np.array([True]*len(test))
reco_frac, reco_frac_err = {}, {}
for composition in comp_list+['total']:
comp_mask = MC_comp_mask[composition]
# Get number of MC comp in each reco energy bin
num_MC_energy = np.histogram(test.log_energy[comp_mask],
bins=energybins.log_energy_bins)[0]
num_MC_energy_err = np.sqrt(num_MC_energy)
# Get number of correctly identified comp in each reco energy bin
num_reco_energy = np.histogram(test.log_energy[comp_mask & correctly_identified_mask],
bins=energybins.log_energy_bins)[0]
num_reco_energy_err = np.sqrt(num_reco_energy)
# Calculate correctly identified fractions as a function of MC energy
reco_frac[composition], reco_frac_err[composition] = comp.ratio_error(
num_reco_energy, num_reco_energy_err,
num_MC_energy, num_MC_energy_err)
return reco_frac, reco_frac_err
In [14]:
# Split training data into CV training and testing folds
kf = KFold(n_splits=10)
frac_correct_folds = defaultdict(list)
fold_num = 0
print('Fold ', end='')
for train_index, test_index in kf.split(sim_train.X):
fold_num += 1
print('{}...'.format(fold_num), end='')
reco_frac, reco_frac_err = get_frac_correct(sim_train[train_index],
sim_train[test_index],
pipeline, comp_list)
for composition in comp_list:
frac_correct_folds[composition].append(reco_frac[composition])
frac_correct_folds['total'].append(reco_frac['total'])
frac_correct_gen_err = {key: np.std(frac_correct_folds[key], axis=0) for key in frac_correct_folds}
# scores = np.array(frac_correct_folds['total'])
# score = scores.mean(axis=1).mean()
# score_std = scores.mean(axis=1).std()
In [15]:
avg_frac_correct_data = {'values': np.mean(frac_correct_folds['total'], axis=0), 'errors': np.std(frac_correct_folds['total'], axis=0)}
avg_frac_correct, avg_frac_correct_err = comp.analysis.averaging_error(**avg_frac_correct_data)
In [16]:
reco_frac, reco_frac_stat_err = get_frac_correct(sim_train, sim_test, pipeline, comp_list)
In [17]:
# Plot fraction of events correctlt classified vs energy
fig, ax = plt.subplots()
for composition in comp_list + ['total']:
err = np.sqrt(frac_correct_gen_err[composition]**2 + reco_frac_stat_err[composition]**2)
plotting.plot_steps(energybins.log_energy_midpoints, reco_frac[composition], err, ax,
color_dict[composition], composition)
plt.xlabel('$\log_{10}(E_{\mathrm{reco}}/\mathrm{GeV})$')
ax.set_ylabel('Fraction correctly identified')
ax.set_ylim([0.0, 1.0])
ax.set_xlim([energybins.log_energy_min, energybins.log_energy_max])
ax.grid()
leg = plt.legend(loc='upper center', frameon=False,
bbox_to_anchor=(0.5, # horizontal
1.1),# vertical
ncol=len(comp_list)+1, fancybox=False)
# set the linewidth of each legend object
for legobj in leg.legendHandles:
legobj.set_linewidth(3.0)
cv_str = 'Accuracy: {:0.2f}\% (+/- {:0.1f}\%)'.format(avg_frac_correct*100, avg_frac_correct_err*100)
ax.text(7.4, 0.2, cv_str,
ha="center", va="center", size=10,
bbox=dict(boxstyle='round', fc="white", ec="gray", lw=0.8))
plt.savefig('/home/jbourbeau/public_html/figures/frac-correct-{}.png'.format(pipeline_str))
plt.show()
In [18]:
# Plot the two-class decision scores
classifier_score = pipeline.decision_function(sim_train.X)
light_mask = sim_train.le.inverse_transform(sim_train.y) == 'light'
heavy_mask = sim_train.le.inverse_transform(sim_train.y) == 'heavy'
fig, ax = plt.subplots()
score_bins = np.linspace(-1, 1, 50)
ax.hist(classifier_score[light_mask], bins=score_bins, label='light', alpha=0.75)
ax.hist(classifier_score[heavy_mask], bins=score_bins, label='heavy', alpha=0.75)
ax.grid()
ax.legend()
plt.show()
In [13]:
import multiprocessing as mp
kf = KFold(n_splits=10)
frac_correct_folds = defaultdict(list)
# Define an output queue
output = mp.Queue()
# define a example function
def rand_string(length, output):
""" Generates a random string of numbers, lower- and uppercase chars. """
rand_str = ''.join(random.choice(
string.ascii_lowercase
+ string.ascii_uppercase
+ string.digits)
for i in range(length))
output.put(rand_str)
# Setup a list of processes that we want to run
processes = [mp.Process(target=get_frac_correct,
args=(sim_train[train_index],
sim_train[test_index],
pipeline, comp_list)) for train_index, test_index in kf.split(sim_train.X)]
# Run processes
for p in processes:
p.start()
# Exit the completed processes
for p in processes:
p.join()
# Get process results from the output queue
results = [output.get() for p in processes]
print(results)
[ back to top ]
In [18]:
def get_num_comp_reco(train, test, pipeline, comp_list):
assert isinstance(train, comp.analysis.DataSet), 'train dataset must be a DataSet'
assert isinstance(test, comp.analysis.DataSet), 'test dataset must be a DataSet'
assert train.y is not None, 'train must have true y values'
pipeline.fit(train.X, train.y)
test_predictions = pipeline.predict(test.X)
# Get number of correctly identified comp in each reco energy bin
num_reco_energy, num_reco_energy_err = {}, {}
for composition in comp_list:
# print('composition = {}'.format(composition))
comp_mask = train.le.inverse_transform(test_predictions) == composition
# print('sum(comp_mask) = {}'.format(np.sum(comp_mask)))
print(test.log_energy[comp_mask])
num_reco_energy[composition] = np.histogram(test.log_energy[comp_mask],
bins=energybins.log_energy_bins)[0]
num_reco_energy_err[composition] = np.sqrt(num_reco_energy[composition])
num_reco_energy['total'] = np.histogram(test.log_energy, bins=energybins.log_energy_bins)[0]
num_reco_energy_err['total'] = np.sqrt(num_reco_energy['total'])
return num_reco_energy, num_reco_energy_err
In [19]:
df_sim = comp.load_dataframe(datatype='sim', config='IC79')
In [20]:
df_sim[['log_dEdX', 'num_millipede_particles']].corr()
Out[20]:
In [21]:
max_zenith_rad = df_sim['lap_zenith'].max()
In [22]:
# Get number of events per energy bin
num_reco_energy, num_reco_energy_err = get_num_comp_reco(sim_train, data, pipeline, comp_list)
import pprint
pprint.pprint(num_reco_energy)
pprint.pprint(num_reco_energy_err)
# Solid angle
solid_angle = 2*np.pi*(1-np.cos(max_zenith_rad))
In [23]:
print(num_reco_energy['light'].sum())
print(num_reco_energy['heavy'].sum())
frac_light = num_reco_energy['light'].sum()/num_reco_energy['total'].sum()
print(frac_light)
In [24]:
# Live-time information
goodrunlist = pd.read_table('/data/ana/CosmicRay/IceTop_GRL/IC79_2010_GoodRunInfo_4IceTop.txt', skiprows=[0, 3])
goodrunlist.head()
Out[24]:
In [25]:
livetimes = goodrunlist['LiveTime(s)']
livetime = np.sum(livetimes[goodrunlist['Good_it_L2'] == 1])
print('livetime (seconds) = {}'.format(livetime))
print('livetime (days) = {}'.format(livetime/(24*60*60)))
In [26]:
fig, ax = plt.subplots()
for composition in comp_list + ['total']:
# Calculate dN/dE
y = num_reco_energy[composition]
y_err = num_reco_energy_err[composition]
plotting.plot_steps(energybins.log_energy_midpoints, y, y_err,
ax, color_dict[composition], composition)
ax.set_yscale("log", nonposy='clip')
plt.xlabel('$\log_{10}(E_{\mathrm{reco}}/\mathrm{GeV})$')
ax.set_ylabel('Counts')
# ax.set_xlim([6.3, 8.0])
# ax.set_ylim([10**-6, 10**-1])
ax.grid(linestyle=':')
leg = plt.legend(loc='upper center', frameon=False,
bbox_to_anchor=(0.5, # horizontal
1.1),# vertical
ncol=len(comp_list)+1, fancybox=False)
# set the linewidth of each legend object
for legobj in leg.legendHandles:
legobj.set_linewidth(3.0)
plt.savefig('/home/jbourbeau/public_html/figures/rate.png')
plt.show()
In [27]:
fig, ax = plt.subplots()
for composition in comp_list + ['total']:
# Calculate dN/dE
y = num_reco_energy[composition]
y_err = num_reco_energy_err[composition]
# Add time duration
# y = y / livetime
# y_err = y / livetime
y, y_err = comp.analysis.ratio_error(y, y_err, livetime, 0.005*livetime)
plotting.plot_steps(energybins.log_energy_midpoints, y, y_err,
ax, color_dict[composition], composition)
ax.set_yscale("log", nonposy='clip')
plt.xlabel('$\log_{10}(E_{\mathrm{reco}}/\mathrm{GeV})$')
ax.set_ylabel('Rate [s$^{-1}$]')
# ax.set_xlim([6.3, 8.0])
# ax.set_ylim([10**-6, 10**-1])
ax.grid(linestyle=':')
leg = plt.legend(loc='upper center', frameon=False,
bbox_to_anchor=(0.5, # horizontal
1.1),# vertical
ncol=len(comp_list)+1, fancybox=False)
# set the linewidth of each legend object
for legobj in leg.legendHandles:
legobj.set_linewidth(3.0)
plt.savefig('/home/jbourbeau/public_html/figures/rate.png')
plt.show()
In [28]:
df_sim, cut_dict_sim = comp.load_dataframe(datatype='sim', config='IC79', return_cut_dict=True)
selection_mask = np.array([True] * len(df_sim))
standard_cut_keys = ['IceTopQualityCuts', 'lap_InIce_containment',
'num_hits_1_60',
# 'num_hits_1_60', 'max_qfrac_1_60',
'InIceQualityCuts']
for key in standard_cut_keys:
selection_mask *= cut_dict_sim[key]
df_sim = df_sim[selection_mask]
In [23]:
def get_energy_res(df_sim, energy_bins):
reco_log_energy = df_sim['lap_log_energy'].values
MC_log_energy = df_sim['MC_log_energy'].values
energy_res = reco_log_energy - MC_log_energy
bin_centers, bin_medians, energy_err = comp.analysis.data_functions.get_medians(reco_log_energy,
energy_res,
energy_bins)
return np.abs(bin_medians)
In [37]:
def counts_to_flux(counts, counts_err, eff_area=156390.673059, livetime=1):
# Calculate dN/dE
y = counts/energybins.energy_bin_widths
y_err = counts_err/energybins.energy_bin_widths
# Add effective area
eff_area = np.array([eff_area]*len(y))
eff_area_error = np.array([0.01 * eff_area]*len(y_err))
y, y_err = comp.analysis.ratio_error(y, y_err, eff_area, eff_area_error)
# Add solid angle
y = y / solid_angle
y_err = y_err / solid_angle
# Add time duration
# y = y / livetime
# y_err = y / livetime
livetime = np.array([livetime]*len(y))
flux, flux_err = comp.analysis.ratio_error(y, y_err, livetime, 0.01*livetime)
# Add energy scaling
scaled_flux = energybins.energy_midpoints**2.7 * flux
scaled_flux_err = energybins.energy_midpoints**2.7 * flux_err
return scaled_flux, scaled_flux_err
In [38]:
# Plot fraction of events vs energy
# fig, ax = plt.subplots(figsize=(8, 6))
fig = plt.figure()
ax = plt.gca()
for composition in comp_list + ['total']:
y, y_err = counts_to_flux(num_reco_energy[composition], num_reco_energy_err[composition], livetime=livetime)
plotting.plot_steps(energybins.log_energy_midpoints, y, y_err, ax, color_dict[composition], composition)
ax.set_yscale("log", nonposy='clip')
plt.xlabel('$\log_{10}(E_{\mathrm{reco}}/\mathrm{GeV})$')
# ax.set_ylabel('$\mathrm{E}^{2.7} \\frac{\mathrm{dN}}{\mathrm{dE dA d\Omega dt}} \ [\mathrm{GeV}^{1.7} \mathrm{m}^{-2} \mathrm{sr}^{-1} \mathrm{s}^{-1}]$')
ax.set_ylabel('$\mathrm{E}^{2.7} \ J(E) \ [\mathrm{GeV}^{1.7} \mathrm{m}^{-2} \mathrm{sr}^{-1} \mathrm{s}^{-1}]$')
ax.set_xlim([6.4, 9.0])
ax.set_ylim([10**2, 10**5])
ax.grid(linestyle='dotted', which="both")
# Add 3-year scraped flux
df_proton = pd.read_csv('3yearscraped/proton', sep='\t', header=None, names=['energy', 'flux'])
df_helium = pd.read_csv('3yearscraped/helium', sep='\t', header=None, names=['energy', 'flux'])
df_light = pd.DataFrame.from_dict({'energy': df_proton.energy,
'flux': df_proton.flux + df_helium.flux})
df_oxygen = pd.read_csv('3yearscraped/oxygen', sep='\t', header=None, names=['energy', 'flux'])
df_iron = pd.read_csv('3yearscraped/iron', sep='\t', header=None, names=['energy', 'flux'])
df_heavy = pd.DataFrame.from_dict({'energy': df_oxygen.energy,
'flux': df_oxygen.flux + df_iron.flux})
# if comp_class:
# ax.plot(np.log10(df_light.energy), df_light.flux, label='3 yr light',
# marker='.', ls=':')
# ax.plot(np.log10(df_heavy.energy), df_heavy.flux, label='3 yr heavy',
# marker='.', ls=':')
# ax.plot(np.log10(df_heavy.energy), df_heavy.flux+df_light.flux, label='3 yr total',
# marker='.', ls=':')
# else:
# ax.plot(np.log10(df_proton.energy), df_proton.flux, label='3 yr proton',
# marker='.', ls=':')
# ax.plot(np.log10(df_helium.energy), df_helium.flux, label='3 yr helium',
# marker='.', ls=':', color=color_dict['He'])
# ax.plot(np.log10(df_oxygen.energy), df_oxygen.flux, label='3 yr oxygen',
# marker='.', ls=':', color=color_dict['O'])
# ax.plot(np.log10(df_iron.energy), df_iron.flux, label='3 yr iron',
# marker='.', ls=':', color=color_dict['Fe'])
# ax.plot(np.log10(df_iron.energy), df_proton.flux+df_helium.flux+df_oxygen.flux+df_iron.flux, label='3 yr total',
# marker='.', ls=':', color='C2')
leg = plt.legend(loc='upper center', frameon=False,
bbox_to_anchor=(0.5, # horizontal
1.15),# vertical
ncol=len(comp_list)+1, fancybox=False)
# set the linewidth of each legend object
for legobj in leg.legendHandles:
legobj.set_linewidth(3.0)
plt.savefig('/home/jbourbeau/public_html/figures/spectrum.png')
plt.show()
In [24]:
if not comp_class:
# Add 3-year scraped flux
df_proton = pd.read_csv('3yearscraped/proton', sep='\t', header=None, names=['energy', 'flux'])
df_helium = pd.read_csv('3yearscraped/helium', sep='\t', header=None, names=['energy', 'flux'])
df_oxygen = pd.read_csv('3yearscraped/oxygen', sep='\t', header=None, names=['energy', 'flux'])
df_iron = pd.read_csv('3yearscraped/iron', sep='\t', header=None, names=['energy', 'flux'])
# Plot fraction of events vs energy
fig, axarr = plt.subplots(2, 2, figsize=(8, 6))
for composition, ax in zip(comp_list + ['total'], axarr.flatten()):
# Calculate dN/dE
y = num_reco_energy[composition]/energybins.energy_bin_widths
y_err = num_reco_energy_err[composition]/energybins.energy_bin_widths
# Add effective area
y, y_err = comp.analysis.ratio_error(y, y_err, eff_area, eff_area_error)
# Add solid angle
y = y / solid_angle
y_err = y_err / solid_angle
# Add time duration
y = y / livetime
y_err = y / livetime
y = energybins.energy_midpoints**2.7 * y
y_err = energybins.energy_midpoints**2.7 * y_err
plotting.plot_steps(energybins.log_energy_midpoints, y, y_err, ax, color_dict[composition], composition)
# Load 3-year flux
df_3yr = pd.read_csv('3yearscraped/{}'.format(composition), sep='\t',
header=None, names=['energy', 'flux'])
ax.plot(np.log10(df_3yr.energy), df_3yr.flux, label='3 yr {}'.format(composition),
marker='.', ls=':', color=color_dict[composition])
ax.set_yscale("log", nonposy='clip')
# ax.set_xscale("log", nonposy='clip')
ax.set_xlabel('$\log_{10}(E_{\mathrm{reco}}/\mathrm{GeV})$')
ax.set_ylabel('$\mathrm{E}^{2.7} \\frac{\mathrm{dN}}{\mathrm{dE dA d\Omega dt}} \ [\mathrm{GeV}^{1.7} \mathrm{m}^{-2} \mathrm{sr}^{-1} \mathrm{s}^{-1}]$')
ax.set_xlim([6.3, 8])
ax.set_ylim([10**3, 10**5])
ax.grid(linestyle='dotted', which="both")
ax.legend()
plt.savefig('/home/jbourbeau/public_html/figures/spectrum.png')
plt.show()
[ back to top ]
In [56]:
bin_midpoints, _, counts, counts_err = comp.get1d('/home/jbourbeau/PyUnfold/unfolded_output_h3a.root', 'NC', 'Unf_ks_ACM/bin0')
In [57]:
light_counts = counts[::2]
heavy_counts = counts[1::2]
light_counts, heavy_counts
Out[57]:
In [58]:
fig, ax = plt.subplots()
for composition in comp_list + ['total']:
y, y_err = counts_to_flux(num_reco_energy[composition], num_reco_energy_err[composition], livetime=livetime)
plotting.plot_steps(energybins.log_energy_midpoints, y, y_err, ax, color_dict[composition], composition)
h3a_light_flux, h3a_flux_err = counts_to_flux(light_counts, np.sqrt(light_counts), livetime=livetime)
h3a_heavy_flux, h3a_flux_err = counts_to_flux(heavy_counts, np.sqrt(heavy_counts), livetime=livetime)
ax.plot(energybins.log_energy_midpoints, h3a_light_flux, ls=':', label='h3a light unfolded')
ax.plot(energybins.log_energy_midpoints, h3a_heavy_flux, ls=':', label='h3a heavy unfolded')
ax.set_yscale("log", nonposy='clip')
plt.xlabel('$\log_{10}(E_{\mathrm{reco}}/\mathrm{GeV})$')
# ax.set_ylabel('$\mathrm{E}^{2.7} \\frac{\mathrm{dN}}{\mathrm{dE dA d\Omega dt}} \ [\mathrm{GeV}^{1.7} \mathrm{m}^{-2} \mathrm{sr}^{-1} \mathrm{s}^{-1}]$')
ax.set_ylabel('$\mathrm{E}^{2.7} \ J(E) \ [\mathrm{GeV}^{1.7} \mathrm{m}^{-2} \mathrm{sr}^{-1} \mathrm{s}^{-1}]$')
ax.set_xlim([6.4, 9.0])
ax.set_ylim([10**2, 10**5])
ax.grid(linestyle='dotted', which="both")
leg = plt.legend(loc='upper center', frameon=False,
bbox_to_anchor=(0.5, # horizontal
1.15),# vertical
ncol=len(comp_list)+1, fancybox=False)
# set the linewidth of each legend object
for legobj in leg.legendHandles:
legobj.set_linewidth(3.0)
plt.savefig('/home/jbourbeau/public_html/figures/spectrum-unfolded.png')
plt.show()
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: