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

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

49 

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

51 preferred_event = event['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 = event['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 = client.files(gracedb_id, voevents[-1]['filename']).read() 

94 root = lxml.etree.fromstring(voevent_text) 

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

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

97 classifications = {} 

98 source_classification = {} 

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

100 p_astro_pipeline = None 

101 if p_astros: 

102 for p_astro in p_astros.getchildren(): 

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

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

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

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

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

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

109 for message in reversed(logs): 

110 filename = message['filename'] 

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

112 filename != 'p_astro.json': 

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

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

115 for key in classifications.keys()): 

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

117 break 

118 if p_astro_pipeline == 'rapidpe_rift': 

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

120 if em_brights: 

121 for em_bright in em_brights.getchildren(): 

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

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

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

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

126 

127 skymaps = {} 

128 for voevent in voevents: 

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

130 root = lxml.etree.fromstring(voevent_text) 

131 alert_type = root.find( 

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

133 if alert_type == 'earlywarning': 

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

135 early_warning_alert = True 

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

137 if url is None: 

138 continue 

139 url = url.attrib['value'] 

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

141 skyloc_pipeline = guess_skyloc_pipeline(filename) 

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

143 if filename not in skymaps: 

144 skymaps[filename] = dict( 

145 alert_type=alert_type, 

146 pipeline=skyloc_pipeline, 

147 filename=filename, 

148 latency=issued_time-event_time.gps) 

149 if skyloc_pipeline.lower() not in citation_index: 

150 citation_index[skyloc_pipeline.lower()] = \ 

151 max(citation_index.values()) + 1 

152 skymaps = list(skymaps.values()) 

153 

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

155 

156 kwargs = dict( 

157 subject='Identification', 

158 gracedb_id=gracedb_id, 

159 gracedb_service_url=urllib.parse.urlunsplit( 

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

161 group=preferred_event['group'], 

162 pipeline=preferred_pipeline, 

163 all_pipelines=pipelines, 

164 early_warning_alert=early_warning_alert, 

165 early_warning_pipelines=early_warning_pipelines, 

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

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

168 far=preferred_event['far'], 

169 utctime=event_time.iso, 

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

171 skymaps=skymaps, 

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

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

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

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

176 include_ellipse=None, 

177 classifications=classifications, 

178 p_astro_pipeline=p_astro_pipeline, 

179 citation_index=citation_index) 

180 

181 if skymaps: 

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

183 cls = [50, 90] 

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

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

186 kwargs.update( 

187 preferred_skymap=preferred_skymap, 

188 cl=cls[-1], 

189 include_ellipse=include_ellipse, 

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

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

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

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

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

195 ellipse_area=area, 

196 greedy_area=greedy_area) 

197 try: 

198 distmu, distsig = get_distances_skymap_gracedb(gracedb_id, 

199 preferred_skymap, 

200 client) 

201 kwargs.update( 

202 distmu=distmu, 

203 distsig=distsig) 

204 except TypeError: 

205 pass 

206 

207 return kwargs 

208 

209 

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

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

212 """Compose GCN Circular draft""" 

213 

214 if client is None: 

215 client = rest.GraceDb(service) 

216 

217 kwargs = main_dict(gracedb_id, client=client) 

218 kwargs.update(authors=authors) 

219 kwargs.update(change_significance_statement=False) 

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

221 

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

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

224 

225 if mailto: 

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

227 url = pattern.format( 

228 urllib.parse.quote(subject), 

229 urllib.parse.quote(body)) 

230 webbrowser.open(url) 

231 else: 

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

233 

234 

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

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

237 """Compose EM_COINC RAVEN GCN Circular draft""" 

238 

239 if client is None: 

240 client = rest.GraceDb(service) 

241 

242 kwargs = dict() 

243 kwargs = _update_raven_parameters(gracedb_id, kwargs, client) 

244 kwargs.update(main_dict(gracedb_id, client=client)) 

245 kwargs.update(update_alert=False) 

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

247 # Add RAVEN citation 

248 citation_index = kwargs['citation_index'] 

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

250 kwargs['citation_index'] = citation_index 

251 

252 subject = ( 

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

254 .strip()) 

255 body = ( 

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

257 .strip()) 

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

259 

260 

261def compose_llama( 

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

263 icecube_alert=None, remove_text_wrap=False, 

264 client=None): 

265 """Compose GRB LLAMA GCN Circular draft. 

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

267 

268 if client is None: 

269 client = rest.GraceDb(service) 

270 

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

272 

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

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

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

276 tl_datetime = str(astropy.time.Time( 

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

278 th_datetime = str(astropy.time.Time( 

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

280 

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

282 kwargs = dict( 

283 gracedb_service_url=urllib.parse.urlunsplit( 

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

285 gracedb_id=gracedb_id, 

286 llama=True, 

287 icecube_alert=icecube_alert, 

288 event_time=event_time, 

289 tl_datetime=tl_datetime, 

290 th_datetime=th_datetime, 

291 authors=authors) 

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

293 

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

295 kwargs.update(citation_index=citation_index) 

296 

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

298 

299 llama_stat_filename = 'significance_subthreshold_lvc-i3.json' 

300 if llama_stat_filename in files: 

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

302 llama_fap = llama_stat_file["p_value"] 

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

304 lines = [] 

305 for neutrino in neutrinos: 

306 # Build aligned string 

307 vals = [] 

308 dt = (event_time - 

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

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

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

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

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

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

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

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

317 

318 kwargs.update(llama_fap=llama_fap, 

319 neutrinos=lines) 

320 

321 subject = ( 

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

323 .strip()) 

324 if icecube_alert: 

325 body = ( 

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

327 .strip()) 

328 else: 

329 body = ( 

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

331 .strip()) 

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

333 

334 

335def compose_grb_medium_latency( 

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

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

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

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

340 

341 if client is None: 

342 client = rest.GraceDb(service) 

343 

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

345 search = event['search'] 

346 

347 if search != 'GRB': 

348 return 

349 

350 group = event['group'] 

351 pipeline = event['pipeline'] 

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

353 

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

355 kwargs = dict( 

356 gracedb_service_url=urllib.parse.urlunsplit( 

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

358 gracedb_id=gracedb_id, 

359 group=group, 

360 grb=True, 

361 pipeline=pipeline, 

362 external_trigger=external_trigger, 

363 exclusions=[], 

364 detections=[]) 

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

366 

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

368 

369 citation_index = {} 

370 index = 1 

371 xpipeline_fap_file = 'false_alarm_probability_x.json' 

372 if xpipeline_fap_file in files: 

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

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

375 # Create detection/exclusion circular based on given argument 

376 # Use default cutoff if not provided 

377 xpipeline_detection = (use_detection_template if use_detection_template 

378 is not None else online_xpipeline_fap < 0.001) 

379 if xpipeline_detection: 

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

381 kwargs.update(online_xpipeline_fap=online_xpipeline_fap) 

382 else: 

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

384 xpipeline_distances_file = 'distances_x.json' 

385 xpipeline_distances = client.files(gracedb_id, 

386 xpipeline_distances_file).json() 

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

388 kwargs.update(burst_exclusion=burst_exclusion) 

389 citation_index['xpipeline'] = index 

390 index += 1 

391 

392 pygrb_fap_file = 'false_alarm_probability_pygrb.json' 

393 if pygrb_fap_file in files: 

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

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

396 # Create detection/exclusion circular based on given argument 

397 # Use default cutoff if not provided 

398 pygrb_detection = (use_detection_template if use_detection_template 

399 is not None else online_pygrb_fap < 0.01) 

400 if pygrb_detection: 

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

402 kwargs.update(online_pygrb_fap=online_pygrb_fap) 

403 else: 

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

405 pygrb_distances_file = 'distances_pygrb.json' 

406 pygrb_distances = client.files(gracedb_id, 

407 pygrb_distances_file).json() 

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

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

410 kwargs.update(nsbh_exclusion=nsbh_exclusion, 

411 nsns_exclusion=nsns_exclusion) 

412 citation_index['pygrb'] = index 

413 

414 kwargs.update(citation_index=citation_index) 

415 

416 subject = ( 

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

418 .strip()) 

419 body = ( 

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

421 .strip()) 

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

423 

424 

425def compose_update(gracedb_id, authors=(), 

426 service=rest.DEFAULT_SERVICE_URL, 

427 update_types=['sky_localization', 'p_astro', 

428 'em_bright', 'raven'], 

429 remove_text_wrap=False, 

430 client=None): 

431 """Compose GCN update circular""" 

432 if client is None: 

433 client = rest.GraceDb(service) 

434 

435 kwargs = main_dict(gracedb_id, client=client) 

436 kwargs.pop('citation_index', None) 

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

438 

439 if isinstance(update_types, str): 

440 update_types = update_types.split(',') 

441 

442 # Adjust files for update type alert 

443 citation_index = {} 

444 skymaps = [] 

445 index = 1 

446 

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

448 kwargs.update(updated_skymap=updated_skymap) 

449 skymaps = [updated_skymap] 

450 if 'sky_localization' in update_types: 

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

452 index += 1 

453 if 'p_astro' in update_types and \ 

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

455 citation_index['rapidpe_rift'] = index 

456 index += 1 

457 if 'em_bright' in update_types: 

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

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

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

461 index += 1 

462 citation_index['em_bright'] = index 

463 index += 1 

464 if 'raven' in update_types: 

465 citation_index['raven'] = index 

466 

467 kwargs.update(skymaps=skymaps, 

468 citation_index=citation_index, 

469 all_pipelines=[], 

470 update_alert=True) 

471 

472 if 'raven' in update_types: 

473 kwargs = _update_raven_parameters(gracedb_id, kwargs, client) 

474 

475 kwargs.update(authors=authors) 

476 kwargs.update(change_significance_statement=False) 

477 kwargs.update(subject='Update') 

478 kwargs.update(update_types=update_types) 

479 

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

481 body = env.get_template( 

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

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

484 

485 

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

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

488 """Compose GCN retraction circular""" 

489 if client is None: 

490 client = rest.GraceDb(service) 

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

492 preferred_event = event['preferred_event_data'] 

493 labels = event['labels'] 

494 earlywarning = \ 

495 ('EARLY_WARNING' in labels and 

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

497 

498 kwargs = dict( 

499 subject='Retraction', 

500 gracedb_id=gracedb_id, 

501 group=preferred_event['group'], 

502 earlywarning=earlywarning, 

503 authors=authors 

504 ) 

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

506 

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

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

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

510 

511 

512def read_map_gracedb(graceid, filename, client): 

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

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

515 shutil.copyfileobj(remotefile, localfile) 

516 localfile.flush() 

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

518 

519 

520def get_distances_skymap_gracedb(graceid, filename, client): 

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

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

523 shutil.copyfileobj(remotefile, localfile) 

524 localfile.flush() 

525 header = getheader(localfile.name, 1) 

526 try: 

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

528 except KeyError: 

529 pass 

530 

531 

532def read_map_from_path(path, client): 

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

534 

535 

536def align_number_string(nums, positions): 

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

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

539 for i, val in enumerate(nums)) 

540 return ''.join(gen) 

541 

542 

543def mask_cl(p, level=90): 

544 pflat = p.ravel() 

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

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

547 cls = np.empty_like(pflat) 

548 cls[i] = cs 

549 cls = cls.reshape(p.shape) 

550 return cls <= 1e-2 * level 

551 

552 

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

554 """Produce table of sky map overlaps""" 

555 if client is None: 

556 client = rest.GraceDb(service) 

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

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

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

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

561 nside = hp.npix2nside(npix) 

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

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

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

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

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

567 

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

569 

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

571 

572 

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

574 ratio_ellipse_cl_areas=1.35): 

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

576 

577 Parameters 

578 ---------- 

579 graceid: str 

580 ID of the trigger used by GraceDB 

581 filename: str 

582 File name of sky map 

583 client: class 

584 REST API client for HTTP connection 

585 cls: array-like 

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

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

588 ratio_ellipse_cl_areas: float 

589 Ratio between ellipse area and minimal credible area from cl 

590 """ 

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

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

593 try: 

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

595 skymap = read_map_gracedb(graceid, new_filename, client) 

596 except (IOError, rest.HTTPError): 

597 skymap = read_map_gracedb(graceid, filename, client) 

598 else: 

599 skymap = read_map_gracedb(graceid, filename, client) 

600 

601 # Convert to an array if necessary 

602 if np.isscalar(cls): 

603 cls = [cls] 

604 cls = np.asarray(cls) 

605 

606 # Pass array of contour inteverals to get areas 

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

608 greedy_areas = np.asarray(result.contour_areas) 

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

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

611 

612 # Only use ellipse if every confidence interval passes 

613 use_ellipse = \ 

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

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

616 greedy_areas[-1]) 

617 

618 

619def _update_raven_parameters(gracedb_id, kwargs, client): 

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

621 

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

623 

624 if 'EM_COINC' not in event['labels']: 

625 raise ValueError("No EM_COINC label for {}".format(gracedb_id)) 

626 

627 preferred_event = event['preferred_event_data'] 

628 group = preferred_event['group'] 

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

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

631 

632 em_event_id = event['em_type'] 

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

634 em_event_gpstime = float(em_event['gpstime']) 

635 external_pipeline = em_event['pipeline'] 

636 # Get all other pipelines, removing duplicates 

637 other_ext_pipelines = \ 

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

639 in event['em_events'])] 

