GW Example

v2 – Refactor with localization

v3 – Refactor with PATH

v3.1 – Recalculate with proper prior normalization

[1]:

%matplotlib notebook

[2]:

# imports
from importlib import reload
import os
from importlib.resources import files as resource_files
import numpy as np

import healpy as hp
import pandas

from astropy.io import fits
from astropy.table import Table
from astropy.coordinates import SkyCoord

from astropath import bayesian
from astropath import localization
from astropath import path

LIGO example – GW170817

[3]:

lfile = os.path.join(str(resource_files('astropath').joinpath('data')), 'gw_examples',
                     'GW170817_skymap.fits.gz')

[4]:

gw170817 = hp.read_map(lfile)

[5]:

header = fits.open(lfile)[1].header
header

[5]:

XTENSION= 'BINTABLE'           / binary table extension
BITPIX  =                    8 / array data type
NAXIS   =                    2 / number of array dimensions
NAXIS1  =                   32 / length of dimension 1
NAXIS2  =             12582912 / length of dimension 2
PCOUNT  =                    0 / number of group parameters
GCOUNT  =                    1 / number of groups
TFIELDS =                    4 / number of table fields
TTYPE1  = 'PROB    '
TFORM1  = 'D       '
TUNIT1  = 'pix-1   '
TTYPE2  = 'DISTMU  '
TFORM2  = 'D       '
TUNIT2  = 'Mpc     '
TTYPE3  = 'DISTSIGMA'
TFORM3  = 'D       '
TUNIT3  = 'Mpc     '
TTYPE4  = 'DISTNORM'
TFORM4  = 'D       '
TUNIT4  = 'Mpc-2   '
PIXTYPE = 'HEALPIX '           / HEALPIX pixelisation
ORDERING= 'NESTED  '           / Pixel ordering scheme: RING, NESTED, or NUNIQ
COORDSYS= 'C       '           / Ecliptic, Galactic or Celestial (equatorial)
NSIDE   =                 1024 / Resolution parameter of HEALPIX
INDXSCHM= 'IMPLICIT'           / Indexing: IMPLICIT or EXPLICIT
OBJECT  = 'GW170817'           / Unique identifier for this event
DATE-OBS= '2017-08-17T12:41:04.429464' / UTC date of the observation
MJD-OBS =    57982.52852348908 / modified Julian date of the observation
DATE    = '2019-05-08T21:51:56.613850' / UTC date of file creation
CREATOR = 'ligo-skymap-from-samples' / Program that created this file
ORIGIN  = 'LIGO/Virgo'         / Organization responsible for this FITS file
DISTMEAN=    38.03408225450563 / Posterior mean distance (Mpc)
DISTSTD =    7.499686641911211 / Posterior standard deviation of distance (Mpc)
VCSVERS = 'ligo.skymap 0.0.15' / Software version
VCSREV  = '65fc6500a1e117fec2e27550a8e9b10c9792ffca' / Software revision (Git)
DATE-BLD= '2018-09-04T14:19:20' / Software build date
HISTORY
HISTORY Generated by running the following script:
HISTORY ligo-skymap-from-samples --seed 170817 --outdir . --objid GW170817 --fit
HISTORY soutname MCMC_TF2_LowSpin_AllSky.fits.gz MCMC_TF2_LowSpin_AllSky_post.da
HISTORY t

[6]:

hp.mollview(
    gw170817,
    coord=["C"],
    title="GW170817",
    #unit="mK",
    #norm="hist",
    #min=-1,
    #max=1,
)
hp.graticule()

0.0 180.0 -180.0 180.0

Galaxies

https://vizier.u-strasbg.fr/viz-bin/VizieR?-source=VII/275

[7]:

galfile = os.path.join(str(resource_files('astropath').joinpath('data')), 'gw_examples',
                     'GW170817_galaxies.csv')

[8]:

cut_galaxies = pandas.read_csv(galfile, index_col=0)

[9]:

cut_galaxies.head()

[9]:

	RAJ2000	DEJ2000	Dist	Bmag	BmagHyp	ImagHyp	modzHyp	mod0Hyp	logd25Hyp	logr25Hyp	...	BmagGWGC	maj	min	PAGWGC	DistGWGC	Kmag2MPZ	Bmag2MPZ	zsp2MPZ	zph2MPZ
