Coverage for ligo/followup_advocate/__init__.py: 91%

1import math as mth 

2import os 

3import shutil 

4import tempfile 

5import urllib 

6import webbrowser 

7 

8import astropy.coordinates as coord 

9import astropy.time 

10import astropy.units as u 

11import healpy as hp 

12import lxml.etree 

13import numpy as np 

14from astropy.io.fits import getheader 

15from ligo.gracedb import rest 

16from ligo.skymap.io.fits import read_sky_map 

17from ligo.skymap.postprocess.ellipse import find_ellipse 

18from ligo.skymap.postprocess.crossmatch import crossmatch 

19 

20from .jinja import env 

21import importlib.metadata 

22__version__ = importlib.metadata.version(__name__) 

23 

24 

25def authors(authors, service=rest.DEFAULT_SERVICE_URL): 

26 """Write GCN Circular author list""" 

27 return env.get_template('authors.jinja2').render(authors=authors).strip() 

28 

29 

30def guess_skyloc_pipeline(filename): 

31 keys = ['cWB', 'BAYESTAR', 'Bilby', 'LIB', 'LALInference', 

32 'oLIB', 'MLy', 'UNKNOWN'] 

33 skyloc_pipelines_dict = dict(zip([x.lower() for x in keys], keys)) 

34 skyloc_pipelines_dict['rapidpe_rift'] = 'RapidPE-RIFT' 

35 try: 

36 return skyloc_pipelines_dict[filename.split('.')[0].lower()] 

37 except KeyError: 

38 return filename.split('.')[0] 

39 

40 

41def text_width(remove_text_wrap): 

42 """Return width of text wrap based on whether we wish to wrap the lines or 

43 not.""" 

44 return 9999 if remove_text_wrap else 79 

45 

46 

47def main_dict(gracedb_id, client, raven_coinc=False): 

48 """Create general dictionary to pass to compose circular""" 

49 

50 superevent = client.superevent(gracedb_id).json() 

51 preferred_event = superevent['preferred_event_data'] 

52 preferred_pipeline = preferred_event['pipeline'] 

53 if preferred_pipeline.lower() == 'cwb' and \ 

54 preferred_event['search'].lower() == 'bbh': 

55 preferred_pipeline = 'cwb_bbh' 

56 early_warning_pipelines = [] 

57 pipelines = [] 

58 gw_events = superevent['gw_events'] 

59 early_warning_alert = False 

60 

61 for gw_event in gw_events: 

62 gw_event_dict = client.event(gw_event).json() 

63 pipeline = gw_event_dict['pipeline'] 

64 search = gw_event_dict['search'] 

65 if pipeline.lower() == 'cwb' and search.lower() == 'bbh': 

66 pipeline = 'cwb_bbh' 

67 if pipeline not in pipelines: 

68 pipelines.append(pipeline) 

69 if pipeline not in early_warning_pipelines and \ 

70 search == 'EarlyWarning': 

71 early_warning_pipelines.append(pipeline) 

72 # Sort to get alphabetical order 

73 pipelines.sort(key=str.lower) 

74 early_warning_pipelines.sort(key=str.lower) 

75 

76 voevents = client.voevents(gracedb_id).json()['voevents'] 

77 if not voevents: 

78 raise ValueError( 

79 "{} has no VOEvent to generate circulars from.".format( 

80 gracedb_id)) 

81 

82 citation_index = {pipeline.lower(): pipelines.index(pipeline) + 1 for 

83 pipeline in pipelines} 

84 for pipeline in early_warning_pipelines: 

85 if pipeline.lower() != 'mbta': 

86 citation_index[pipeline.lower() + '_earlywarning'] = \ 

87 max(citation_index.values()) + 1 

88 

89 gpstime = float(preferred_event['gpstime']) 

90 event_time = astropy.time.Time(gpstime, format='gps').utc 

91 

92 # Grab latest p_astro and em_bright values from lastest VOEvent 

93 voevent_text_latest = \ 

