KDE PDFs Via GenericPDFRatio#
This notebook will walk you through using the generic PDF implementation in csky to implement the KDE spatial PDFs that were used in the most recent NorthernTracks time integrated analysis. In addition to reproducing the time integrated results for NGC 1068, we’ll additionally show how this can be combined with the untriggered flare likelihood, allowing the KDE spatial PDF description to be used for transient analyses
[1]:
#Imports we'll need
import numpy as np
import csky as cy
import os
import photospline as psp
from csky import cext, hyp
import time
We’ll create our function that calculates the generic PDF ratio, as well as some helper functions to evaluate the provided KDEs and mask out events that are far away (to save computation time). We also need to load up the spline tables that the improved point source people have generated
[2]:
base_path='/data/ana/analyses/northern_tracks/version-005-p00/KDE_PDFs_v007/IC86_pass2/'
spatial_pdf = psp.SplineTable(os.path.join(base_path, 'sig_E_psi_photospline_v006_4D.fits'))
energy_pdf = psp.SplineTable(os.path.join(base_path, 'E_dec_photospline_v006_3D.fits'))
bkg_pdf = psp.SplineTable(os.path.join(base_path, 'bg_2d_photospline.fits'))
def GreatCircleDistance(ra_1, dec_1, ra_2, dec_2):
'''Compute the great circle distance between two events'''
'''All coordinates must be given in radians'''
delta_dec = np.abs(dec_1 - dec_2)
delta_ra = np.abs(ra_1 - ra_2)
x = (np.sin(delta_dec / 2.))**2. + np.cos(dec_1) *\
np.cos(dec_2) * (np.sin(delta_ra / 2.))**2.
return 2. * np.arcsin(np.sqrt(x))
def spatial_kde(log_psi, log_ang_err, log_energy, gamma):
'''Evaluates the spatial kde for an event and spectral index hypothesis'''
return spatial_pdf.evaluate_simple([log_ang_err, log_energy, log_psi, np.full(len(log_psi), gamma)], 0)
def energy_kde(log_energy, src_dec, gamma):
'''Evaluates the energy kde for an event and spectral index hypothesis'''
sin_src_dec = np.sin(src_dec)
return energy_pdf.evaluate_simple([log_energy, np.full(len(log_energy), sin_src_dec), np.full(len(log_energy), gamma)], 0)
def calc_sb_ratio(ev_ra, ev_dec, ev_angErr, ev_logE, src, gamma, **kwargs):
'''Calculates the (space x energy) S/B ratio using spatial and energy kde pdfs'''
'''There are 2 cases, one where the mask is provided by the trial runner args, and one
where no mask has been provided'''
'''This is the function that we will be pointing our trial runner to to use as a "generic" pdf ratio'''
if '_mask' in kwargs:
mask = kwargs['_mask']
antimask = mask[mask==False]
src_ra = src['ra']
src_dec = src['dec']
ev_ra_mask = ev_ra[mask]
ev_dec_mask = ev_dec[mask]
ev_angErr_mask = ev_angErr[mask]
ev_logE_mask = ev_logE[mask]
psi = GreatCircleDistance(ev_ra_mask, ev_dec_mask, src_ra, src_dec)
log_psi = np.log10(psi)
log_ang_err = np.log10(ev_angErr_mask)
log_energy = ev_logE_mask
masked_ev_dec = ev_dec_mask
spdfs = spatial_kde(log_psi, log_ang_err, log_energy, gamma)
spatial_norm = np.log(10) * psi * np.sin(psi)
spdfs/= spatial_norm
epdfs = energy_kde(log_energy, src_dec, gamma)
bg_pdfs = bkg_pdf.evaluate_simple([log_energy, np.sin(masked_ev_dec)])
sb_ratio_close = spdfs*epdfs/bg_pdfs
far_evts_sb = np.zeros(len(antimask))
sb_ratio = np.zeros(len(ev_ra))
sb_ratio[mask] = sb_ratio_close
else:
src_ra = src['ra']
src_dec = src['dec']
psi = GreatCircleDistance(ev_ra, ev_dec, src_ra, src_dec)
log_psi = np.log10(psi)
log_ang_err = np.log10(ev_angErr)
log_energy = ev_logE
masked_ev_dec = ev_dec
spdfs = spatial_kde(log_psi, log_ang_err, log_energy, gamma)
spatial_norm = np.log(10) * psi* np.sin(psi)
spdfs/= spatial_norm
epdfs = energy_kde(log_energy, src_dec, gamma)
bg_pdfs = bkg_pdf.evaluate_simple([log_energy, np.sin(masked_ev_dec)])
sb_ratio = spdfs*epdfs/bg_pdfs
return sb_ratio
def dist_mask(arr, src, cut_deg):
'''A simple distance mask that masks out events (arr) farther than cut_deg away from src'''
out = np.zeros(len(arr), dtype=bool)
ras = arr['ra']
decs = arr['dec']
psi = GreatCircleDistance(ras, decs, src.ra, src.dec)
out = psi < np.radians(cut_deg)
return out
Do the usual setup of our analysis object:
[3]:
#Set up our analysis objects to use NTv5
ana_dir = cy.utils.ensure_dir('/data/ana/analyses/')
ana = cy.get_analysis(cy.selections.repo, 'version-005-p00', cy.selections.NTDataSpecs.nt_v5[0:1], dir=ana_dir, min_sigma=0.0)
Setting up Analysis for:
NT86v5
Setting up NT86v5...
Reading /data/ana/analyses/northern_tracks/version-005-p00/IC86_pass2_MC.npy ...
Reading /data/ana/analyses/northern_tracks/version-005-p00/IC86_2011_exp.npy ...
Reading /data/ana/analyses/northern_tracks/version-005-p00/IC86_2012_exp.npy ...
Reading /data/ana/analyses/northern_tracks/version-005-p00/IC86_2013_exp.npy ...
Reading /data/ana/analyses/northern_tracks/version-005-p00/IC86_2014_exp.npy ...
Reading /data/ana/analyses/northern_tracks/version-005-p00/IC86_2015_exp.npy ...
Reading /data/ana/analyses/northern_tracks/version-005-p00/IC86_2016_exp.npy ...
Reading /data/ana/analyses/northern_tracks/version-005-p00/IC86_2017_exp.npy ...
Reading /data/ana/analyses/northern_tracks/version-005-p00/IC86_2018_exp.npy ...
Reading /data/ana/analyses/northern_tracks/version-005-p00/IC86_2019_exp.npy ...
Reading /data/ana/analyses/northern_tracks/version-005-p00/GRL/IC86_2011_exp.npy ...
Reading /data/ana/analyses/northern_tracks/version-005-p00/GRL/IC86_2012_exp.npy ...
Reading /data/ana/analyses/northern_tracks/version-005-p00/GRL/IC86_2013_exp.npy ...
Reading /data/ana/analyses/northern_tracks/version-005-p00/GRL/IC86_2014_exp.npy ...
Reading /data/ana/analyses/northern_tracks/version-005-p00/GRL/IC86_2015_exp.npy ...
Reading /data/ana/analyses/northern_tracks/version-005-p00/GRL/IC86_2016_exp.npy ...
Reading /data/ana/analyses/northern_tracks/version-005-p00/GRL/IC86_2017_exp.npy ...
Reading /data/ana/analyses/northern_tracks/version-005-p00/GRL/IC86_2018_exp.npy ...
Reading /data/ana/analyses/northern_tracks/version-005-p00/GRL/IC86_2019_exp.npy ...
Energy PDF Ratio Model...
* gamma = 4.0000 ...
Signal Acceptance Model...
* gamma = 4.0000 ...
Done.
Now we can set up our first trial runner. This trial runner will be used to generate trials, and we’ll create another one to actually evaluate the likelihood. We do this to ensure that
I’ll be setting this up to point to NGC 1068 (ra=40.667, dec =-0.0069) so we can check against the results that have already been obtained using SkyLLH. In our trial runner setup, we choose our space pdf ratio to be 'generic', and point the function to be used to calc_sb_ratio. features is a dictionary of the args that calc_sb_ratio will need, and fits are the parameters we’d like the likelihood to fit for (in addition to ns). Note that we turn energy=False here,
since our pdf ratio function calc_sb_ratio will be handling both the space and energy components of the likelihood.
[4]:
srcs = cy.sources(np.radians(40.667), np.radians(-0.0069), name=0)
gamma = 2.0
tr = cy.get_trial_runner(ana=ana, src=srcs, space='generic', func=calc_sb_ratio,
features={'ev_ra':'ra', 'ev_dec':'dec', 'ev_angErr':'sigma', 'ev_logE':'log10energy'},
fits={'gamma':(1.0,2.0,3.0,4.0)}, extra_keep=['dec', 'sigma', 'sindec', 'event'],
energy=False,
cut_n_sigma=np.inf,
concat_evs=True,
flux=hyp.PowerLawFlux(gamma))
Now we can generate and evaluate a trial. Note that here we introduce a distance mask so we don’t have to evaluate the KDE splines for events farther than 15 degrees. This significantly speeds up the calculation.
[5]:
trialseed = 0
trial = tr.get_one_trial(n_sig=0, TRUTH=True, seed=trialseed)
for src in srcs:
t1 = time.time()
trial_mask = dist_mask(trial[0][0][0], srcs, 15.)
fit = tr.get_one_fit_from_trial(trial,flat=False,_masks=[trial_mask], TRUTH=True)
print(fit)
t2 = time.time()
print("fit took", t2-t1, "seconds")
(29.204403027486393, {'gamma': 3.1586663282907503, 'ns': 76.96023713945011})
fit took 1.4317007064819336 seconds
Time Dependent fits with KDE PDFs#
We can combine the new KDE spatial/energy PDFs with the untriggered flare likelihood by using a MultiflareTrialRunner, set up like so:
[8]:
mtr = cy.conf.get_multiflare_trial_runner(ana=ana, src=srcs, threshold=1000., space='generic', func=calc_sb_ratio,
features={'ev_ra':'ra', 'ev_dec':'dec', 'ev_angErr':'sigma', 'ev_logE':'log10energy'},
fits={'gamma':(1.0,2.0,3.0,4.0)}, extra_keep=['dec', 'sigma', 'event'],
energy=False,
concat_evs=True,
flux=hyp.PowerLawFlux(gamma),
muonflag=True)
[ ]:
trialseed = 0
mtr_trial = mtr.get_one_trial(seed=trialseed, TRUTH=False)
mtr_fit = mtr.get_one_fit_from_trial(mtr_trial, flat=False, _fmin_epsilon=1e-3)
print("mtr_fit is", mtr_fit)
finished making trial
type of evaluator is <class 'csky.pdf.SpaceTimePDFRatioEvaluator'>
this evaluator does not have separate space and energy components
/cvmfs/icecube.opensciencegrid.org/py3-v4.1.1/RHEL_7_x86_64/lib/python3.7/site-packages/ipykernel_launcher.py:37: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
/cvmfs/icecube.opensciencegrid.org/py3-v4.1.1/RHEL_7_x86_64/lib/python3.7/site-packages/ipykernel_launcher.py:38: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
/cvmfs/icecube.opensciencegrid.org/py3-v4.1.1/RHEL_7_x86_64/lib/python3.7/site-packages/ipykernel_launcher.py:39: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
/cvmfs/icecube.opensciencegrid.org/py3-v4.1.1/RHEL_7_x86_64/lib/python3.7/site-packages/ipykernel_launcher.py:40: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
/cvmfs/icecube.opensciencegrid.org/py3-v4.1.1/RHEL_7_x86_64/lib/python3.7/site-packages/ipykernel_launcher.py:60: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.
There are 590 windows
[ ]: