notebook.community

Edit and run 





In [1]:



    


%load_ext watermark
%watermark -a 'Author: James Bourbeau' -u -d -v -p numpy,matplotlib,scipy,pandas,sklearn,mlxtend



    


















    






Author: James Bourbeau 
last updated: 2018-08-13 

CPython 2.7.13
IPython 5.7.0

numpy 1.14.5
matplotlib 2.2.2
scipy 1.1.0
pandas 0.23.1
sklearn 0.19.1
mlxtend 0.12.0

Cosmic-ray composition unfolded spectrum

Table of contents

  1. Define analysis free parameters
  2. Data preprocessing
  3. Fitting random forest
  4. Fraction correctly identified
  5. Spectrum
  6. Unfolding




In [2]:



    


from __future__ import division, print_function
import os
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
from matplotlib.animation import FuncAnimation
import seaborn as sns
import matplotlib as mpl
import pyprind

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 comptools as comp

sns.set_context(context='paper', font_scale=1.5)
# color_dict allows for a consistent color-coding for each composition
color_dict = comp.get_color_dict()

%matplotlib inline



    


















    






/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/seaborn/apionly.py:6: UserWarning: As seaborn no longer sets a default style on import, the seaborn.apionly module is deprecated. It will be removed in a future version.
  warnings.warn(msg, UserWarning)

Define analysis free parameters

[ back to top ]

Whether or not to train on 'light' and 'heavy' composition classes, or the individual compositions





In [3]:



    


# config = 'IC79.2010'
config = 'IC86.2012'
num_groups = 2
comp_list = comp.get_comp_list(num_groups=num_groups)



    











In [4]:



    


comp_list



    


















    
Out[4]:








['light', 'heavy']

Get composition classifier pipeline

Define energy binning for this analysis





In [5]:



    


energybins = comp.get_energybins(config)



    











In [6]:



    


log_energy_min = energybins.log_energy_min
log_energy_max = energybins.log_energy_max

Data preprocessing

[ back to top ]

  1. Load simulation/data dataframe and apply specified quality cuts
  2. Extract desired features from dataframe
  3. Get separate testing and training datasets




In [7]:



    


df_sim_train, df_sim_test = comp.load_sim(config=config,
                                          log_energy_min=log_energy_min,
                                          log_energy_max=log_energy_max)
feature_list, feature_labels = comp.get_training_features()



    


















    