94 client.files(gracedb_id, voevents[-1]['filename']).read() 

95 root = lxml.etree.fromstring(voevent_text_latest) 

96 p_astros = root.find('./What/Group[@name="Classification"]') 

97 em_brights = root.find('./What/Group[@name="Properties"]') 

98 classifications = {} 

99 source_classification = {} 

100 # Only try to load if present to prevent errors with .getchildren() 

101 p_astro_pipeline = None 

102 

103 # Grab logs for either RapidPE checks or chirp mass bins 

104 if p_astros or em_brights: 

105 logs = client.logs(gracedb_id).json()['log'] 

106 

107 if p_astros: 

108 for p_astro in p_astros.getchildren(): 

109 if p_astro.attrib.get('value'): 

110 classifications[p_astro.attrib['name']] = \ 

111 float(p_astro.attrib['value']) * 100 

112 # Figure out which pipeline uploaded p_astro, usually the first one 

113 # FIXME: Replace with more efficient method in the future 

114 for message in reversed(logs): 

115 filename = message['filename'] 

116 if filename and '.p_astro.json' in filename and \ 

117 filename != 'p_astro.json': 

118 p_astro = client.files(gracedb_id, filename).json() 

119 if all(mth.isclose(p_astro[key] * 100, classifications[key]) 

120 for key in classifications.keys()): 

121 p_astro_pipeline = filename.split('.')[0].lower() 

122 break 

123 if p_astro_pipeline == 'rapidpe_rift': 

124 citation_index['rapidpe_rift'] = max(citation_index.values()) + 1 

125 if em_brights: 

126 for em_bright in em_brights.getchildren(): 

127 if em_bright.attrib.get('value'): 

128 source_classification[em_bright.attrib['name']] = \ 

129 float(em_bright.attrib['value']) * 100 

130 if len(source_classification) > 0: 

131 citation_index['em_bright'] = max(citation_index.values()) + 1 

132 

133 skymaps = {} 

134 voevents_text = [] 

135 for voevent in voevents: 

136 # Don't load latest voevent since already loaded from before 

137 if voevent == voevents[-1]: 

138 voevent_text = voevent_text_latest 

139 # If earlier voevent, load 

140 else: 

141 voevent_text = client.files(gracedb_id, voevent['filename']).read() 

142 root = lxml.etree.fromstring(voevent_text) 

143 alert_type = root.find( 

144 './What/Param[@name="AlertType"]').attrib['value'].lower() 

145 if alert_type == 'earlywarning': 

146 # Add text for early warning detection if one early warning alert 

147 early_warning_alert = True 

148 url = root.find('./What/Group/Param[@name="skymap_fits"]') 

149 if url is None: 

150 continue 

151 url = url.attrib['value'] 

152 _, filename = os.path.split(url) 

153 skyloc_pipeline = guess_skyloc_pipeline(filename) 

154 issued_time = astropy.time.Time(root.find('./Who/Date').text).gps 

155 if filename not in skymaps: 

156 skymaps[filename] = dict( 

157 alert_type=alert_type, 

158 pipeline=skyloc_pipeline, 

159 filename=filename, 

160 latency=issued_time-event_time.gps) 

161 if skyloc_pipeline.lower() not in citation_index: 

162 citation_index[skyloc_pipeline.lower()] = \ 

163 max(citation_index.values()) + 1 

164 voevents_text.append(voevent_text) 

165 skymaps = list(skymaps.values()) 

166 

167 o = urllib.parse.urlparse(client.service_url) 

168 

