Coverage for /opt/anaconda/lib/python3.7/site-packages/llama/files/coinc_significance/subthreshold.py : 19%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""
2``FileHandler`` class for calculating significance of joint events with
3subthreshold gravitational wave triggers.
5NOTE: equations after 8 increment up by 1. equations after 22 increment up by
62.
8@author: dvesk, stefco
9"""
11import logging
12from math import pi
13from numpy import (
14 zeros,
15 dot,
16)
17from llama import detectors
18from llama.filehandler import JSONFile
19from llama.files.i3 import IceCubeNeutrinoList
20from llama.files.lvc_skymap import (
21 LvcSkymapHdf5, LvcDistancesJson
22)
23from llama.files.skymap_info import SkymapInfo
24from llama.files.gracedb import LVCGraceDbEventData, PAstro
25from llama.files.healpix import (
26 HEALPixSkyMap,
27)
28from llama.utils import (
29 RemoteFileCacher,
30)
31from llama.files.slack import SlackReceiptLlama
32from llama.files.coinc_significance.utils import (
33 # ODDS_RATIO_KEYS,
34 r0,
35 load_neutrino_background,
36 p_value,
37 Pempgw,
38 t_overlap,
39 angle_list,
40 Aeff_skymap,
41 search_parameters,
42)
43from llama.files.lvalert_json import LVAlertJSON
45LOGGER = logging.getLogger(__name__)
47# post O3a SNR update
48GW_BG_RHOS = RemoteFileCacher(
49 'https://nyc3.digitaloceanspaces.com/llama/llama/objects/'
50 '239c1ce1ce0fe263653a84ecb219e1483f830829e1cd87361a7e6ddec22ba6bb'
51)
53def emptyskymap(val, skymap):
54 """A ``HEALPixSkyMap`` initialized with a uniform value.
56 Parameters
57 ----------
58 val : float or int
59 The value to assign each pixel.
60 skymap : HEALPixSkyMap
61 A skymap whose shape should be used as the template.
63 Returns
64 -------
65 skymap : HEALPixSkyMap
66 A skymap shaped like ``skymap`` but with all pixel values set to
67 ``val``.
68 """
69 return HEALPixSkyMap(val,nest=skymap.nest)
72def five(tgw, gw_skymap, ρgw, dist,neutrinolist, search_params, Eₙ,pempnu, cwb_Boolean):
73 """Calculate the signal likelihood (Eq. 5).
75 Parameters
76 ----------
77 tgw : float
78 GPS time of the GW trigger
79 gw_skymap : llama.files.healpix.skymap.HEALPixSkyMap
80 A HEALPix skymap for the GW specifying pixel values and locations.
81 ρgw : float
82 The signal-to-noise ratio (SNR) of the GW
83 dist : float
84 The probability-weighted average reconstructed distance of the GW event
85 neutrino : llama.files.i3.Neutrino
86 The neutrino candidate that is being compared to the GW skymap
87 search_params : IceCubeLvcSearchParameters
88 A collection of constant parameters used for this search.
89 Eₙ : float
90 The neutrino energy. This should be capped at the maximum energy of a
91 nearby background neutrino. This way the signal likelihood does not
92 drop (as it does for increasing energy) while the background stays the
93 same (since the background calculation always assumes at least one
94 neutrino).
95 cwb_Boolean : bool
96 Whether this is a CWB event or not.
98 Returns
99 -------
100 P_H_S1 : float
101 The signal hypothesis likelihood for a single neutrino:
103 AllskyIntegral(P(x_gw, x_nu | θ, H_s)*P(θ | H_s))
105 """
106 from scipy import integrate
107 from scipy.stats import poisson
108 count=-1
109 if(not(cwb_Boolean)):
110 def innerfive(Eₙ):
111 return poisson.pmf(1,expnu(Eₙ,dist,search_params)*angle_list()[int(neutrino.zenith/5)])*PEₙ(Eₙ,search_params)
112 summednuskymap=emptyskymap(0,gw_skymap)
113 for neutrino in neutrinolist:
114 count+=1
115 a=emptyskymap(t_overlap(tgw, neutrino.gps, search_params)/pempnu[count]*integrate.quad(innerfive,search_params.Eₙmin,search_params.Eₙmax)[0],gw_skymap)
116 #print('pempnu',pempnu[count])
117 summednuskymap += a*Aeff_skymap(gw_skymap,Eₙ[count])*gw_skymap.psf_gaussian(neutrino.ra, neutrino.dec, neutrino.sigma,normalize=False)
118 numerator = PEgw((dist*ρgw/search_params.k0)**2,search_params)*dist*((search_params.Egwmax)-(search_params.Egwmin))*search_params.k0**2/ρgw**2
119 else:
120 def innerfive(Eₙ,d):
121 return PEgw((d*ρgw/search_params.k0)**2,search_params)*d**2*poisson.pmf(1,expnu(Eₙ,d,search_params)*angle_list()[int(neutrino.zenith/5)])*PEₙ(Eₙ,search_params)
122 summednuskymap=emptyskymap(0,gw_skymap)
123 for neutrino in neutrinolist:
124 count+=1
125 a=emptyskymap(t_overlap(tgw, neutrino.gps, search_params)/pempnu[count]*integrate.nquad(innerfive,[[search_params.Eₙmin,search_params.Eₙmax],[0,r0(search_params.Egwmax, search_params)]])[0],gw_skymap)
126 #print('pempnu',pempnu[count])
127 summednuskymap += a*Aeff_skymap(gw_skymap,Eₙ[count])*gw_skymap.psf_gaussian(neutrino.ra, neutrino.dec, neutrino.sigma,normalize=False)
128 numerator = 1
129 denominator = ((search_params.tnuplus-search_params.tnuminus)*(search_params.tgwplus-search_params.tgwminus))
130 integral = sum(
131 fA(90-gw_skymap.dec, gw_skymap.ra, tgw) *
132 (gw_skymap *
133 summednuskymap
134 ).pixels
135 )
136 return (integral * numerator / denominator)
139def five2(tgw, gw_skymap, ρgw,dist, neutrino, neutrino_i, search_params, Eₙ,
140 Eₙ_i):
141 """Calculate the signal likelihood for 2 coincident neutrinos (modified Eq.
142 5)
144 Parameters
145 ----------
146 tgw : float
147 GPS time of the GW trigger
148 gw_skymap : llama.files.healpix.skymap.HEALPixSkyMap
149 A HEALPix skymap for the GW specifying pixel values and locations.
150 ρgw : float
151 The signal-to-noise ratio (SNR) of the GW
152 dist : float
153 The probability-weighted average reconstructed distance of the GW event
154 neutrino : llama.files.i3.Neutrino
155 The neutrino candidate that is being compared to the GW skymap as well
156 as ``neutrino_i``.
157 neutrino_i : llama.files.i3.Neutrino
158 The ith neutrino candidate that is being compared to the GW skymap as
159 well as ``neutrino``.
160 search_params : IceCubeLvcSearchParameters
161 A collection of constant parameters used for this search.
162 Eₙ : float
163 The neutrino energy. This should be capped at the maximum energy of a
164 nearby background neutrino. This way the signal likelihood does not
165 drop (as it does for increasing energy) while the background stays the
166 same (since the background calculation always assumes at least one
167 neutrino).
168 Eₙ_i : float
169 Same as Eₙ, but for the ith neutrino.
171 Returns
172 -------
173 P_H_S2 : float
174 The signal hypothesis likelihood for a double-neutrino detection:
176 AllskyIntegral(P(x_gw, x_nu, x_nu_i | θ, H_s)*P(θ | H_s))
178 """
179 from scipy import integrate
180 from scipy.stats import poisson
181 def innerfive2(Eₙ):
182 return poisson.pmf(2,expnu(Eₙ,dist,search_params)*angle_list()[int((neutrino.zenith+neutrino_i.zenith)/10)])*PEₙ(Eₙ,search_params)
183 numerator = PEgw((dist*ρgw/search_params.k0)**2,search_params)*dist**2*integrate.quad(innerfive2,search_params.Eₙmin,search_params.Eₙmax)[0]
184 denominator = ( (search_params.tnuplus-search_params.tnuminus)*(search_params.tgwplus-search_params.tgwminus))
185 integral = sum(
186 fA(90-gw_skymap.dec, gw_skymap.ra, tgw) *
187 (
188 gw_skymap* Aeff_skymap(gw_skymap,Eₙ)*Aeff_skymap(gw_skymap,Eₙ_i) * gw_skymap.psf_gaussian(neutrino.ra, neutrino.dec, neutrino.sigma, normalize=False)*gw_skymap.psf_gaussian(neutrino_i.ra, neutrino_i.dec, neutrino_i.sigma, normalize=False)
189 ).pixels)
190 return (integral *t_overlap(max(tgw, neutrino.gps, neutrino_i.gps),min(tgw, neutrino.gps, neutrino_i.gps),search_params) *
191 numerator / denominator)
194def fivewogw(gw_skymap,neutrinolist, search_params, Eₙ,pempnu):
195 """Calculate the coincidence likelihood (Eq. 58).
197 Parameters
198 ----------
200 neutrino : llama.files.i3.Neutrino
201 The neutrino candidate that is being compared to the GW skymap
202 search_params : IceCubeLvcSearchParameters
203 A collection of constant parameters used for this search.
204 Eₙ : float
205 The neutrino energy. This should be capped at the maximum energy of a
206 nearby background neutrino. This way the signal likelihood does not
207 drop (as it does for increasing energy) while the background stays the
208 same (since the background calculation always assumes at least one
209 neutrino).
211 Returns
212 -------
213 P_H_c1 : float
214 The coincidence hypothesis likelihood for a single neutrino:
216 AllskyIntegral(P(x_nu | θ, H_c)*P(θ | H_c))
218 """
219 from scipy import integrate
220 from scipy.stats import poisson
221 def innerfivewogw(Eₙ,dist):
222 return poisson.pmf(1,expnu(Eₙ,dist,search_params)*angle_list()[int(neutrino.zenith/5)])*PEₙ(Eₙ,search_params)*dist**2
223 count=-1
224 summednuskymap=emptyskymap(0,gw_skymap)
225 for neutrino in neutrinolist:
226 count+=1
227 a=emptyskymap(1/pempnu[count]*integrate.nquad(innerfivewogw,[[search_params.Eₙmin,search_params.Eₙmax],[0,r0(search_params.Egwmax, search_params)]])[0],gw_skymap)
228 summednuskymap += a *Aeff_skymap(gw_skymap,Eₙ[count])*gw_skymap.psf_gaussian(neutrino.ra, neutrino.dec, neutrino.sigma,normalize=True)
229 numerator = 1#integrate.nquad(innerfivewogw,[[search_params.Eₙmin,search_params.Eₙmax],[0,r0(search_params.Egwmax, search_params)]])[0]
230 #denominator = ((search_params.tnuplus-search_params.tnuminus)*(search_params.tgwplus-search_params.tgwminus))
231 integral = sum(summednuskymap.pixels)
232 return integral*numerator#(integral * t_overlap(tgw, neutrino.gps, search_params) *
233 #numerator / denominator)
236def five2wogw(neutrino, neutrino_i, search_params, Eₙ, Eₙ_i):
237 """NOT USED Calculate the coincidence likelihood for 2 coincident neutrinos (modified Eq.
238 58)
240 Parameters
241 ----------
243 neutrino : llama.files.i3.Neutrino
244 The neutrino candidate that is being compared to the GW skymap as well
245 as ``neutrino_i``.
246 neutrino_i : llama.files.i3.Neutrino
247 The ith neutrino candidate that is being compared to the GW skymap as
248 well as ``neutrino``.
249 search_params : IceCubeLvcSearchParameters
250 A collection of constant parameters used for this search.
251 Eₙ : float
252 The neutrino energy. This should be capped at the maximum energy of a
253 nearby background neutrino. This way the signal likelihood does not
254 drop (as it does for increasing energy) while the background stays the
255 same (since the background calculation always assumes at least one
256 neutrino).
257 Eₙ_i : float
258 Same as Eₙ, but for the ith neutrino.
260 Returns
261 -------
262 P_H_c2 : float
263 The coincidence hypothesis likelihood for a double-neutrino detection:
265 AllskyIntegral(P(x_nu, x_nu_i | θ, H_c)*P(θ | H_c))
267 """
268 from scipy import integrate
269 from scipy.stats import poisson
270 def innerfive2wogw(Eₙ,dist):
271 return poisson.pmf(2,expnu(Eₙ,dist,search_params)*angle_list()[int((neutrino.zenith+neutrino_i.zenith)/10)])*PEₙ(Eₙ,search_params)*dist**2
272 numerator = integrate.nquad(innerfive2wogw,[[search_params.Eₙmin,search_params.Eₙmax],[0,r0(search_params.Egwmax, search_params)]])[0]
273 denominator = ( (search_params.tnuplus-search_params.tnuminus)**2)
274 integral = sum(
275 (
276 Aeff_skymap(gw_skymap,Eₙ)*Aeff_skymap(gw_skymap,Eₙ_i) * gw_skymap.psf_gaussian(neutrino.ra, neutrino.dec, neutrino.sigma, normalize=False)*gw_skymap.psf_gaussian(neutrino_i.ra, neutrino_i.dec, neutrino_i.sigma, normalize=True)
277 ).pixels)
278 return (integral *t_overlap(neutrino.gps, neutrino_i.gps,search_params) *
279 numerator / denominator)
282def PEgw(Egw, search_params):
283 """Probability of isotropic-equivalent GW emission energy Egw given the
284 signal hypothesis, P(E_{GW}).
286 Parameters
287 ----------
288 Egw : float
289 The assumed isotropic GW emission energy. We will marginalize over this
290 value to get Rdet and the signal likelihoods.
291 search_params : IceCubeLvcSearchParameters
292 A collection of constant parameters used for this search.
294 Returns
295 -------
296 p_e_gw : float
297 P(E_{GW}).
298 """
299 import numpy as np
300 return (Egw*np.log(search_params.Egwmax/search_params.Egwmin))**-1
303def PEₙ(Eₙ, search_params):
304 """Probability of isotropic-equivalent neutrino emission energy Eₙ given
305 the signal hypothesis, P(E_{nu}).
307 Parameters
308 ----------
309 Eₙ : float
310 The assumed isotropic neutrino emission energy. We will marginalize
311 over this value to get Rdet and the signal likelihoods.
312 search_params : IceCubeLvcSearchParameters
313 A collection of constant parameters used for this search.
315 Returns
316 -------
317 p_e_nu : float
318 P(E_{nu}).
319 """
320 import numpy as np
321 return (Eₙ*np.log(search_params.Eₙmax/search_params.Eₙmin))**-1
324def fA(θ, phi, t):
325 """The antenna pattern at direction θ, phi and time t."""
326 # TODO calculate actual antenna factor.
327 return 1
330def iA(search_params):
331 """Integrated antenna factor for detector network given by
332 ``detectors``."""
333 # TODO implement this and fA.
334 return 1
337def expnu(Eₙ, r, search_params):
338 """Expected number of neutrinos for a given emission energy and distance,
339 `<n_ν(E_ν, r)>`
341 Parameters
342 ----------
343 Eₙ : float
344 Neutrino emission energy [kg Mpc**2 / sec**2].
345 r : float
346 Distance to the hypothesized neutrino source. [Mpc]
347 search_params : IceCubeLvcSearchParameters
348 A collection of constant parameters used for this search.
350 Returns
351 -------
352 exp_nu : float
353 `<n_ν(E_ν, r)>`
354 """
355 return search_params.nu_51_100*(Eₙ*10)*(r/100)**-2
358def odds_ratio(tgw, gw_skymap, ρgw, dist, search_params,
359 neutrinolist, p_ter, cwb_Boolean):
360 """Calculate the odds ratio for the signal vs. the null+chance coincidence
361 hypothesis. This is the main measure of significance for the pipeline. Also
362 returns the factors that went into calculating the odds ratio. **Calculates
363 the odds ratio for an entire ensemble of neutrino candidates** in contrast
364 with the OPA version (which calculates p-value on a per-neutrino basis).
366 Parameters
367 ----------
368 tgw : float
369 GPS time of the GW trigger
370 gw_skymap : llama.files.healpix.skymap.HEALPixSkyMap
371 A HEALPix skymap for the GW specifying pixel values and locations.
372 ρgw : float
373 The signal-to-noise ratio (SNR) of the GW
374 dist : float
375 The probability-weighted average reconstructed distance of the GW event
376 search_params : IceCubeLvcSearchParameters
377 A collection of constant parameters used for this search.
378 neutrinolist : list
379 A list of ``llama.files.i3.Neutrino``
380 instances describing the other neutrinos in this search (used for the
381 2-neutrino search, ``five2``).
382 cwb_Boolean : bool
383 Whether this is a CWB event or not.
385 Returns
386 -------
387 odds : float
388 The odds ratio:
389 P(H_s|x_gw, x_nu)
390 -------------------------------------
391 P(H_0|x_gw, x_nu) + P(H_c|x_gw, x_nu)
392 pₕˢ1 : float
393 The signal probability (times the marginal likelihood, cancels
394 everywhere) for the single-neutrino detection case:
395 P(H_s_1|x_gw, x_nu) * P(x_gw, x_nu)
396 pₕˢ2 : float
397 The signal probability (times the marginal likelihood, cancels
398 everywhere) for the double-neutrino detection case:
399 P(H_s_2|x_gw, x_nu1, x_nu2) * P(x_gw, x_nu1, x_nu2)
400 pₕ0 : float
401 The null hypothesis probability (times the marginal likelihood, cancels
402 everywhere):
403 P(H_0|x_gw, x_nu) * P(x_gw, x_nu)
404 pₕᶜ : float
405 The chance coincidence hypothesis probability (times the marginal
406 likelihood, cancels everywhere):
407 P(H_0|x_gw, x_nu) * P(x_gw, x_nu)
408 p_value : float
409 The p-value or false-alarm probability (FAP) of the given
410 ``odds_ratio`` assuming the given source population and detector
411 configuration.
412 dump : dict
413 Local variables with scalar values from the function (for debugging).
414 The contents of this object are not guaranteed to equal anything in
415 particular and should not be used for anything science-related.
416 """
417 import numpy as np
418 from scipy import integrate
419 from scipy.stats import poisson
420 # load search_params into local variables
421 # locals().update(dict(search_params))
422 p = search_params
423 nu_51_100 = p.nu_51_100
424 ndotgwnu = p.ndotgwnu
425 fb = p.fb
426 k0 = p.k0
427 tgwplus = p.tgwplus
428 tgwminus = p.tgwminus
429 tnuplus = p.tnuplus
430 tnuminus = p.tnuminus
431 ratebggw = p.ratebggw
432 ratebgnu = p.ratebgnu
433 Egwmax = p.Egwmax
434 Egwmin = p.Egwmin
435 Eₙmax = p.Eₙmax
436 Eₙmin = p.Eₙmin
437 ϵmin = p.ϵmin
438 ϵmax = p.ϵmax
439 δω = p.δω
440 ΔE = p.ΔE
441 ρmin = p.ρmin
442 #rnumax = p.rnumax
443 #tnu = neutrinolist.gps # GPS time of the neutrino
444 #θₙ = neutrinolist.zenith #zenith deg
445 #phinu = neutrinolist.azimuth #azimuth deg
446 #sigmanu = neutrinolist.sigma #deg2, neutrino angular uncertainty
447 #Eₙ = neutrinolist.energy #GeV
450 θbg, EGev = load_neutrino_background()
451 p_xgw_given_h0 = 1.0*(Pempgw(ρgw, GW_BG_RHOS)*ratebggw+3*ρmin**3/ρgw**4*4/3*pi*r0(Egwmax,search_params)**3*ndotgwnu)*p_ter/(ratebggw+4/3*pi*r0(Egwmax,search_params)**3*ndotgwnu)
452 iₙ=-1
453 # initialize the max and min neutrino energies; we will adjust these to
454 # match the actual max and min values from the neutrino population.
455 maxEₙ=zeros(len(neutrinolist))+0.0
456 minEₙ=zeros(len(neutrinolist))+1000000.0
457 Pₙᵉᵐᵖ=zeros(len(neutrinolist))
458 Eₙ=zeros(len(neutrinolist))
459 θₙ=zeros(len(neutrinolist))
460 count=np.zeros((len(neutrinolist), len(neutrinolist)))
461 for i in range(len(neutrinolist)):
462 count[i][i]=i+1
463 for neutrino in neutrinolist:
464 iₙ +=1
465 Eₙ[iₙ] = neutrino.energy #GeV
466 θₙ[iₙ] = neutrino.zenith
468 # icecube-centric coordinates (zenith, azimuth correspond to θbg, phibg
469 # for the background neutrinos; θ-90 give declination
471 # for i in range (1, len(θbg)):
472 # if(abs(θₙ[iₙ]-θbg[i]) < δω and EGev[i]>maxenu[iₙ]):
473 # maxenu[iₙ] = EGev[i]
474 # if (enu[iₙ] > maxenu[iₙ]):
475 # enu[iₙ] = maxenu[iₙ]
476 # for i in range (1, len(θbg)):
477 # if (abs(θₙ[iₙ] - θbg[i]) < δω and EGev[i] < minenu[iₙ]):
478 # minenu[iₙ] = EGev[i]
480 maxEₙ[iₙ] = max(EGev[(abs(θₙ[iₙ]-θbg)<δω)])
481 minEₙ[iₙ] = min(EGev[(abs(θₙ[iₙ]-θbg)<δω)])
483 if (Eₙ[iₙ] < minEₙ[iₙ]):
484 Eₙ[iₙ] = minEₙ[iₙ]
485 # equations 31 and 33 in the paper at this git revision
486 while (Pₙᵉᵐᵖ[iₙ]<5):
487 # for i in range (1,len(θbg)):
488 # if ( abs(θₙ[iₙ]-θbg[i])<δω and abs(enu[iₙ]-EGev[i])<enu[iₙ]*ΔE):
489 # Pₙᵉᵐᵖ[iₙ]=Pₙᵉᵐᵖ[iₙ]+1
490 Pₙᵉᵐᵖ[iₙ] = dot((abs(θₙ[iₙ]-θbg)<δω)*1,
491 (abs(Eₙ[iₙ]-EGev)<Eₙ[iₙ]*ΔE)*1)
492 # make sure at least five reasonable background neutrinos are predicted
493 # to avoid calculating 0 probability of being in the background.
494 if(Pₙᵉᵐᵖ[iₙ]<5):
495 ΔE *= 1.5
496 δω *= 1.5
497 Pₙᵉᵐᵖ[iₙ] = 0
498 if(Eₙ[iₙ] < np.median(EGev)):
499 Eₙ[iₙ] = Eₙ[iₙ]*10
500 else:
501 Eₙ[iₙ] = Eₙ[iₙ]/10
502 energy_factor = 1.0
503 if(θₙ[iₙ]>80 and neutrino.energy>Eₙ[iₙ]):
504 energy_factor=(neutrino.energy/Eₙ[iₙ])**3.7
505 Eₙ[iₙ]=neutrino.energy
508 Pₙᵉᵐᵖ[iₙ]=Pₙᵉᵐᵖ[iₙ]/(len(θbg)*2*pi*abs(np.cos(pi/180.0*(θₙ[iₙ]-δω))-np.cos(pi/180.0*(θₙ[iₙ]+δω)))*2*ΔE*Eₙ[iₙ])/energy_factor
510 #fiveresult2=0
511 #counthere=0
512 #fiveresult2wogw=0
513 #for neutrino_i in neutrinolist:
514 # counthere=counthere+1
515 # if(iₙ>counthere-1 and not(counthere in count[:][iₙ]) and ((θₙ[iₙ]-neutrino_i.zenith)**2+(phinu[iₙ]-neutrino_i.azimuth)**2)**0.5<5*(sigmanu[iₙ]+neutrino_i.sigma)):
516 # count[iₙ][counthere-1]=counthere
517 # count[counthere-1][iₙ]=iₙ+1
518 # fiveresult2= fiveresult2 + five2(tgw, gw_skymap, ρgw,dist, neutrino, neutrino_i, search_params, enu[iₙ], enu[counthere-1])/Pₙᵉᵐᵖ[counthere-1]/Pₙᵉᵐᵖ[iₙ]
519 # Eq. 36 in this revision's copy of the paper. Time overlap cancels out.
521 pₕˢ1 = (five(tgw, gw_skymap, ρgw,dist, neutrinolist, search_params, Eₙ,Pₙᵉᵐᵖ, cwb_Boolean))*ndotgwnu*poisson.pmf(max(0,len(neutrinolist)-1),ratebgnu*500)
523 p_xnu_given_h0 = poisson.pmf(len(neutrinolist),ratebgnu*500)
525 p_h0_all = float(
526 p_xgw_given_h0 *
527 poisson.pmf(len(neutrinolist), ratebgnu*500)
528 ) / (4/3*pi*r0(Egwmax,search_params)**3)
530 # Eq. 38 in the paper
531 if(not(cwb_Boolean)):
532 p_xgw_xnu_given_hgw_c = (dist)*PEgw((ρgw*dist/k0)**2,search_params)*p_xnu_given_h0*((search_params.Egwmax)-(search_params.Egwmin))*search_params.k0**2/ρgw**2
533 else:
534 def coincidence_integral(d):
535 return (d)**2*PEgw((ρgw*d/k0)**2,search_params)
536 p_xgw_xnu_given_hgw_c = integrate.quad(coincidence_integral,0,r0(search_params.Egwmax, search_params))[0]*p_xnu_given_h0
538 # we will sum over the pₕˢ2 values; initialize to 0
539 #pₕˢ2 = fiveresult2*ndotgwnu*poisson.pmf(max(0,len(neutrinolist)-2),ratebgnu*500)
541 #pₕ0 = 0
542 pₕᶜ = p_xgw_xnu_given_hgw_c *ndotgwnu+ p_xgw_given_h0*(fivewogw(gw_skymap, neutrinolist, search_params, Eₙ,Pₙᵉᵐᵖ))*ndotgwnu*poisson.pmf(max(0,len(neutrinolist)-1),ratebgnu*500)# + fourtyone*fourtyeight # P(H_c|x_gw, x_nu)
543 # make sure we don't count clustered neutrinos in the single neutrino odds
544 single_neutrino_odds = (pₕˢ1 ) / (pₕᶜ + p_h0_all)
545 pval = p_value(
546 single_neutrino_odds,
547 'subthreshold',
548 search_params.population,
549 search_params.sensitivity,
550 search_params.instruments
551 )
552 dump = {} # {k: v for k, v in locals().items() if isinstance(v, Number)}
554 return (single_neutrino_odds, pₕˢ1, 0, p_h0_all, pₕᶜ, pval, dump)
557@SlackReceiptLlama.upload_this()
558@JSONFile.set_class_attributes
559class CoincSignificanceSubthresholdI3Lvc(JSONFile):
560 """Calculates the significance (Bayes Factor) of a joint GW+HEN event."""
562 TIMEOUT = 600
563 FILENAME = "significance_subthresh_lvc-i3.json"
564 DEPENDENCIES = (
565 LVAlertJSON,
566 SkymapInfo,
567 IceCubeNeutrinoList,
568 LvcSkymapHdf5,
569 LVCGraceDbEventData,
570 LvcDistancesJson,
571 PAstro,
572 )
573 DETECTORS = (detectors.IceCube, detectors.LVC,)
575 @property
576 def combined_p_value(self):
577 """Return the combined p-value for the whole event."""
578 return self.read_json()["outputs"]["neutrino_gw_combined"]["p_value"]
580 def _generate(self):
581 import numpy as np
582 # GW event parameters
583 # TODO get ρ from skymap_info
584 # TODO add method for getting ρ from skymap_info
585 tgw = SkymapInfo(self).event_time_gps
586 gracedb = LVCGraceDbEventData(self)
587 ρgw = gracedb.snr
588 instruments = gracedb.instruments
589 population = PAstro(self).most_likely_population
590 p_ter = PAstro(self).p_terr
591 cwb_Boolean = SkymapInfo(self).pipeline == 'CWB'
592 if(cwb_Boolean):
593 dist=1#just a placeholder in order to avoid errors from functions' inputs, will not used in functions if cwb_Boolean is True
594 else:
595 dist = LvcDistancesJson(self).distmean
596 gw_skymap = LvcSkymapHdf5(self).get_healpix()[0]
597 # external parameters mostly related to the run
598 # TODO put real values in for the search parameters here.
599 LOGGER.debug('Initializing search parameters.')
600 search_params = search_parameters(population, instruments)
602 neutrinos = IceCubeNeutrinoList(self).neutrinos
603 LOGGER.debug("BEGIN significance calculation (%s neutrinos)",
604 len(neutrinos))
605 odds, pₕˢ1, pₕˢ2, pₕ0, pₕᶜ, pval, dump = odds_ratio(
606 tgw,
607 gw_skymap,
608 ρgw,
609 dist,
610 #neutrino,
611 search_params,
612 neutrinos,
613 #count,
614 #count2
615 p_ter,
616 cwb_Boolean,
617 )
618 if (not np.isfinite(odds) and not np.isposinf(odds)) or odds < 0:
619 LOGGER.error(("Got a bad odds value {}, setting to 0.0, "
620 "for GW event {} with params:\n"
621 "{}").format(odds, SkymapInfo(self).graceid,
622 dump))
623 odds = 0.
624 if odds == -0.0:
625 LOGGER.debug(("Got -0.0 for odds value, setting to 0.0, for "
626 "GW event {} with params:\n"
627 "{}").format(SkymapInfo(self).graceid, dump))
628 odds = 0.
630 neutrino_gw_combined = {
631 "odds_ratio": odds,
632 "p_h_s1": pₕˢ1,
633 "p_h_s2": pₕˢ2,
634 "p_h_0": pₕ0,
635 "p_h_c": pₕᶜ,
636 "p_value": pval,
637 "dump": dump,
638 }
639 # make a JSON file with the full analysis as output
640 result_dict = {
641 "inputs": {
642 "gw_skymap_info": SkymapInfo(self).read_json(),
643 "neutrino_info": IceCubeNeutrinoList(self).read_json(),
644 "search_parameters": dict(search_params),
645 },
646 "outputs": {
647 "neutrino_gw_combined": neutrino_gw_combined,
648 },
649 }
650 self._write_json(result_dict)