Coverage for /opt/anaconda/lib/python3.7/site-packages/llama/files/coinc_significance/utils.py : 29%
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"""
2Utility functions for calculating significance of joint events.
4NOTE: equations after 8 increment up by 1. equations after 22 increment up by
52.
7@author: dvesk, stefco
8"""
10import logging
11from llama.files.healpix import (
12 HEALPixSkyMap,
13)
14from llama.utils import (
15 parameter_factory,
16 memoize,
17 RemoteFileCacher,
18)
19from llama.com.s3 import PrivateFileCacher
21BACKGROUND_NEUTRINOS = RemoteFileCacher(
22 'https://nyc3.digitaloceanspaces.com/llama/llama/objects/'
23 '1e43fa4add9b0e2bff5c856898db5b78c9aacf53bea6e950f1fbc637d373f634'
24)
25NEUTRINO_POPULATIONS = PrivateFileCacher(
26 'llama/objects/2ae02042fd06fb977765cdd69440b3dc6cda0a13a27b3ef5cf2a3dd894e3e387',
27 'llama',
28)
29NEUTRINO_POPULATIONS_TARBALL = PrivateFileCacher(
30 'llama/objects/080b0a6c4ea33b878921d9bbf64ed95a589eab262cff7248edc1506118d8f284'
31 'llama',
32)
33LOGGER = logging.getLogger(__name__)
34ODDS_RATIO_KEYS = [
35 'bayes_factor', # due to a typo in early version of the pipeline
36 'odds_ratio',
37]
39def r0(Egw, search_params):
40 """Maximum GW detection distance for isotropic-equivalent emission energy
41 Egw."""
42 return search_params.k0 * Egw**0.5 / search_params.ρmin
45@memoize()
46def load_neutrino_background():
47 """Load neutrino background data.
49 Returns
50 -------
51 thetabg : array-like
52 Array of background neutrino zenith angles [deg] as reconstructed by
53 IceCube.
54 EGev : array-like
55 Array of background neutrino energies [GeV] as reconstructed by
56 IceCube.
57 """
58 import pandas
59 neutrinos = pandas.read_hdf(BACKGROUND_NEUTRINOS.get())
60 return neutrinos['thetabg'], neutrinos['EGev']
63def p_value(odds_ratio, search_type, population, sensitivity, instruments):
64 """Calculate the p-value for a given odds ratio.
66 Parameters
67 ----------
68 odds_ratio : float
69 The odds ratio for signal vs. background/null; is here being
70 interpreted as a test statistic to make the analysis robust against
71 constant factors in the priors.
72 search_type : str
73 The type of search performed. This decides which background population
74 will be used. Distinguishes between subthreshold, open public alert,
75 and in the future, potentially other simulated backgrounds.
76 population : str
77 The name of the source population against which this ``odds_ratio``
78 will be compared, e.g. 'bns' for binary neutron stars. *NOTE*: This
79 just has to be set so that it references an existing group in the
80 ``NEUTRINO_POPULATIONS`` HDF5 file; there are currently no explict
81 standards for naming the ``population``. See
82 ``populations-hdf5-collater`` and the contents of the directory
83 containing ``NEUTRINO_POPULATIONS`` for available options.
84 sensitivity : str
85 The detector configuration and sensitivity used, e.g. 'design' for HLV
86 at design sensitivity. See the *NOTE* on ``population`` for naming
87 information.
88 instruments : list
89 The instruments participating in this search as a list of detector
90 names, e.g. ``[H1, L1, V1]`` for all three O3 GW detectors. If a
91 background population is available for this set of instruments, it will
92 be used; otherwise, at attempt will be made to use a population
93 averaging over all sets of participating detectors (for populations
94 where full separate cases could not be feasibly generated).
96 Returns
97 -------
98 p_value : float or NoneType
99 The p-value or false-alarm probability (FAP) of the given
100 ``odds_ratio`` assuming the given source population and detector
101 configuration. Returns ``None`` if a matching background does not yet
102 exist for this search configuration.
103 """
104 import h5py
105 det = ''.join(sorted(instruments))
106 with h5py.File(NEUTRINO_POPULATIONS.get(), 'r') as bg_odds_file:
107 try:
108 bg_odds = bg_odds_file[search_type][population][sensitivity][det][()]
109 except (ValueError, KeyError):
110 LOGGER.error("No background population found for search_type: %s, "
111 "population: %s, sensitivity: %s, detector: %s ",
112 search_type, population, sensitivity, det)
113 return None
114 return 1 - bg_odds.searchsorted(odds_ratio) / len(bg_odds)
117def Pempgw(ρgw, gw_bg_ρs):
118 """The empirical neutrino background for a given GW SNR ``ρgw``.
120 Parameters
121 ----------
122 ρgw : float
123 The SNR of the current GW trigger.
124 gw_bg_ρs : RemoteFileCacher
125 A file with a list of historical background SNR values.
126 """
127 import numpy as np
128 binnumber=50
129 with gw_bg_ρs.get().open() as nubgfile:
130 cont = nubgfile.readlines()
131 cont2 = np.zeros((len(cont), 1))
132 for i in range(0, len(cont)):
133 cont2[i] = float(cont[i])
134 maxcont2=max(cont2)[0]
135 mincont2=min(cont2)[0]
136 hist, binedges = np.histogram(cont2, bins=binnumber, range=(mincont2, maxcont2))
137 i = 1
138 while ρgw > binedges[i]:
139 i = i+1
140 if (i == len(binedges)):
141 i = i-1
142 break
143 # conservatively find at least one matching background event even if the
144 # bin is empty (since the tails are sparsely populated).
145 return max(hist[i-1],1)*1.0/hist.sum()/(maxcont2-mincont2)*binnumber
148# pylint: disable=invalid-name
149IceCubeLvcSearchParameters = parameter_factory(
150 "IceCubeLvcSearchParameters",
151 """A collection of search parameters for IceCube/LVC joint searches for the
152 LLAMA 2.0 pipeline. Energies should be in ergs, distances in Mpc.
154 https://arxiv.org/abs/1810.11467
155 """,
156 # TODO add units
157 nu_51_100="Detector-specific parameter for neutrino observatories (Eq 10)",
158 fb="Neutrino beaming factor",
159 expected_multimessenger_event_detection_rate_density="""Expected
160 rate-density of GW+HEN events detected [1/(m**3 t)]""",
161 k0="""A constant that relates ρ, event distance, and emission energy for
162 typical GWs via the equation r = k0 * Egw**0.5 / ρgw.""",
163 k1="""A constant that relates the lowest admissible ρ, the GW emission
164 energy, and the maximum possible source distance for those values via
165 the equation r_max = k1 * Egw**0.5 / ρmin.""",
166 tgwminus="""Max number of seconds BEFORE event time that GW emission can
167 occur in the model.""",
168 tgwplus="""Max number of seconds AFTER event time that GW emission can
169 occur in the model.""",
170 tnuminus="""Max number of seconds BEFORE event time that neutrino emission
171 can occur in the model.""",
172 tnuplus="""Max number of seconds AFTER event time that neutrino emission
173 can occur in the model.""",
174 ratebggw="Background rates for GW channel (Eq. 30) in Hz",
175 ndotgwnu="""a constant such that ndotgwnu*P(Egw,Enu|Hs) corresponds to the
176 differential rate density of GW+neutrino events (Eq. 11)""",
177 ratebgnu="Background rates for neutrino channel (Eq. 30) in Hz",
178 Egwmax="""Maximum isotropic-equivalent GW energy emission accounted for by
179 the model [GeV].""",
180 Egwmin="""Minimum isotropic-equivalent GW energy emission accounted for by
181 the model [GeV].""",
182 # TODO add units.
183 Eₙmax="""Maximum isotropic-equivalent neutrino energy emission accounted
184 for by the model.""",
185 # TODO add units.
186 Eₙmin="""Minimum isotropic-equivalent neutrino energy emission accounted
187 for by the model.""",
188 ϵmin="Minimum detectable neutrino energy (Eq. 22)",
189 ϵmax="""Maximum neutrino energy (above which the emission rate is
190 negligible); serves as an upper bound for the integral in Eq. 22.""",
191 ρmin="""Threshold GW SNR (ρgw) under which GW events are not analyzed
192 by the pipeline. Affects event candidate selection and constrains the
193 maximum detected GW source distance. (Eq 18, 47)""",
194 δω="""Max angular distance between neutrino reconstructed direction
195 and background neutrino reconstructed directions to use when
196 calculating chance coincidence hypothesis.""",
197 ΔE="""Max energy difference between neutrino reconstructed energy and
198 background neutrino reconstructed energy to use when calculating chance
199 coincidence hypothesis.""",
200 population="""The name of the source population to use when calculating the
201 p-value (see note in ``p_value`` on which names are valid), e.g. 'bns'
202 for binary neutron starts.""",
203 sensitivity="""The detector configuration and sensitivity to
204 assume when calculating the p-value (see note in ``p_value`` on which
205 names to use and which configurations are available), e.g. 'design' for
206 HLV at design sensitivity.""",
207 instruments="""The instruments that participated in the discovery of a
208 particular trigger in list form, e.g. ``["H1", "L1", "V1"]`` for all
209 LIGO/Virgo GW detectors.""",
210 #rnumax=rnumax
211)
214def t_overlap(tgw, tnu, search_params):
215 """Integral of the source-time-dependent part of the likelihood integral
216 (Eq. 5) finding the probability of getting ``tgw`` and ``tnu`` (the
217 detection times of the gravitational wave trigger and neutrino trigger,
218 respectively) given the specified ``search_params``. Factors out from the
219 rest of Eq. 5."""
220 return max(
221 0,
222 (min(tgw+search_params.tgwplus, tnu+search_params.tnuplus) -
223 max(tgw+search_params.tgwminus, tnu+search_params.tnuminus))
224 )
227# TODO define which exact angle we are using in our lookup table; dec? zenith?
228def aeff_lookup_table():
229 """
230 Get a dictionary mapping from strings describing neutrino
231 energy bounds to tuples of the form (factor, exponent) where ``factor`` is a
232 lookup table of the constant factor multiplying the effective area (with
233 indices representing 10-degree increments in the angle theta) and
234 ``exponent`` is the exponent to which the neutrino energy ``e`` is raised
235 (since we assume a rough power law for Aeff, again with indices representing
236 10-degree increments in the angle theta). The final value will be
237 ``factor*e**exponent``.
239 NOTE that we need to include an item for [180, 190) degrees because
240 int(180./10) will give us 18, indexing into the 19th entry in the lookup
241 table, and we want to include this edge case.
242 """
243 import numpy as np
244 return {
245 4: ( # e < 10**2 [GeV]
246 np.array([
247 -15. , -15. , -15. , -15. ,
248 -15. , -15. , -15. , -15. ,
249 -4.90235479, -4.38396659, -4.40832267, -4.7313624 ,
250 -4.58396307, -4.73739424, -5.00846592, -4.78305924,
251 -4.47899964, -4.18126993, -4.18
252 ]),
253 np.array([
254 0. , 0. , 0. , 0. , 0. ,
255 0. , 0. , 0. , 0.34436205, 0.12582712,
256 0.17628614, 0.4954696 , 0.37370245, 0.57628784, 0.90438348,
257 0.94177117, 0.90763264, 0.74932548, 0.75
258 ])
259 ),
260 5: ( # e < 10**2.5 [GeV]
261 np.array([
262 -15. , -15. , -15. , -15. ,
263 -15. , -15. , -15. , -15. ,
264 -4.55799274, -4.25813947, -4.23203653, -4.2358928 ,
265 -4.21026062, -4.1611064 , -4.10408244, -3.84128807,
266 -3.57136699, -3.43194445, -3.43
267 ]),
268 np.array([
269 0. , 0. , 0. , 0. , 0. ,
270 0. , 0. , 0. , 1.20473396, 0.98717625,
271 0.97966312, 1.00738257, 0.99134297, 0.94732103, 0.88415571,
272 0.62595829, 0.31468422, 0.11159207, 0.11
273 ])
274 ),
275 6: ( # e < 10**3 [GeV]
276 np.array([
277 -15. , -15. , -15. , -15. ,
278 -15. , -15. , -15. , -15. ,
279 -3.35325878, -3.27096322, -3.25237341, -3.22851022,
280 -3.21891764, -3.21378537, -3.21992673, -3.21532979,
281 -3.25668277, -3.32035238, -3.32
282 ]),
283 np.array([
284 0. , 0. , 0. , 0. , 0. ,
285 0. , 9.48859873, 10.16760375, -0.6849539 , -0.80286694,
286 -0.85374556, -0.88919918, -0.89323734, -0.9175427 , -0.95471555,
287 -1.02075346, -1.01410299, -0.93018747, -0.93
288 ])
289 ),
290 7: ( # e < 10**3.5 [GeV]
291 np.array([
292 -15. , -15. , -15. , -15. ,
293 -15. , -15. , -5.51140127, -4.83239625,
294 -4.03821268, -4.07383017, -4.10611897, -4.1177094 ,
295 -4.11215499, -4.13132807, -4.17464228, -4.23608325,
296 -4.27078576, -4.25053985, -4.25
297 ]),
298 np.array([
299 9.89429222, 9.7875918 , 9.31902768, 9.09135011, 9.53211507,
300 9.89219328, 0.16757057, 0.41436273, -0.84342435, -0.96672298,
301 -1.05037202, -0.99455757, -1.06779287, -1.09967731, -1.08900755,
302 -1.0719317 , -1.1455688 , -1.2305211 , -1.23
303 ])
304 ),
305 8: ( # e < 10**4 [GeV]
306 np.array([
307 -5.10570778, -5.2124082 , -5.68097232, -5.90864989, -5.46788493,
308 -5.10780672, -5.3438307 , -4.41803352, -4.88163703, -5.04055314,
309 -5.15649099, -5.11226698, -5.17994785, -5.23100538, -5.26364983,
310 -5.30801495, -5.41635455, -5.48106095, -5.48
311 ]),
312 np.array([
313 0.00979211, 0.11743399, 0.69796885, 0.8364719 , 0.36332801,
314 0.13797984, 0.42222146, -0.54389681, -1.05252801, -1.19674539,
315 -1.08660649, -1.18304152, -1.19964391, -1.22403293, -1.30984807,
316 -1.76417344, -1.22666052, -1.31624722, -1.32
317 ])
318 ),
319 9: ( # e < 10**4.5 [GeV]
320 np.array([
321 -5.09591566, -5.09497421, -4.98300347, -5.07217799, -5.10455692,
322 -4.96982688, -4.92160924, -4.96193033, -5.93416504, -6.23729853,
323 -6.24309749, -6.2953085 , -6.37959176, -6.4550383 , -6.5734979 ,
324 -7.07218839, -6.64301508, -6.79730817, -6.80
325 ]),
326 np.array([
327 -0.25713648, -0.40441414, -0.38103129, -0.27439524, -0.28844569,
328 -0.46017066, -0.62468374, -0.92680854, -1.22303767, -1.15094019,
329 -1.31501771, -1.5002799 , -1.22590789, -1.85832137, -1.76168438,
330 -0.72358103, -1.36116491, -5.68343383, -5.68
331 ])
332 ),
333 10: ( # e < 10**5 [GeV]
334 np.array([
335 -5.35305214, -5.49938835, -5.36403475, -5.34657322,
336 -5.39300262, -5.42999754, -5.54629298, -5.88873887,
337 -7.15720271, -7.38823872, -7.5581152 , -7.7955884 ,
338 -7.60549965, -8.31335968, -8.33518228, -7.79576942,
339 -8.00417999, -12.480742, -12.48
340 ]),
341 np.array([
342 -1.17141307, -1.21879427, -0.97977989, -0.85173484, -1.09922416,
343 -1.23275651, -1.17944999, -1.19841619, -1.16663105, -1.18244783,
344 -1.66127277, -1.47143556, -1.55130616, -6.68664032, -1.75273715,
345 -3.44202581, -6.99582001, -2.519258 , -2.52
346 ])
347 ),
348 11: ( # e < 10**5.5 [GeV]
349 np.array([
350 -6.52446521, -6.71818261, -6.34381465, -6.19830806,
351 -6.49222677, -6.66275405, -6.72574296, -7.08715507,
352 -8.32383376, -8.57068655, -9.21938797, -9.26702396,
353 -9.15680581, -15. , -10.08791943, -11.23779523,
354 -15. , -15. , -15.
355 ]),
356 np.array([
357 -8.47553479, -1.46568724, -1.21868672, -1.4595885 , -1.18975634,
358 -0.99450981, -1.22154436, -1.01422582, -1.24586562, -1.01925748,
359 -1.28972563, -1.17564258, -2.35048237, 4.1899023 , -4.91208057,
360 -3.76220477, 0. , 0. ,0.
361 ])
362 ),
363 12: ( # e < 10**6 [GeV]
364 np.array([
365 -15. , -8.18386986, -7.56250137, -7.65789656,
366 -7.68198311, -7.65726385, -7.94728732, -8.10138089,
367 -9.56969938, -9.58994404, -10.50911359, -10.44266654,
368 -11.50728818, -10.8100977 , -15. , -15. ,
369 -15. , -15. , -15.
370 ]),
371 np.array([
372 0. , 0.00655032, -0.60124358, -0.58206732, -0.60303634,
373 -1.01475895, -0.99433117, -0.44069475, -5.43030062, -0.53937906,
374 -4.49088641, -4.55733346, -3.49271182, -4.1899023 , 0. ,
375 0. , 0. , 0. , 0.
376 ])
377 ),
378 13: ( # e < 10**6.5 [GeV]
379 np.array([
380 -15. , -8.17731954, -8.16374495, -8.23996388,
381 -8.28501945, -8.6720228 , -8.94161849, -8.54207563,
382 -15. , -10.12932309, -15. , -15. ,
383 -15. , -15. , -15. , -15. ,
384 -15. , -15. , -15.
385 ]),
386 np.array([
387 6.4773239 , -6.82268046, -6.83625505, -0.81428677, -1.05563951,
388 -0.77352849, -0.4472135 , -6.45792437, 4.27147988, -4.87067691,
389 0. , 0. , 0. , 0. , 0. ,
390 0. , 0. , 0. , 0.
391 ])
392 ),
393 14: ( # e < 10**7 [GeV]
394 np.array([
395 -8.5226761 , -15. , -15. , -9.05425065,
396 -9.34065896, -9.44555129, -9.388832 , -15. ,
397 -10.72852012, -15. , -15. , -15. ,
398 -15. , -15. , -15. , -15. ,
399 -15. , -15. , -15.]),
400 np.array([
401 -6.4773239 , 0. , 0. , -5.94574935, -5.65934104,
402 -0.47052734, -5.611168 , 4.23464539, -4.27147988, 0. ,
403 0. , 0. , 0. , 0. , 0. ,
404 0. , 0. , 0. , 0.
405 ])
406 ),
407 15: ( # e < 10**7.5 [GeV]
408 np.array([
409 -15. , -15. , -15. , -15. ,
410 -15. , -9.91607863, -15. , -10.76535461,
411 -15. , -15. , -15. , -15. ,
412 -15. , -15. , -15. , -15. ,
413 -15. , -15. ,-15
414 ]),
415 np.array([
416 0. , 0. , 0. , 0. , 0. ,
417 -5.08392137, 0. , -4.23464539, 0. , 0. ,
418 0. , 0. , 0. , 0. , 0. ,
419 0. , 0. , 0. ,0.
420 ])
421 ),
422 'high_e': ( # e >= 10**8 [GeV]
423 np.array( # Constant Factor
424 1 * [-15.0] + # [0, 10) [deg]
425 18 * [-15.0] # [10, 190) [deg]
426 ),
427 np.array( # Exponent applied to the neutrino energy
428 1 * [0.] + # [0, 10) [deg]
429 18 * [0.] # [10, 190) [deg]
430 )
431 ),
432 }
435# relative occurence probability of a signal neutrino in 5deg zenith bands
436# starting from 0
437def angle_list():
438 import numpy as np
439 return np.array([0.064912 , 0.06930208, 0.07447563, 0.0793188 ,
440 0.08591584, 0.09197046, 0.10022791, 0.10723858,
441 0.11913539, 0.12817956, 0.14293096, 0.17679627,
442 0.21033154, 0.25710645, 0.35197546, 0.54138516,
443 1.33562294, 2.59808902, 2.66636028, 2.44928372,
444 2.25773049, 2.11718214, 2.01375073, 1.91195977,
445 1.82192818, 1.72508975, 1.63648361, 1.57418033,
446 1.5086993, 1.42200218, 1.28744103, 1.21574383,
447 1.15708913, 1.10138159, 1.04943515, 0.92742034])
450def Aeff(e, theta):
451 import numpy as np
452 log_e = np.log10(e)
453 # angle-dependent from here on out; use the angles to index into
454 # aeff_lookup_table.
455 angle_index = np.array((theta)/10, dtype=int)
456 # default factor, exponent for Aeff at highest energies
457 factor, exponent = aeff_lookup_table()['high_e']
458 # see if we need factor, exponent for Aeff at lower energies
459 for log_e_upper in range (5,17):
460 if log_e < log_e_upper/2.0:
461 factor, exponent = aeff_lookup_table()[log_e_upper-1]
462 break
463 a_eff = 10**(factor[angle_index]+(log_e-log_e_upper/2.0+0.5)*2.0*exponent[angle_index])
464 # if input was a scalar, still return a scalar
465 if a_eff.shape == ():
466 return float(a_eff)
467 return a_eff
470def Aeff_skymap(skymap, enu):
471 """Get the energy- and direction-dependent effective areas of the neutrino
472 detector for every direction provided in ``skymap`` (Eq. 21, 22)
474 Parameters
475 ----------
476 skymap : llama.files.healpix.HEALPixSkyMap
477 A ``HEALPixSkyMap`` instance with defined right-ascensions and
478 declinations.
479 enu : float
480 The reconstructed neutrino energy.
481 """
482 return HEALPixSkyMap(Aeff(enu, 90+skymap.dec), nest=skymap.nest)
485# TODO get IceCubeLvcSearchParameters into auto documentation at some point,
486# since it has descriptions of all of our parameters
487def search_parameters(population, instruments):
488 """Return an ``IceCubeLvcSearchParameters`` instance with all default
489 values and the given ``population`` and ``instruments``.
491 population : str
492 The name of the source population against which this ``odds_ratio``
493 will be compared, e.g. 'bns' for binary neutron stars. *NOTE*: This
494 just has to be set so that it references an existing group in the
495 ``NEUTRINO_POPULATIONS`` HDF5 file; there are currently no explict
496 standards for naming the ``population``. See
497 ``populations-hdf5-collater`` and the contents of the directory
498 containing ``NEUTRINO_POPULATIONS`` for available options.
499 instruments : list
500 The GW detectors participating in this search as a list of detector
501 names, e.g. ``[H1, L1, V1]`` for all three O3 GW detectors. If a
502 background population is available for this set of instruments, it will
503 be used; otherwise, at attempt will be made to use a population
504 averaging over all sets of participating instruments (for populations
505 where full separate cases could not be feasibly generated).
507 Returns
508 -------
509 search_params : IceCubeLvcSearchParameters
510 Default search parameter values along with the provided ``population``
511 and ``instruments``.
512 """
513 return IceCubeLvcSearchParameters(
514 nu_51_100=1.1, #*1.05*10**-49
515 ndotgwnu=1500./(365*24*3600)*10**-9, #(1/(MpC^3s))
516 fb=10.0,
517 k0=400.8, #(s^1.5/kg^0.5)
518 k1=400.8,
519 # TODO unused it seems; remove
520 expected_multimessenger_event_detection_rate_density=11,
521 tgwplus=250.0, #s
522 tgwminus=-250.0,
523 tnuplus=250.0,
524 tnuminus=-250.0,
525 ratebggw=1./(3600*24), #1/s
526 ratebgnu=6.4/1000, # 1/s
527 Egwmax=18.7, #kg MpC^2/s^2
528 Egwmin=18.7*10**-4,
529 Eₙmax=0.1, #kg MpC^2/s^2
530 Eₙmin=10**-6,
531 ϵmin=10.0, # GeV
532 ϵmax=10**9, # GeV
533 ρmin=5.59,
534 δω=2.5,
535 ΔE=0.3,
536 population=population,
537 instruments=instruments,
538 sensitivity='o3',
539 )