169 kwargs = dict( 

170 subject='Identification', 

171 gracedb_id=gracedb_id, 

172 gracedb_service_url=urllib.parse.urlunsplit( 

173 (o.scheme, o.netloc, '/superevents/', '', '')), 

174 group=preferred_event['group'], 

175 pipeline=preferred_pipeline, 

176 all_pipelines=pipelines, 

177 early_warning_alert=early_warning_alert, 

178 early_warning_pipelines=early_warning_pipelines, 

179 gpstime='{0:.03f}'.format(round(float(preferred_event['gpstime']), 3)), 

180 search=preferred_event.get('search', ''), 

181 far=preferred_event['far'], 

182 utctime=event_time.iso, 

183 instruments=preferred_event['instruments'].split(','), 

184 skymaps=skymaps, 

185 prob_has_ns=source_classification.get('HasNS'), 

186 prob_has_remnant=source_classification.get('HasRemnant'), 

187 prob_has_massgap=source_classification.get('HasMassGap'), 

188 prob_has_ssm=source_classification.get('HasSSM'), 

189 source_classification=source_classification, 

190 include_ellipse=None, 

191 classifications=classifications, 

192 p_astro_pipeline=p_astro_pipeline, 

193 citation_index=citation_index) 

194 

195 if skymaps: 

196 preferred_skymap = skymaps[-1]['filename'] 

197 cls = [50, 90] 

198 include_ellipse, ra, dec, a, b, pa, area, greedy_area = \ 

199 uncertainty_ellipse(gracedb_id, preferred_skymap, client, cls=cls) 

200 kwargs.update( 

201 preferred_skymap=preferred_skymap, 

202 cl=cls[-1], 

203 include_ellipse=include_ellipse, 

204 ra=coord.Longitude(ra*u.deg), 

205 dec=coord.Latitude(dec*u.deg), 

206 a=coord.Angle(a*u.deg), 

207 b=coord.Angle(b*u.deg), 

208 pa=coord.Angle(pa*u.deg), 

209 ellipse_area=area, 

210 greedy_area=greedy_area) 

211 try: 

212 distmu, distsig = get_distances_skymap_gracedb(gracedb_id, 

213 preferred_skymap, 

214 client) 

215 kwargs.update( 

216 distmu=distmu, 

217 distsig=distsig) 

218 except TypeError: 

219 pass 

220 

221 if raven_coinc: 

222 kwargs = _update_raven_parameters(superevent, kwargs, client, 

223 voevents_text) 

224 return kwargs 

225 

226 

227def compose(gracedb_id, authors=(), mailto=False, remove_text_wrap=False, 

228 service=rest.DEFAULT_SERVICE_URL, client=None): 

229 """Compose GCN Circular draft""" 

230 

231 if client is None: 

232 client = rest.GraceDb(service) 

233 

234 kwargs = main_dict(gracedb_id, client=client) 

235 kwargs.update(authors=authors) 

236 kwargs.update(change_significance_statement=False) 

237 kwargs.update(text_width=text_width(remove_text_wrap)) 

238 

239 subject = env.get_template('subject.jinja2').render(**kwargs).strip() 

240 body = env.get_template('initial_circular.jinja2').render(**kwargs).strip() 

241 

242 if mailto: 

243 pattern = 'mailto:emfollow@ligo.org,dac@ligo.org?subject={0}&body={1}' 

244 url = pattern.format( 

245 urllib.parse.quote(subject), 

246 urllib.parse.quote(body)) 

247 webbrowser.open(url) 

248 else: 

249 return '{0}\n\n{1}'.format(subject, body) 

250 

251 

252def compose_raven(gracedb_id, authors=(), remove_text_wrap=False, 

253 service=rest.DEFAULT_SERVICE_URL, client=None): 

254 """Compose EM_COINC RAVEN GCN Circular draft""" 

255 

256 if client is None: 

257 client = rest.GraceDb(service) 

258 

259 kwargs = dict() 

260 kwargs.update(main_dict(gracedb_id, client=client, raven_coinc=True)) 

261 kwargs.update(update_alert=False) 

262 kwargs.update(text_width=text_width(remove_text_wrap)) 

263 # Add RAVEN citation 

264 citation_index = kwargs['citation_index'] 

265 citation_index['raven'] = max(citation_index.values()) + 1 

266 kwargs['citation_index'] = citation_index 

267 

268 subject = ( 

269 env.get_template('RAVEN_subject.jinja2').render(**kwargs) 

270 .strip()) 

