Coverage for /opt/anaconda/lib/python3.7/site-packages/llama/files/coinc_significance/opa.py : 16%
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
3high-confidence gravitational wave triggers from Open Public Alerts (OPAs).
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 llama import detectors
14from llama.filehandler import JSONFile
15from llama.files.i3 import IceCubeNeutrinoList
16from llama.files.lvc_skymap import (
17 LvcSkymapHdf5, LvcDistancesJson
18)
19from llama.files.skymap_info import SkymapInfo
20from llama.files.gracedb import LVCGraceDbEventData, PAstro
21from llama.utils import (
22 RemoteFileCacher,
23)
24from llama.files.slack import SlackReceiptLlama
25from llama.files.coinc_significance.utils import (
26 ODDS_RATIO_KEYS,
27 # r0,
28 load_neutrino_background,
29 p_value,
30 Pempgw,
31 t_overlap,
32 angle_list,
33 Aeff_skymap,
34 search_parameters,
35)
37LOGGER = logging.getLogger(__name__)
39# original pre-O3 SNR list
40GW_BG_RHOS = RemoteFileCacher(
41 'https://nyc3.digitaloceanspaces.com/llama/llama/objects/'
42 '88722142e6310591f918f524a901c4376f41d22f9c8cd22c9b598d5c66c0010c'
43)
45def five(tgw, gw_skymap, ρgw, dist,neutrino, search_params, enu, cwb_Boolean):
46 """Calculate the signal likelihood (Eq. 5).
48 Parameters
49 ----------
50 tgw : float
51 GPS time of the GW trigger
52 gw_skymap : llama.files.healpix.skymap.HEALPixSkyMap
53 A HEALPix skymap for the GW specifying pixel values and locations.
54 ρgw : float
55 The signal-to-noise ratio (SNR) of the GW
56 dist : float
57 The probability-weighted average reconstructed distance of the GW event
58 neutrino : llama.files.i3.Neutrino
59 The neutrino candidate that is being compared to the GW skymap
60 search_params : IceCubeLvcSearchParameters
61 A collection of constant parameters used for this search.
62 enu : float
63 The neutrino energy. This should be capped at the maximum energy of a
64 nearby background neutrino. This way the signal likelihood does not
65 drop (as it does for increasing energy) while the background stays the
66 same (since the background calculation always assumes at least one
67 neutrino).
68 cwb_Boolean : bool
69 Whether this is a CWB event or not.
71 Returns
72 -------
73 P_H_S1 : float
74 The signal hypothesis likelihood for a single neutrino:
76 AllskyIntegral(P(x_gw, x_nu | theta, H_s)*P(theta | H_s))
78 """
79 from scipy import integrate
80 from scipy.stats import poisson
81 if(not(cwb_Boolean)):
82 def innerfive(Enu):
83 return poisson.pmf(1,expnu(Enu,dist,search_params)*angle_list()[int(neutrino.zenith/5)])*PEnu(Enu,search_params)
84 numerator = PEgw((dist*ρgw/search_params.k0)**2,search_params)*dist**2*integrate.quad(innerfive,search_params.Eₙmin,search_params.Eₙmax)[0]
85 else:
86 def innerfive(Enu,d):
87 return PEgw((d*ρgw/search_params.k0)**2,search_params)*d**2*poisson.pmf(1,expnu(Enu,d,search_params)*angle_list()[int(neutrino.zenith/5)])*PEnu(Enu,search_params)
88 numerator = integrate.nquad(innerfive,[[search_params.Eₙmin,search_params.Eₙmax],[0,r0(search_params.Egwmax, search_params)]])[0]
89 denominator = ((search_params.tnuplus-search_params.tnuminus)*(search_params.tgwplus-search_params.tgwminus))
90 integral = sum(
91 fA(90-gw_skymap.dec, gw_skymap.ra, tgw) *
92 (gw_skymap *
93 Aeff_skymap(gw_skymap,enu) *
95 gw_skymap.psf_gaussian(neutrino.ra, neutrino.dec, neutrino.sigma,
96 normalize=False)
97 ).pixels
98 )
99 return (integral * t_overlap(tgw, neutrino.gps, search_params) *
100 numerator / denominator)
102def five2(tgw, gw_skymap, ρgw,dist, neutrino, neutrino_i, search_params, enu,
103 enu_i, cwb_Boolean):
104 """Calculate the signal likelihood for 2 coincident neutrinos (modified Eq.
105 5)
107 Parameters
108 ----------
109 tgw : float
110 GPS time of the GW trigger
111 gw_skymap : llama.files.healpix.skymap.HEALPixSkyMap
112 A HEALPix skymap for the GW specifying pixel values and locations.
113 ρgw : float
114 The signal-to-noise ratio (SNR) of the GW
115 dist : float
116 The probability-weighted average reconstructed distance of the GW event
117 neutrino : llama.files.i3.Neutrino
118 The neutrino candidate that is being compared to the GW skymap as well
119 as ``neutrino_i``.
120 neutrino_i : llama.files.i3.Neutrino
121 The ith neutrino candidate that is being compared to the GW skymap as
122 well as ``neutrino``.
123 search_params : IceCubeLvcSearchParameters
124 A collection of constant parameters used for this search.
125 enu : float
126 The neutrino energy. This should be capped at the maximum energy of a
127 nearby background neutrino. This way the signal likelihood does not
128 drop (as it does for increasing energy) while the background stays the
129 same (since the background calculation always assumes at least one
130 neutrino).
131 enu_i : float
132 Same as enu, but for the ith neutrino.
133 cwb_Boolean : bool
134 Whether this is a CWB event or not.
136 Returns
137 -------
138 P_H_S2 : float
139 The signal hypothesis likelihood for a double-neutrino detection:
141 AllskyIntegral(P(x_gw, x_nu, x_nu_i | theta, H_s)*P(theta | H_s))
143 """
144 from scipy import integrate
145 from scipy.stats import poisson
146 if(not(cwb_Boolean)):
147 def innerfive2(Enu):
148 return poisson.pmf(2,expnu(Enu,dist,search_params)*angle_list()[int((neutrino.zenith+neutrino_i.zenith)/10)])*PEnu(Enu,search_params)
149 numerator = PEgw((dist*ρgw/search_params.k0)**2,search_params)*dist**2*integrate.quad(innerfive2,search_params.Eₙmin,search_params.Eₙmax)[0]
150 else:
151 def innerfive2(Enu,d):
152 return PEgw((d*ρgw/search_params.k0)**2,search_params)*poisson.pmf(2,expnu(Enu,d,search_params)*angle_list()[int((neutrino.zenith+neutrino_i.zenith)/10)])*PEnu(Enu,search_params)*d**2
153 numerator = integrate.nquad(innerfive2,[[search_params.Eₙmin,search_params.Eₙmax],[0,r0(search_params.Egwmax, search_params)]])[0]
154 denominator = ( (search_params.tnuplus-search_params.tnuminus)*(search_params.tgwplus-search_params.tgwminus))
155 integral = sum(
156 fA(90-gw_skymap.dec, gw_skymap.ra, tgw) *
157 (
158 gw_skymap* Aeff_skymap(gw_skymap,enu)*Aeff_skymap(gw_skymap,enu_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)
159 ).pixels)
160 return (integral *t_overlap(max(tgw, neutrino.gps, neutrino_i.gps),min(tgw, neutrino.gps, neutrino_i.gps),search_params) *
161 numerator / denominator)
164def PEgw(Egw, search_params):
165 """Probability of isotropic-equivalent GW emission energy Egw given the
166 signal hypothesis, P(E_{GW}).
168 Parameters
169 ----------
170 Egw : float
171 The assumed isotropic GW emission energy. We will marginalize over this
172 value to get Rdet and the signal likelihoods.
173 search_params : IceCubeLvcSearchParameters
174 A collection of constant parameters used for this search.
176 Returns
177 -------
178 p_e_gw : float
179 P(E_{GW}).
180 """
181 import numpy as np
182 return (Egw*np.log(search_params.Egwmax/search_params.Egwmin))**-1
185def PEnu(Enu, search_params):
186 """Probability of isotropic-equivalent neutrino emission energy Enu given
187 the signal hypothesis, P(E_{nu}).
189 Parameters
190 ----------
191 Enu : float
192 The assumed isotropic neutrino emission energy. We will marginalize
193 over this value to get Rdet and the signal likelihoods.
194 search_params : IceCubeLvcSearchParameters
195 A collection of constant parameters used for this search.
197 Returns
198 -------
199 p_e_nu : float
200 P(E_{nu}).
201 """
202 import numpy as np
203 return (Enu*np.log(search_params.Eₙmax/search_params.Eₙmin))**-1
206def fA(theta, phi, t):
207 """The antenna pattern at direction theta, phi and time t."""
208 # TODO calculate actual antenna factor.
209 return 1
212def iA(search_params):
213 """Integrated antenna factor for detector network given by
214 ``detectors``."""
215 # TODO implement this and fA.
216 return 1
219def expnu(Enu, r, search_params):
220 """Expected number of neutrinos for a given emission energy and distance,
221 < n_{\nu}(E_{\nu}, r) >
223 Parameters
224 ----------
225 Enu : float
226 Neutrino emission energy [kg Mpc**2 / sec**2].
227 r : float
228 Distance to the hypothesized neutrino source. [Mpc]
229 search_params : IceCubeLvcSearchParameters
230 A collection of constant parameters used for this search.
232 Returns
233 -------
234 exp_nu : float
235 < n_{\nu}(E_{\nu}, r) >
236 """
237 return search_params.nu_51_100*(Enu*10)*(r/100)**-2
240def odds_ratio(tgw, gw_skymap, ρgw, dist, neutrino, search_params,
241 neutrinolist, count, count2, cwb_Boolean):
242 """Calculate the odds ratio for the signal vs. the null+chance coincidence
243 hypothesis. This is the main measure of significance for the pipeline. Also
244 returns the factors that went into calculating the odds ratio.
246 Parameters
247 ----------
248 tgw : float
249 GPS time of the GW trigger
250 gw_skymap : llama.files.healpix.skymap.HEALPixSkyMap
251 A HEALPix skymap for the GW specifying pixel values and locations.
252 ρgw : float
253 The signal-to-noise ratio (SNR) of the GW
254 dist : float
255 The probability-weighted average reconstructed distance of the GW event
256 neutrino : llama.files.i3.Neutrino
257 The neutrino candidate that is being compared to the GW skymap
258 search_params : IceCubeLvcSearchParameters
259 A collection of constant parameters used for this search.
260 neutrinolist : list
261 A list of ``llama.files.i3.Neutrino``
262 instances describing the other neutrinos in this search (used for the
263 2-neutrino search, ``five2``).
264 count : np.array
265 An array tracking whether two neutrinos have been compared. If an entry
266 is 1, then they are considered to have been compared (the diagonal
267 entries should be set to 1 at the start).
268 count2 : int
269 The index of the current primary neutrino in ``neutrinolist``.
270 cwb_Boolean : bool
271 Whether this is a CWB event or not.
273 Returns
274 -------
275 odds : float
276 The odds ratio:
277 P(H_s|x_gw, x_nu)
278 -------------------------------------
279 P(H_0|x_gw, x_nu) + P(H_c|x_gw, x_nu)
280 p_h_s1 : float
281 The signal probability (times the marginal likelihood, cancels
282 everywhere) for the single-neutrino detection case:
283 P(H_s_1|x_gw, x_nu) * P(x_gw, x_nu)
284 p_h_s2 : float
285 The signal probability (times the marginal likelihood, cancels
286 everywhere) for the double-neutrino detection case:
287 P(H_s_2|x_gw, x_nu1, x_nu2) * P(x_gw, x_nu1, x_nu2)
288 p_h_0 : float
289 The null hypothesis probability (times the marginal likelihood, cancels
290 everywhere):
291 P(H_0|x_gw, x_nu) * P(x_gw, x_nu)
292 p_h_c : float
293 The chance coincidence hypothesis probability (times the marginal
294 likelihood, cancels everywhere):
295 P(H_0|x_gw, x_nu) * P(x_gw, x_nu)
296 p_value : float
297 The p-value or false-alarm probability (FAP) of the given
298 ``odds_ratio`` assuming the given source population and detector
299 configuration.
300 dump : dict
301 Local variables with scalar values from the function (for debugging).
302 The contents of this object are not guaranteed to equal anything in
303 particular and should not be used for anything science-related.
304 """
305 import numpy as np
306 from scipy import integrate
307 from scipy.stats import poisson
308 # load search_params into local variables
309 # locals().update(dict(search_params))
310 p = search_params
311 nu_51_100 = p.nu_51_100
312 ndotgwnu = p.ndotgwnu
313 fb = p.fb
314 k0 = p.k0
315 tgwplus = p.tgwplus
316 tgwminus = p.tgwminus
317 tnuplus = p.tnuplus
318 tnuminus = p.tnuminus
319 ratebggw = p.ratebggw
320 ratebgnu = p.ratebgnu
321 Egwmax = p.Egwmax
322 Egwmin = p.Egwmin
323 Eₙmax = p.Eₙmax
324 Eₙmin = p.Eₙmin
325 ϵmin = p.ϵmin
326 ϵmax = p.ϵmax
327 δω = p.δω
328 ΔE = p.ΔE
329 ρmin = p.ρmin
330 #rnumax = p.rnumax
332 tnu = neutrino.gps # GPS time of the neutrino
333 thetanu = neutrino.zenith #zenith deg
334 phinu = neutrino.azimuth #azimuth deg
335 sigmanu = neutrino.sigma #deg2, neutrino angular uncertainty
336 enu = neutrino.energy #GeV
341 # initialize the max and min neutrino energies; we will adjust these to
342 # match the actual max and min values from the neutrino population.
344 # REMOVE THIS BLOCK for speed
345 # maxenu=0.0
346 # minenu=1000000.0
347 # for i in range (1, len(thetabg)):
348 # if(abs(thetanu-thetabg[i]) < δω and EGev[i]>maxenu):
349 # maxenu = EGev[i]
350 # if (enu > maxenu):
351 # enu = maxenu
352 # for i in range (1, len(thetabg)):
353 # if (abs(thetanu - thetabg[i]) < δω and EGev[i] < minenu):
354 # minenu = EGev[i]
356 thetabg, EGev = load_neutrino_background()
358 #replace with these two lines
359 maxenu = max(EGev[(abs(thetanu-thetabg)<δω)])
360 minenu = min(EGev[(abs(thetanu-thetabg)<δω)])
362 # icecube-centric coordinates (zenith, azimuth correspond to thetabg, phibg
363 # for the background neutrinos; theta-90 give declination
364 if (enu < minenu):
365 enu = minenu
367 # equations 31 and 33 in the paper at this git revision
368 p_xgw_given_h0 = 1.0*Pempgw(ρgw, GW_BG_RHOS)
370 Pempnu=0
372 while (Pempnu<5):
373 #remove these three lines:
374 # for i in range (1,len(thetabg)):
375 # if ( abs(thetanu-thetabg[i])<δω and abs(enu-EGev[i])<enu*ΔE):
376 # Pempnu=Pempnu+1
377 #replace with this line:
378 Pempnu=np.dot((abs(thetanu-thetabg)<δω)*1,(abs(enu-EGev)<enu*ΔE)*1)
380 # make sure at least five reasonable background neutrinos are predicted
381 # to avoid calculating 0 probability of being in the background.
382 if(Pempnu<5):
383 ΔE = ΔE*1.5
384 δω = δω*1.5
385 Pempnu = 0
386 if(enu < np.median(EGev)):
387 enu = enu*10
388 else:
389 enu = enu/10
390 energy_factor = 1.0
391 if(thetanu>80 and neutrino.energy>enu):
392 energy_factor=(neutrino.energy/enu)**3.7
393 enu=neutrino.energy
396 Pempnu=Pempnu/(len(thetabg)*2*pi*abs(np.cos(pi/180.0*(thetanu-δω))-np.cos(pi/180.0*(thetanu+δω)))*2*ΔE*enu)/energy_factor
397 # Eq. 36 in this revision's copy of the paper. Time overlap cancels out.
398 p_xnu_given_h0 = float(Pempnu)
399 # background coincidence rate
400 Rbg = ratebggw*ratebgnu
401 # Eq. 26 in the paper
402 p_h0_given_xgw_xnu = float(p_xgw_given_h0 * p_xnu_given_h0)
403 #def fourtyfour(r,Egw):
404 # return r**2*PEgw(Egw, search_params)*(Tobs(tgw)*N44)**-1
405 #def fourtyseven(r,Enu):
406 # return r**2*poisson.pmf(1,expnu(Enu, r, search_params))*(Tobs(tgw)*N47)**-1
408 #def innerfourtyfive(Egw):
409 # return PEgw(Egw, search_params)*r0(Egw, search_params)**3/3
410 # Eq. 45 in the paper, but with the integral cancelled elsewhere:
411 p_h_gw_c = ratebgnu*ndotgwnu#*4*pi*R41*integrate.quad(innerfourtyfive,Egwmin,Egwmax)[0]
412 # Eq. 38 in the paper
413 if(not(cwb_Boolean)):
414 p_xgw_xnu_given_hgw_c = (dist)**2*PEgw((ρgw*dist/k0)**2,search_params)*p_xnu_given_h0
415 else:
416 def coincidence_integral(d):
417 return (d)**2*PEgw((ρgw*d/k0)**2,search_params)
418 p_xgw_xnu_given_hgw_c = integrate.quad(coincidence_integral,0,r0(search_params.Egwmax, search_params))[0]*p_xnu_given_h0
419 #def innerfourtyeight(r,Enu):
420 # return r**2*(poisson.pmf(1, expnu(Enu, r, search_params)))*PEnu(Enu, search_params)
421 #TODO
422 #fourtyeight=2.0*ratebggw*Tobs(tgw)*ndotgwnu*4*pi*#integrate.nquad(innerfourtyeight,[[0,rnumax],[Eₙmin,Eₙmax]])[0]
423 #fourtyeight_2=2.0*ratebggw*Tobs(tgw)*ndotgwnu*4*pi*
425 #def innerthirtyeight(Egw):
426 # return p_xnu_given_h0*fourtyfour(k0/ρgw*Egw**0.5,Egw)
427 #p_xgw_xnu_given_hgw_c=fourtyfour(dist,(dist*ρgw/k0)**2)*p_xnu_given_h0#integrate.nquad(innerthirtyeight,[[Egwmin,Egwmax]])[0]
428 #def innerfourtyone(r,Enu):
429 # return fourtyseven()
430 #integral = sum(
431 # fA(90-gw_skymap.dec, gw_skymap.ra, tgw) *
432 # (
433 # Aeff_skymap(gw_skymap,enu) *
434 # gw_skymap.psf_gaussian(neutrino.ra, neutrino.dec, neutrino.sigma, normalize=True)
435 # ).pixels
436 #)
437 #fourtyone=p_xgw_given_h0*integral/(4*pi*Ne*enu**2)
438 # we will sum over the p_h_s2 values; initialize to 0
439 fiveresult2=0
440 counthere=0
442 for neutrino_i in neutrinolist:
443 δω=p.δω
444 ΔE=p.ΔE
445 counthere=counthere+1
446 if(not(counthere in count[:][count2-1]) and ((thetanu-neutrino_i.zenith)**2+(phinu-neutrino_i.azimuth)**2)**0.5<5*(sigmanu+neutrino_i.sigma)):
447 count[count2-1][counthere-1]=counthere
448 count[counthere-1][count2-1]=count2
449 enu_i=neutrino_i.energy
450 if(enu_i>maxenu):
451 enu_i=maxenu
452 if(enu_i<minenu):
453 enu_i=minenu
454 Pempnu2=0
455 # make sure at least five background neutrinos are predicted to
456 # avoid calculating 0 probability of being in the background.
457 while (Pempnu2<5):
458 Pempnu2 = np.dot(
459 (abs(neutrino_i.zenith-thetabg)<δω)*1,
460 (abs(enu_i-EGev)<enu_i*ΔE)*1
461 )
462 if(Pempnu<5):
463 ΔE=ΔE*1.5
464 δω=δω*1.5
465 Pempnu=0
466 if(enu_i<np.median(EGev)):
467 enu_i=enu_i*10
468 else:
469 enu_i=enu_i/10
470 energy_factor2=1.0
471 if(neutrino_i.zenith>80 and neutrino_i.energy>enu_i):
472 energy_factor2=(neutrino_i.energy/enu_i)**3.7
473 enu_i=neutrino_i.energy
474 Pempnu2=Pempnu2/(len(thetabg)*2*pi*abs(np.cos(pi/180.0*(neutrino_i.zenith-δω))-np.cos(pi/180.0*(neutrino_i.zenith+δω)))*2*ΔE*enu_i)/energy_factor2
475 fiveresult2= fiveresult2 + five2(tgw, gw_skymap, ρgw,dist, neutrino, neutrino_i, search_params, enu, enu_i, cwb_Boolean)/Pempnu2
477 p_h_s1 = (five(tgw, gw_skymap, ρgw,dist, neutrino, search_params, enu, cwb_Boolean))*ndotgwnu/fb # P(H_s|x_gw, x_nu)
479 p_h_s2 = fiveresult2*ndotgwnu/fb
480 #p_h_0 = Rbg * p_h0_given_xgw_xnu # P(H_0|x_gw, x_nu)
481 p_h_0 = 0
482 p_h_c = p_xgw_xnu_given_hgw_c * p_h_gw_c# + fourtyone*fourtyeight # P(H_c|x_gw, x_nu)
483 # make sure we don't count clustered neutrinos in the single neutrino odds
484 single_neutrino_odds = (p_h_s1 + p_h_s2) / (p_h_c + p_h_0)
485 pval = p_value(
486 single_neutrino_odds,
487 'opa',
488 search_params.population,
489 search_params.sensitivity,
490 search_params.instruments
491 )
492 dump = {} # {k: v for k, v in locals().items() if isinstance(v, Number)}
494 return (single_neutrino_odds, p_h_s1, p_h_s2, p_h_0, p_h_c, pval, dump)
497@SlackReceiptLlama.upload_this()
498@JSONFile.set_class_attributes
499class CoincSignificanceI3Lvc(JSONFile):
500 """Calculates the significance (Bayes Factor) of a joint GW+HEN event."""
502 TIMEOUT = 600
503 FILENAME = "significance_lvc-i3.json"
504 DEPENDENCIES = (
505 SkymapInfo,
506 IceCubeNeutrinoList,
507 LvcSkymapHdf5,
508 LVCGraceDbEventData,
509 LvcDistancesJson,
510 PAstro,
511 )
512 DETECTORS = (detectors.IceCube, detectors.LVC,)
514 def odds(self):
515 """Return a list of the odds ratios calculated. Ordering corresponds to
516 the ordering of neutrinos in the ``IceCubeNeutrinoList`` file used to
517 calculate the odds ratios."""
518 triggers = self.read_json()['outputs']
519 # an early version of the analysis code included a typo that called the
520 # odds ratio the bayes factor, so try both possibilities.
521 for odds_ratio_key in ODDS_RATIO_KEYS:
522 try:
523 return [n[odds_ratio_key]
524 for n in triggers['neutrino_gw_triggers']]
525 except KeyError:
526 pass
527 raise ValueError(("Could not parse the odds ratio values in "
528 "the provided list of trigger significance "
529 "calculations:\n{}").format(triggers))
531 @property
532 def p_values(self):
533 """Return a list of the p-values calculated. Ordering corresponds to
534 the ordering of neutrinos in the ``IceCubeNeutrinoList`` file used to
535 calculate the odds ratios."""
536 triggers = self.read_json()['outputs']
537 return [n['p_value'] for n in triggers['neutrino_gw_triggers']]
539 @property
540 def combined_p_value(self):
541 """Return the combined p-value for the whole event."""
542 return self.read_json()["outputs"]["neutrino_gw_combined"]["p_value"]
544 def _generate(self):
545 import numpy as np
546 # GW event parameters
547 # TODO get ρ from skymap_info
548 # TODO add method for getting ρ from skymap_info
549 tgw = SkymapInfo(self).event_time_gps
550 gracedb = LVCGraceDbEventData(self)
551 ρgw = gracedb.snr
552 instruments = gracedb.instruments
553 population = PAstro(self).most_likely_population
554 cwb_Boolean = SkymapInfo(self).pipeline == 'CWB'
555 if(cwb_Boolean):
556 dist=1#just a placeholder in order to avoid errors from functions' inputs, will not used in functions if cwb_Boolean is True
557 else:
558 dist = LvcDistancesJson(self).distmean
559 gw_skymap = LvcSkymapHdf5(self).get_healpix()[0]
560 # external parameters mostly related to the run
561 # TODO put real values in for the search parameters here.
562 LOGGER.debug('Initializing search parameters.')
563 search_params = search_parameters(population, instruments)
565 # calculate odds ratio for each neutrino/gw combination and format
566 # output joint trigger dictionaries with RA/Dec units, energy, sigma,
567 # and odds ratio for improved readability
568 neutrino_gw_triggers = list()
569 neutrinos = IceCubeNeutrinoList(self).neutrinos
570 count=np.zeros((len(neutrinos), len(neutrinos)))
571 for i in range(len(neutrinos)):
572 count[i][i]=i+1
573 count2=0
574 LOGGER.debug("BEGIN significance calculation (%s neutrinos)",
575 len(neutrinos))
576 for neutrino in neutrinos:
577 count2=count2+1
578 LOGGER.debug("Starting neutrino %s/%s", count2, len(neutrinos))
579 odds, p_h_s1, p_h_s2, p_h_0, p_h_c, pval, dump = odds_ratio(
580 tgw,
581 gw_skymap,
582 ρgw,
583 dist,
584 neutrino,
585 search_params,
586 neutrinos,
587 count,
588 count2,
589 cwb_Boolean,
590 )
591 if (not np.isfinite(odds) and not np.isposinf(odds)) or odds < 0:
592 LOGGER.error(("Got a bad odds value {}, setting to 0.0, "
593 "for GW event {} with params:\n"
594 "{}").format(odds, SkymapInfo(self).graceid,
595 dump))
596 odds = 0.
597 if odds == -0.0:
598 LOGGER.debug(("Got -0.0 for odds value, setting to 0.0, for "
599 "GW event {} with params:\n"
600 "{}").format(SkymapInfo(self).graceid, dump))
601 odds = 0.
602 neutrino_gw_triggers.append(
603 {
604 'ra': neutrino.ra,
605 'dec': neutrino.dec,
606 'sigma': neutrino.sigma,
607 'energy': neutrino.energy,
608 'dt': neutrino.gps - tgw,
609 'odds_ratio': odds,
610 'p_h_s1': p_h_s1,
611 'p_h_s2': p_h_s2,
612 'p_h_0': p_h_0,
613 'p_h_c': p_h_c,
614 'p_value': pval,
615 'dump': dump,
616 }
617 )
619 # get the combined p-value
620 odds_sum = sum(n['odds_ratio'] for n in neutrino_gw_triggers)
621 neutrino_gw_combined = {
622 "odds_ratio": odds_sum,
623 "p_value": p_value(
624 odds_sum,
625 'opa',
626 search_params.population,
627 search_params.sensitivity,
628 search_params.instruments,
629 ),
630 }
631 # make a JSON file with the full analysis as output
632 result_dict = {
633 "inputs": {
634 "gw_skymap_info": SkymapInfo(self).read_json(),
635 "neutrino_info": IceCubeNeutrinoList(self).read_json(),
636 "search_parameters": dict(search_params),
637 },
638 "outputs": {
639 "neutrino_gw_triggers": neutrino_gw_triggers,
640 "neutrino_gw_combined": neutrino_gw_combined,
641 },
642 }
643 self._write_json(result_dict)