In [1]:
%load_ext watermark
%watermark -u -d -v -p numpy,matplotlib,scipy,pandas,sklearn,mlxtend


last updated: 2017-08-22 

CPython 2.7.13
IPython 5.3.0

numpy 1.12.1
matplotlib 2.0.2
scipy 0.19.0
pandas 0.20.1
sklearn 0.18.1
mlxtend 0.7.0

In [2]:
%matplotlib inline
from __future__ import division, print_function
from collections import defaultdict
import os
import numpy as np
from scipy import optimize
from scipy.stats import chisquare
import pandas as pd
import matplotlib.pyplot as plt
import seaborn.apionly as sns

import comptools as comp

color_dict = comp.analysis.get_color_dict()


/home/jbourbeau/.virtualenvs/composition/lib/python2.7/site-packages/sklearn/cross_validation.py:44: DeprecationWarning: This module was deprecated in version 0.18 in favor of the model_selection module into which all the refactored classes and functions are moved. Also note that the interface of the new CV iterators are different from that of this module. This module will be removed in 0.20.
  "This module will be removed in 0.20.", DeprecationWarning)

Load simulation DataFrame and apply quality cuts

[ back to top ]


In [3]:
# config = 'IC79'
config = 'IC86.2012'
df_sim = comp.load_sim(config=config, test_size=0)

In [5]:
df_sim