640 # Remove preferred pipeline 

641 other_ext_pipelines.remove(external_pipeline) 

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

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

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

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

646 and not snews) 

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

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

649 far_grb = em_event['far'] 

650 

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

652 kwargs.update( 

653 gracedb_service_url=urllib.parse.urlunsplit( 

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

655 gracedb_id=gracedb_id, 

656 group=group, 

657 external_pipeline=external_pipeline, 

658 external_trigger=external_trigger_id, 

659 snews=snews, 

660 grb=grb, 

661 subthreshold=subthreshold, 

662 subthreshold_targeted=subthreshold_targeted, 

663 other_ext_pipelines=sorted(other_ext_pipelines), 

664 far_grb=far_grb, 

665 latency=abs(round(em_event_gpstime-gpstime, 1)), 

666 beforeafter='before' if gpstime > em_event_gpstime else 'after') 

667 

668 if grb: 

669 # Grab GRB coincidence FARs 

670 time_coinc_far = event['time_coinc_far'] 

671 space_time_coinc_far = event['space_coinc_far'] 

672 kwargs.update( 

673 time_coinc_far=time_coinc_far, 

674 space_time_coinc_far=space_time_coinc_far, 

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

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

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

678 

679 # Find combined sky maps for GRB 

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

681 combined_skymaps = {} 

682 if not voevents: 

683 raise ValueError( 

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

685 gracedb_id)) 

686 for i, voevent in enumerate(voevents): 

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

688 root = lxml.etree.fromstring(voevent_text) 

689 alert_type = root.find( 

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

691 # Check if significant GW alert already queued or sent 

692 change_significance_statement = \ 

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

694 event['labels']) 

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

696 if url is None: 

697 continue 

698 url = url.attrib['value'] 

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

700 issued_time = astropy.time.Time( 

701 root.find('./Who/Date').text).gps 

702 if filename not in combined_skymaps: 

703 combined_skymaps[filename] = dict(