In [1]:
%load_ext watermark
%watermark -u -d -v -p numpy,matplotlib,scipy,pandas,sklearn,mlxtend
In [47]:
%matplotlib inline
from __future__ import division, print_function
from collections import defaultdict
import os
import glob
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import healpy as hp
from scipy.stats import ks_2samp
import multiprocessing as mp
from scipy.special import erfcinv
import pyprind
import dask
from dask import delayed, multiprocessing
from dask.diagnostics import ProgressBar
from icecube import astro
import comptools as comp
import comptools.analysis.plotting as plotting
import comptools.anisotropy.anisotropy as anisotropy
color_dict = comp.analysis.get_color_dict()
In [70]:
config = ['IC86.2011', 'IC86.2012', 'IC86.2013', 'IC86.2014','IC86.2015']
n_side = 64
smooth_rad = 5.0
scale = 3
decmax = -55
years_str = '2011-2015'
composition='all'
low_energy = True
n_bins = 36
In [71]:
df = anisotropy.get_test_stats(config=config, low_energy=low_energy, smooth=smooth_rad, n_bins=n_bins)
In [72]:
df['proj_RI_red_chi2'][df['proj_RI_red_chi2'].isnull()].index
Out[72]:
In [73]:
df['proj_RI_red_chi2'].describe()
Out[73]:
In [74]:
def calc_chi2(ri_0, ri_1, ri_err_0, ri_err_1):
chi2 = np.sum((ri_0-ri_1)**2/(ri_err_0**2+ri_err_1**2)) / ri_0.shape[0]
return chi2
In [75]:
def get_chi2_data(config):
kwargs_relint = {'smooth': smooth_rad, 'scale': None, 'decmax': -55, 'low_energy': True}
relint_light = anisotropy.get_map(config=config, name='relint', composition='light', **kwargs_relint)
relint_heavy = anisotropy.get_map(config=config, name='relint', composition='heavy', **kwargs_relint)
relerr_light = anisotropy.get_map(config=config, name='relerr', composition='light', **kwargs_relint)
relerr_heavy = anisotropy.get_map(config=config, name='relerr', composition='heavy', **kwargs_relint)
ri_light, ri_err_light, ra, _ = anisotropy.get_proj_relint(relint_light, relerr_light, n_bins=n_bins)
ri_heavy, ri_err_heavy, ra, _ = anisotropy.get_proj_relint(relint_heavy, relerr_heavy, n_bins=n_bins)
chi2_data = calc_chi2(ri_light, ri_heavy, ri_err_light, ri_err_heavy)
return chi2_data
In [76]:
def calc_sig_ks_pval(config):
kwargs_relint = {'smooth': smooth_rad, 'scale': None, 'decmax': -55, 'low_energy': True}
# Get 2D Li-Ma significance maps
sig_light = anisotropy.get_map(config=config, name='sig', composition='light', **kwargs_relint)
sig_heavy = anisotropy.get_map(config=config, name='sig', composition='heavy', **kwargs_relint)
# Construct masks
is_good_mask_light = (sig_light != hp.UNSEEN) & ~np.isnan(sig_light)
is_good_mask_heavy = (sig_heavy != hp.UNSEEN) & ~np.isnan(sig_heavy)
# Calculate ks test statistic and corresponding p-value
ks_statistic, pval = ks_2samp(sig_light[is_good_mask_light], sig_heavy[is_good_mask_heavy])
# Calculate ~area between cumulative distributions
bins = np.linspace(0, 5, 50)
# counts_light = np.histogram(sig_light[is_good_mask_light], bins=bins)[0]
# counts_heavy = np.histogram(sig_heavy[is_good_mask_heavy], bins=bins)[0]
counts_light = np.histogram(np.abs(sig_light[is_good_mask_light]), bins=bins)[0]
counts_heavy = np.histogram(np.abs(sig_heavy[is_good_mask_heavy]), bins=bins)[0]
sig_cumsum_diff_area = np.abs(np.cumsum(counts_light) - np.cumsum(counts_heavy)).sum()
return ks_statistic, pval, sig_cumsum_diff_area
In [77]:
plt.errorbar(ra, ri_light, yerr=ri_err_light)
plt.errorbar(ra, ri_heavy, yerr=ri_err_heavy)
plt.grid()
plt.show()
In [78]:
chi2_data = get_chi2_data(config)
sig_ks_dist_data, sig_pval_data, sig_cumsum_diff_area = calc_sig_ks_pval(config)
In [79]:
sig_ks_dist_data, sig_pval_data, sig_cumsum_diff_area
Out[79]:
In [80]:
test_stat = df['proj_RI_red_chi2']
pval = np.sum(test_stat > chi2_data)/len(test_stat)
print(pval)
significance = erfcinv(2*pval)*np.sqrt(2)
print(significance)
In [81]:
fig, ax = plt.subplots()
chi2_max = 4
chi2_bins = np.linspace(0, chi2_max, 100)
counts = np.histogram(df['proj_RI_red_chi2'], bins=chi2_bins)[0]
# ax = plotting.plot_steps(chi2_bins, counts, yerr=np.sqrt(counts))
ax = plotting.plot_steps(chi2_bins, counts, yerr=np.sqrt(counts), label='Random trials')
ax.axvline(chi2_data, marker='None', ls='-.', color='C2', lw=1.5,
label='Light-heavy split \n ($\mathrm{\chi^2_{red}} = '+'${:0.2f})'.format(chi2_data))
ax.set_xlabel('$\mathrm{\chi^2_{red}}$')
ax.set_ylabel('Counts')
ax.set_title('IC86 2011-2015 random trials')
# ax.text(0.3, 4, 'p-value $= {:g}$'.format(pval) + \
# '\n significance $= {:0.2f}\sigma$'.format(significance))
ax.set_ylim(1e-1)
ax.set_xlim(0, chi2_max)
ax.set_yscale("log", nonposy='clip')
ax.legend()
ax.grid()
plt.savefig(os.path.join(comp.paths.figures_dir, 'random-trials-chi2_IC86.2011-2015.png'))
plt.show()
In [16]:
df['sig_ks_dist'].plot(kind='hist', bins=50)
Out[16]:
In [17]:
df['sig_cumsum_diff_area'].plot(kind='hist', bins=50)
Out[17]:
In [18]:
np.log10(df['sig_pval']).plot(kind='hist', bins=50, logy=True)
Out[18]:
In [113]:
fig, ax = plt.subplots()
sig_pval_bins = np.linspace(-300, 0, 100)
counts = np.histogram(np.log10(df['sig_pval']), bins=sig_pval_bins)[0]
# ax = plotting.plot_steps(chi2_bins, counts, yerr=np.sqrt(counts))
ax = plotting.plot_steps(sig_pval_bins, counts, yerr=np.sqrt(counts), label='Random trials')
ax.axvline(np.log10(sig_pval_data), marker='None', ls='-.', color='C2', lw=1.5,
)
# label='Light-heavy split \n ($\mathrm{\chi^2_{red}} = '+'${:0.2f})'.format(chi2_data))
ax.set_xlabel('KS p-value')
ax.set_ylabel('Counts')
ax.set_title('IC86 2011-2015 random trials')
# ax.text(1, 5, 'p-value $= {:g}$'.format(pval) + \
# '\n significance $= {:0.2f}\sigma$'.format(significance))
# ax.set_ylim((1e-1, 1e2))
# ax.set_xlim(0, chi2_max)
# ax.set_yscale("log", nonposy='clip')
ax.legend()
ax.grid()
# plt.savefig(os.path.join(comp.paths.figures_dir, 'random-trials-chi2_IC86.2012-2014.png'))
plt.show()
In [35]:
df['sig_ks_dist'].describe()
Out[35]:
In [47]:
fig, ax = plt.subplots()
sig_ks_bins = np.linspace(0, 0.4, 75)
counts = np.histogram(df['sig_ks_dist'], bins=sig_ks_bins)[0]
# ax = plotting.plot_steps(chi2_bins, counts, yerr=np.sqrt(counts))
ax = plotting.plot_steps(sig_ks_bins, counts, yerr=np.sqrt(counts), label='Random trials')
ax.axvline(sig_ks_dist_data, marker='None', ls='-.', color='C2', lw=1.5,
label='Light-heavy split')
# label='Light-heavy split \n ($\mathrm{\chi^2_{red}} = '+'${:0.2f})'.format(chi2_data))
ax.set_xlabel('KS distance')
ax.set_ylabel('Counts')
ax.set_title('IC86 2011-2015 random trials')
test_stat = df['sig_ks_dist']
pval = np.sum(test_stat > sig_ks_dist_data)/len(test_stat)
significance = erfcinv(2*pval)*np.sqrt(2)
ax.text(0.06, 2, 'p-value $= {:g}$'.format(pval) + \
'\n significance $= {:0.2f}\sigma$'.format(significance))
ax.set_ylim(0)
# ax.set_xlim(0, chi2_max)
ax.set_yscale("log", nonposy='clip')
ax.set_ylim(1e-1)
ax.legend()
ax.grid()
plt.savefig(os.path.join(comp.paths.figures_dir, 'random-trials-ks-dist_IC86.2011-2015.png'))
plt.show()
In [115]:
test_stat = df['sig_ks_dist']
pval = np.sum(test_stat > sig_ks_dist_data)/len(test_stat)
print(pval)
significance = erfcinv(2*pval)*np.sqrt(2)
print(significance)
In [125]:
fig, ax = plt.subplots()
sig_cumsum_diff_bins = np.linspace(0, 70000, 100)
counts = np.histogram(df['sig_cumsum_diff_area'], bins=sig_cumsum_diff_bins)[0]
# ax = plotting.plot_steps(chi2_bins, counts, yerr=np.sqrt(counts))
ax = plotting.plot_steps(sig_cumsum_diff_bins, counts, yerr=np.sqrt(counts), label='Random trials')
ax.axvline(sig_cumsum_diff_area, marker='None', ls='-.', color='C2', lw=1.5,
)
# label='Light-heavy split \n ($\mathrm{\chi^2_{red}} = '+'${:0.2f})'.format(chi2_data))
ax.set_xlabel('Cumulative sum area proxy')
ax.set_ylabel('Counts')
ax.set_title('IC86 2011-2015 random trials')
# ax.text(1, 5, 'p-value $= {:g}$'.format(pval) + \
# '\n significance $= {:0.2f}\sigma$'.format(significance))
# ax.set_ylim((1e-1, 1e2))
# ax.set_xlim(0, chi2_max)
# ax.set_yscale("log", nonposy='clip')
ax.legend()
ax.grid()
# plt.savefig(os.path.join(comp.paths.figures_dir, 'random-trials-chi2_IC86.2012-2014.png'))
plt.show()
In [126]:
test_stat = df['sig_cumsum_diff_area']
pval = np.sum(test_stat > sig_cumsum_diff_area)/len(test_stat)
print(pval)
significance = erfcinv(2*pval)*np.sqrt(2)
print(significance)
In [40]:
fig, ax = plt.subplots()
sig_ks_bins = np.linspace(0, 0.5, 100)
configs = [['IC86.2012'], ['IC86.2012', 'IC86.2013', 'IC86.2014']]
# configs = [['IC86.2012'], ['IC86.2012', 'IC86.2013'], ['IC86.2012', 'IC86.2013', 'IC86.2014']]
for config, color in zip(configs, ['C0', 'C1', 'C2']):
print('On {}...'.format(config))
df = anisotropy.get_test_stats(config=config, low_energy=True)
counts = np.histogram(df['sig_ks_dist'], bins=sig_ks_bins)[0]
ax = plotting.plot_steps(sig_ks_bins, counts, yerr=np.sqrt(counts), label='{} years'.format(len(config)), color=color, ax=ax)
# chi2_data = get_chi2_data(config)
# ax.axvline(chi2_data, marker='None', ls='-.', color=color, lw=1.5)
# # label='Light-heavy split \n ($\mathrm{\chi^2_{red}} = '+'${:0.2f})'.format(chi2_data))
# ax.axvline(6.1, marker='None', ls='-.', color=color, lw=1.5)
ax.set_xlabel('KS distance')
ax.set_ylabel('Counts')
ax.set_title('1000 random trials')
# ax.text(1, 5, 'p-value $= {:g}$'.format(pval) + \
# '\n significance $= {:0.2f}\sigma$'.format(significance))
# ax.set_ylim((1e-1, 1e2))
# ax.set_xlim(0, chi2_max)
# ax.set_yscale("log", nonposy='clip')
ax.legend()
ax.grid()
# fig.text(0.5, 0.0, 'RA [ $^{\circ}$]', ha='center', fontsize=16)
# fig.text(0.0, 0.5, '$\mathrm{\langle RI \\rangle }$', va='center', rotation='vertical', fontsize=16)
plt.tight_layout()
# plt.savefig(os.path.join(comp.paths.figures_dir, 'trials_compare-running-sum.png'))
plt.show()
In [41]:
fig, ax = plt.subplots()
sig_pval_bins = np.linspace(-300, 0, 100)
configs = [['IC86.2012'], ['IC86.2012', 'IC86.2013', 'IC86.2014']]
# configs = [['IC86.2012'], ['IC86.2012', 'IC86.2013'], ['IC86.2012', 'IC86.2013', 'IC86.2014']]
for config, color in zip(configs, ['C0', 'C1', 'C2']):
print('On {}...'.format(config))
df = anisotropy.get_test_stats(config=config, low_energy=True)
counts = np.histogram(np.log10(df['sig_pval']), bins=sig_pval_bins)[0]
ax = plotting.plot_steps(sig_pval_bins, counts, yerr=np.sqrt(counts), label='{} years'.format(len(config)), color=color, ax=ax)
# chi2_data = get_chi2_data(config)
# ax.axvline(chi2_data, marker='None', ls='-.', color=color, lw=1.5)
# # label='Light-heavy split \n ($\mathrm{\chi^2_{red}} = '+'${:0.2f})'.format(chi2_data))
# ax.axvline(6.1, marker='None', ls='-.', color=color, lw=1.5)
ax.set_xlabel('KS p-value')
ax.set_ylabel('Counts')
ax.set_title('1000 random trials')
# ax.text(1, 5, 'p-value $= {:g}$'.format(pval) + \
# '\n significance $= {:0.2f}\sigma$'.format(significance))
# ax.set_ylim((1e-1, 1e2))
# ax.set_xlim(0, chi2_max)
# ax.set_yscale("log", nonposy='clip')
ax.legend()
ax.grid()
# fig.text(0.5, 0.0, 'RA [ $^{\circ}$]', ha='center', fontsize=16)
# fig.text(0.0, 0.5, '$\mathrm{\langle RI \\rangle }$', va='center', rotation='vertical', fontsize=16)
plt.tight_layout()
# plt.savefig(os.path.join(comp.paths.figures_dir, 'trials_compare-running-sum.png'))
plt.show()
In [4]:
fig, ax = plt.subplots()
chi2_max = 25
chi2_bins = np.linspace(0, 30, 100)
for config, color in zip(['IC86.2012', 'IC86.2013', 'IC86.2014'],
['C0', 'C1', 'C2']):
df = anisotropy.get_test_stats(config=config, low_energy=True)
counts = np.histogram(df.chi2, bins=chi2_bins)[0]
ax = plotting.plot_steps(chi2_bins, counts, yerr=np.sqrt(counts), label=config, color=color)
# ax.axvline(chi2_data, marker='None', ls='-.', color='C2', lw=1.5,
# label='Light-heavy split \n ($\mathrm{\chi^2_{red}} = '+'${:0.2f})'.format(chi2_data))
ax.set_xlabel('$\mathrm{\chi^2_{red}}$')
ax.set_ylabel('Counts')
ax.set_title('1000 random trials')
# ax.text(1, 5, 'p-value $= {:g}$'.format(pval) + \
# '\n significance $= {:0.2f}\sigma$'.format(significance))
ax.set_ylim((1e-1, 1e2))
ax.set_xlim(0, chi2_max)
ax.set_yscale("log", nonposy='clip')
ax.legend()
ax.grid()
plt.savefig(os.path.join(comp.paths.figures_dir, 'random-trials-chi2_compare-years.png'))
plt.show()
In [11]:
fig, axarr = plt.subplots(2, 3, figsize=(10, 8), sharex=True, sharey=True)
for config, ax in zip(['IC86.2011', 'IC86.2012', 'IC86.2013', 'IC86.2014', 'IC86.2015'],
axarr.flatten()):
print('On {}...'.format(config))
kwargs_relint = {'smooth': 20, 'scale': None, 'decmax': -55, 'low_energy': True}
relint_light = anisotropy.get_map(config=config, name='relint', composition='light', **kwargs_relint)
relint_heavy = anisotropy.get_map(config=config, name='relint', composition='heavy', **kwargs_relint)
relerr_light = anisotropy.get_map(config=config, name='relerr', composition='light', **kwargs_relint)
relerr_heavy = anisotropy.get_map(config=config, name='relerr', composition='heavy', **kwargs_relint)
ri_light, ri_err_light, ra, _ = anisotropy.get_proj_relint(relint_light, relerr_light, n_bins=24)
ri_heavy, ri_err_heavy, ra, _ = anisotropy.get_proj_relint(relint_heavy, relerr_heavy, n_bins=24)
ax.errorbar(ra, ri_light, yerr=ri_err_light, label='light', marker='.', ls=':')
ax.errorbar(ra, ri_heavy, yerr=ri_err_heavy, label='heavy', marker='.', ls=':')
ax.axhline(0, marker='None', ls='-.', c='k')
chi2_data = calc_chi2(ri_light, ri_heavy, ri_err_light, ri_err_heavy)
ax.set_title('{}: '.format(config)+'$\mathrm{\chi^2_{red}} =$'+ ' {:0.1f}'.format(chi2_data))
for i, row in enumerate(axarr):
for j, cell in enumerate(row):
if i == len(axarr) - 1:
cell.set_xlabel('RA [ $^{\circ}$]')
if j == 0:
cell.set_ylabel('$\mathrm{\langle RI \\rangle }$')
ax.grid()
ax.legend()
# fig.text(0.5, 0.0, 'RA [ $^{\circ}$]', ha='center', fontsize=16)
# fig.text(0.0, 0.5, '$\mathrm{\langle RI \\rangle }$', va='center', rotation='vertical', fontsize=16)
plt.tight_layout()
plt.savefig(os.path.join(comp.paths.figures_dir, 'light-heavy-split_compare-years.png'))
plt.show()
In [24]:
fig, axarr = plt.subplots(2, 3, figsize=(10, 8), sharex=True, sharey=True)
configs = [['IC86.2011'], ['IC86.2011', 'IC86.2012'], ['IC86.2011', 'IC86.2012', 'IC86.2013'],
['IC86.2011', 'IC86.2012', 'IC86.2013', 'IC86.2014'],
['IC86.2011', 'IC86.2012', 'IC86.2013', 'IC86.2014', 'IC86.2015']]
for config, ax in zip(configs, axarr.flatten()):
print('On {}...'.format(config))
kwargs_relint = {'smooth': 20, 'scale': None, 'decmax': -55, 'low_energy': True}
relint_light = anisotropy.get_map(config=config, name='relint', composition='light', **kwargs_relint)
relint_heavy = anisotropy.get_map(config=config, name='relint', composition='heavy', **kwargs_relint)
relerr_light = anisotropy.get_map(config=config, name='relerr', composition='light', **kwargs_relint)
relerr_heavy = anisotropy.get_map(config=config, name='relerr', composition='heavy', **kwargs_relint)
ri_light, ri_err_light, ra, _ = anisotropy.get_proj_relint(relint_light, relerr_light, n_bins=24)
ri_heavy, ri_err_heavy, ra, _ = anisotropy.get_proj_relint(relint_heavy, relerr_heavy, n_bins=24)
ax.errorbar(ra, ri_light, yerr=ri_err_light, label='light', marker='.', ls=':')
ax.errorbar(ra, ri_heavy, yerr=ri_err_heavy, label='heavy', marker='.', ls=':')
ax.axhline(0, marker='None', ls='-.', c='k')
popt_light, perr, _ = anisotropy.get_proj_fit_params(np.deg2rad(ra), ri_light, sigmay=ri_err_light)
popt_heavy, perr, _ = anisotropy.get_proj_fit_params(np.deg2rad(ra), ri_heavy, sigmay=ri_err_heavy)
ax.plot(ra, anisotropy.cos_fit_func(np.deg2rad(ra), *popt_light[:3]), marker='None', color=color_dict['light'])
ax.plot(ra, anisotropy.cos_fit_func(np.deg2rad(ra), *popt_heavy[:3]), marker='None', color=color_dict['heavy'])
chi2_data = calc_chi2(ri_light, ri_heavy, ri_err_light, ri_err_heavy)
ax.set_title('{} years: '.format(len(config))+'$\mathrm{\chi^2_{red}} =$'+ ' {:0.1f}'.format(chi2_data))
for i, row in enumerate(axarr):
for j, cell in enumerate(row):
if i == len(axarr) - 1:
cell.set_xlabel('RA [ $^{\circ}$]')
if j == 0:
cell.set_ylabel('$\mathrm{\langle RI \\rangle }$')
ax.grid()
ax.legend()
# fig.text(0.5, 0.0, 'RA [ $^{\circ}$]', ha='center', fontsize=16)
# fig.text(0.0, 0.5, '$\mathrm{\langle RI \\rangle }$', va='center', rotation='vertical', fontsize=16)
plt.tight_layout()
plt.savefig(os.path.join(comp.paths.figures_dir, 'light-heavy-split_compare-running-sum.png'))
plt.show()
In [14]:
fig, ax = plt.subplots()
chi2_max = 25
chi2_bins = np.linspace(0, 30, 100)
configs = [['IC86.2012'], ['IC86.2012', 'IC86.2013'], ['IC86.2012', 'IC86.2013', 'IC86.2014']]
for config, color in zip(configs, ['C0', 'C1', 'C2']):
print('On {}...'.format(config))
df = anisotropy.get_test_stats(config=config, low_energy=True)
counts = np.histogram(df.chi2, bins=chi2_bins)[0]
ax = plotting.plot_steps(chi2_bins, counts, yerr=np.sqrt(counts), label='{} years'.format(len(config)), color=color, ax=ax)
chi2_data = get_chi2_data(config)
ax.axvline(chi2_data, marker='None', ls='-.', color=color, lw=1.5)
# label='Light-heavy split \n ($\mathrm{\chi^2_{red}} = '+'${:0.2f})'.format(chi2_data))
ax.axvline(6.1, marker='None', ls='-.', color=color, lw=1.5)
ax.set_xlabel('$\mathrm{\chi^2_{red}}$')
ax.set_ylabel('Counts')
ax.set_title('1000 random trials')
# ax.text(1, 5, 'p-value $= {:g}$'.format(pval) + \
# '\n significance $= {:0.2f}\sigma$'.format(significance))
ax.set_ylim((1e-1, 1e2))
ax.set_xlim(0, chi2_max)
ax.set_yscale("log", nonposy='clip')
ax.legend()
ax.grid()
# fig.text(0.5, 0.0, 'RA [ $^{\circ}$]', ha='center', fontsize=16)
# fig.text(0.0, 0.5, '$\mathrm{\langle RI \\rangle }$', va='center', rotation='vertical', fontsize=16)
plt.tight_layout()
plt.savefig(os.path.join(comp.paths.figures_dir, 'trials_compare-running-sum.png'))
plt.show()
In [31]:
plt.plot(counts.cumsum())
plt.xlabel('')
plt.grid()
In [31]:
n_bins
Out[31]:
In [33]:
np.sum(test_stat > chi2_data)
Out[33]:
In [8]:
map_dir = '/data/user/jbourbeau/composition/IC86.2013_data/anisotropy/random_trials/'
sample_0_file_pattern = os.path.join(map_dir, 'random_split_0_trial-*.fits')
sample_1_file_pattern = os.path.join(map_dir, 'random_split_1_trial-*.fits')
infiles_sample_0 = sorted(glob.glob(sample_0_file_pattern))
infiles_sample_1 = sorted(glob.glob(sample_1_file_pattern))
# kwargs_relint = {'smooth': 20, 'scale': None, 'decmax': -55}
n_trials = 2
In [66]:
def get_proj_trial_dict(config, trial_num):
proj_dict = {}
proj_dict['trial'] = trial_num
trial_files_0 = []
trial_files_1 = []
for c in config:
map_dir = os.path.join(comp.paths.comp_data_dir, c + '_data',
'anisotropy/random_trials')
sample_0_file = os.path.join(map_dir, 'random_split_0_trial-{}.fits'.format(trial_num))
trial_files_0.append(sample_0_file)
sample_1_file = os.path.join(map_dir, 'random_split_1_trial-{}.fits'.format(trial_num))
trial_files_1.append(sample_1_file)
kwargs_relint = {'smooth': smooth_rad, 'scale': None, 'decmax': -55}
relint_0 = anisotropy.get_map(files=trial_files_0, name='relint', **kwargs_relint)
relint_1 = anisotropy.get_map(files=trial_files_1, name='relint', **kwargs_relint)
relerr_0 = anisotropy.get_map(files=trial_files_0, name='relerr', **kwargs_relint)
relerr_1 = anisotropy.get_map(files=trial_files_1, name='relerr', **kwargs_relint)
ri_0, ri_0_err, ra, _ = anisotropy.get_proj_relint(relint_0, relerr_0, n_bins=n_bins)
ri_1, ri_1_err, ra, _ = anisotropy.get_proj_relint(relint_1, relerr_1, n_bins=n_bins)
proj_dict['ri_0'] = ri_0
proj_dict['ri_0_err'] = ri_0_err
proj_dict['ri_1'] = ri_1
proj_dict['ri_1_err'] = ri_1_err
chi2 = calc_chi2(ri_0, ri_1, ri_0_err, ri_1_err)
proj_dict['chi2'] = chi2
return proj_dict
In [67]:
n_trials = 16
proj_dicts = [delayed(get_proj_trial_dict)(config, trial_num) for trial_num in range(n_trials)]
proj_df = delayed(pd.DataFrame.from_records)(proj_dicts, index='trial')
In [68]:
with ProgressBar() as bar:
# with dask.set_options(pool=mp.Pool(10)) as _, ProgressBar() as bar:
proj_df = proj_df.compute(get=multiprocessing.get, num_workers=10)
proj_df.head()
Out[68]:
In [69]:
# n_side = np.sqrt(n_trials).astype(int)
n_side = 4
# n_bins = 72
# fig, ax = plt.subplots(n, n, figsize=(10, 8), sharex=True, sharey=True)
fig, axarr = plt.subplots(n_side, n_side, figsize=(10, 8), sharex=True, sharey=True)
for trial, ax in zip(proj_df.index, axarr.flatten()):
proj_trial = proj_df.loc[trial]
ra_bins = np.linspace(0, 360, n_bins + 1)
ax = plotting.plot_steps(ra_bins, proj_trial['ri_0'], yerr=proj_trial['ri_0_err'], color='C3', ax=ax)
ax = plotting.plot_steps(ra_bins, proj_trial['ri_1'], yerr=proj_trial['ri_1_err'], color='C4', ax=ax)
ax.set_title('$\mathrm{\chi^2_{red}}$' + ' = {:0.1f}'.format(proj_trial['chi2']))
ax.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
ax.set_xlim(ra_bins.min(), ra_bins.max())
# ax.set_ylim(-1e-2, 1e-2)
ax.invert_xaxis()
ax.grid()
for i, row in enumerate(axarr):
for j, cell in enumerate(row):
if i == len(axarr) - 1:
cell.set_xlabel('RA [ $^{\circ}$]')
if j == 0:
cell.set_ylabel('$\mathrm{\langle RI \\rangle }$')
plt.tight_layout()
plt.savefig(os.path.join(comp.paths.figures_dir, 'random-trials-grid.png'))
plt.show()
In [33]:
n_trials = 16
sig_trials_0 = []
sig_trials_1 = []
for trial_num in pyprind.prog_bar(range(n_trials)):
trial_files_0 = []
trial_files_1 = []
for c in config:
map_dir = os.path.join(comp.paths.comp_data_dir, c + '_data',
'anisotropy/random_trials')
sample_0_file = os.path.join(map_dir, 'random_split_0_trial-{}.fits'.format(trial_num))
trial_files_0.append(sample_0_file)
sample_1_file = os.path.join(map_dir, 'random_split_1_trial-{}.fits'.format(trial_num))
trial_files_1.append(sample_1_file)
kwargs_relint = {'smooth': smooth_rad, 'scale': None, 'decmax': -55}
# kwargs_relint = {'smooth': 20, 'scale': None, 'decmax': -55}
# Get 2D Li-Ma significance maps
sig_0 = anisotropy.get_map(files=trial_files_0, name='sig', **kwargs_relint)
sig_1 = anisotropy.get_map(files=trial_files_1, name='sig', **kwargs_relint)
sig_trials_0.append(sig_0)
sig_trials_1.append(sig_1)
In [34]:
n_side = 4
# fig, ax = plt.subplots(n, n, figsize=(10, 8), sharex=True, sharey=True)
fig, axarr = plt.subplots(n_side, n_side, figsize=(10, 8), sharex=True, sharey=False)
for trial, ax in enumerate(axarr.flatten()):
trial = int(trial)
print('On trial {}...'.format(trial))
# Get 2D Li-Ma significance maps
sig_0 = sig_trials_0[trial]
sig_1 = sig_trials_1[trial]
# sig_0 = np.abs(sig_trials_0[trial])
# sig_1 = np.abs(sig_trials_1[trial])
# Construct masks
is_good_mask_0 = (sig_0 != hp.UNSEEN) & ~np.isnan(sig_0)
is_good_mask_1 = (sig_1 != hp.UNSEEN) & ~np.isnan(sig_1)
bins = np.linspace(-5, 5, 50)
# bins = np.linspace(0, 5, 50)
counts_0, _, _ = ax.hist(sig_0[is_good_mask_0], bins=bins, alpha=0.6, color='C3', label='Split 0')
counts_1, _, _ = ax.hist(sig_1[is_good_mask_1], bins=bins, alpha=0.6, color='C4', label='Split 1')
ks_statistic, pval = ks_2samp(sig_0[is_good_mask_0], sig_1[is_good_mask_1])
ax.set_title('KS p-value = {:0.1e}'.format(pval))
# ax.set_title('$\mathrm{\chi^2_{red}}$'+' = {:0.2f}'.format(chi2_trials[trial]))
ax.grid()
# ax.set_xlim(0, 5)
# ax.legend()
for i, row in enumerate(axarr):
for j, cell in enumerate(row):
if i == len(axarr) - 1:
cell.set_xlabel('$\mathrm{\sigma_{LM}}$')
if j == 0:
cell.set_ylabel('Counts')
plt.tight_layout()
# fig.text(0.5, 0.00, 'RA [ $^{\circ}$]', ha='center')
# fig.text(0.00, 0.5, '$\mathrm{\langle RI \\rangle }$', va='center', rotation='vertical')
plt.savefig(os.path.join(comp.paths.figures_dir, 'random-trials-grid.png'))
plt.show()
In [141]:
n_side = 4
# fig, ax = plt.subplots(n, n, figsize=(10, 8), sharex=True, sharey=True)
fig, axarr = plt.subplots(n_side, n_side, figsize=(10, 8), sharex=True, sharey=False)
for trial, ax in enumerate(axarr.flatten()):
trial = int(trial)
print('On trial {}...'.format(trial))
# Get 2D Li-Ma significance maps
sig_0 = sig_trials_0[trial]
sig_1 = sig_trials_1[trial]
# sig_0 = np.abs(sig_trials_0[trial])
# sig_1 = np.abs(sig_trials_1[trial])
# Construct masks
is_good_mask_0 = (sig_0 != hp.UNSEEN) & ~np.isnan(sig_0)
is_good_mask_1 = (sig_1 != hp.UNSEEN) & ~np.isnan(sig_1)
# bins = np.linspace(0, 5, 50)
bins = np.linspace(-5, 5, 50)
counts_0 = np.histogram(sig_0[is_good_mask_0], bins=bins)[0]
counts_1 = np.histogram(sig_1[is_good_mask_1], bins=bins)[0]
ax.plot(np.cumsum(counts_0), alpha=0.7, color='C2', label='Split 0')
ax.plot(np.cumsum(counts_1), alpha=0.7, color='C3', label='Split 1')
ax.set_title(np.abs(np.cumsum(counts_0) - np.cumsum(counts_1)).sum())
# ax.set_title('$\mathrm{\chi^2_{red}}$'+' = {:0.2f}'.format(chi2_trials[trial]))
ax.grid()
# ax.set_xlim(0, 5)
# ax.legend()
for i, row in enumerate(axarr):
for j, cell in enumerate(row):
if i == len(axarr) - 1:
cell.set_xlabel('$\mathrm{\sigma_{LM}}$')
if j == 0:
cell.set_ylabel('Counts')
plt.tight_layout()
# fig.text(0.5, 0.00, 'RA [ $^{\circ}$]', ha='center')
# fig.text(0.00, 0.5, '$\mathrm{\langle RI \\rangle }$', va='center', rotation='vertical')
# plt.savefig(os.path.join(comp.paths.figures_dir, 'random-trials-grid.png'))
plt.show()
In [77]:
np.cumsum(counts_0)
Out[77]:
In [9]:
kwargs_relint = {'smooth': 20, 'scale': None, 'decmax': -55}
file_0 = infiles_sample_0[501]
print(file_0)
file_1 = infiles_sample_1[501]
print(file_1)
relint_0 = anisotropy.get_map(files=file_0, name='relint', **kwargs_relint)
relint_1 = anisotropy.get_map(files=file_1, name='relint', **kwargs_relint)
relerr_0 = anisotropy.get_map(files=file_0, name='relerr', **kwargs_relint)
relerr_1 = anisotropy.get_map(files=file_1, name='relerr', **kwargs_relint)
ri_0, ri_err_0, ra, ra_err = anisotropy.get_proj_relint(relint_0, relerr_0, n_bins=24)
ri_1, ri_err_1, ra, ra_err = anisotropy.get_proj_relint(relint_1, relerr_1, n_bins=24)
plt.errorbar(ra, ri_0, yerr=ri_err_0)
plt.errorbar(ra, ri_1, yerr=ri_err_1)
plt.grid()
plt.show()
In [4]:
kwargs_data = {'config': config, 'low_energy': low_energy, 'smooth': False, 'scale': False, 'decmax': decmax}
In [14]:
n_total = anisotropy.get_num_events(config=config, composition='all', decmax=decmax)
n_light = anisotropy.get_num_events(config=config, composition='light', decmax=decmax)
n_heavy = anisotropy.get_num_events(config=config, composition='heavy', decmax=decmax)
In [15]:
n_light/n_total, n_heavy/n_total
Out[15]:
In [69]:
data_heavy = anisotropy.get_map('data', composition='heavy', **kwargs_data)
data_heavy[ data_heavy < 0 ] = 0
'{:g}'.format(data_heavy.sum())
Out[69]:
In [36]:
kwargs_relint = {'config': config, 'low_energy': low_energy, 'smooth': smooth_rad, 'scale': scale, 'decmax': decmax}
In [37]:
relint_all = anisotropy.get_map('relint', composition='all', **kwargs_relint)
relint_light = anisotropy.get_map('relint', composition='light', **kwargs_relint)
relint_heavy = anisotropy.get_map('relint', composition='heavy', **kwargs_relint)
In [38]:
print(rel_int_all.max())
print(rel_int_all[rel_int_all != hp.UNSEEN].min())
title = 'Relative Intensity [$\\times \ 10^{}$]'.format('{'+str(-scale)+'}')
fig, ax = anisotropy.plot_skymap(relint_all, smooth=smooth_rad, decmax=decmax, scale=scale,
color_palette='RdBu_r', symmetric=True, llabel='IC86 2012-2015',
cbar_title=title, cbar_min=-2.5, cbar_max=2.5, polar=True)
outfile = 'IC86-{}_relint_{}_nside-{}_smooth-{:0.1f}.png'.format(years_str, 'all', n_side, smooth_rad)
plt.savefig(os.path.join(comp.paths.figures_dir, outfile))
In [39]:
print(rel_int_light.max())
print(rel_int_light[rel_int_light != hp.UNSEEN].min())
title = 'Relative Intensity [$\\times \ 10^{}$]'.format('{'+str(-scale)+'}')
fig, ax = anisotropy.plot_skymap(relint_light, smooth=smooth_rad, decmax=decmax, scale=scale,
color_palette='RdBu_r', symmetric=True, llabel='IC86 2012-2015',
cbar_title=title, cbar_min=-2.5, cbar_max=2.5, polar=True)
outfile = 'IC86-{}_relint_{}_nside-{}_smooth-{:0.1f}.png'.format(years_str, 'light', n_side, smooth_rad)
plt.savefig(os.path.join(comp.paths.figures_dir, outfile))
In [40]:
print(rel_int_heavy.max())
print(rel_int_heavy[rel_int_heavy != hp.UNSEEN].min())
title = 'Relative Intensity [$\\times \ 10^{}$]'.format('{'+str(-scale)+'}')
fig, ax = anisotropy.plot_skymap(relint_heavy, smooth=smooth_rad, decmax=decmax, scale=scale,
color_palette='RdBu_r', symmetric=True, llabel='IC86 2012-2015',
cbar_title=title, cbar_min=-2.5, cbar_max=2.5, polar=True)
outfile = 'IC86-{}_relint_{}_nside-{}_smooth-{:0.1f}.png'.format(years_str, 'heavy', n_side, smooth_rad)
plt.savefig(os.path.join(comp.paths.figures_dir, outfile))
In [46]:
kwargs_sig = {'config': config, 'low_energy': low_energy, 'smooth': smooth_rad, 'scale': None, 'decmax': decmax}
In [47]:
sig_all = anisotropy.get_map('sig', composition='all', **kwargs_sig)
sig_light = anisotropy.get_map('sig', composition='light', **kwargs_sig)
sig_heavy = anisotropy.get_map('sig', composition='heavy', **kwargs_sig)
In [48]:
print(sig_all.max())
print(sig_all[sig_all != hp.UNSEEN].min())
title = 'Significance [$\mathrm{\sigma_{LM}}$]'
fig, ax = anisotropy.plot_skymap(sig_all, smooth=smooth_rad, decmax=decmax, scale=scale,
color_palette='RdBu_r', symmetric=True, llabel='IC86 2012-2015',
cbar_title=title, cbar_min=-4.1, cbar_max=4.1, polar=True)
outfile = 'IC86-{}_sig_{}_nside-{}_smooth-{:0.1f}.png'.format(years_str, 'all', n_side, smooth_rad)
plt.savefig(os.path.join(comp.paths.figures_dir, outfile))
In [49]:
print(sig_light.max())
print(sig_light[sig_light != hp.UNSEEN].min())
title = 'Significance [$\mathrm{\sigma_{LM}}$]'
fig, ax = anisotropy.plot_skymap(sig_light, smooth=smooth_rad, decmax=decmax, scale=scale,
color_palette='RdBu_r', symmetric=True, llabel='IC86 2012-2015',
cbar_title=title, cbar_min=-4.1, cbar_max=4.1, polar=True)
outfile = 'IC86-{}_sig_{}_nside-{}_smooth-{:0.1f}.png'.format(years_str, 'light', n_side, smooth_rad)
plt.savefig(os.path.join(comp.paths.figures_dir, outfile))
In [50]:
print(sig_heavy.max())
print(sig_heavy[sig_heavy != hp.UNSEEN].min())
title = 'Significance [$\mathrm{\sigma_{LM}}$]'
fig, ax = anisotropy.plot_skymap(sig_heavy, smooth=smooth_rad, decmax=decmax, scale=scale,
color_palette='RdBu_r', symmetric=True, llabel='IC86 2012-2015',
cbar_title=title, cbar_min=-4.1, cbar_max=4.1, polar=True)
outfile = 'IC86-{}_sig_{}_nside-{}_smooth-{:0.1f}.png'.format(years_str, 'heavy', n_side, smooth_rad)
plt.savefig(os.path.join(comp.paths.figures_dir, outfile))
In [30]:
rel_int_diff = anisotropy.get_relint_diff(**kwargs_relint)
print(rel_int_diff.max())
print(rel_int_diff[rel_int_diff != hp.UNSEEN].min())
In [31]:
title = '$\mathrm{\Delta \ RI}$'
fig, ax = anisotropy.plot_skymap(rel_int_diff, smooth=smooth_rad, decmax=decmax, scale=scale,
color_palette='BrBG', symmetric=True, llabel='IC86 2012-2013',
cbar_title=title, cbar_min=-2.5, cbar_max=2.5, polar=True)
outfile = 'IC86-{}_RI-diff_nside-{}_smooth-{:0.1f}.png'.format(years_str, n_side, smooth_rad)
plt.savefig(os.path.join(comp.paths.figures_dir, outfile))
[ back to top ]
In [4]:
kwargs_relint = {'config': config, 'low_energy': low_energy, 'smooth': smooth_rad, 'scale': None, 'decmax': decmax}
In [5]:
relint_all = anisotropy.get_map('relint', composition='all', **kwargs_relint)
relint_all_err = anisotropy.get_map('relerr', composition='all', **kwargs_relint)
In [80]:
ri_all, ri_all_err, ra, ra_err = anisotropy.get_proj_relint(relint_all, relint_all_err, n_bins=30)
In [81]:
fig, ax = plt.subplots()
ax.errorbar(ra, ri_all, yerr=ri_all_err, marker='.', ls='None', c='C2', label='all')
ax.axhline(0, marker='None', ls='-.', c='k')
ax.set_ylabel('$\mathrm{\langle RI \\rangle }$')
ax.set_xlabel('RA [ $^{\circ}$]')
ax.grid()
ax.set_ylim(-2e-3, 2e-3)
ax.invert_xaxis()
ax.legend()
outfile = 'IC86-{}_proj-RI-all_nside-{}_smooth-{:0.1f}.png'.format(years_str, n_side, smooth_rad)
plt.savefig(os.path.join(comp.paths.figures_dir, outfile))
plt.show()
In [82]:
relint_light = anisotropy.get_map('relint', composition='light', **kwargs_relint)
relint_light_err = anisotropy.get_map('relerr', composition='light', **kwargs_relint)
ri_light, ri_light_err, ra, ra_err = anisotropy.get_proj_relint(relint_light, relint_light_err, n_bins=30)
In [83]:
relint_heavy = anisotropy.get_map('relint', composition='heavy', **kwargs_relint)
relint_heavy_err = anisotropy.get_map('relerr', composition='heavy', **kwargs_relint)
ri_heavy, ri_heavy_err, ra, ra_err = anisotropy.get_proj_relint(relint_heavy, relint_heavy_err, n_bins=30)
In [84]:
popt_light, perr_light, chi2 = anisotropy.get_proj_fit_params(np.deg2rad(ra), ri_light, sigmay=ri_light_err, l=3)
popt_heavy, perr_heavy, chi2 = anisotropy.get_proj_fit_params(np.deg2rad(ra), ri_heavy, sigmay=ri_heavy_err, l=3)
In [85]:
fig, ax = plt.subplots()
ax.errorbar(ra, ri_light, yerr=ri_light_err, marker='.', ls='None', label='light')
ax.errorbar(ra, ri_heavy, yerr=ri_heavy_err, marker='.', ls='None', label='heavy')
ax.axhline(0, marker='None', ls='-.', c='k')
ax.set_ylabel('$\mathrm{\langle RI \\rangle }$')
ax.set_xlabel('RA [ $^{\circ}$]')
ax.legend()
ax.grid()
ax.invert_xaxis()
ax.set_ylim(-2e-3, 2e-3)
outfile = 'IC86-{}_proj-RI-comps_nside-{}_smooth-{:0.1f}.png'.format(years_str, n_side, smooth_rad)
plt.savefig(os.path.join(comp.paths.figures_dir, outfile))
plt.show()
In [86]:
fig, ax = plt.subplots()
ax.errorbar(ra, ri_light, yerr=ri_light_err, marker='.', ls='None', label='light')
ax.errorbar(ra, ri_heavy, yerr=ri_heavy_err, marker='.', ls='None', label='heavy')
ax.errorbar(ra, ri_all, yerr=ri_all_err, marker='.', ls='None', label='all')
ax.axhline(0, marker='None', ls='-.', c='k')
ax.set_ylabel('$\mathrm{\langle RI \\rangle }$')
ax.set_xlabel('RA [ $^{\circ}$]')
ax.legend()
ax.grid()
ax.invert_xaxis()
ax.set_ylim(-2e-3, 2e-3)
outfile = 'IC86-{}_proj-RI_nside-{}_smooth-{:0.1f}.png'.format(years_str, n_side, smooth_rad)
plt.savefig(os.path.join(comp.paths.figures_dir, outfile))
plt.show()
In [87]:
fig, ax = plt.subplots()
ax.errorbar(ra, ri_light, yerr=ri_light_err, marker='.', ls='None', label='light')
ax.errorbar(ra, ri_heavy, yerr=ri_heavy_err, marker='.', ls='None', label='heavy')
ax.errorbar(ra, ri_all, yerr=ri_all_err, marker='.', ls='None', label='all')
ax.axhline(0, marker='None', ls='-.', c='k')
ax.plot(ra, anisotropy.cos_fit_func(np.deg2rad(ra), *popt_light), color='C0', marker='None')
ax.plot(ra, anisotropy.cos_fit_func(np.deg2rad(ra), *popt_heavy), color='C1', marker='None')
ax.set_ylabel('$\mathrm{\langle RI \\rangle }$')
ax.set_xlabel('RA [ $^{\circ}$]')
ax.legend()
ax.grid()
ax.invert_xaxis()
ax.set_ylim(-2e-3, 2e-3)
outfile = 'IC86-{}_proj-RI-fit_nside-{}_smooth-{:0.1f}.png'.format(years_str, n_side, smooth_rad)
plt.savefig(os.path.join(comp.paths.figures_dir, outfile))
plt.show()
In [88]:
amp_light = popt_light[1]
amp_heavy = popt_heavy[1]
amp_err_light = perr_light[1]
amp_err_heavy = perr_heavy[1]
phase_light = np.rad2deg(popt_light[2])
phase_heavy = np.rad2deg(popt_heavy[2])
phase_err_light = np.rad2deg(perr_light[2])
phase_err_heavy = np.rad2deg(perr_heavy[2])
In [89]:
fig, ax = plt.subplots()
ax.errorbar(ra, ri_light, yerr=ri_light_err, marker='.', ls='None', label='light')
ax.errorbar(ra, ri_heavy, yerr=ri_heavy_err, marker='.', ls='None', label='heavy')
ax.errorbar(ra, ri_all, yerr=ri_all_err, marker='.', ls='None', label='all')
ax.axhline(0, marker='None', ls='-.', c='k')
ax.plot(ra, anisotropy.cos_fit_func(np.deg2rad(ra), *popt_light[:3]), color='C0', marker='None')
ax.plot(ra, anisotropy.cos_fit_func(np.deg2rad(ra), *popt_heavy[:3]), color='C1', marker='None')
ax.set_ylabel('$\mathrm{\langle RI \\rangle }$')
ax.set_xlabel('RA [ $^{\circ}$]')
ax.legend()
ax.grid()
ax.invert_xaxis()
ax.set_ylim(-2e-3, 2e-3)
light_amp_str = 'light amp = {:0.1g} +/- {:0.1g}'.format(amp_light, amp_err_light)
light_phase_str = 'light phase = {:0.2f} {} +/- {:0.2f} {}'.format(phase_light, '$^{\circ}$', phase_err_light, '$^{\circ}$')
ax.text(350, 0.00225, light_amp_str + '\n' + light_phase_str)
heavy_amp_str = 'heavy amp = {:0.1g} +/- {:0.1g}'.format(amp_heavy, amp_err_heavy)
heavy_phase_str = 'heavy phase = {:0.2f} {} +/- {:0.2f} {}'.format(phase_heavy, '$^{\circ}$', phase_err_heavy, '$^{\circ}$')
ax.text(175, 0.00225, heavy_amp_str + '\n' + heavy_phase_str)
outfile = 'IC86-{}_proj-RI-dipolefit_nside-{}_smooth-{:0.1f}.png'.format(years_str, n_side, smooth_rad)
plt.savefig(os.path.join(comp.paths.figures_dir, outfile))
plt.show()
[ back to top ]
In [66]:
config = ['IC86.2012', 'IC86.2013', 'IC86.2014', 'IC86.2015']
n_side = 64
smooth_rad = 20.0
scale = 3
decmax = -55
figures_dir = '/home/jbourbeau/public_html/figures/composition-anisotropy/cross-check-random-split'
years_str = '2012-2015'
In [67]:
rel_int_0 = anisotropy.get_map('relint', config=config, composition='random_0',
smooth=smooth_rad, scale=scale, decmax=decmax)
rel_int_1 = anisotropy.get_map('relint', config=config, composition='random_1',
smooth=smooth_rad, scale=scale, decmax=decmax)
print(rel_int_0.max())
print(rel_int_0[rel_int_0 != hp.UNSEEN].min())
print(rel_int_1.max())
print(rel_int_1[rel_int_1 != hp.UNSEEN].min())
In [68]:
title = 'Relative Intensity [$\\times \ 10^{}$]'.format('{'+str(-scale)+'}')
fig, ax = anisotropy.plot_skymap(rel_int_0, smooth=smooth_rad, decmax=decmax, scale=scale,
color_palette='RdBu_r', symmetric=True, llabel='IC86 2012-2013',
cbar_title=title, cbar_min=-2.5, cbar_max=2.5, polar=True)
outfile = 'IC86-{}_relint_{}_nside-{}_smooth-{:0.1f}.png'.format(years_str, 'random_0', n_side, smooth_rad)
plt.savefig(os.path.join(figures_dir, outfile))
In [69]:
title = 'Relative Intensity [$\\times \ 10^{}$]'.format('{'+str(-scale)+'}')
fig, ax = anisotropy.plot_skymap(rel_int_1, smooth=smooth_rad, decmax=decmax, scale=scale,
color_palette='RdBu_r', symmetric=True, llabel='IC86 2012-2013',
cbar_title=title, cbar_min=-2.5, cbar_max=2.5, polar=True)
outfile = 'IC86-{}_relint_{}_nside-{}_smooth-{:0.1f}.png'.format(years_str, 'random_1', n_side, smooth_rad)
plt.savefig(os.path.join(figures_dir, outfile))
In [70]:
sig_0 = anisotropy.get_map('sig', composition='random_0', config=config, smooth=smooth_rad, decmax=decmax)
sig_1 = anisotropy.get_map('sig', composition='random_1', config=config, smooth=smooth_rad, decmax=decmax)
In [71]:
print(sig_0.max())
print(sig_0[sig_0 != hp.UNSEEN].min())
print(sig_1.max())
print(sig_1[sig_1 != hp.UNSEEN].min())
In [72]:
title = 'Significance [$\mathrm{\sigma_{LM}}$]'
fig, ax = anisotropy.plot_skymap(sig_0, smooth=smooth_rad, decmax=decmax, scale=scale,
color_palette='RdBu_r', symmetric=True, llabel='IC86 2012-2013',
cbar_title=title, cbar_min=-3.5, cbar_max=3.5, polar=True)
outfile = 'IC86-{}_sig_{}_nside-{}_smooth-{:0.1f}.png'.format(years_str, 'random_0', n_side, smooth_rad)
plt.savefig(os.path.join(figures_dir, outfile))
In [73]:
title = 'Significance [$\mathrm{\sigma_{LM}}$]'
fig, ax = anisotropy.plot_skymap(sig_1, smooth=smooth_rad, decmax=decmax, scale=scale,
color_palette='RdBu_r', symmetric=True, llabel='IC86 2012-2013',
cbar_title=title, cbar_min=-3.5, cbar_max=3.5, polar=True)
outfile = 'IC86-{}_sig_{}_nside-{}_smooth-{:0.1f}.png'.format(years_str, 'random_1', n_side, smooth_rad)
plt.savefig(os.path.join(figures_dir, outfile))
In [74]:
relint = anisotropy.get_map('relint', config=config, composition='all',
smooth=smooth_rad, scale=None, decmax=decmax)
relint_err = anisotropy.get_map('relerr', config=config, composition='all',
smooth=smooth_rad, scale=None, decmax=decmax)
In [75]:
ri_all, ri_all_err, ra, ra_err = anisotropy.get_proj_relint(relint, relint_err)
In [76]:
fig, ax = plt.subplots()
ax.errorbar(ra, ri_all, yerr=ri_all_err, marker='.', ls='None')
ax.axhline(0, marker='None', ls='-.', c='k')
ax.set_ylabel('$\mathrm{\langle RI \\rangle }$')
ax.set_xlabel('RA [ $^{\circ}$]')
ax.grid()
ax.invert_xaxis()
plt.show()
In [91]:
relint_0 = anisotropy.get_map('relint', config=config, composition='random_0',
smooth=smooth_rad, scale=None, decmax=decmax)
relint_err_0 = anisotropy.get_map('relerr', config=config, composition='random_0',
smooth=smooth_rad, scale=None, decmax=decmax)
ri_0, ri_err_0, ra, ra_err = anisotropy.get_proj_relint(relint_0, relint_err_0)
In [92]:
relint_1 = anisotropy.get_map('relint', config=config, composition='random_1',
smooth=smooth_rad, scale=None, decmax=decmax)
relint_err_1 = anisotropy.get_map('relerr', config=config, composition='random_1',
smooth=smooth_rad, scale=None, decmax=decmax)
ri_1, ri_err_1, ra, ra_err = anisotropy.get_proj_relint(relint_1, relint_err_1)
In [93]:
popt_0, perr, chi2 = anisotropy.get_proj_fit_params(np.deg2rad(ra), ri_0, sigmay=ri_err_0)
popt_1, perr, chi2 = anisotropy.get_proj_fit_params(np.deg2rad(ra), ri_1, sigmay=ri_err_1)
In [94]:
fig, ax = plt.subplots()
ax.errorbar(ra, ri_0, yerr=ri_err_0, marker='.', ls='None', label='Subset 0')
ax.errorbar(ra, ri_1, yerr=ri_err_1, marker='.', ls='None', label='Subset 1')
ax.axhline(0, marker='None', ls='-.', c='k')
ax.set_ylabel('$\mathrm{\langle RI \\rangle }$')
ax.set_xlabel('RA [ $^{\circ}$]')
ax.legend()
ax.grid()
ax.invert_xaxis()
plt.show()
In [95]:
fig, ax = plt.subplots()
ax.errorbar(ra, ri_0, yerr=ri_err_0, marker='.', ls='None', label='Subset 0')
ax.errorbar(ra, ri_1, yerr=ri_err_1, marker='.', ls='None', label='Subset 1')
ax.errorbar(ra, ri_all, yerr=ri_all_err, marker='.', ls='None', label='all')
ax.axhline(0, marker='None', ls='-.', c='k')
ax.set_ylabel('$\mathrm{\langle RI \\rangle }$')
ax.set_xlabel('RA [ $^{\circ}$]')
ax.legend()
ax.grid()
ax.invert_xaxis()
plt.show()
In [96]:
fig, ax = plt.subplots()
ax.errorbar(ra, ri_0, yerr=ri_err_0, marker='.', ls='None', label='Subset 0')
ax.errorbar(ra, ri_1, yerr=ri_err_1, marker='.', ls='None', label='Subset 1')
ax.errorbar(ra, ri_all, yerr=ri_all_err, marker='.', ls='None', label='all')
ax.axhline(0, marker='None', ls='-.', c='k')
ax.plot(ra, anisotropy.cos_fit_func(np.deg2rad(ra), *popt_0), color='C0', marker='None')
ax.plot(ra, anisotropy.cos_fit_func(np.deg2rad(ra), *popt_1), color='C1', marker='None')
ax.set_ylabel('$\mathrm{\langle RI \\rangle }$')
ax.set_xlabel('RA [ $^{\circ}$]')
ax.legend()
ax.grid()
ax.invert_xaxis()
plt.show()
In [97]:
fig, ax = plt.subplots()
ax.errorbar(ra, ri_0, yerr=ri_err_0, marker='.', ls='None', label='Subset 0')
ax.errorbar(ra, ri_1, yerr=ri_err_1, marker='.', ls='None', label='Subset 1')
ax.errorbar(ra, ri_all, yerr=ri_all_err, marker='.', ls='None', label='all')
ax.axhline(0, marker='None', ls='-.', c='k')
ax.plot(ra, anisotropy.cos_fit_func(np.deg2rad(ra), *popt_0[:3]), color='C0', marker='None')
ax.plot(ra, anisotropy.cos_fit_func(np.deg2rad(ra), *popt_1[:3]), color='C1', marker='None')
ax.set_ylabel('$\mathrm{\langle RI \\rangle }$')
ax.set_xlabel('RA [ $^{\circ}$]')
ax.legend()
ax.grid()
ax.invert_xaxis()
ax.set_ylim(-2e-3, 2e-3)
outfile = 'IC86-{}_proj-RI-dipolefit_random_nside-{}_smooth-{:0.1f}.png'.format(years_str, n_side, smooth_rad)
plt.savefig(os.path.join(figures_dir, outfile))
plt.show()
In [98]:
'{:g}'.format(popt_0[1]), '{:g}'.format(popt_1[1])
Out[98]:
In [99]:
np.rad2deg(popt_0[2]), np.rad2deg(popt_1[2])
Out[99]:
In [ ]:
In [3]:
# Load full DataFrame for config
mypaths = comp.get_paths()
df_file = os.path.join(mypaths.comp_data_dir, 'IC86.2013' + '_data',
'anisotropy_dataframe.hdf')
with pd.HDFStore(df_file) as store:
data_df = store['dataframe']
In [6]:
data_df.head()
Out[6]:
In [5]:
n_side = 64
n_pix = hp.nside2npix(n_side)
bdt_score_skymap = np.zeros(n_pix)
for idx, row in data_df.iterrows():
local_zenith = row['lap_zenith']
local_azimuth = row['lap_azimuth']
local_time = row['start_time_mjd']
ra, dec = astro.dir_to_equa(local_zenith, local_azimuth, local_time)
hp_theta, hp_phi = anisotropy.equatorial_to_healpy(ra, dec)
pix = hp.ang2pix(n_side, hp_theta, hp_phi)
bdt_score_skymap[pix] += row['score']
In [80]:
fig, axarr = plt.subplots(1, 2)
In [82]:
ax = axarr[0]
In [ ]:
ax.s
In [ ]: