Coverage for ligo/followup_advocate/__init__.py: 91%
393 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-12-11 21:40 +0000
1import math as mth
2import os
3import shutil
4import tempfile
5import urllib
6import webbrowser
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
20from .jinja import env
21from . import _version
22__version__ = _version.get_versions()['version']
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()
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]
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
47def main_dict(gracedb_id, client):
48 """Create general dictionary to pass to compose circular"""
50 event = client.superevent(gracedb_id).json()
51 preferred_event = event['preferred_event_data']
52 preferred_pipeline = preferred_event['pipeline']
53 early_warning_pipelines = []
54 pipelines = []
55 gw_events = event['gw_events']
56 early_warning_alert = False
58 for gw_event in gw_events:
59 gw_event_dict = client.event(gw_event).json()
60 pipeline = gw_event_dict['pipeline']
61 search = gw_event_dict['search']
62 if pipeline not in pipelines:
63 pipelines.append(pipeline)
64 if pipeline not in early_warning_pipelines and \
65 search == 'EarlyWarning':
66 early_warning_pipelines.append(pipeline)
67 # Sort to get alphabetical order
68 pipelines.sort(key=str.lower)
69 early_warning_pipelines.sort(key=str.lower)
71 voevents = client.voevents(gracedb_id).json()['voevents']
72 if not voevents:
73 raise ValueError(
74 "{} has no VOEvent to generate circulars from.".format(
75 gracedb_id))
77 citation_index = {pipeline.lower(): pipelines.index(pipeline) + 1 for
78 pipeline in pipelines}
79 for pipeline in early_warning_pipelines:
80 if pipeline.lower() != 'mbta':
81 citation_index[pipeline.lower() + '_earlywarning'] = \
82 max(citation_index.values()) + 1
84 gpstime = float(preferred_event['gpstime'])
85 event_time = astropy.time.Time(gpstime, format='gps').utc
87 # Grab latest p_astro and em_bright values from lastest VOEvent
88 voevent_text = client.files(gracedb_id, voevents[-1]['filename']).read()
89 root = lxml.etree.fromstring(voevent_text)
90 p_astros = root.find('./What/Group[@name="Classification"]')
91 em_brights = root.find('./What/Group[@name="Properties"]')
92 classifications = {}
93 source_classification = {}
94 # Only try to load if present to prevent errors with .getchildren()
95 p_astro_pipeline = None
96 if p_astros:
97 for p_astro in p_astros.getchildren():
98 if p_astro.attrib.get('value'):
99 classifications[p_astro.attrib['name']] = \
100 float(p_astro.attrib['value']) * 100
101 # Figure out which pipeline uploaded p_astro, usually the first one
102 # FIXME: Replace with more efficient method in the future
103 logs = client.logs(gracedb_id).json()['log']
104 for message in reversed(logs):
105 filename = message['filename']
106 if filename and '.p_astro.json' in filename and \
107 filename != 'p_astro.json':
108 p_astro = client.files(gracedb_id, filename).json()
109 if all(mth.isclose(p_astro[key] * 100, classifications[key])
110 for key in classifications.keys()):
111 p_astro_pipeline = filename.split('.')[0].lower()
112 break
113 if p_astro_pipeline == 'rapidpe_rift':
114 citation_index['rapidpe_rift'] = max(citation_index.values()) + 1
115 if em_brights:
116 for em_bright in em_brights.getchildren():
117 if em_bright.attrib.get('value'):
118 source_classification[em_bright.attrib['name']] = \
119 float(em_bright.attrib['value']) * 100
120 citation_index['em_bright'] = max(citation_index.values()) + 1
122 skymaps = {}
123 for voevent in voevents:
124 voevent_text = client.files(gracedb_id, voevent['filename']).read()
125 root = lxml.etree.fromstring(voevent_text)
126 alert_type = root.find(
127 './What/Param[@name="AlertType"]').attrib['value'].lower()
128 if alert_type == 'earlywarning':
129 # Add text for early warning detection if one early warning alert
130 early_warning_alert = True
131 url = root.find('./What/Group/Param[@name="skymap_fits"]')
132 if url is None:
133 continue
134 url = url.attrib['value']
135 _, filename = os.path.split(url)
136 skyloc_pipeline = guess_skyloc_pipeline(filename)
137 issued_time = astropy.time.Time(root.find('./Who/Date').text).gps
138 if filename not in skymaps:
139 skymaps[filename] = dict(
140 alert_type=alert_type,
141 pipeline=skyloc_pipeline,
142 filename=filename,
143 latency=issued_time-event_time.gps)
144 if skyloc_pipeline.lower() not in citation_index:
145 citation_index[skyloc_pipeline.lower()] = \
146 max(citation_index.values()) + 1
147 skymaps = list(skymaps.values())
149 o = urllib.parse.urlparse(client.service_url)
151 kwargs = dict(
152 subject='Identification',
153 gracedb_id=gracedb_id,
154 gracedb_service_url=urllib.parse.urlunsplit(
155 (o.scheme, o.netloc, '/superevents/', '', '')),
156 group=preferred_event['group'],
157 pipeline=preferred_pipeline,
158 all_pipelines=pipelines,
159 early_warning_alert=early_warning_alert,
160 early_warning_pipelines=early_warning_pipelines,
161 gpstime='{0:.03f}'.format(round(float(preferred_event['gpstime']), 3)),
162 search=preferred_event.get('search', ''),
163 far=preferred_event['far'],
164 utctime=event_time.iso,
165 instruments=preferred_event['instruments'].split(','),
166 skymaps=skymaps,
167 prob_has_ns=source_classification.get('HasNS'),
168 prob_has_remnant=source_classification.get('HasRemnant'),
169 prob_has_massgap=source_classification.get('HasMassGap'),
170 include_ellipse=None,
171 classifications=classifications,
172 p_astro_pipeline=p_astro_pipeline,
173 citation_index=citation_index)
175 if skymaps:
176 preferred_skymap = skymaps[-1]['filename']
177 cls = [50, 90]
178 include_ellipse, ra, dec, a, b, pa, area, greedy_area = \
179 uncertainty_ellipse(gracedb_id, preferred_skymap, client, cls=cls)
180 kwargs.update(
181 preferred_skymap=preferred_skymap,
182 cl=cls[-1],
183 include_ellipse=include_ellipse,
184 ra=coord.Longitude(ra*u.deg),
185 dec=coord.Latitude(dec*u.deg),
186 a=coord.Angle(a*u.deg),
187 b=coord.Angle(b*u.deg),
188 pa=coord.Angle(pa*u.deg),
189 ellipse_area=area,
190 greedy_area=greedy_area)
191 try:
192 distmu, distsig = get_distances_skymap_gracedb(gracedb_id,
193 preferred_skymap,
194 client)
195 kwargs.update(
196 distmu=distmu,
197 distsig=distsig)
198 except TypeError:
199 pass
201 return kwargs
204def compose(gracedb_id, authors=(), mailto=False, remove_text_wrap=False,
205 service=rest.DEFAULT_SERVICE_URL, client=None):
206 """Compose GCN Circular draft"""
208 if client is None:
209 client = rest.GraceDb(service)
211 kwargs = main_dict(gracedb_id, client=client)
212 kwargs.update(authors=authors)
213 kwargs.update(change_significance_statement=False)
214 kwargs.update(text_width=text_width(remove_text_wrap))
216 subject = env.get_template('subject.jinja2').render(**kwargs).strip()
217 body = env.get_template('initial_circular.jinja2').render(**kwargs).strip()
219 if mailto:
220 pattern = 'mailto:emfollow@ligo.org,dac@ligo.org?subject={0}&body={1}'
221 url = pattern.format(
222 urllib.parse.quote(subject),
223 urllib.parse.quote(body))
224 webbrowser.open(url)
225 else:
226 return '{0}\n\n{1}'.format(subject, body)
229def compose_raven(gracedb_id, authors=(), remove_text_wrap=False,
230 service=rest.DEFAULT_SERVICE_URL, client=None):
231 """Compose EM_COINC RAVEN GCN Circular draft"""
233 if client is None:
234 client = rest.GraceDb(service)
236 kwargs = dict()
237 kwargs = _update_raven_parameters(gracedb_id, kwargs, client)
238 kwargs.update(main_dict(gracedb_id, client=client))
239 kwargs.update(update_alert=False)
240 kwargs.update(text_width=text_width(remove_text_wrap))
241 # Add RAVEN citation
242 citation_index = kwargs['citation_index']
243 citation_index['raven'] = max(citation_index.values()) + 1
244 kwargs['citation_index'] = citation_index
246 subject = (
247 env.get_template('RAVEN_subject.jinja2').render(**kwargs)
248 .strip())
249 body = (
250 env.get_template('RAVEN_circular.jinja2').render(**kwargs)
251 .strip())
252 return '{0}\n\n{1}'.format(subject, body)
255def compose_llama(
256 gracedb_id, authors=(), service=rest.DEFAULT_SERVICE_URL,
257 icecube_alert=None, remove_text_wrap=False,
258 client=None):
259 """Compose GRB LLAMA GCN Circular draft.
260 Here, gracedb_id will be a GRB superevent ID in GraceDb."""
262 if client is None:
263 client = rest.GraceDb(service)
265 superevent = client.superevent(gracedb_id).json()
267 gpstime = float(superevent['t_0'])
268 tl, th = gpstime - 500, gpstime + 500
269 event_time = astropy.time.Time(gpstime, format='gps').utc
270 tl_datetime = str(astropy.time.Time(
271 tl, format='gps').isot).replace('T', ' ')
272 th_datetime = str(astropy.time.Time(
273 th, format='gps').isot).replace('T', ' ')
275 o = urllib.parse.urlparse(client.service_url)
276 kwargs = dict(
277 gracedb_service_url=urllib.parse.urlunsplit(
278 (o.scheme, o.netloc, '/superevents/', '', '')),
279 gracedb_id=gracedb_id,
280 llama=True,
281 icecube_alert=icecube_alert,
282 event_time=event_time,
283 tl_datetime=tl_datetime,
284 th_datetime=th_datetime,
285 authors=authors)
286 kwargs.update(text_width=text_width(remove_text_wrap))
288 citation_index = {'llama': 1, 'llama_stat': 2}
289 kwargs.update(citation_index=citation_index)
291 files = client.files(gracedb_id).json()
293 llama_stat_filename = 'significance_subthreshold_lvc-i3.json'
294 if llama_stat_filename in files:
295 llama_stat_file = client.files(gracedb_id, llama_stat_filename).json()
296 llama_fap = llama_stat_file["p_value"]
297 neutrinos = llama_stat_file["inputs"]["neutrino_info"]
298 lines = []
299 for neutrino in neutrinos:
300 # Build aligned string
301 vals = []
302 dt = (event_time -
303 astropy.time.Time(neutrino['mjd'],
304 format='mjd')).to(u.s).value
305 vals.append('{:.2f}'.format(dt))
306 vals.append('{:.2f}'.format(neutrino['ra']))
307 vals.append('{:.2f}'.format(neutrino['dec']))
308 vals.append('{:.2f}'.format(neutrino['sigma']))
309 vals.append('{:.4f}'.format(llama_fap))
310 lines.append(align_number_string(vals, [0, 11, 21, 40, 59]))
312 kwargs.update(llama_fap=llama_fap,
313 neutrinos=lines)
315 subject = (
316 env.get_template('llama_subject.jinja2').render(**kwargs)
317 .strip())
318 if icecube_alert:
319 body = (
320 env.get_template('llama_alert_circular.jinja2').render(**kwargs)
321 .strip())
322 else:
323 body = (
324 env.get_template('llama_track_circular.jinja2').render(**kwargs)
325 .strip())
326 return '{0}\n\n{1}'.format(subject, body)
329def compose_grb_medium_latency(
330 gracedb_id, authors=(), service=rest.DEFAULT_SERVICE_URL,
331 use_detection_template=None, remove_text_wrap=False, client=None):
332 """Compose GRB Medium Latency GCN Circular draft.
333 Here, gracedb_id will be a GRB external trigger ID in GraceDb."""
335 if client is None:
336 client = rest.GraceDb(service)
338 event = client.event(gracedb_id).json()
339 search = event['search']
341 if search != 'GRB':
342 return
344 group = event['group']
345 pipeline = event['pipeline']
346 external_trigger = event['extra_attributes']['GRB']['trigger_id']
348 o = urllib.parse.urlparse(client.service_url)
349 kwargs = dict(
350 gracedb_service_url=urllib.parse.urlunsplit(
351 (o.scheme, o.netloc, '/events/', '', '')),
352 gracedb_id=gracedb_id,
353 group=group,
354 grb=True,
355 pipeline=pipeline,
356 external_trigger=external_trigger,
357 exclusions=[],
358 detections=[])
359 kwargs.update(text_width=text_width(remove_text_wrap))
361 files = client.files(gracedb_id).json()
363 citation_index = {}
364 index = 1
365 xpipeline_fap_file = 'false_alarm_probability_x.json'
366 if xpipeline_fap_file in files:
367 xpipeline_fap = client.files(gracedb_id, xpipeline_fap_file).json()
368 online_xpipeline_fap = xpipeline_fap.get('Online Xpipeline')
369 # Create detection/exclusion circular based on given argument
370 # Use default cutoff if not provided
371 xpipeline_detection = (use_detection_template if use_detection_template
372 is not None else online_xpipeline_fap < 0.001)
373 if xpipeline_detection:
374 kwargs['detections'].append('xpipeline')
375 kwargs.update(online_xpipeline_fap=online_xpipeline_fap)
376 else:
377 kwargs['exclusions'].append('xpipeline')
378 xpipeline_distances_file = 'distances_x.json'
379 xpipeline_distances = client.files(gracedb_id,
380 xpipeline_distances_file).json()
381 burst_exclusion = xpipeline_distances.get('Burst Exclusion')
382 kwargs.update(burst_exclusion=burst_exclusion)
383 citation_index['xpipeline'] = index
384 index += 1
386 pygrb_fap_file = 'false_alarm_probability_pygrb.json'
387 if pygrb_fap_file in files:
388 pygrb_fap = client.files(gracedb_id, pygrb_fap_file).json()
389 online_pygrb_fap = pygrb_fap.get('Online PyGRB')
390 # Create detection/exclusion circular based on given argument
391 # Use default cutoff if not provided
392 pygrb_detection = (use_detection_template if use_detection_template
393 is not None else online_pygrb_fap < 0.01)
394 if pygrb_detection:
395 kwargs['detections'].append('pygrb')
396 kwargs.update(online_pygrb_fap=online_pygrb_fap)
397 else:
398 kwargs['exclusions'].append('pygrb')
399 pygrb_distances_file = 'distances_pygrb.json'
400 pygrb_distances = client.files(gracedb_id,
401 pygrb_distances_file).json()
402 nsns_exclusion = pygrb_distances.get('NSNS Exclusion')
403 nsbh_exclusion = pygrb_distances.get('NSBH Exclusion')
404 kwargs.update(nsbh_exclusion=nsbh_exclusion,
405 nsns_exclusion=nsns_exclusion)
406 citation_index['pygrb'] = index
408 kwargs.update(citation_index=citation_index)
410 subject = (
411 env.get_template('medium_latency_grb_subject.jinja2').render(**kwargs)
412 .strip())
413 body = (
414 env.get_template('medium_latency_grb_circular.jinja2').render(**kwargs)
415 .strip())
416 return '{0}\n\n{1}'.format(subject, body)
419def compose_update(gracedb_id, authors=(),
420 service=rest.DEFAULT_SERVICE_URL,
421 update_types=['sky_localization', 'p_astro',
422 'em_bright', 'raven'],
423 remove_text_wrap=False,
424 client=None):
425 """Compose GCN update circular"""
426 if client is None:
427 client = rest.GraceDb(service)
429 kwargs = main_dict(gracedb_id, client=client)
430 kwargs.pop('citation_index', None)
431 kwargs.update(text_width=text_width(remove_text_wrap))
433 if isinstance(update_types, str):
434 update_types = update_types.split(',')
436 # Adjust files for update type alert
437 citation_index = {}
438 skymaps = []
439 index = 1
441 updated_skymap = kwargs.get('skymaps')[-1]
442 kwargs.update(updated_skymap=updated_skymap)
443 skymaps = [updated_skymap]
444 if 'sky_localization' in update_types:
445 citation_index[updated_skymap['pipeline'].lower()] = index
446 index += 1
447 if 'p_astro' in update_types and \
448 kwargs.get('p_astro_pipeline') == 'rapidpe_rift':
449 citation_index['rapidpe_rift'] = index
450 index += 1
451 if 'em_bright' in update_types:
452 # If not already cited, cite sky map pipeline for em_bright
453 if updated_skymap['pipeline'].lower() not in citation_index.keys():
454 citation_index[updated_skymap['pipeline'].lower()] = index
455 index += 1
456 citation_index['em_bright'] = index
457 index += 1
458 if 'raven' in update_types:
459 citation_index['raven'] = index
461 kwargs.update(skymaps=skymaps,
462 citation_index=citation_index,
463 all_pipelines=[],
464 update_alert=True)
466 if 'raven' in update_types:
467 kwargs = _update_raven_parameters(gracedb_id, kwargs, client)
469 kwargs.update(authors=authors)
470 kwargs.update(change_significance_statement=False)
471 kwargs.update(subject='Update')
472 kwargs.update(update_types=update_types)
474 subject = env.get_template('subject.jinja2').render(**kwargs).strip()
475 body = env.get_template(
476 'update_circular.jinja2').render(**kwargs).strip()
477 return '{0}\n\n{1}'.format(subject, body)
480def compose_retraction(gracedb_id, authors=(), remove_text_wrap=False,
481 service=rest.DEFAULT_SERVICE_URL, client=None):
482 """Compose GCN retraction circular"""
483 if client is None:
484 client = rest.GraceDb(service)
485 event = client.superevent(gracedb_id).json()
486 preferred_event = event['preferred_event_data']
487 labels = event['labels']
488 earlywarning = \
489 ('EARLY_WARNING' in labels and
490 {'EM_SelectedConfident', 'SIGNIF_LOCKED'}.isdisjoint(labels))
492 kwargs = dict(
493 subject='Retraction',
494 gracedb_id=gracedb_id,
495 group=preferred_event['group'],
496 earlywarning=earlywarning,
497 authors=authors
498 )
499 kwargs.update(text_width=text_width(remove_text_wrap))
501 subject = env.get_template('subject.jinja2').render(**kwargs).strip()
502 body = env.get_template('retraction.jinja2').render(**kwargs).strip()
503 return '{0}\n\n{1}'.format(subject, body)
506def read_map_gracedb(graceid, filename, client):
507 with tempfile.NamedTemporaryFile(mode='w+b') as localfile:
508 remotefile = client.files(graceid, filename, raw=True)
509 shutil.copyfileobj(remotefile, localfile)
510 localfile.flush()
511 return read_sky_map(localfile.name, moc=True)
514def get_distances_skymap_gracedb(graceid, filename, client):
515 with tempfile.NamedTemporaryFile(mode='w+b') as localfile:
516 remotefile = client.files(graceid, filename, raw=True)
517 shutil.copyfileobj(remotefile, localfile)
518 localfile.flush()
519 header = getheader(localfile.name, 1)
520 try:
521 return header['distmean'], header['diststd']
522 except KeyError:
523 pass
526def read_map_from_path(path, client):
527 return read_map_gracedb(*path.split('/'), client)[0]
530def align_number_string(nums, positions):
531 positions.append(len(nums[-1]))
532 gen = (val + ' ' * (positions[i+1]-positions[i]-len(val))
533 for i, val in enumerate(nums))
534 return ''.join(gen)
537def mask_cl(p, level=90):
538 pflat = p.ravel()
539 i = np.flipud(np.argsort(p))
540 cs = np.cumsum(pflat[i])
541 cls = np.empty_like(pflat)
542 cls[i] = cs
543 cls = cls.reshape(p.shape)
544 return cls <= 1e-2 * level
547def compare_skymaps(paths, service=rest.DEFAULT_SERVICE_URL, client=None):
548 """Produce table of sky map overlaps"""
549 if client is None:
550 client = rest.GraceDb(service)
551 filenames = [path.split('/')[1] for path in paths]
552 pipelines = [guess_skyloc_pipeline(filename) for filename in filenames]
553 probs = [read_map_from_path(path, client) for path in paths]
554 npix = max(len(prob) for prob in probs)
555 nside = hp.npix2nside(npix)
556 deg2perpix = hp.nside2pixarea(nside, degrees=True)
557 probs = [hp.ud_grade(prob, nside, power=-2) for prob in probs]
558 masks = [mask_cl(prob) for prob in probs]
559 areas = [mask.sum() * deg2perpix for mask in masks]
560 joint_areas = [(mask & masks[-1]).sum() * deg2perpix for mask in masks]
562 kwargs = dict(params=zip(filenames, pipelines, areas, joint_areas))
564 return env.get_template('compare_skymaps.jinja2').render(**kwargs)
567def uncertainty_ellipse(graceid, filename, client, cls=[50, 90],
568 ratio_ellipse_cl_areas=1.35):
569 """Compute uncertainty ellipses for a given sky map
571 Parameters
572 ----------
573 graceid: str
574 ID of the trigger used by GraceDB
575 filename: str
576 File name of sky map
577 client: class
578 REST API client for HTTP connection
579 cls: array-like
580 List of percentage of minimal credible area used to check whether the
581 areas are close to an ellipse, returning the values of the final item
582 ratio_ellipse_cl_areas: float
583 Ratio between ellipse area and minimal credible area from cl
584 """
585 if filename.endswith('.gz'):
586 # Try using the multi-res sky map if it exists
587 try:
588 new_filename = filename.replace('.fits.gz', '.multiorder.fits')
589 skymap = read_map_gracedb(graceid, new_filename, client)
590 except (IOError, rest.HTTPError):
591 skymap = read_map_gracedb(graceid, filename, client)
592 else:
593 skymap = read_map_gracedb(graceid, filename, client)
595 # Convert to an array if necessary
596 if np.isscalar(cls):
597 cls = [cls]
598 cls = np.asarray(cls)
600 # Pass array of contour inteverals to get areas
601 result = crossmatch(skymap, contours=cls / 100)
602 greedy_areas = np.asarray(result.contour_areas)
603 ra, dec, a, b, pa, ellipse_areas = find_ellipse(skymap, cl=cls)
604 a, b = np.asarray(a), np.asarray(b)
606 # Only use ellipse if every confidence interval passes
607 use_ellipse = \
608 np.all(ellipse_areas <= ratio_ellipse_cl_areas * greedy_areas)
609 return (use_ellipse, ra, dec, a[-1], b[-1], pa, ellipse_areas[-1],
610 greedy_areas[-1])
613def _update_raven_parameters(gracedb_id, kwargs, client):
614 """Update kwargs with parameters for RAVEN coincidence"""
616 event = client.superevent(gracedb_id).json()
618 if 'EM_COINC' not in event['labels']:
619 raise ValueError("No EM_COINC label for {}".format(gracedb_id))
621 preferred_event = event['preferred_event_data']
622 group = preferred_event['group']
623 gpstime = float(preferred_event['gpstime'])
624 event_time = astropy.time.Time(gpstime, format='gps').utc
626 em_event_id = event['em_type']
627 em_event = client.event(em_event_id).json()
628 em_event_gpstime = float(em_event['gpstime'])
629 external_pipeline = em_event['pipeline']
630 # Get all other pipelines, removing duplicates
631 other_ext_pipelines = \
632 [*set(client.event(id).json()['pipeline'] for id
633 in event['em_events'])]
634 # Remove preferred pipeline
635 other_ext_pipelines.remove(external_pipeline)
636 # FIXME in GraceDb: Even SNEWS triggers have an extra attribute GRB.
637 external_trigger_id = em_event['extra_attributes']['GRB']['trigger_id']
638 snews = (em_event['pipeline'] == 'SNEWS')
639 grb = (em_event['search'] in ['GRB', 'SubGRB', 'SubGRBTargeted', 'MDC']
640 and not snews)
641 subthreshold = em_event['search'] in ['SubGRB', 'SubGRBTargeted']
642 subthreshold_targeted = em_event['search'] == 'SubGRBTargeted'
643 far_grb = em_event['far']
645 o = urllib.parse.urlparse(client.service_url)
646 kwargs.update(
647 gracedb_service_url=urllib.parse.urlunsplit(
648 (o.scheme, o.netloc, '/superevents/', '', '')),
649 gracedb_id=gracedb_id,
650 group=group,
651 external_pipeline=external_pipeline,
652 external_trigger=external_trigger_id,
653 snews=snews,
654 grb=grb,
655 subthreshold=subthreshold,
656 subthreshold_targeted=subthreshold_targeted,
657 other_ext_pipelines=other_ext_pipelines,
658 far_grb=far_grb,
659 latency=abs(round(em_event_gpstime-gpstime, 1)),
660 beforeafter='before' if gpstime > em_event_gpstime else 'after')
662 if grb:
663 # Grab GRB coincidence FARs
664 time_coinc_far = event['time_coinc_far']
665 space_time_coinc_far = event['space_coinc_far']
666 kwargs.update(
667 time_coinc_far=time_coinc_far,
668 space_time_coinc_far=space_time_coinc_far,
669 ext_ra=em_event['extra_attributes']['GRB']['ra'],
670 ext_dec=em_event['extra_attributes']['GRB']['dec'],
671 ext_error=em_event['extra_attributes']['GRB']['error_radius'])
673 # Find combined sky maps for GRB
674 voevents = client.voevents(gracedb_id).json()['voevents']
675 combined_skymaps = {}
676 if not voevents:
677 raise ValueError(
678 "{} has no VOEvent to generate circulars from.".format(
679 gracedb_id))
680 for i, voevent in enumerate(voevents):
681 voevent_text = client.files(gracedb_id, voevent['filename']).read()
682 root = lxml.etree.fromstring(voevent_text)
683 alert_type = root.find(
684 './What/Param[@name="AlertType"]').attrib['value'].lower()
685 # Check if significant GW alert already queued or sent
686 change_significance_statement = \
687 {'EM_SelectedConfident', 'SIGNIF_LOCKED'}.isdisjoint(
688 event['labels'])
689 url = root.find('./What/Group/Param[@name="joint_skymap_fits"]')
690 if url is None:
691 continue
692 url = url.attrib['value']
693 _, filename = os.path.split(url)
694 issued_time = astropy.time.Time(
695 root.find('./Who/Date').text).gps
696 if filename not in combined_skymaps:
697 combined_skymaps[filename] = dict(
698 alert_type=alert_type,
699 filename=filename,
700 latency=issued_time-event_time.gps)
702 if combined_skymaps:
703 combined_skymaps = list(combined_skymaps.values())
704 combined_skymap = combined_skymaps[-1]['filename']
705 cls = [50, 90]
706 include_ellipse, ra, dec, a, b, pa, area, greedy_area = \
707 uncertainty_ellipse(gracedb_id, combined_skymap, client,
708 cls=cls)
709 kwargs.update(
710 combined_skymap=combined_skymap,
711 combined_skymaps=combined_skymaps,
712 cl=cls[-1],
713 combined_skymap_include_ellipse=include_ellipse,
714 combined_skymap_ra=coord.Longitude(ra*u.deg),
715 combined_skymap_dec=coord.Latitude(dec*u.deg),
716 combined_skymap_a=coord.Angle(a*u.deg),
717 combined_skymap_b=coord.Angle(b*u.deg),
718 combined_skymap_pa=coord.Angle(pa*u.deg),
719 combined_skymap_ellipse_area=area,
720 combined_skymap_greedy_area=greedy_area,
721 change_significance_statement=change_significance_statement)
723 return kwargs