Out[5]:
FractionContainment_Laputop_IceTop FractionContainment_Laputop_InIce FractionContainment_MCPrimary_IceTop FractionContainment_MCPrimary_InIce IceTopMaxSignal IceTopMaxSignalInEdge IceTopMaxSignalString IceTopNeighbourMaxSignal InIce_charge_1_60 MC_azimuth ... lap_cos_zenith log_s50 log_s80 log_s125 log_s180 log_s250 log_s500 log_dEdX target is_thinned
0 0.723129 0.975129 0.730482 0.983024 18.432690 0 10 14.470111 245.350001 2.912543 ... 0.970416 1.089618 0.417746 -0.243566 -0.800776 -1.315759 -2.442802 0.527069 0 False
1 0.197626 0.803790 0.216988 0.827441 24.921978 0 79 19.548155 122.400000 2.912543 ... 0.969793 1.133598 0.410655 -0.299151 -0.895983 -1.446662 -2.649023 0.528187 0 False
2 0.758493 0.741260 0.778203 0.776185 12.654794 0 15 4.653794 115.975000 2.912543 ... 0.968690 0.941354 0.351498 -0.231935 -0.725515 -1.183174 -2.189262 0.376604 0 False
3 0.307622 0.593746 0.295832 0.578842 34.789658 0 54 8.977652 139.450000 2.912543 ... 0.967562 0.894099 0.349628 -0.190710 -0.649078 -1.075016 -2.014172 0.300163 0 False
4 0.584972 0.929127 0.571597 0.949283 16.521660 0 17 15.387228 30.375000 2.912543 ... 0.972129 0.996699 0.410424 -0.169610 -0.660412 -1.115568 -2.116376 0.444314 0 False
5 0.599457 0.825440 0.603860 0.834725 915.565613 0 64 266.637909 72.200000 2.912543 ... 0.970588 0.965718 0.472972 -0.018253 -0.436491 -0.826276 -1.689150 0.669241 0 False
6 0.518328 0.687292 0.482674 0.619027 13.730797 0 24 6.145182 53.975000 2.912543 ... 0.969633 0.875428 0.315497 -0.239522 -0.709884 -1.146627 -2.108583 0.291250 0 False
7 0.707230 0.219441 0.738094 0.231852 618.074707 0 52 397.986542 35.950000 2.912543 ... 0.972208 1.340168 0.541520 -0.240171 -0.895737 -1.499329 -2.813338 0.148202 0 False
8 0.704877 0.789085 0.701785 0.782489 8.787918 0 15 4.478899 115.325000 2.912543 ... 0.969743 0.988941 0.374217 -0.232830 -0.745702 -1.220742 -2.263505 0.494186 0 False
9 0.672916 0.596048 0.675481 0.649017 28.609421 0 29 12.178488 180.050000 5.673565 ... 0.932002 1.058991 0.515847 -0.023232 -0.480570 -0.905580 -1.842779 0.733548 0 False
10 0.537961 0.760892 0.543832 0.771945 101.481819 0 38 68.172859 244.275001 5.673565 ... 0.925519 1.140504 0.560367 -0.013838 -0.499876 -0.950742 -1.942496 0.758869 0 False
11 0.587308 0.706469 0.556582 0.741011 36.809093 0 19 21.983028 112.250000 5.673565 ... 0.927107 0.935091 0.408557 -0.114750 -0.559202 -0.972602 -1.885305 0.574857 0 False
12 0.676321 0.597996 0.674722 0.655864 20.867756 0 20 10.946673 165.950000 5.673565 ... 0.930699 1.007560 0.502415 -0.000583 -0.428441 -0.826892 -1.708051 0.363861 0 False
13 0.255498 0.933506 0.248732 0.951277 37.628693 0 80 21.158876 27.025000 5.673565 ... 0.926979 1.077422 0.507438 -0.057127 -0.535289 -0.979059 -1.955840 0.663036 0 False
14 0.440063 0.895023 0.457791 0.920354 45.083260 0 18 25.149473 47.575000 5.673565 ... 0.926911 1.384774 0.500941 -0.361636 -1.083291 -1.746422 -3.186059 0.399265 0 False
15 0.582348 0.771855 0.594584 0.761230 9.815289 0 19 7.491657 185.700000 5.673565 ... 0.924079 1.210193 0.524700 -0.149544 -0.717321 -1.241824 -2.388954 0.744836 0 False
16 0.860122 0.885281 0.856127 0.930048 35.029404 0 10 19.274426 160.075000 5.673565 ... 0.926482 1.246090 0.565818 -0.103470 -0.667197 -1.188052 -2.327484 1.227303 0 False
17 0.430628 0.817969 0.438451 0.821270 24.654524 0 27 7.365587 67.475000 5.673565 ... 0.926303 1.174117 0.551162 -0.063700 -0.582959 -1.063752 -2.118654 0.634211 0 False
18 0.389801 0.887318 0.378421 0.896916 19.827286 0 27 6.790532 59.850000 5.673565 ... 0.927895 1.099503 0.533604 -0.027082 -0.502074 -0.942989 -1.913746 0.510647 0 False
19 0.609966 0.591132 0.585553 0.570004 726.153137 0 64 211.173141 40.025000 1.778023 ... 0.946838 1.605879 0.861154 0.130667 -0.483064 -1.048967 -2.283451 -0.028871 0 False
20 0.830225 0.403236 0.845329 0.348111 24.713930 0 70 23.905867 30.175000 1.778023 ... 0.947264 1.205546 0.599397 0.000491 -0.505729 -0.974776 -2.004894 -0.132718 0 False
21 0.527111 0.587919 0.542497 0.578386 18.544727 0 63 18.182594 52.900000 1.778023 ... 0.951781 1.430346 0.783446 0.145847 -0.391988 -0.889517 -1.979732 0.157762 0 False
22 0.540229 0.848432 0.568934 0.800526 89.009895 0 43 45.976791 31.975000 1.778023 ... 0.947796 1.484707 0.856598 0.236841 -0.286416 -0.770812 -1.833314 0.075981 0 False
23 0.872215 0.295916 0.873934 0.281987 138.764206 0 61 87.493927 1585.625003 2.836055 ... 0.968237 2.424690 1.789795 1.163596 0.635075 0.145936 -0.926573 1.531194 0 False
24 0.780010 0.870553 0.770539 0.852982 361.826416 0 71 181.568573 1315.250002 2.836055 ... 0.969415 2.472788 1.807009 1.151484 0.599002 0.088278 -1.029779 1.519787 0 False
25 0.520361 0.366598 0.540827 0.435061 102.640343 0 28 66.634789 63.025000 6.080910 ... 0.999282 1.012832 0.416685 -0.172722 -0.671182 -1.133238 -2.148605 0.232347 0 False
26 0.601685 0.669965 0.584634 0.657250 26.275906 0 53 10.203528 50.400000 6.080910 ... 0.999663 1.160530 0.527211 -0.097493 -0.624792 -1.112829 -2.183014 0.180237 0 False
27 0.590722 0.538213 0.596413 0.578119 7.266950 0 16 6.468921 226.825000 6.080910 ... 0.999701 1.280508 0.606102 -0.057616 -0.616791 -1.133546 -2.264325 0.501576 0 False
28 0.146048 0.150852 0.154011 0.156084 61.693287 0 46 33.477386 97.975000 6.080910 ... 0.999546 1.312978 0.568528 -0.161699 -0.775217 -1.340928 -2.575007 0.054224 0 False
29 0.286743 0.349657 0.300086 0.343006 10.924444 0 34 6.257939 81.125000 6.080910 ... 0.999377 1.074479 0.441134 -0.183594 -0.710913 -1.198968 -2.269192 0.091126 0 False
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
82332 0.542809 0.561129 0.530729 0.529757 57.596844 0 24 27.887409 937.075002 2.627922 ... 0.987102 1.832492 1.249378 0.672346 0.183997 -0.268949 -1.265094 1.182762 1 False
82333 0.294916 0.623198 0.302544 0.615514 239.760452 0 25 188.959457 678.250001 2.627922 ... 0.988149 1.928603 1.321527 0.721741 0.214801 -0.254895 -1.286380 1.275560 1 False
82334 0.738466 0.703938 0.739964 0.737874 4072.527832 0 65 2644.014648 761.125001 2.627922 ... 0.988996 1.874953 1.269561 0.671375 0.165743 -0.302775 -1.331775 1.215308 1 False
82335 0.207036 0.512296 0.210388 0.495595 86.559837 0 35 65.736382 1209.775000 2.627922 ... 0.986751 1.903834 1.226126 0.559274 -0.002464 -0.521526 -1.657175 1.283511 1 False
82336 0.453064 0.456914 0.440695 0.442785 44.511902 0 63 32.275497 798.575001 2.627922 ... 0.987127 1.845008 1.276345 0.713035 0.235898 -0.206949 -1.181782 1.267335 1 False
82337 0.121348 0.470837 0.128776 0.454453 492.699646 0 35 221.639572 1149.700001 2.627922 ... 0.987238 1.892744 1.263684 0.643024 0.119029 -0.366032 -1.429937 1.377935 1 False
82338 0.635091 0.061544 0.629203 0.022597 75.334190 0 53 56.082336 1092.975000 2.627922 ... 0.987379 1.938709 1.270489 0.612644 0.058267 -0.454165 -1.575823 1.196107 1 False
82339 0.738779 0.656895 0.731486 0.648937 39.639874 0 23 30.376133 651.150000 2.627922 ... 0.987808 1.914205 1.305516 0.704198 0.196007 -0.274815 -1.308679 1.241338 1 False
82340 0.686438 0.600983 0.697426 0.594575 29.673231 0 23 21.280340 818.050001 2.627922 ... 0.987653 1.981304 1.327809 0.683946 0.140994 -0.361146 -1.461088 1.238614 1 False
82341 0.448052 0.423898 0.444762 0.447377 33.037098 0 63 16.491507 1164.200002 2.627922 ... 0.988608 1.818292 1.200671 0.590873 0.075753 -0.401313 -1.448348 1.306535 1 False
82342 0.240652 0.364578 0.238602 0.359069 243.007614 0 44 213.703964 1707.475002 2.627922 ... 0.987126 1.761964 1.172454 0.589349 0.096038 -0.361379 -1.366957 1.475373 1 False
82343 0.800449 0.529363 0.778693 0.479889 209.052917 0 29 130.362976 1349.874999 5.014167 ... 0.990901 2.207728 1.546145 0.894604 0.345377 -0.162415 -1.274284 1.550695 1 False
82344 0.574775 0.888929 0.605562 0.914142 158.047302 0 64 108.113846 2354.025004 5.014167 ... 0.991212 2.162995 1.545272 0.935377 0.420178 -0.056958 -1.104144 1.658689 1 False
82345 0.493224 0.581653 0.499767 0.547518 296.404510 0 24 196.540741 1590.375000 5.014167 ... 0.990528 2.100234 1.503943 0.914399 0.415827 -0.046330 -1.061908 1.655536 1 False
82346 0.851119 0.298224 0.850215 0.289337 53.923611 0 11 39.480839 3059.825000 5.014167 ... 0.991553 2.031591 1.460165 0.894231 0.414951 -0.029827 -1.008735 1.652591 1 False
82347 0.365299 0.695959 0.378649 0.740548 3224.828125 0 54 1236.548950 1216.550001 5.014167 ... 0.991348 2.095459 1.479481 0.871244 0.357398 -0.118518 -1.163130 1.513332 1 False
82348 0.716641 0.510959 0.716226 0.492493 129.779099 0 8 103.443375 1148.350000 5.014167 ... 0.990631 2.166213 1.554291 0.949904 0.439205 -0.033878 -1.072509 1.470731 1 False
82349 0.546523 0.919652 0.547509 0.924929 36.919479 0 62 34.910767 1792.774999 5.014167 ... 0.991208 2.028903 1.461241 0.898881 0.422521 -0.019626 -0.992984 1.527920 1 False
82350 0.532956 0.769086 0.525053 0.766322 1031.610107 0 43 416.934906 2814.575003 5.014167 ... 0.991250 1.984331 1.440527 0.900822 0.442972 0.017501 -0.920670 1.647303 1 False
82351 0.727196 0.368021 0.748807 0.381469 1071.053955 0 12 676.296936 3323.875003 5.014167 ... 0.992074 2.326395 1.661596 1.007000 0.455278 -0.054762 -1.171375 1.636114 1 False
82352 0.576680 0.838821 0.548199 0.817465 256.502533 0 53 205.766678 2176.700004 5.014167 ... 0.990802 2.157948 1.520484 0.891845 0.361330 -0.129604 -1.205903 1.709313 1 False
82353 0.265742 0.538603 0.268912 0.507577 192.648438 0 37 114.862701 1992.825001 5.014167 ... 0.990975 2.213551 1.583482 0.961865 0.437087 -0.048678 -1.114071 1.600897 1 False
82354 0.619260 0.872802 0.614582 0.854302 61.503860 0 52 46.346455 1240.900000 5.014167 ... 0.990761 2.079346 1.457895 0.844460 0.326368 -0.153375 -1.206059 1.504653 1 False
82355 0.117316 0.390192 0.124162 0.383958 704.127075 0 80 322.663452 2558.725001 5.014167 ... 0.991164 2.048976 1.477931 0.912359 0.433374 -0.011137 -0.989483 1.512289 1 False
82356 0.408937 0.121274 0.414314 0.129645 96.984573 0 18 89.660645 1697.975000 5.014167 ... 0.991376 2.061345 1.469949 0.885052 0.390278 -0.068458 -1.076817 1.541928 1 False
82357 0.387460 0.715748 0.405245 0.743389 159.834549 0 56 115.355331 1442.100002 5.014167 ... 0.991605 2.050513 1.451326 0.859031 0.358212 -0.105969 -1.125819 1.654524 1 False
82358 0.729396 0.963783 0.724129 0.953653 1500.063843 0 52 1416.501099 3454.899995 5.014167 ... 0.991167 2.280382 1.619550 0.968722 0.420077 -0.087190 -1.197952 1.678574 1 False
82359 0.574530 0.898915 0.579627 0.891866 101.676804 0 64 52.037830 2256.150000 5.014167 ... 0.990737 1.986260 1.417811 0.854704 0.377733 -0.064964 -1.039483 1.524207 1 False
82360 0.813172 0.850899 0.807022 0.868772 242.377502 0 39 242.377502 1096.275002 5.548864 ... 0.928684 2.010327 1.454495 0.903368 0.436186 0.002307 -0.953603 1.659950 1 False
82361 0.431919 0.829721 0.431718 0.850483 48.837177 0 28 29.609440 941.450001 5.548864 ... 0.929411 2.127366 1.464447 0.811637 0.261373 -0.247352 -1.361192 1.704321 1 False

82362 rows × 81 columns


In [6]:
# df_sim, cut_dict_sim = comp.load_dataframe(datatype='sim', config=config, return_cut_dict=True)
# selection_mask = np.array([True] * len(df_sim))
# # standard_cut_keys = ['IceTopQualityCuts', 'lap_InIce_containment',
# # #                 'num_hits_1_60', 'max_qfrac_1_60',
# #                 'InIceQualityCuts', 'num_hits_1_60']

# standard_cut_keys = ['passed_IceTopQualityCuts', 'FractionContainment_Laputop_InIce',
#                      'passed_InIceQualityCuts', 'num_hits_1_60']
# # for cut in ['MilliNCascAbove2', 'MilliQtotRatio', 'MilliRloglBelow2', 'StochRecoSucceeded']:
# #     standard_cut_keys += ['InIceQualityCuts_{}'.format(cut)]
    
# for key in standard_cut_keys:
#     selection_mask *= cut_dict_sim[key]
#     print(key, np.sum(selection_mask))

# df_sim = df_sim[selection_mask]

Define energy binning for this analysis


In [7]:
log_energy_bins = np.arange(5.0, 9.51, 0.05)
# log_energy_bins = np.arange(5.0, 9.51, 0.1)
energy_bins = 10**log_energy_bins
energy_midpoints = (energy_bins[1:] + energy_bins[:-1]) / 2