271 body = ( 

272 env.get_template('RAVEN_circular.jinja2').render(**kwargs) 

273 .strip()) 

274 return '{0}\n\n{1}'.format(subject, body) 

275 

276 

277def compose_llama( 

278 gracedb_id, authors=(), service=rest.DEFAULT_SERVICE_URL, 

279 icecube_alert=None, remove_text_wrap=False, 

280 client=None): 

281 """Compose GRB LLAMA GCN Circular draft. 

282 Here, gracedb_id will be a GRB superevent ID in GraceDb.""" 

283 

284 if client is None: 

285 client = rest.GraceDb(service) 

286 

287 superevent = client.superevent(gracedb_id).json() 

288 

289 gpstime = float(superevent['t_0']) 

290 tl, th = gpstime - 500, gpstime + 500 

291 event_time = astropy.time.Time(gpstime, format='gps').utc 

292 tl_datetime = str(astropy.time.Time( 

293 tl, format='gps').isot).replace('T', ' ') 

294 th_datetime = str(astropy.time.Time( 

295 th, format='gps').isot).replace('T', ' ') 

296 

297 o = urllib.parse.urlparse(client.service_url) 

298 kwargs = dict( 

299 gracedb_service_url=urllib.parse.urlunsplit( 

300 (o.scheme, o.netloc, '/superevents/', '', '')), 

301 gracedb_id=gracedb_id, 

302 llama=True, 

303 icecube_alert=icecube_alert, 

304 event_time=event_time, 

305 tl_datetime=tl_datetime, 

306 th_datetime=th_datetime, 

307 authors=authors) 

308 kwargs.update(text_width=text_width(remove_text_wrap)) 

309 

310 citation_index = {'llama': 1, 'llama_stat': 2} 

311 kwargs.update(citation_index=citation_index) 

312 

313 files = client.files(gracedb_id).json() 

314 

315 llama_stat_filename = 'significance_subthreshold_lvc-i3.json' 

316 if llama_stat_filename in files: 

317 llama_stat_file = client.files(gracedb_id, llama_stat_filename).json() 

318 llama_fap = llama_stat_file["p_value"] 

319 neutrinos = llama_stat_file["inputs"]["neutrino_info"] 

320 lines = [] 

321 for neutrino in neutrinos: 

322 # Build aligned string 

323 vals = [] 

324 dt = (event_time - 

325 astropy.time.Time(neutrino['mjd'], 

326 format='mjd')).to(u.s).value 

327 vals.append('{:.2f}'.format(dt)) 

328 vals.append('{:.2f}'.format(neutrino['ra'])) 

329 vals.append('{:.2f}'.format(neutrino['dec'])) 

330 vals.append('{:.2f}'.format(neutrino['sigma'])) 

331 vals.append('{:.4f}'.format(llama_fap)) 

332 lines.append(align_number_string(vals, [0, 11, 21, 40, 59])) 

333 

334 kwargs.update(llama_fap=llama_fap, 

335 neutrinos=lines) 

336 

337 subject = ( 

338 env.get_template('llama_subject.jinja2').render(**kwargs) 

339 .strip()) 

340 if icecube_alert: 

341 body = ( 

342 env.get_template('llama_alert_circular.jinja2').render(**kwargs) 

343 .strip()) 

344 else: 

345 body = ( 

346 env.get_template('llama_track_circular.jinja2').render(**kwargs) 

347 .strip()) 

348 return '{0}\n\n{1}'.format(subject, body) 

349 

350 

351def compose_grb_medium_latency( 

352 gracedb_id, authors=(), service=rest.DEFAULT_SERVICE_URL, 

353 use_detection_template=None, remove_text_wrap=False, client=None): 

354 """Compose GRB Medium Latency GCN Circular draft. 

355 Here, gracedb_id will be a GRB external trigger ID in GraceDb.""" 

356 

357 if client is None: 

358 client = rest.GraceDb(service) 

359 