667013	357.169220	-2.362345	92.028	15.33	15.81	14.02	34.92	NaN	0.85	0.19	...	15.33	0.708	0.457	136.2	92.028	12.761	15.622	0.02219	0.031973
667018	150.132019	-31.244976	31.583	12.60	12.95	10.66	32.62	NaN	1.31	0.15	...	12.60	2.042	1.446	179.9	31.583	8.792	12.490	0.00813	0.008867
667030	76.325119	-9.147248	44.347	12.53	13.03	10.86	33.33	NaN	1.43	0.35	...	12.53	1.816	0.890	112.0	44.347	8.818	12.735	0.01090	0.010972
667044	338.249908	-20.094194	80.264	14.13	14.43	13.35	34.57	NaN	1.05	0.03	...	14.13	0.830	0.772	NaN	80.264	12.792	14.575	0.01905	0.011343
667076	138.796539	29.252758	88.167	16.19	17.13	NaN	34.83	NaN	0.72	0.46	...	16.19	0.525	0.182	2.4	88.167	13.127	16.320	0.02080	0.062994

5 rows × 31 columns

Coordinates

[10]:

cut_gal_coord = SkyCoord(ra=cut_galaxies.RAJ2000, dec=cut_galaxies.DEJ2000, unit='deg')

[11]:

ngc_4993 = SkyCoord('13h09m47.706s -23d23m01.79s', frame='icrs')

[12]:

np.min(ngc_4993.separation(cut_gal_coord).to('arcmin'))

[12]:

$6.96812e-05\mathrm{{}^{\prime}}$

[13]:

np.argmin(ngc_4993.separation(cut_gal_coord))

[13]:

[14]:

cut_galaxies.iloc[11]

[14]:

RAJ2000      197.448776
DEJ2000      -23.383831
Dist          33.806000
Bmag          12.870000
BmagHyp       13.450000
ImagHyp       11.310000
modzHyp       33.130000
mod0Hyp             NaN
logd25Hyp      1.190000
logr25Hyp      0.060000
logdcHyp       1.260000
PAHyp        173.200000
BmagHypC      12.870000
ImagHypC      11.080000
U-BHypC             NaN
B-VHypC             NaN
Jmag2         10.292000
Hmag2          9.599000
Kmag2          9.330000
a_b            0.900000
PAK          -20.000000
BmagGWGC      12.870000
maj            0.932000
min            0.792000
PAGWGC       173.200000
DistGWGC      33.806000
Kmag2MPZ       9.285000
Bmag2MPZ      13.160000
zsp2MPZ        0.009680
zph2MPZ        0.012307
Flag           0.000000
Name: 667145, dtype: float64

PATH time

Instantiate

[15]:

Path = path.PATH()

Candidates

[16]:

Path.init_candidates(cut_galaxies.RAJ2000.values,
                     cut_galaxies.DEJ2000.values,
                     cut_galaxies.maj.values,
                     mag=cut_galaxies.Bmag.values)

Priors

[18]:

# Candidates
Path.init_cand_prior('inverse', P_U=0.)
# Offsets
Path.init_theta_prior('exp', 6., scale=1.)

Localization

[19]:

Path.init_localization('healpix',
               healpix_data=gw170817,
               healpix_nside=header['NSIDE'],
               healpix_ordering='NESTED',
               healpix_coord='C')

Calculate Priors

[20]:

P_O = Path.calc_priors()

Calculate Posteriors

[21]:

P_Ox, P_Ux = Path.calc_posteriors('local', box_hwidth=30.)

[22]:

Path.candidates.sort_values('P_Ox', ascending=False)

[22]:

	ra	dec	ang_size	mag	P_O	p_xO	P_Ox
11	197.448776	-23.383831	0.932	12.87	0.006478	2.798955e-04	9.999998e-01
144	193.861099	-20.892786	0.776	15.73	0.000128	3.154919e-09	2.219717e-07
64	194.469193	-22.638481	1.047	15.52	0.000167	3.191201e-11	2.942908e-09
347	206.735901	-39.884365	0.757	14.88	0.000388	2.935595e-16	6.282478e-14
297	192.973251	-12.604387	0.468	16.27	0.000064	1.317743e-17	4.682630e-16
...	...	...	...	...	...	...	...
132	180.701767	21.145529	0.871	16.00	0.000090	0.000000e+00	0.000000e+00
131	110.026970	53.010151	0.757	14.28	0.000875	0.000000e+00	0.000000e+00
129	140.950027	2.112878	0.630	15.07	0.000301	0.000000e+00	0.000000e+00
128	21.381145	32.136272	1.582	13.63	0.002167	0.000000e+00	0.000000e+00
389	171.504303	1.984145	0.775	14.04	0.001219	0.000000e+00	0.000000e+00

390 rows × 7 columns

F’ing A

Multi order

https://emfollow.docs.ligo.org/userguide/tutorial/multiorder_skymaps.html

[28]:

mofile = os.path.join(str(resource_files('astropath').joinpath('data')), 'gw_examples',
                     'GW190814_PublicationSamples.multiorder.fits')