energy_min_fit, energy_max_fit = 5.8, 7.0
midpoints_fitmask = (energy_midpoints >= 10**energy_min_fit) & (energy_midpoints <= 10**energy_max_fit)

In [8]:
log_energy_bins


Out[8]:
array([ 5.  ,  5.05,  5.1 ,  5.15,  5.2 ,  5.25,  5.3 ,  5.35,  5.4 ,
        5.45,  5.5 ,  5.55,  5.6 ,  5.65,  5.7 ,  5.75,  5.8 ,  5.85,
        5.9 ,  5.95,  6.  ,  6.05,  6.1 ,  6.15,  6.2 ,  6.25,  6.3 ,
        6.35,  6.4 ,  6.45,  6.5 ,  6.55,  6.6 ,  6.65,  6.7 ,  6.75,
        6.8 ,  6.85,  6.9 ,  6.95,  7.  ,  7.05,  7.1 ,  7.15,  7.2 ,
        7.25,  7.3 ,  7.35,  7.4 ,  7.45,  7.5 ,  7.55,  7.6 ,  7.65,
        7.7 ,  7.75,  7.8 ,  7.85,  7.9 ,  7.95,  8.  ,  8.05,  8.1 ,
        8.15,  8.2 ,  8.25,  8.3 ,  8.35,  8.4 ,  8.45,  8.5 ,  8.55,
        8.6 ,  8.65,  8.7 ,  8.75,  8.8 ,  8.85,  8.9 ,  8.95,  9.  ,
        9.05,  9.1 ,  9.15,  9.2 ,  9.25,  9.3 ,  9.35,  9.4 ,  9.45,  9.5 ])