360 event = client.event(gracedb_id).json() 

361 search = event['search'] 

362 

363 if search != 'GRB': 

364 return 

365 

366 group = event['group'] 

367 pipeline = event['pipeline'] 

368 external_trigger = event['extra_attributes']['GRB']['trigger_id'] 

369 

370 o = urllib.parse.urlparse(client.service_url) 

371 kwargs = dict( 

372 gracedb_service_url=urllib.parse.urlunsplit( 

373 (o.scheme, o.netloc, '/events/', '', '')), 

374 gracedb_id=gracedb_id, 

375 group=group, 

376 grb=True, 

377 pipeline=pipeline, 

378 external_trigger=external_trigger, 

379 exclusions=[], 

380 detections=[]) 

381 kwargs.update(text_width=text_width(remove_text_wrap)) 

382 

383 files = client.files(gracedb_id).json() 

384 

385 citation_index = {} 

386 index = 1 

387 xpipeline_fap_file = 'false_alarm_probability_x.json' 

388 if xpipeline_fap_file in files: 

389 xpipeline_fap = client.files(gracedb_id, xpipeline_fap_file).json() 

390 online_xpipeline_fap = xpipeline_fap.get('Online Xpipeline') 

391 # Create detection/exclusion circular based on given argument 

392 # Use default cutoff if not provided 

393 xpipeline_detection = (use_detection_template if use_detection_template 

394 is not None else online_xpipeline_fap < 0.001) 

395 if xpipeline_detection: 

396 kwargs['detections'].append('xpipeline') 

397 kwargs.update(online_xpipeline_fap=online_xpipeline_fap) 

398 else: 

399 kwargs['exclusions'].append('xpipeline') 

400 xpipeline_distances_file = 'distances_x.json' 

401 xpipeline_distances = client.files(gracedb_id, 

402 xpipeline_distances_file).json() 

403 burst_exclusion = xpipeline_distances.get('Burst Exclusion') 

404 kwargs.update(burst_exclusion=burst_exclusion) 

405 citation_index['xpipeline'] = index 

406 index += 1 

407 

408 pygrb_fap_file = 'false_alarm_probability_pygrb.json' 

409 if pygrb_fap_file in files: 

410 pygrb_fap = client.files(gracedb_id, pygrb_fap_file).json() 

411 online_pygrb_fap = pygrb_fap.get('Online PyGRB') 

412 # Create detection/exclusion circular based on given argument 

413 # Use default cutoff if not provided 

414 pygrb_detection = (use_detection_template if use_detection_template 

415 is not None else online_pygrb_fap < 0.01) 

416 if pygrb_detection: 

417 kwargs['detections'].append('pygrb') 

418 kwargs.update(online_pygrb_fap=online_pygrb_fap) 

419 else: 

420 kwargs['exclusions'].append('pygrb') 

421 pygrb_distances_file = 'distances_pygrb.json' 

422 pygrb_distances = client.files(gracedb_id, 

423 pygrb_distances_file).json() 

424 nsns_exclusion = pygrb_distances.get('NSNS Exclusion') 

425 nsbh_exclusion = pygrb_distances.get('NSBH Exclusion') 

426 kwargs.update(nsbh_exclusion=nsbh_exclusion, 

427 nsns_exclusion=nsns_exclusion) 

428 citation_index['pygrb'] = index 

429 

430 kwargs.update(citation_index=citation_index) 

431 

432 subject = ( 

433 env.get_template('medium_latency_grb_subject.jinja2').render(**kwargs) 

434 .strip()) 

435 body = ( 

436 env.get_template('medium_latency_grb_circular.jinja2').render(**kwargs) 

437 .strip()) 

438 return '{0}\n\n{1}'.format(subject, body) 

439 

440 

441def compose_update(gracedb_id, authors=(), 

442 service=rest.DEFAULT_SERVICE_URL, 

443 update_types=['sky_localization', 'p_astro', 

444 'em_bright', 'raven'], 

445 remove_text_wrap=False, 

446 client=None): 