/home/jbourbeau/cr-composition/.virtualenv/lib/python2.7/site-packages/dask/base.py:835: UserWarning: The get= keyword has been deprecated. Please use the scheduler= keyword instead with the name of the desired scheduler like 'threads' or 'processes'
  warnings.warn("The get= keyword has been deprecated. "

















In [8]:



    


feature_list



    


















    
Out[8]:








['lap_cos_zenith', 'log_s125', 'log_dEdX']

Spectrum

[ back to top ]





In [9]:



    


unfolding_dir  = os.path.join(comp.paths.comp_data_dir, config, 'unfolding')
print(unfolding_dir)



    


















    






/data/user/jbourbeau/composition/IC86.2012/unfolding

















In [10]:



    


# Solid angle
theta_max = 40 if config == 'IC79.2010' else 65
# solid_angle = 2*np.pi*(np.cos(df_sim_train['lap_zenith'].min())-np.cos(df_sim_train['lap_zenith'].max()))
solid_angle = np.pi*np.sin(np.deg2rad(theta_max))**2
solid_angle



    


















    
Out[10]:








2.580484742999784

















In [11]:



    


livetime, livetime_err = comp.get_detector_livetime(config=config)
livetime, livetime_err



    


















    
Out[11]:








(28442248.920955755, 2699.132817731916)

















In [12]:



    


# Get simulation thrown areas for each energy bin
thrown_radii = comp.simfunctions.get_sim_thrown_radius(energybins.log_energy_midpoints)
thrown_areas = np.pi * thrown_radii**2
thrown_area = thrown_areas.max()
print('thrown_area = {}'.format(thrown_area))



    


















    






thrown_area = 9079202.76887

















In [13]:



    


geom_factor = (np.cos(np.deg2rad(theta_max)) + df_sim_train.lap_cos_zenith.min()) / 2
geom_factor



    


















    
Out[13]:








0.6173882793393718

Plot unfolded flux at each iteration





In [14]:



    


# model_name = 'Jeffreys'
model_name = 'H4a'
# model_name = 'Polygonato'



    











In [15]:



    


df_file = os.path.join(unfolding_dir,
                       'pyunfold_output_{}-groups.hdf'.format(num_groups))
df_unfolding = pd.read_hdf(df_file, model_name, mode='r')



    











In [16]:



    


df_unfolding



    


















    
Out[16]:











  
    

      
      stat_err
      sys_err
      ts_iter
      ts_stopping
      unfolded
    

  
  
    

      0
      [119755.348903795, 77650.72880662314, 58050.19...
      [18066643.41515722, 10918005.246765772, 741786...
      0.116809
      0.005
      [220613064.0903719, 126759425.37446691, 115415...
    

    

      1
      [144457.249105327, 109375.82728253765, 67938.3...
      [22236272.524407927, 15297361.72887334, 863039...
      0.018127
      0.005
      [213823943.32869592, 128581892.9095077, 108272...
    

    

      2
      [167639.80108754037, 138996.986581983, 80084.0...
      [25953555.089152317, 19618213.72094939, 100464...
      0.004339
      0.005
      [212886685.3884521, 133397361.82888678, 104073...
    

  




















In [17]:



    


fig, ax = plt.subplots()

def update(i):
    ax.clear()
    idx_iter, df_iter = i
    color_dict_iter = {}
    color_dict_iter['light'] = sns.color_palette('Blues', len(df_unfolding.index)).as_hex()[idx_iter]
    color_dict_iter['heavy'] = sns.color_palette('Oranges', len(df_unfolding.index)).as_hex()[idx_iter]
    color_dict_iter['total'] = sns.color_palette('Greens', len(df_unfolding.index)).as_hex()[idx_iter]

    counts, counts_sys_err, counts_stat_err = {}, {}, {}
    for idx, composition in enumerate(comp_list):
        counts[composition] = df_iter['n_c'][idx::num_groups]
        counts_sys_err[composition] = df_iter['sys_err'][idx::num_groups]
        counts_stat_err[composition] = df_iter['stat_err'][idx::num_groups]
        
    for idx, composition in enumerate(comp_list):
        if idx == 0:
            counts['total'] = np.zeros_like(counts[composition])
            counts_sys_err['total'] = np.zeros_like(counts[composition])
            counts_stat_err['total'] = np.zeros_like(counts[composition])
        counts['total'] += counts[composition]
        counts_sys_err['total'] += counts_sys_err[composition]**2
        counts_stat_err['total'] += counts_stat_err[composition]**2
    counts_sys_err['total'] = np.sqrt(counts_sys_err['total'])
    counts_stat_err['total'] = np.sqrt(counts_stat_err['total'])    
    
    for composition in comp_list + ['total']:    
        flux, flux_err_sys = comp.analysis.get_flux(counts[composition], counts_sys_err[composition],
                                                     energybins=energybins.energy_bins,
                                                     eff_area=thrown_area,
                                                     livetime=livetime, livetime_err=livetime_err, 
                                                     solid_angle=solid_angle)
        flux, flux_err_stat = comp.analysis.get_flux(counts[composition], counts_stat_err[composition],
                                                     energybins=energybins.energy_bins,
                                                     eff_area=thrown_area,
                                                     livetime=livetime, livetime_err=livetime_err, 
                                                     solid_angle=solid_angle)
        
        plotting.plot_steps(energybins.log_energy_bins, flux, yerr=flux_err_sys,
                            ax=ax, alpha=0.4, fillalpha=0.4,  
                            color=color_dict[composition])
        ax.errorbar(energybins.log_energy_midpoints, flux, yerr=flux_err_stat,  
                    color=color_dict[composition], ls='None', marker='.', 
                    label=composition)
        

    ax.set_yscale("log", nonposy='clip')
    ax.set_xlabel('$\mathrm{\log_{10}(E_{reco}/GeV)}$')
    ax.set_ylabel('$\mathrm{ E^{2.7} \ J(E) \ [GeV^{1.7} m^{-2} sr^{-1} s^{-1}]}$')
    ax.set_title('Iteration: {}'.format(idx_iter+1))
    ax.set_xlim(6.4, 7.8)
    ax.set_ylim([1e3, 7e4])
    ax.grid(linestyle='dotted', which="both")
    ax.legend()
    
    return ax

anim = FuncAnimation(fig, update, frames=list(df_unfolding.iterrows()),
                     interval=1250)
iter_unfold_outfile = os.path.join(comp.paths.figures_dir, 'unfolding', config, 'iterations', 
                                   model_name, '{}-groups'.format(num_groups),
                                   'flux_iter_{}-prior.gif'.format(model_name))
comp.check_output_dir(iter_unfold_outfile)
anim.save(iter_unfold_outfile, dpi=300, writer='imagemagick')



    


















    






/home/jbourbeau/.virtualenvs/composition/lib/python2.7/site-packages/matplotlib/figure.py:1743: UserWarning: This figure includes Axes that are not compatible with tight_layout, so its results might be incorrect.
  warnings.warn("This figure includes Axes that are not "










    
























In [32]:



    


# for idx_iter, df_iter in df_unfolding.iterrows():

#     fig, ax = plt.subplots()
#     color_dict_iter = {}
#     color_dict_iter['light'] = sns.color_palette('Blues', len(df_unfolding.index)).as_hex()[idx_iter]
#     color_dict_iter['heavy'] = sns.color_palette('Oranges', len(df_unfolding.index)).as_hex()[idx_iter]
#     color_dict_iter['total'] = sns.color_palette('Greens', len(df_unfolding.index)).as_hex()[idx_iter]

#     counts, counts_sys_err, counts_stat_err = {}, {}, {}
#     for idx, composition in enumerate(comp_list):
#         counts[composition] = df_iter['n_c'][idx::num_groups]
#         counts_sys_err[composition] = df_iter['sys_err'][idx::num_groups]
#         counts_stat_err[composition] = df_iter['stat_err'][idx::num_groups]
        
#     for idx, composition in enumerate(comp_list):
#         if idx == 0:
#             counts['total'] = np.zeros_like(counts[composition])
#             counts_sys_err['total'] = np.zeros_like(counts[composition])
#             counts_stat_err['total'] = np.zeros_like(counts[composition])
#         counts['total'] += counts[composition]
#         counts_sys_err['total'] += counts_sys_err[composition]**2
#         counts_stat_err['total'] += counts_stat_err[composition]**2
#     counts_sys_err['total'] = np.sqrt(counts_sys_err['total'])
#     counts_stat_err['total'] = np.sqrt(counts_stat_err['total'])
    
# #     counts['total'] = counts['light'] + counts['heavy']
# #     counts_sys_err['total'] = np.sqrt(counts_sys_err['light']**2 + counts_sys_err['heavy']**2)
# #     counts_stat_err['total'] = np.sqrt(counts_stat_err['light']**2 + counts_stat_err['heavy']**2)
    
    
#     for composition in comp_list + ['total']:    
#         flux, flux_err_sys = comp.analysis.get_flux(counts[composition], counts_sys_err[composition],
#                                                      energybins=energybins.energy_bins,
#                                                      eff_area=thrown_area,
#                                                      livetime=livetime, livetime_err=livetime_err, 
#                                                      solid_angle=solid_angle)
#         flux, flux_err_stat = comp.analysis.get_flux(counts[composition], counts_stat_err[composition],
#                                                      energybins=energybins.energy_bins,
#                                                      eff_area=thrown_area,
#                                                      livetime=livetime, livetime_err=livetime_err, 
#                                                      solid_angle=solid_angle)
        
#         plotting.plot_steps(energybins.log_energy_bins, flux, yerr=flux_err_sys,
#                             ax=ax, alpha=0.4, fillalpha=0.4,  
#                             color=color_dict[composition])
#         ax.errorbar(energybins.log_energy_midpoints, flux, yerr=flux_err_stat,  
#                     color=color_dict[composition], ls='None', marker='.', 
#                     label=composition)
        

#     ax.set_yscale("log", nonposy='clip')
#     ax.set_xlabel('$\mathrm{\log_{10}(E_{reco}/GeV)}$')
#     ax.set_ylabel('$\mathrm{ E^{2.7} \ J(E) \ [GeV^{1.7} m^{-2} sr^{-1} s^{-1}]}$')
#     ax.set_title('Iteration: {}'.format(idx_iter+1))
# #     stopping_cond_str = 'Stopping condition: KS p-val $< {}$'.format('0.{}'.format(ks_pval_str.split('p')[-1]))
# #     bbox_props = dict(boxstyle='round', ec='gray', lw=0.8, fc='white', alpha=0.5)
# #     ax.text(7.2, 7e4, stopping_cond_str, bbox=bbox_props)

#     ax.set_xlim(6.4, 7.8)
#     ax.set_ylim([1e3, 7e4])
#     ax.grid(linestyle='dotted', which="both")
#     ax.legend()
# #     ax.legend(loc='center left', bbox_to_anchor=(1, 0.5), frameon=False)
#     iter_unfold_outfile = os.path.join(comp.paths.figures_dir, 'unfolding', config, 'iterations', 
#                                        model_name, '{}-groups'.format(num_groups),
#                                        'flux_iter_{}.png'.format(idx_iter))
#     comp.check_output_dir(iter_unfold_outfile)
#     plt.savefig(iter_unfold_outfile)
#     plt.show()

Plot comparison of final unfolded flux for different priors





In [28]:



    


# model_names = ['Jeffreys']
model_names = ['H4a']
# model_names = ['Jeffreys', 'H3a', 'H4a']
# model_names = ['Jeffreys', 'h3a', 'Hoerandel5']

Load 3-year analysis flux





In [29]:



    


flux_data_3yr_array = np.load(os.path.join(comp.paths.comp_data_dir, 'fraction_spectrum_fit_data.npy'))



    











In [30]:



    


flux_data_3yr_format = ['log_energy_midpoints',
                        'PPlus_flux', 'PPlus_err_down', 'PPlus_err_up', 'PPlus_err_down2', 'PPlus_err_up2',
                        'He4Nucleus_flux', 'He4Nucleus_err_down', 'He4Nucleus_err_up', 'He4Nucleus_err_down2', 'He4Nucleus_err_up2',
                        'O16Nucleus_flux', 'O16Nucleus_err_down', 'O16Nucleus_err_up', 'O16Nucleus_err_down2', 'O16Nucleus_err_up2',
                        'Fe56Nucleus_flux', 'Fe56Nucleus_err_down', 'Fe56Nucleus_err_up', 'Fe56Nucleus_err_down2', 'Fe56Nucleus_err_up2']



    











In [31]:



    


flux_data_3yr = {key: value.astype(float) for key, value in zip(flux_data_3yr_format, flux_data_3yr_array)}
flux_comp_3yr = {}
MC_comp_list = ['PPlus', 'He4Nucleus', 'O16Nucleus', 'Fe56Nucleus']
for MC_comp in MC_comp_list:
    comp_group = comp.composition_encoding.composition_group_labels([MC_comp], num_groups=num_groups)[0]
    print(comp_group)
    if '{}_flux'.format(comp_group) not in flux_comp_3yr:
        flux_comp_3yr['{}_flux'.format(comp_group)] = flux_data_3yr['{}_flux'.format(MC_comp)]
        flux_comp_3yr['{}_err_down'.format(comp_group)] = flux_data_3yr['{}_err_down'.format(MC_comp)]**2
        flux_comp_3yr['{}_err_up'.format(comp_group)] = flux_data_3yr['{}_err_up'.format(MC_comp)]**2
    else:
        flux_comp_3yr['{}_flux'.format(comp_group)] += flux_data_3yr['{}_flux'.format(MC_comp)]
        flux_comp_3yr['{}_err_down'.format(comp_group)] += flux_data_3yr['{}_err_down'.format(MC_comp)]**2
        flux_comp_3yr['{}_err_up'.format(comp_group)] += flux_data_3yr['{}_err_up'.format(MC_comp)]**2

for composition in comp_list:
    flux_comp_3yr['{}_err_down'.format(composition)] = np.sqrt(flux_comp_3yr['{}_err_down'.format(composition)])
    flux_comp_3yr['{}_err_up'.format(composition)] = np.sqrt(flux_comp_3yr['{}_err_up'.format(composition)])

flux_comp_3yr['total_flux'] = np.sum((flux_comp_3yr['{}_flux'.format(composition)] for composition in comp_list),
                                     axis=0)
flux_comp_3yr['total_err_down'] = np.sqrt(np.sum((flux_comp_3yr['{}_err_down'.format(composition)]**2 for composition in comp_list),
                                     axis=0))
flux_comp_3yr['total_err_up'] = np.sqrt(np.sum((flux_comp_3yr['{}_err_up'.format(composition)]**2 for composition in comp_list),
                                     axis=0))



    


















    






light
light
heavy
heavy

















In [32]:



    


plot_3yr = False



    











In [34]:



    


fig, ax = plt.subplots()
marker_dict = {model: marker for model, marker in zip(model_names, '.^*o')}
color_dict_model = {}
color_dict_model['light'] = sns.color_palette('Blues', len(model_names)+1).as_hex()[::-1]
color_dict_model['intermediate'] = sns.color_palette('Reds', len(model_names)+1).as_hex()[::-1]
color_dict_model['heavy'] = sns.color_palette('Oranges', len(model_names)+1).as_hex()[::-1]

color_dict_model['PPlus'] = sns.color_palette('Blues', len(model_names)+1).as_hex()[::-1]
color_dict_model['O16Nucleus'] = sns.color_palette('Reds', len(model_names)+1).as_hex()[::-1]
color_dict_model['He4Nucleus'] = sns.color_palette('Purples', len(model_names)+1).as_hex()[::-1]
color_dict_model['Fe56Nucleus'] = sns.color_palette('Oranges', len(model_names)+1).as_hex()[::-1]

color_dict_model['total'] = sns.color_palette('Greens', len(model_names)+1).as_hex()[::-1]
for idx_model, model_name in enumerate(model_names):
    
    df_file = os.path.join(unfolding_dir, 
                           'pyunfold_output_{}-groups.hdf'.format(num_groups))
    df_unfolding = pd.read_hdf(df_file, model_name, mode='r')
    df_final_iter = df_unfolding.iloc[-1]
    
    counts, counts_sys_err, counts_stat_err = {}, {}, {}
    for idx, composition in enumerate(comp_list):
        counts[composition] = df_final_iter['unfolded'][idx::num_groups]
        counts_sys_err[composition] = df_final_iter['sys_err'][idx::num_groups]
        counts_stat_err[composition] = df_final_iter['stat_err'][idx::num_groups]
        
    for idx, composition in enumerate(comp_list):
        if idx == 0:
            counts['total'] = np.zeros_like(counts[composition])
            counts_sys_err['total'] = np.zeros_like(counts[composition])
            counts_stat_err['total'] = np.zeros_like(counts[composition])
        counts['total'] += counts[composition]
        counts_sys_err['total'] += counts_sys_err[composition]**2
        counts_stat_err['total'] += counts_stat_err[composition]**2
    counts_sys_err['total'] = np.sqrt(counts_sys_err['total'])
    counts_stat_err['total'] = np.sqrt(counts_stat_err['total'])
    
    for composition in comp_list + ['total']:    
        flux, flux_err_sys = comp.get_flux(counts[composition], counts_sys_err[composition],
                                                     energybins=energybins.energy_bins,
                                                     eff_area=thrown_area,
#                                                      eff_area=thrown_area*geom_factor,
                                                     livetime=livetime, livetime_err=livetime_err, 
                                                     solid_angle=solid_angle)
        flux, flux_err_stat = comp.get_flux(counts[composition], counts_stat_err[composition],
                                                     energybins=energybins.energy_bins,
                                                     eff_area=thrown_area,
#                                                      eff_area=thrown_area*geom_factor, 
                                                     livetime=livetime, livetime_err=livetime_err, 
                                                     solid_angle=solid_angle)
        
        comp.plot_steps(energybins.log_energy_bins, flux, yerr=flux_err_sys,
                            ax=ax, alpha=0.7, fillalpha=0.5,  
                            color=color_dict_model[composition][idx_model])

        ax.errorbar(energybins.log_energy_midpoints, flux, yerr=flux_err_stat,  
                    color=color_dict_model[composition][idx_model], ls='None', marker=marker_dict[model_name], 
                    label=composition, alpha=0.8)
#                     label=composition + ' ({})'.format(model_name), alpha=0.8)


# Add 3-year composition analysis flux for comparison 
if plot_3yr:
    for composition in comp_list + ['total']:
        ax.errorbar(flux_data_3yr['log_energy_midpoints'], 
                    flux_comp_3yr['{}_flux'.format(composition)],
                    yerr=[flux_comp_3yr['{}_err_down'.format(composition)],
                          flux_comp_3yr['{}_err_up'.format(composition)]],  
                    ls='None', marker='*', color=color_dict[composition],
    #                 label='3-year ({})'.format(composition),
                    alpha=0.75)
    
    
ax.set_yscale("log", nonposy='clip')
ax.set_xlabel('$\mathrm{\log_{10}(E_{reco}/GeV)}$')
ax.set_ylabel('$\mathrm{ E^{2.7} \ J(E) \ [GeV^{1.7} m^{-2} sr^{-1} s^{-1}]}$')

# ax.set_xlim(6.4, 9.0)
ax.set_xlim(6.4, 7.8)
ax.set_ylim([1e3, 7e4])

ax.grid(linestyle='dotted', which="both")
ax.legend(loc='lower left', ncol=len(model_names), fontsize=10)
# leg = plt.legend(loc='upper center', frameon=False,
#           bbox_to_anchor=(0.5,  # horizontal
#                           1.15),# vertical 
#           ncol=len(comp_list)+1, fancybox=False)

if plot_3yr:
    flux_outfile = os.path.join(comp.paths.figures_dir, 'unfolding', config, 
                                'flux_{}-groups-3yr-comparison.png'.format(num_groups))
else:
    flux_outfile = os.path.join(comp.paths.figures_dir, 'unfolding', config, 
                                'flux_{}-groups.png'.format(num_groups))
comp.check_output_dir(flux_outfile)
plt.savefig(flux_outfile)
plt.show()



    


















    
























In [21]:



    


model_flux = comp.model_flux('H4a', energy=energybins.energy_midpoints,
                                                            num_groups=num_groups)
model_flux.head()



    


















    
Out[21]:











  
    

      
      flux_PPlus
      flux_He4Nucleus
      flux_O16Nucleus
      flux_Fe56Nucleus
      flux_total
    

  
  
    

      0
      2.952013e-13
      4.381007e-13
      2.391451e-13
      1.631829e-13
      1.135630e-12
    

    

      1
      1.495900e-13
      2.325688e-13
      1.303989e-13
      8.947391e-14
      6.020317e-13
    

    

      2
      7.458560e-14
      1.221825e-13
      7.098235e-14
      4.903990e-14
      3.167904e-13
    

    

      3
      3.649089e-14
      6.336501e-14
      3.855464e-14
      2.686341e-14
      1.652739e-13
    

    

      4
      1.747554e-14
      3.234035e-14
      2.088294e-14
      1.470413e-14
      8.540295e-14
    

  




















In [38]:



    


if num_groups == 4:
    fig, axarr = plt.subplots(2, 2, sharex=True, sharey=True)
    models = ['Jeffreys',
#               'H4a',
             ]
    for model_name, marker, ax in zip(models, '.^*X', axarr.flat):
    #     model_name = 'H4a'
    #     model_name = 'Jeffreys'
        df_file = os.path.join(unfolding_dir, 
                               'pyunfold_output_{}-groups.hdf'.format(num_groups))
        df_unfolding = pd.read_hdf(df_file, model_name, mode='r')
    #     print(df_unfolding)
        df_final_iter = df_unfolding.iloc[-1]

        counts, counts_sys_err, counts_stat_err = {}, {}, {}
        for idx, composition in enumerate(comp_list):
            counts[composition] = df_final_iter['n_c'][idx::num_groups]
            counts_sys_err[composition] = df_final_iter['sys_err'][idx::num_groups]
            counts_stat_err[composition] = df_final_iter['stat_err'][idx::num_groups]

        for idx, composition in enumerate(comp_list):
            if idx == 0:
                counts['total'] = np.zeros_like(counts[composition])
                counts_sys_err['total'] = np.zeros_like(counts[composition])
                counts_stat_err['total'] = np.zeros_like(counts[composition])
            counts['total'] += counts[composition]
            counts_sys_err['total'] += counts_sys_err[composition]**2
            counts_stat_err['total'] += counts_stat_err[composition]**2
        counts_sys_err['total'] = np.sqrt(counts_sys_err['total'])
        counts_stat_err['total'] = np.sqrt(counts_stat_err['total'])


        flux_total, flux_total_err_sys = comp.get_flux(counts['total'], counts_sys_err['total'],
                                                         energybins=energybins.energy_bins,
                                                         eff_area=thrown_area,
                                                         livetime=livetime, livetime_err=livetime_err, 
                                                         solid_angle=solid_angle)
        flux_total, flux_total_err_stat = comp.get_flux(counts['total'], counts_stat_err['total'],
                                                     energybins=energybins.energy_bins,
                                                     eff_area=thrown_area,
                                                     livetime=livetime, livetime_err=livetime_err, 
                                                     solid_angle=solid_angle)


        for composition, ax in zip(comp_list, axarr.flatten()):

            flux, flux_err_sys = comp.get_flux(counts[composition], counts_sys_err[composition],
                                                         energybins=energybins.energy_bins,
                                                         eff_area=thrown_area,
                                                         livetime=livetime, livetime_err=livetime_err, 
                                                         solid_angle=solid_angle)
            flux, flux_err_stat = comp.get_flux(counts[composition], counts_stat_err[composition],
                                                         energybins=energybins.energy_bins,
                                                         eff_area=thrown_area,
                                                         livetime=livetime, livetime_err=livetime_err, 
                                                         solid_angle=solid_angle)

            comp.plot_steps(energybins.log_energy_bins, flux, yerr=flux_err_sys,
                                ax=ax, alpha=0.4, fillalpha=0.4, lw=0.8, 
                                color=color_dict[composition])

            ax.errorbar(energybins.log_energy_midpoints, flux, yerr=flux_err_stat,  
                        color=color_dict[composition], ls='None', marker=marker, 
                        label=composition + ' ({})'.format(model_name), alpha=0.8)

            # Plot total flux
            comp.plot_steps(energybins.log_energy_bins, flux_total, yerr=flux_total_err_sys,
                                ax=ax, alpha=0.4, fillalpha=0.4,  
                                color=color_dict['total'])

            ax.errorbar(energybins.log_energy_midpoints, flux_total, yerr=flux_total_err_stat,  
                        color=color_dict['total'], ls='None', marker=marker, 
    #                     label='total',
                        alpha=0.8)



        # Add 3-year composition analysis flux for comparison 
            if plot_3yr:
                ax.errorbar(flux_data_3yr['log_energy_midpoints'], 
                            flux_comp_3yr['{}_flux'.format(composition)],
                            yerr=[flux_comp_3yr['{}_err_down'.format(composition)],
                                  flux_comp_3yr['{}_err_up'.format(composition)]],  
                            ls='None', marker='*', color=color_dict[composition],
                            label='{} (3-year)'.format(composition),
                            alpha=0.75)
                ax.errorbar(flux_data_3yr['log_energy_midpoints'], 
                            flux_comp_3yr['total_flux'],
                            yerr=[flux_comp_3yr['total_err_down'],
                                  flux_comp_3yr['total_err_up']],  
                            ls='None', marker='*', color=color_dict['total'],
                            alpha=0.75)


            ax.set_yscale("log", nonposy='clip')
            ax.set_xlim(6.4, 7.8)
            ax.set_ylim([1e3, 7e4])

            ax.grid(linestyle='dotted', which="both")
            ax.legend(loc='lower left', ncol=1, fontsize=10)
        
    fig.text(0.5, -0.01, '$\mathrm{\log_{10}(E/GeV)}$', ha='center')
    fig.text(-0.01, 0.5, '$\mathrm{ E^{2.7} \ J(E) \ [GeV^{1.7} m^{-2} sr^{-1} s^{-1}]}$',
             va='center', rotation='vertical')
    
#     flux_outfile = os.path.join(comp.paths.figures_dir, 'unfolding', config, 
#                                 'flux_grid_{}-groups_{}-prior.png'.format(num_groups, model_name))
    flux_outfile = os.path.join(comp.paths.figures_dir, 'unfolding', config, 
                                'flux_grid_{}-groups.png'.format(num_groups))
    comp.check_output_dir(flux_outfile)
    plt.savefig(flux_outfile, dpi=1000)
    plt.show()



    


















    






/home/jbourbeau/.virtualenvs/composition/lib/python2.7/site-packages/matplotlib/scale.py:111: RuntimeWarning: invalid value encountered in less_equal
  out[a <= 0] = -1000










    
























In [23]:



    


fig, ax = plt.subplots()
# Add 3-year composition analysis flux for comparison 
for composition in comp_list + ['total']:
    ax.errorbar(flux_data_3yr['log_energy_midpoints'], 
                flux_comp_3yr['{}_flux'.format(composition)],
                yerr=[flux_comp_3yr['{}_err_down'.format(composition)],
                      flux_comp_3yr['{}_err_up'.format(composition)]],  
                ls='None', marker='*', color=color_dict[composition],
                label='3-year ({})'.format(composition),
                alpha=0.75)
    
ax.set_yscale("log", nonposy='clip')
ax.set_xlabel('$\mathrm{\log_{10}(E_{reco}/GeV)}$')
ax.set_ylabel('$\mathrm{ E^{2.7} \ J(E) \ [GeV^{1.7} m^{-2} sr^{-1} s^{-1}]}$')

# ax.set_xlim(6.4, 9.0)
ax.set_xlim(6.4, 7.8)
ax.set_ylim([1e3, 7e4])

ax.grid(linestyle='dotted', which="both")
ax.legend(loc='lower left', ncol=2, fontsize=10)
# ax.legend(loc='center left', bbox_to_anchor=(1, 0.5), frameon=False)cd 

# leg = plt.legend(loc='upper center', frameon=False,
#           bbox_to_anchor=(0.5,  # horizontal
#                           1.15),# vertical 
#           ncol=len(comp_list)+1, fancybox=False)
flux_3yr_outfile = os.path.join(comp.paths.figures_dir, 'unfolding', config, 
                                'flux_3yr-analysis_{}-groups.png'.format(num_groups))
comp.check_output_dir(flux_3yr_outfile)
# plt.savefig(flux_3yr_outfile)
plt.show()



    


















    
























In [ ]:



    






    











In [ ]:

Content source: jrbourbeau/cr-composition

Similar notebooks:

notebook.community | gallery | about