In [9]:
np.log10(energy_midpoints[midpoints_fitmask])


Out[9]:
array([ 5.82571916,  5.87571916,  5.92571916,  5.97571916,  6.02571916,
        6.07571916,  6.12571916,  6.17571916,  6.22571916,  6.27571916,
        6.32571916,  6.37571916,  6.42571916,  6.47571916,  6.52571916,
        6.57571916,  6.62571916,  6.67571916,  6.72571916,  6.77571916,
        6.82571916,  6.87571916,  6.92571916,  6.97571916])

Define functions to be fit to effective area


In [10]:
def constant(energy, c):
    return c

def linefit(energy, m, b):
    return m*np.log10(energy) + b

def sigmoid_flat(energy, p0, p1, p2):
    return p0 / (1 + np.exp(-p1*np.log10(energy) + p2))

def sigmoid_slant(energy, p0, p1, p2, p3):
    return (p0 + p3*np.log10(energy)) / (1 + np.exp(-p1*np.log10(energy) + p2))

In [11]:
def red_chisquared(obs, fit, sigma, n_params):
    zero_mask = sigma != 0
    return np.nansum(((obs[zero_mask] - fit[zero_mask])/sigma[zero_mask]) ** 2) / (len(obs[zero_mask]) - n_params)
#     return np.sum(((obs - fit)/sigma) ** 2) / (len(obs) - 1 - n_params)

In [12]:
np.sum(midpoints_fitmask)-3


Out[12]:
21

Calculate effective areas


In [13]:
eff_area, eff_area_error, _ = comp.calculate_effective_area_vs_energy(df_sim, energy_bins)
eff_area_light, eff_area_error_light, _ = comp.calculate_effective_area_vs_energy(df_sim[df_sim.MC_comp_class == 'light'], energy_bins)
eff_area_heavy, eff_area_error_heavy, _ = comp.calculate_effective_area_vs_energy(df_sim[df_sim.MC_comp_class == 'heavy'], energy_bins)


Calculating effective area...
Simulation set 12360: 19781 files
Simulation set 12362: 19808 files
Simulation set 12630: 19787 files
Simulation set 12631: 19834 files
Calculating effective area...
Simulation set 12360: 19781 files
Simulation set 12630: 19787 files
Calculating effective area...
Simulation set 12362: 19808 files
Simulation set 12631: 19834 files