447 """Compose GCN update circular""" 

448 if client is None: 

449 client = rest.GraceDb(service) 

450 

451 kwargs = main_dict(gracedb_id, client=client, 

452 raven_coinc='raven' in update_types) 

453 kwargs.pop('citation_index', None) 

454 kwargs.update(text_width=text_width(remove_text_wrap)) 

455 

456 if isinstance(update_types, str): 

457 update_types = update_types.split(',') 

458 

459 # Adjust files for update type alert 

460 citation_index = {} 

461 skymaps = [] 

462 index = 1 

463 

464 updated_skymap = kwargs.get('skymaps')[-1] 

465 kwargs.update(updated_skymap=updated_skymap) 

466 skymaps = [updated_skymap] 

467 if 'sky_localization' in update_types: 

468 citation_index[updated_skymap['pipeline'].lower()] = index 

469 index += 1 

470 if 'p_astro' in update_types and \ 

471 kwargs.get('p_astro_pipeline') == 'rapidpe_rift': 

472 citation_index['rapidpe_rift'] = index 

473 index += 1 

474 if 'em_bright' in update_types: 

475 # If not already cited, cite sky map pipeline for em_bright 

476 if updated_skymap['pipeline'].lower() not in citation_index.keys(): 

477 citation_index[updated_skymap['pipeline'].lower()] = index 

478 index += 1 

479 citation_index['em_bright'] = index 

480 index += 1 

481 if 'raven' in update_types: 

482 citation_index['raven'] = index 

483 

484 kwargs.update(skymaps=skymaps, 

485 citation_index=citation_index, 

486 all_pipelines=[], 

487 update_alert=True) 

488 

489 kwargs.update(authors=authors) 

490 kwargs.update(change_significance_statement=False) 

491 kwargs.update(subject='Update') 

492 kwargs.update(update_types=update_types) 

493 

494 subject = env.get_template('subject.jinja2').render(**kwargs).strip() 

495 body = env.get_template( 

496 'update_circular.jinja2').render(**kwargs).strip() 

497 return '{0}\n\n{1}'.format(subject, body) 

498 

499 

500def compose_retraction(gracedb_id, authors=(), remove_text_wrap=False, 

501 service=rest.DEFAULT_SERVICE_URL, client=None): 

502 """Compose GCN retraction circular""" 

503 if client is None: 

504 client = rest.GraceDb(service) 

505 event = client.superevent(gracedb_id).json() 

506 preferred_event = event['preferred_event_data'] 

507 labels = event['labels'] 

508 earlywarning = \ 

509 ('EARLY_WARNING' in labels and 

510 {'EM_SelectedConfident', 'SIGNIF_LOCKED'}.isdisjoint(labels)) 

511 

512 kwargs = dict( 

513 subject='Retraction', 

514 gracedb_id=gracedb_id, 

515 group=preferred_event['group'], 

516 earlywarning=earlywarning, 

517 authors=authors 

518 ) 

519 kwargs.update(text_width=text_width(remove_text_wrap)) 

520 

521 subject = env.get_template('subject.jinja2').render(**kwargs).strip() 

522 body = env.get_template('retraction.jinja2').render(**kwargs).strip() 

523 return '{0}\n\n{1}'.format(subject, body) 

524 

525 

526def read_map_gracedb(graceid, filename, client): 

527 with tempfile.NamedTemporaryFile(mode='w+b') as localfile: 

528 remotefile = client.files(graceid, filename, raw=True) 

529 shutil.copyfileobj(remotefile, localfile) 

530 localfile.flush() 

531 return read_sky_map(localfile.name, moc=True) 

532 

533 

534def get_distances_skymap_gracedb(graceid, filename, client): 

535 with tempfile.NamedTemporaryFile(mode='w+b') as localfile: 

536 remotefile = client.files(graceid, filename, raw=True) 

537 shutil.copyfileobj(remotefile, localfile) 

538 localfile.flush() 

