Hide keyboard shortcuts

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). 

4 

5NOTE: equations after 8 increment up by 1. equations after 22 increment up by 

62. 

7 

8@author: dvesk, stefco 

9""" 

10 

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) 

36 

37LOGGER = logging.getLogger(__name__) 

38 

39# original pre-O3 SNR list 

40GW_BG_RHOS = RemoteFileCacher( 

41 'https://nyc3.digitaloceanspaces.com/llama/llama/objects/' 

42 '88722142e6310591f918f524a901c4376f41d22f9c8cd22c9b598d5c66c0010c' 

43) 

44 

45def five(tgw, gw_skymap, ρgw, dist,neutrino, search_params, enu, cwb_Boolean): 

46 """Calculate the signal likelihood (Eq. 5). 

47 

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. 

70 

71 Returns 

72 ------- 

73 P_H_S1 : float 

74 The signal hypothesis likelihood for a single neutrino: 

75 

76 AllskyIntegral(P(x_gw, x_nu | theta, H_s)*P(theta | H_s)) 

77 

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) * 

94 

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) 

101 

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) 

106  

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. 

135 

136 Returns 

137 ------- 

138 P_H_S2 : float 

139 The signal hypothesis likelihood for a double-neutrino detection: 

140 

141 AllskyIntegral(P(x_gw, x_nu, x_nu_i | theta, H_s)*P(theta | H_s)) 

142 

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) 

162 

163 

164def PEgw(Egw, search_params): 

165 """Probability of isotropic-equivalent GW emission energy Egw given the 

166 signal hypothesis, P(E_{GW}). 

167 

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. 

175 

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 

183 

184 

185def PEnu(Enu, search_params): 

186 """Probability of isotropic-equivalent neutrino emission energy Enu given 

187 the signal hypothesis, P(E_{nu}). 

188 

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. 

196 

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 

204 

205 

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 

210 

211 

212def iA(search_params): 

213 """Integrated antenna factor for detector network given by 

214 ``detectors``.""" 

215 # TODO implement this and fA. 

216 return 1 

217 

218 

219def expnu(Enu, r, search_params): 

220 """Expected number of neutrinos for a given emission energy and distance, 

221 < n_{\nu}(E_{\nu}, r) > 

222 

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. 

231 

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 

238 

239 

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. 

245  

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. 

272 

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 

331 

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 

337 

338 

339 

340 

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. 

343 

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] 

355 

356 thetabg, EGev = load_neutrino_background() 

357 

358 #replace with these two lines 

359 maxenu = max(EGev[(abs(thetanu-thetabg)<δω)]) 

360 minenu = min(EGev[(abs(thetanu-thetabg)<δω)]) 

361 

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 

366 

367 # equations 31 and 33 in the paper at this git revision 

368 p_xgw_given_h0 = 1.0*Pempgw(ρgw, GW_BG_RHOS) 

369 

370 Pempnu=0 

371 

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) 

379 

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 

394 

395 

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 

407 

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* 

424 

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 

441 

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 

476 

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) 

478 

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)} 

493 

494 return (single_neutrino_odds, p_h_s1, p_h_s2, p_h_0, p_h_c, pval, dump) 

495 

496 

497@SlackReceiptLlama.upload_this() 

498@JSONFile.set_class_attributes 

499class CoincSignificanceI3Lvc(JSONFile): 

500 """Calculates the significance (Bayes Factor) of a joint GW+HEN event.""" 

501 

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,) 

513 

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)) 

530 

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']] 

538 

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"] 

543 

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) 

564 

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 ) 

618 

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)