In [71]:
eff_area, eff_area_error, _ = comp.analysis.get_effective_area(df_sim,
                                                                energy_bins, energy='MC')
eff_area_light, eff_area_error_light, _ = comp.analysis.get_effective_area(
                                                                df_sim[df_sim.MC_comp_class == 'light'],
                                                                energy_bins, energy='MC')
eff_area_heavy, eff_area_error_heavy, _ = comp.analysis.get_effective_area(
                                                                df_sim[df_sim.MC_comp_class == 'heavy'],
                                                                energy_bins, energy='MC')


Calculating effective area...
simlist = [12360 12362 12630 12631]
Simulation set 12360: 19781 files
type(sim) = <type 'numpy.int64'>
sim = 12360
Simulation set 12362: 19808 files
type(sim) = <type 'numpy.int64'>
sim = 12362
Simulation set 12630: 19787 files
type(sim) = <type 'numpy.int64'>
sim = 12630
Simulation set 12631: 19834 files
type(sim) = <type 'numpy.int64'>
sim = 12631
num_ptypes = 4
Calculating effective area...
simlist = [12360 12630]
Simulation set 12360: 19781 files
type(sim) = <type 'numpy.int64'>
sim = 12360
Simulation set 12630: 19787 files
type(sim) = <type 'numpy.int64'>
sim = 12630
num_ptypes = 2
Calculating effective area...
simlist = [12362 12631]
Simulation set 12362: 19808 files
type(sim) = <type 'numpy.int64'>
sim = 12362
Simulation set 12631: 19834 files
type(sim) = <type 'numpy.int64'>
sim = 12631
num_ptypes = 2

In [14]:
eff_area_light


Out[14]:
array([  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
         5.24026972e+01,   5.34145001e+01,   3.26247872e+02,
         1.31356318e+03,   8.64194357e+02,   2.88359479e+03,
         4.32422452e+03,   8.96141417e+03,   1.21499869e+04,
         2.19251818e+04,   3.35217211e+04,   4.60892729e+04,
         6.00292849e+04,   8.13773734e+04,   7.84626788e+04,
         9.71846210e+04,   1.07211232e+05,   1.23782486e+05,
         1.32868673e+05,   1.21280851e+05,   1.31204956e+05,
         1.39347403e+05,   1.21366758e+05,   1.30505274e+05,
         1.33588715e+05,   1.41885917e+05,   1.47158581e+05,
         1.55710566e+05,   1.44095655e+05,   1.53291849e+05,
         1.39841590e+05,   1.49307117e+05,   1.33264659e+05,
         1.43413024e+05,   1.33854688e+05,   1.36742008e+05,
         1.41580651e+05,   1.54560314e+05,   1.34595585e+05,
         1.30231468e+05,   1.44190774e+05,   1.51709532e+05,
         1.21735053e+05,   1.36834143e+05,   1.39435180e+05,
         1.61788813e+05,   1.29339411e+05,   1.48903444e+05,
         1.42750787e+05,   1.39730662e+05,   1.61179526e+05,
         1.77414987e+05,   1.25415581e+05,   1.59301968e+05,
         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   0.00000000e+00,   0.00000000e+00])

Fit functions to effective area data


In [15]:
p0 = [1.5e5, 8.0, 50.0]
popt_light, pcov_light = optimize.curve_fit(sigmoid_flat, energy_midpoints[midpoints_fitmask],
                                            eff_area_light[midpoints_fitmask], p0=p0,
                                            sigma=eff_area_error_light[midpoints_fitmask])
popt_heavy, pcov_heavy = optimize.curve_fit(sigmoid_flat, energy_midpoints[midpoints_fitmask],
                                            eff_area_heavy[midpoints_fitmask], p0=p0,
                                            sigma=eff_area_error_heavy[midpoints_fitmask])

In [16]:
print(popt_light)
print(popt_heavy)


[  1.41746363e+05   7.77283311e+00   4.63861302e+01]
[  1.36838018e+05   9.38458915e+00   5.66947059e+01]

In [17]:
perr_light = np.sqrt(np.diag(pcov_light))
print(perr_light)
perr_heavy = np.sqrt(np.diag(pcov_heavy))
print(perr_heavy)