539 header = getheader(localfile.name, 1) 

540 try: 

541 return header['distmean'], header['diststd'] 

542 except KeyError: 

543 pass 

544 

545 

546def read_map_from_path(path, client): 

547 return read_map_gracedb(*path.split('/'), client)[0] 

548 

549 

550def align_number_string(nums, positions): 

551 positions.append(len(nums[-1])) 

552 gen = (val + ' ' * (positions[i+1]-positions[i]-len(val)) 

553 for i, val in enumerate(nums)) 

554 return ''.join(gen) 

555 

556 

557def mask_cl(p, level=90): 

558 pflat = p.ravel() 

559 i = np.flipud(np.argsort(p)) 

560 cs = np.cumsum(pflat[i]) 

561 cls = np.empty_like(pflat) 

562 cls[i] = cs 

563 cls = cls.reshape(p.shape) 

564 return cls <= 1e-2 * level 

565 

566 

567def compare_skymaps(paths, service=rest.DEFAULT_SERVICE_URL, client=None): 

568 """Produce table of sky map overlaps""" 

569 if client is None: 

570 client = rest.GraceDb(service) 

571 filenames = [path.split('/')[1] for path in paths] 

572 pipelines = [guess_skyloc_pipeline(filename) for filename in filenames] 

573 probs = [read_map_from_path(path, client) for path in paths] 

574 npix = max(len(prob) for prob in probs) 

575 nside = hp.npix2nside(npix) 

576 deg2perpix = hp.nside2pixarea(nside, degrees=True) 

577 probs = [hp.ud_grade(prob, nside, power=-2) for prob in probs] 

578 masks = [mask_cl(prob) for prob in probs] 

579 areas = [mask.sum() * deg2perpix for mask in masks] 

580 joint_areas = [(mask & masks[-1]).sum() * deg2perpix for mask in masks] 

581 

582 kwargs = dict(params=zip(filenames, pipelines, areas, joint_areas)) 

583 

584 return env.get_template('compare_skymaps.jinja2').render(**kwargs) 

585 

586 

587def uncertainty_ellipse(graceid, filename, client, cls=[50, 90], 

588 ratio_ellipse_cl_areas=1.35): 

589 """Compute uncertainty ellipses for a given sky map 

590 

591 Parameters 

592 ---------- 

593 graceid: str 

594 ID of the trigger used by GraceDB 

595 filename: str 

596 File name of sky map 

597 client: class 

598 REST API client for HTTP connection 

599 cls: array-like 

600 List of percentage of minimal credible area used to check whether the 

601 areas are close to an ellipse, returning the values of the final item 

602 ratio_ellipse_cl_areas: float 

603 Ratio between ellipse area and minimal credible area from cl 

604 """ 

605 if filename.endswith('.gz'): 

606 # Try using the multi-res sky map if it exists 

607 try: 

608 new_filename = filename.replace('.fits.gz', '.multiorder.fits') 

609 skymap = read_map_gracedb(graceid, new_filename, client) 

610 except (IOError, rest.HTTPError): 

611 skymap = read_map_gracedb(graceid, filename, client) 

612 else: 

613 skymap = read_map_gracedb(graceid, filename, client) 

614 

615 # Convert to an array if necessary 

616 if np.isscalar(cls): 

617 cls = [cls] 

618 cls = np.asarray(cls) 

619 

620 # Pass array of contour inteverals to get areas 

621 result = crossmatch(skymap, contours=cls / 100) 

622 greedy_areas = np.asarray(result.contour_areas) 

623 ra, dec, a, b, pa, ellipse_areas = find_ellipse(skymap, cl=cls) 

624 a, b = np.asarray(a), np.asarray(b) 

625 

626 # Only use ellipse if every confidence interval passes 

627 use_ellipse = \ 

628 np.all(ellipse_areas <= ratio_ellipse_cl_areas * greedy_areas) 

629 return (use_ellipse, ra, dec, a[-1], b[-1], pa, ellipse_areas[-1], 

630 greedy_areas[-1]) 