[  2.53034467e+03   6.21806613e-01   3.67077785e+00]
[  3.01846146e+03   7.11873552e-01   4.22664184e+00]

In [18]:
avg = (popt_light[0] + popt_heavy[0]) / 2
print('avg eff area = {}'.format(avg))


avg eff area = 139292.190378

In [19]:
eff_area_light


Out[19]:
array([  0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
         5.24026972e+01,   5.34145001e+01,   3.26247872e+02,
         1.31356318e+03,   8.64194357e+02,   2.88359479e+03,
         4.32422452e+03,   8.96141417e+03,   1.21499869e+04,
         2.19251818e+04,   3.35217211e+04,   4.60892729e+04,
         6.00292849e+04,   8.13773734e+04,   7.84626788e+04,
         9.71846210e+04,   1.07211232e+05,   1.23782486e+05,
         1.32868673e+05,   1.21280851e+05,   1.31204956e+05,
         1.39347403e+05,   1.21366758e+05,   1.30505274e+05,
         1.33588715e+05,   1.41885917e+05,   1.47158581e+05,
         1.55710566e+05,   1.44095655e+05,   1.53291849e+05,
         1.39841590e+05,   1.49307117e+05,   1.33264659e+05,
         1.43413024e+05,   1.33854688e+05,   1.36742008e+05,
         1.41580651e+05,   1.54560314e+05,   1.34595585e+05,
         1.30231468e+05,   1.44190774e+05,   1.51709532e+05,
         1.21735053e+05,   1.36834143e+05,   1.39435180e+05,
         1.61788813e+05,   1.29339411e+05,   1.48903444e+05,
         1.42750787e+05,   1.39730662e+05,   1.61179526e+05,
         1.77414987e+05,   1.25415581e+05,   1.59301968e+05,
         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   0.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   0.00000000e+00,   0.00000000e+00])

In [20]:
light_chi2 = red_chisquared(eff_area_light, sigmoid_flat(energy_midpoints, *popt_light),
               eff_area_error_light, len(popt_light))
print(light_chi2)
heavy_chi2 = red_chisquared(eff_area_heavy,
               sigmoid_flat(energy_midpoints, *popt_heavy),
               eff_area_error_heavy, len(popt_heavy))
print(heavy_chi2)


38.0787059901
9.92114463142

Plot result


In [21]:
fig, ax = plt.subplots()
# plot effective area data points with poisson errors
ax.errorbar(np.log10(energy_midpoints), eff_area_light, yerr=eff_area_error_light,
            ls='None', marker='.')
ax.errorbar(np.log10(energy_midpoints), eff_area_heavy, yerr=eff_area_error_heavy,
            ls='None', marker='.')

# plot corresponding sigmoid fits to effective area
x = 10**np.arange(5.0, 9.5, 0.01)
ax.plot(np.log10(x), sigmoid_flat(x, *popt_light),
        color=color_dict['light'], label='light', marker='None', ls='-')
ax.plot(np.log10(x), sigmoid_flat(x, *popt_heavy),
        color=color_dict['heavy'], label='heavy', marker='None')


avg_eff_area = (sigmoid_flat(x, *popt_light) + sigmoid_flat(x, *popt_heavy)) / 2
ax.plot(np.log10(x), avg_eff_area,
        color=color_dict['total'], label='avg', marker='None')
ax.fill_between(np.log10(x),
        avg_eff_area-0.01*avg_eff_area,
        avg_eff_area+0.01*avg_eff_area,
        color=color_dict['total'], alpha=0.5)

ax.axvline(6.4, marker='None', ls='-.', color='k')

ax.set_ylabel('Effective area [m$^2$]')
ax.set_xlabel('$\mathrm{\log_{10}(E_{true}/GeV)}$')
# ax.set_title('$\mathrm{A_{eff} = 143177 \pm 1431.77 \ m^2}$')
ax.grid()
# ax.set_ylim([0, 180000])
ax.set_xlim([5.4, 8.1])
ax.set_title(config)

#set label style        
ax.ticklabel_format(style='sci',axis='y')
ax.yaxis.major.formatter.set_powerlimits((0,0))

leg = plt.legend(title='True composition')
for legobj in leg.legendHandles:
    legobj.set_linewidth(2.0)