631 

632 

633def _update_raven_parameters(superevent, kwargs, client, voevents_text): 

634 """Update kwargs with parameters for RAVEN coincidence""" 

635 

636 gracedb_id = superevent['superevent_id'] 

637 

638 if 'EM_COINC' not in superevent['labels']: 

639 raise ValueError("No EM_COINC label for {}".format( 

640 gracedb_id)) 

641 

642 preferred_event = superevent['preferred_event_data'] 

643 group = preferred_event['group'] 

644 gpstime = float(preferred_event['gpstime']) 

645 event_time = astropy.time.Time(gpstime, format='gps').utc 

646 em_event_id = superevent['em_type'] 

647 

648 # FIXME: Grab more info from the latest VOEvent if deemed suitable 

649 em_event = client.event(em_event_id).json() 

650 external_pipeline = em_event['pipeline'] 

651 # Get all other pipelines, removing duplicates 

652 other_ext_pipelines = \ 

653 [*set(client.event(id).json()['pipeline'] for id 

654 in superevent['em_events'])] 

655 # Remove preferred pipeline 

656 other_ext_pipelines.remove(external_pipeline) 

657 # FIXME in GraceDb: Even SNEWS triggers have an extra attribute GRB. 

658 external_trigger_id = em_event['extra_attributes']['GRB']['trigger_id'] 

659 snews = (em_event['pipeline'] == 'SNEWS') 

660 grb = (em_event['search'] in ['GRB', 'SubGRB', 'SubGRBTargeted', 'MDC'] 

661 and not snews) 

662 subthreshold = em_event['search'] in ['SubGRB', 'SubGRBTargeted'] 

663 subthreshold_targeted = em_event['search'] == 'SubGRBTargeted' 

664 far_grb = em_event['far'] 

665 

666 voevent_text_latest = voevents_text[-1] 

667 root = lxml.etree.fromstring(voevent_text_latest) 

668 time_diff = \ 

669 root.find('./What/Group/Param[@name="Time_Difference"]') 

670 time_diff = float(time_diff.attrib['value']) 

671 

672 o = urllib.parse.urlparse(client.service_url) 

673 kwargs.update( 

674 gracedb_service_url=urllib.parse.urlunsplit( 

675 (o.scheme, o.netloc, '/superevents/', '', '')), 

676 gracedb_id=gracedb_id, 

677 group=group, 

678 external_pipeline=external_pipeline, 

679 external_trigger=external_trigger_id, 

680 snews=snews, 

681 grb=grb, 

682 subthreshold=subthreshold, 

683 subthreshold_targeted=subthreshold_targeted, 

684 other_ext_pipelines=sorted(other_ext_pipelines), 

685 far_grb=far_grb, 

686 latency=abs(round(time_diff, 1)), 

687 beforeafter='before' if time_diff < 0 else 'after') 

688 

689 if grb: 

690 # Grab GRB coincidence FARs 

691 time_coinc_far = superevent['time_coinc_far'] 

692 space_time_coinc_far = superevent['space_coinc_far'] 

693 kwargs.update( 

694 time_coinc_far=time_coinc_far, 

695 space_time_coinc_far=space_time_coinc_far, 

696 ext_ra=em_event['extra_attributes']['GRB']['ra'], 

697 ext_dec=em_event['extra_attributes']['GRB']['dec'], 

698 ext_error=em_event['extra_attributes']['GRB']['error_radius']) 

699 

700 # Find combined sky maps for GRB 

701 combined_skymaps = {} 

702 for i, voevent_text in enumerate(voevents_text): 

703 root = lxml.etree.fromstring(voevent_text) 

704 alert_type = root.find( 

705 './What/Param[@name="AlertType"]').attrib['value'].lower() 

706 # Check if significant GW alert already queued or sent 

707 change_significance_statement = \ 

708 {'EM_SelectedConfident', 'SIGNIF_LOCKED'}.isdisjoint( 

709 superevent['labels'])