# eff_area_outfile = os.path.join(comp.paths.figures_dir, 'effective-area-{}.png'.format(config))
# comp.check_output_dir(eff_area_outfile)
# plt.savefig(eff_area_outfile)
plt.show()


/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 "

Effective area as quality cuts are sequentially applied


In [19]:
df_sim, cut_dict_sim = comp.load_dataframe(datatype='sim', config='IC79', return_cut_dict=True)
standard_cut_keys = ['num_hits_1_60', 'IceTopQualityCuts', 'lap_InIce_containment',
#                 'num_hits_1_60', 'max_qfrac_1_60',
                'InIceQualityCuts']
# for cut in ['MilliNCascAbove2', 'MilliQtotRatio', 'MilliRloglBelow2', 'StochRecoSucceeded']:
#     standard_cut_keys += ['InIceQualityCuts_{}'.format(cut)]

eff_area_dict = {}
eff_area_err_dict = {}
selection_mask = np.array([True] * len(df_sim))
for key in standard_cut_keys:
    selection_mask *= cut_dict_sim[key]
    print(key, np.sum(selection_mask))
    eff_area, eff_area_error, _ = comp.analysis.get_effective_area(df_sim[selection_mask],
                                                                energy_bins, energy='MC')
#     eff_area, eff_area_error = comp.analysis.effective_area.effective_area(df_sim[selection_mask],
#                                                                            np.arange(5.0, 9.51, 0.1))
    eff_area_dict[key] = eff_area
    eff_area_err_dict[key] = eff_area_error


num_hits_1_60 397466
Calculating effective area...
Simulation set 7006: 30000 files
Simulation set 7007: 30000 files
Simulation set 7241: 10000 files
Simulation set 7242: 10000 files
Simulation set 7262: 19999 files
Simulation set 7263: 19998 files
Simulation set 7579: 12000 files
Simulation set 7784: 12000 files
Simulation set 7791: 12000 files
Simulation set 7851: 12000 files
IceTopQualityCuts 397466
Calculating effective area...
Simulation set 7006: 30000 files
Simulation set 7007: 30000 files
Simulation set 7241: 10000 files
Simulation set 7242: 10000 files
Simulation set 7262: 19999 files
Simulation set 7263: 19998 files
Simulation set 7579: 12000 files
Simulation set 7784: 12000 files
Simulation set 7791: 12000 files
Simulation set 7851: 12000 files
lap_InIce_containment 308874
Calculating effective area...
Simulation set 7006: 30000 files
Simulation set 7007: 30000 files
Simulation set 7241: 10000 files
Simulation set 7242: 10000 files
Simulation set 7262: 19999 files
Simulation set 7263: 19998 files
Simulation set 7579: 12000 files
Simulation set 7784: 12000 files
Simulation set 7791: 12000 files
Simulation set 7851: 12000 files
InIceQualityCuts 298466
Calculating effective area...
Simulation set 7006: 30000 files
Simulation set 7007: 30000 files
Simulation set 7241: 10000 files
Simulation set 7242: 10000 files
Simulation set 7262: 19999 files
Simulation set 7263: 19998 files
Simulation set 7579: 12000 files
Simulation set 7784: 12000 files
Simulation set 7791: 12000 files
Simulation set 7851: 12000 files

In [20]:
fig, ax = plt.subplots()
cut_labels = {'num_hits_1_60': 'NStations/NChannels', 'IceTopQualityCuts': 'IceTopQualityCuts',
             'lap_InIce_containment': 'InIce containment', 'InIceQualityCuts': 'InIceQualityCuts'}
for key in standard_cut_keys:
# plot effective area data points with poisson errors
    ax.errorbar(np.log10(energy_midpoints), eff_area_dict[key], yerr=eff_area_err_dict[key],
                ls='None', marker='.', label=cut_labels[key], alpha=0.75)

ax.set_ylabel('Effective area [m$^2$]')
ax.set_xlabel('$\log_{10}(E_{\mathrm{MC}}/\mathrm{GeV})$')
ax.grid()
# ax.set_ylim([0, 180000])
ax.set_xlim([5.4, 9.6])

#set label style        
ax.ticklabel_format(style='sci',axis='y')
ax.yaxis.major.formatter.set_powerlimits((0,0))

leg = plt.legend()

plt.savefig('/home/jbourbeau/public_html/figures/effective-area-cuts.png')
plt.show()



In [ ]:


In [ ]: