Coverage for gwcelery/tasks/external_skymaps.py: 100%
257 statements
« prev ^ index » next coverage.py v7.4.4, created at 2025-01-17 06:48 +0000
1"""Create and upload external sky maps."""
2import io
3import re
4import ssl
5import urllib
6from base64 import b64decode
7from urllib.error import HTTPError
9import astropy_healpix as ah
10import gcn
11import lxml.etree
12# import astropy.utils.data
13import numpy as np
14from astropy import units as u
15from astropy.coordinates import SkyCoord
16from astropy.io.fits import getheader
17from celery import group
18from hpmoc import PartialUniqSkymap
19from hpmoc.utils import reraster, uniq_intersection
20from ligo.skymap.distance import parameters_to_marginal_moments
21from ligo.skymap.io import fits
22from ligo.skymap.moc import bayestar_adaptive_grid
23from ligo.skymap.plot.bayes_factor import plot_bayes_factor
25from .. import _version, app
26from ..util import closing_figures
27from ..util.cmdline import handling_system_exit
28from ..util.tempfile import NamedTemporaryFile
29from . import gracedb, skymaps
31COMBINED_SKYMAP_FILENAME_MULTIORDER = 'combined-ext.multiorder.fits'
32"""Filename of combined sky map in a multiordered format"""
34COMBINED_SKYMAP_FILENAME_PNG = 'combined-ext.png'
35"""Filename of combined sky map plot"""
37FERMI_OFFICIAL_SKYMAP_FILENAME = 'glg_healpix_all_bn_v00'
38"""Filename of sky map from official Fermi GBM analysis"""
40NOTICE_TYPE_DICT = {
41 '53': 'INTEGRAL_WAKEUP',
42 '54': 'INTEGRAL_REFINED',
43 '55': 'INTEGRAL_OFFLINE',
44 '60': 'SWIFT_BAT_GRB_ALERT',
45 '61': 'SWIFT_BAT_GRB_POSITION',
46 '110': 'FERMI_GBM_ALERT',
47 '111': 'FERMI_GBM_FLT_POS',
48 '112': 'FERMI_GBM_GND_POS',
49 '115': 'FERMI_GBM_FINAL_POS',
50 '131': 'FERMI_GBM_SUBTHRESHOLD'
51}
54@app.task(shared=False)
55def create_combined_skymap(se_id, ext_id, preferred_event=None):
56 """Creates and uploads a combined LVK-external skymap, uploading to the
57 external trigger GraceDB page. The filename used for the combined sky map
58 will be 'combined-ext.multiorder.fits' if the external sky map is
59 multi-ordered or 'combined-ext.fits.gz' if the external sky map is not.
60 Will use the GW sky map in from the preferred event if given.
62 Parameters
63 ----------
64 se_id : str
65 Superevent GraceDB ID
66 ext_id : str
67 External event GraceDB ID
68 preferred_event : str
69 Preferred event GraceDB ID. If given, use sky map from preferred event
71 """
72 # gw_id is either superevent or preferred event graceid
73 gw_id = preferred_event if preferred_event else se_id
74 gw_skymap_filename = get_skymap_filename(gw_id, is_gw=True)
75 ext_skymap_filename = get_skymap_filename(ext_id, is_gw=False)
76 # Determine whether external sky map is multiordered or flat
77 ext_moc = '.multiorder.fits' in ext_skymap_filename
79 message = \
80 ('Combined LVK-external sky map using {0} and {1}, with {2} and '
81 '{3}'.format(gw_id, ext_id, gw_skymap_filename, ext_skymap_filename))
82 message_png = (
83 'Mollweide projection of <a href="/api/events/{graceid}/files/'
84 '{filename}">{filename}</a>').format(
85 graceid=ext_id,
86 filename=COMBINED_SKYMAP_FILENAME_MULTIORDER)
88 (
89 _download_skymaps.si(
90 gw_skymap_filename, ext_skymap_filename, gw_id, ext_id
91 )
92 |
93 combine_skymaps.s(ext_moc=ext_moc)
94 |
95 group(
96 gracedb.upload.s(
97 COMBINED_SKYMAP_FILENAME_MULTIORDER,
98 ext_id, message, ['sky_loc', 'ext_coinc']
99 ),
101 skymaps.plot_allsky.s()
102 |
103 gracedb.upload.s(COMBINED_SKYMAP_FILENAME_PNG, ext_id,
104 message_png, ['sky_loc', 'ext_coinc'])
105 )
106 |
107 gracedb.create_label.si('COMBINEDSKYMAP_READY', ext_id)
108 ).delay()
111@app.task(autoretry_for=(ValueError,), retry_backoff=10,
112 retry_backoff_max=600)
113def get_skymap_filename(graceid, is_gw):
114 """Get the skymap FITS filename.
116 If not available, will try again 10 seconds later, then 20, then 40, etc.
117 up to a max 10 minutes retry delay.
119 Parameters
120 ----------
121 graceid : str
122 GraceDB ID
123 is_gw : bool
124 If True, uses method for superevent or preferred event. Otherwise uses
125 method for external event.
127 Returns
128 -------
129 filename : str
130 Filename of latest sky map
132 """
133 gracedb_log = gracedb.get_log(graceid)
134 if is_gw:
135 # Try first to get a multiordered sky map
136 for message in reversed(gracedb_log):
137 filename = message['filename']
138 v = message['file_version']
139 fv = '{},{}'.format(filename, v)
140 if filename.endswith('.multiorder.fits') and \
141 "combined-ext." not in filename:
142 return fv
143 # Try next to get a flattened sky map
144 for message in reversed(gracedb_log):
145 filename = message['filename']
146 v = message['file_version']
147 fv = '{},{}'.format(filename, v)
148 if filename.endswith('.fits.gz') and \
149 "combined-ext." not in filename:
150 return fv
151 else:
152 for message in reversed(gracedb_log):
153 filename = message['filename']
154 v = message['file_version']
155 fv = '{},{}'.format(filename, v)
156 if (filename.endswith('.fits') or filename.endswith('.fit') or
157 filename.endswith('.fits.gz')) and \
158 "combined-ext." not in filename:
159 return fv
160 raise ValueError('No skymap available for {0} yet.'.format(graceid))
163def is_skymap_moc(skymap_bytes):
164 with NamedTemporaryFile(content=skymap_bytes) as skymap_file:
165 # If ordering is EXPLICIT, should be multiordered (MOC)
166 return getheader(skymap_file.name, ext=1)['INDXSCHM'] == 'EXPLICIT'
169@app.task(shared=False)
170def _download_skymaps(gw_filename, ext_filename, gw_id, ext_id):
171 """Download both superevent and external sky map to be combined.
173 Parameters
174 ----------
175 gw_filename : str
176 GW sky map filename
177 ext_filename : str
178 External sky map filename
179 gw_id : str
180 GraceDB ID of GW candidate, either superevent or preferred event
181 ext_id : str
182 GraceDB ID of external candidate
184 Returns
185 -------
186 gw_skymap, ext_skymap : tuple
187 Tuple of gw_skymap and ext_skymap bytes
189 """
190 gw_skymap = gracedb.download(gw_filename, gw_id)
191 ext_skymap = gracedb.download(ext_filename, ext_id)
192 return gw_skymap, ext_skymap
195def combine_skymaps_moc_moc(gw_sky, ext_sky):
196 """This function combines a multi-ordered (MOC) GW sky map with a MOC
197 external skymap.
198 """
199 gw_sky_hpmoc = PartialUniqSkymap(gw_sky["PROBDENSITY"], gw_sky["UNIQ"],
200 name="PROBDENSITY", meta=gw_sky.meta)
201 # Determine the column name in ext_sky and rename it as PROBDENSITY.
202 ext_sky_hpmoc = PartialUniqSkymap(ext_sky["PROBDENSITY"], ext_sky["UNIQ"],
203 name="PROBDENSITY", meta=ext_sky.meta)
205 comb_sky_hpmoc = gw_sky_hpmoc * ext_sky_hpmoc
206 comb_sky_hpmoc /= np.sum(comb_sky_hpmoc.s * comb_sky_hpmoc.area())
207 comb_sky = comb_sky_hpmoc.to_table(name='PROBDENSITY')
209 # Modify GW sky map with new data, ensuring they exist first
210 if 'DISTMU' in gw_sky.keys() and 'DISTSIGMA' in gw_sky.keys():
211 UNIQ = comb_sky['UNIQ']
212 UNIQ_ORIG = gw_sky['UNIQ']
213 intersection = uniq_intersection(UNIQ_ORIG, UNIQ)
214 DIST_MU = reraster(UNIQ_ORIG,
215 gw_sky["DISTMU"],
216 UNIQ,
217 method='copy',
218 intersection=intersection)
219 DIST_SIGMA = reraster(UNIQ_ORIG,
220 gw_sky["DISTSIGMA"],
221 UNIQ,
222 method='copy',
223 intersection=intersection)
224 DIST_NORM = reraster(UNIQ_ORIG,
225 gw_sky["DISTNORM"],
226 UNIQ,
227 method='copy',
228 intersection=intersection)
229 comb_sky.add_columns([DIST_MU, DIST_SIGMA, DIST_NORM],
230 names=['DISTMU', 'DISTSIGMA', 'DISTNORM'])
232 distmean, diststd = parameters_to_marginal_moments(
233 comb_sky['PROBDENSITY'] * comb_sky_hpmoc.area().value,
234 comb_sky['DISTMU'], comb_sky['DISTSIGMA'])
235 comb_sky.meta['distmean'], comb_sky.meta['diststd'] = distmean, diststd
236 if 'instruments' not in ext_sky.meta:
237 ext_sky.meta.update({'instruments': {'external instrument'}})
238 if 'instruments' in comb_sky.meta:
239 comb_sky.meta['instruments'].update(ext_sky.meta['instruments'])
240 if 'HISTORY' in comb_sky.meta:
241 ext_instrument = list(ext_sky.meta['instruments'])[0]
242 comb_sky.meta['HISTORY'].extend([
243 '', 'The values were reweighted by using data from {0}{1}'.format(
244 ('an ' if ext_instrument == 'external instrument'
245 else ''),
246 ext_instrument)])
248 # Remove redundant field
249 if 'ORDERING' in comb_sky.meta:
250 del comb_sky.meta['ORDERING']
251 return comb_sky
254def combine_skymaps_moc_flat(gw_sky, ext_sky, ext_header):
255 """This function combines a multi-ordered (MOC) GW sky map with a flattened
256 external one by re-weighting the MOC sky map using the values of the
257 flattened one.
259 Header info is generally inherited from the GW sky map or recalculated
260 using the combined sky map values.
262 Parameters
263 ----------
264 gw_sky : Table
265 GW sky map astropy Table
266 ext_sky : array
267 External sky map array
268 ext_header : dict
269 Header of external sky map
271 Returns
272 -------
273 comb_sky : Table
274 Table of combined sky map
276 """
277 # Find ra/dec of each GW pixel
278 level, ipix = ah.uniq_to_level_ipix(gw_sky["UNIQ"])
279 nsides = ah.level_to_nside(level)
280 areas = ah.nside_to_pixel_area(nsides)
281 ra_gw, dec_gw = ah.healpix_to_lonlat(ipix, nsides, order='nested')
282 # Find corresponding external sky map indicies
283 ext_nside = ah.npix_to_nside(len(ext_sky))
284 ext_ind = ah.lonlat_to_healpix(
285 ra_gw, dec_gw, ext_nside,
286 order='nested' if ext_header['nest'] else 'ring')
287 # Fix any negative values and normalize
288 ext_sky[ext_sky < 0.] = 0.
289 ext_sky /= np.sum(ext_sky)
290 # Reweight GW prob density by external sky map probabilities
291 gw_sky['PROBDENSITY'] *= ext_sky[ext_ind]
292 gw_sky['PROBDENSITY'] /= \
293 np.sum(gw_sky['PROBDENSITY'] * areas).value
294 # Modify GW sky map with new data, ensuring they exist first
295 if 'DISTMU' in gw_sky.keys() and 'DISTSIGMA' in gw_sky.keys():
296 distmean, diststd = parameters_to_marginal_moments(
297 gw_sky['PROBDENSITY'] * areas.value,
298 gw_sky['DISTMU'], gw_sky['DISTSIGMA'])
299 gw_sky.meta['distmean'], gw_sky.meta['diststd'] = distmean, diststd
300 if 'instruments' not in ext_header:
301 ext_header.update({'instruments': {'external instrument'}})
302 if 'instruments' in gw_sky.meta:
303 gw_sky.meta['instruments'].update(ext_header['instruments'])
304 if 'HISTORY' in gw_sky.meta:
305 ext_instrument = list(ext_header['instruments'])[0]
306 gw_sky.meta['HISTORY'].extend([
307 '', 'The values were reweighted by using data from {0}{1}'.format(
308 ('an ' if ext_instrument == 'external instrument'
309 else ''),
310 ext_instrument)])
311 return gw_sky
314@app.task(shared=False)
315def combine_skymaps(skymapsbytes, ext_moc=True):
316 """This task combines the two input sky maps, in this case the external
317 trigger skymap and the LVK skymap and writes to a temporary output file. It
318 then returns the contents of the file as a byte array.
320 There are separate methods in case the GW sky map is multiordered (we just
321 reweight using the external sky map) or flattened (use standard
322 ligo.skymap.combine method).
324 Parameters
325 ----------
326 skymapbytes : tuple
327 Tuple of gw_skymap and ext_skymap bytes
328 gw_moc : bool
329 If True, assumes the GW sky map is a multi-ordered format
331 Returns
332 -------
333 combinedskymap : bytes
334 Bytes of combined sky map
335 """
336 gw_skymap_bytes, ext_skymap_bytes = skymapsbytes
337 with NamedTemporaryFile(mode='rb', suffix=".fits") as combinedskymap, \
338 NamedTemporaryFile(content=gw_skymap_bytes) as gw_skymap_file, \
339 NamedTemporaryFile(content=ext_skymap_bytes) as ext_skymap_file, \
340 handling_system_exit():
342 gw_skymap = fits.read_sky_map(gw_skymap_file.name, moc=True)
343 # If GW sky map is multiordered, use reweighting method
344 if ext_moc:
345 # Load external sky map
346 ext_skymap = fits.read_sky_map(ext_skymap_file.name, moc=True)
347 # Create and write combined sky map
348 combined_skymap = combine_skymaps_moc_moc(gw_skymap,
349 ext_skymap)
350 # If GW sky map is flattened, use older method
351 else:
352 # Load external sky map
353 ext_skymap, ext_header = fits.read_sky_map(ext_skymap_file.name,
354 moc=False, nest=True)
355 # Create and write combined sky map
356 combined_skymap = combine_skymaps_moc_flat(gw_skymap, ext_skymap,
357 ext_header)
358 fits.write_sky_map(combinedskymap.name, combined_skymap, moc=True)
359 return combinedskymap.read()
362@app.task(shared=False)
363def external_trigger_heasarc(external_id):
364 """Returns the HEASARC FITS file link.
366 Parameters
367 ----------
368 external_id : str
369 GraceDB ID of external event
371 Returns
372 -------
373 heasarc_link : str
374 Guessed HEASARC URL link
376 """
377 gracedb_log = gracedb.get_log(external_id)
378 for message in gracedb_log:
379 if 'Original Data' in message['comment']:
380 filename = message['filename']
381 xmlfile = gracedb.download(urllib.parse.quote(filename),
382 external_id)
383 root = lxml.etree.fromstring(xmlfile)
384 heasarc_url = root.find('./What/Param[@name="LightCurve_URL"]'
385 ).attrib['value']
386 return re.sub(r'quicklook(.*)', 'current/', heasarc_url)
387 raise ValueError('Not able to retrieve HEASARC link for {0}.'.format(
388 external_id))
391@app.task(autoretry_for=(gracedb.RetryableHTTPError,), retry_backoff=30,
392 max_retries=8)
393def get_external_skymap(link, search):
394 """Download the Fermi sky map FITS file and return the contents as a byte
395 array. If GRB, will construct a HEASARC url, while if SubGRB, will use the
396 link directly.
398 If not available, will try again 10 seconds later, then 20, then 40, etc.
399 until up to 15 minutes after initial attempt.
401 Parameters
402 ----------
403 link : str
404 HEASARC URL link
405 search : str
406 Search field of external event
408 Returns
409 -------
410 external_skymap : bytes
411 Bytes of external sky map
413 """
414 if search == 'GRB':
415 # if Fermi GRB, determine final HEASARC link
416 trigger_id = re.sub(r'.*\/(\D+?)(\d+)(\D+)\/.*', r'\2', link)
417 skymap_name = 'glg_healpix_all_bn{0}_v00.fit'.format(trigger_id)
418 skymap_link = link + skymap_name
419 elif search in {'SubGRB', 'FromURL'}:
420 skymap_link = link
421 # FIXME: Under Anaconda on the LIGO Caltech computing cluster, Python
422 # (and curl, for that matter) fail to negotiate TLSv1.2 with
423 # heasarc.gsfc.nasa.gov
424 context = ssl.create_default_context()
425 context.options |= ssl.OP_NO_TLSv1_3
426 # return astropy.utils.data.get_file_contents(
427 # (skymap_link), encoding='binary', cache=False)
428 try:
429 response = urllib.request.urlopen(skymap_link, context=context)
430 return response.read()
431 except HTTPError as e:
432 if e.code == 404:
433 raise gracedb.RetryableHTTPError("Failed to download the sky map."
434 "Retrying...")
435 else:
436 raise
439@app.task(shared=False)
440def read_upload_skymap_from_base64(event, skymap_str):
441 """Decode and upload 64base encoded sky maps from Kafka alerts.
443 Parameters
444 ----------
445 event : dict
446 External event dictionary
447 skymap_str : str
448 Base 64 encoded sky map string
450 """
452 graceid = event['graceid']
454 # Decode base64 encoded string to bytes string
455 skymap_data = b64decode(skymap_str)
457 # Determine filename based on whether is multiordered or flattened
458 skymap_filename = event['pipeline'].lower() + '_skymap.'
459 skymap_filename += ('multiorder.fits' if is_skymap_moc(skymap_data)
460 else 'fits.gz')
462 message = (
463 'Mollweide projection of <a href="/api/events/{graceid}/files/'
464 '{filename}">{filename}</a>').format(
465 graceid=graceid, filename=skymap_filename)
467 (
468 group(
469 gracedb.upload.si(
470 skymap_data,
471 skymap_filename,
472 graceid,
473 'Sky map uploaded from {} via a kafka notice'.format(
474 event['pipeline']),
475 ['sky_loc']),
477 skymaps.plot_allsky.si(skymap_data)
478 |
479 gracedb.upload.s(event['pipeline'].lower() + '_skymap.png',
480 graceid,
481 message,
482 ['sky_loc'])
483 )
484 |
485 gracedb.create_label.si('EXT_SKYMAP_READY', graceid)
486 ).delay()
489@app.task(autoretry_for=(urllib.error.HTTPError, urllib.error.URLError,),
490 retry_backoff=10, retry_backoff_max=1200)
491def get_upload_external_skymap(event, skymap_link=None):
492 """If a Fermi sky map is not uploaded yet, tries to download one and upload
493 to external event. If sky map is not available, passes so that this can be
494 re-run the next time an update GCN notice is received. If GRB, will
495 construct a HEASARC url, while if SubGRB, will use the link directly.
496 If SubGRB or FromURL, downloads a skymap using the provided URL rather
497 than construct one.
499 Parameters
500 ----------
501 event : dict
502 External event dictionary
503 skymap_link : str
504 HEASARC URL link
506 """
507 graceid = event['graceid']
508 search = event['search']
510 if search == 'GRB':
511 skymap_data = \
512 get_external_skymap(external_trigger_heasarc(graceid), search)
513 elif search in {'SubGRB', 'SubGRBTargeted', 'FromURL'}:
514 skymap_data = get_external_skymap(skymap_link, search)
516 skymap_filename_base = \
517 ('external_from_url' if search == 'FromURL'
518 else FERMI_OFFICIAL_SKYMAP_FILENAME)
519 skymap_filename = skymap_filename_base + \
520 ('.multiorder.fits' if is_skymap_moc(skymap_data) else '.fits.gz')
522 fits_message = \
523 ('Downloaded from {}.'.format(skymap_link) if search == 'FromURL'
524 else 'Official sky map from Fermi analysis.')
525 png_message = (
526 'Mollweide projection of <a href="/api/events/{graceid}/files/'
527 '{filename}">{filename}</a>').format(
528 graceid=graceid, filename=skymap_filename)
530 (
531 group(
532 gracedb.upload.si(
533 skymap_data,
534 skymap_filename,
535 graceid,
536 fits_message,
537 ['sky_loc']),
539 skymaps.plot_allsky.si(skymap_data)
540 |
541 gracedb.upload.s(skymap_filename_base + '.png',
542 graceid,
543 png_message,
544 ['sky_loc'])
545 )
546 |
547 gracedb.create_label.si('EXT_SKYMAP_READY', graceid)
548 ).delay()
551def from_cone(pts, ra, dec, error):
552 """
553 Based on the given RA, DEC, and error radius of the center points,
554 it calculates the gaussian pdf.
555 """
556 ras, decs = pts.T
557 center = SkyCoord(ra * u.deg, dec * u.deg)
558 error_radius = error * u.deg
560 pts_loc = SkyCoord(np.rad2deg(ras) * u.deg, np.rad2deg(decs) * u.deg)
562 distance = pts_loc.separation(center)
563 skymap = np.exp(-0.5 * np.square(distance / error_radius).to_value(
564 u.dimensionless_unscaled))
565 return skymap
568def _fisher(distance, error_stat, error_sys):
569 """
570 Calculates the Fisher distribution from Eq. 1 and Eq. 2 of Connaughton
571 et al. 2015 (doi: 10.1088/0067-0049/216/2/32).
572 """
574 error_tot2 = (
575 np.square(np.radians(error_stat)) +
576 np.square(np.radians(error_sys))
577 )
578 kappa = 1 / (0.4356 * error_tot2)
579 return kappa / (2 * np.pi * (np.exp(kappa) - np.exp(-kappa))) \
580 * np.exp(kappa * np.cos(distance))
583def fermi_error_model(pts, ra, dec, error, core, tail, core_weight):
584 """
585 Calculate the Fermi GBM error model from Connaughton et al. 2015 based
586 on the given RA, Dec and error radius of the center point using the model
587 parameters of core radii, tail radii, and core proportion (f in the paper).
588 """
589 ras, decs = pts.T
590 center = SkyCoord(ra * u.deg, dec * u.deg)
592 pts_loc = SkyCoord(np.rad2deg(ras) * u.deg, np.rad2deg(decs) * u.deg)
593 distance = pts_loc.separation(center).rad # Ensure the output is in radian
595 core_component = core_weight * _fisher(distance, error, core)
596 tail_component = (1 - core_weight) * _fisher(distance, error, tail)
598 return core_component + tail_component
601def create_external_skymap(ra, dec, error, pipeline, notice_type=111):
602 """Create a sky map, either a gaussian or a single
603 pixel sky map, given an RA, dec, and error radius.
605 If from Fermi, convolves the sky map with both a core and
606 tail Gaussian and then sums these to account for systematic
607 effects as measured in :doi:`10.1088/0067-0049/216/2/32`
609 If from Swift, converts the error radius from that containing 90% of the
610 credible region to ~68% (see description of Swift error
611 here:`https://gcn.gsfc.nasa.gov/swift.html#tc7`)
613 Parameters
614 ----------
615 ra : float
616 Right ascension in deg
617 dec : float
618 Declination in deg
619 error : float
620 Error radius in deg
621 pipeline : str
622 External trigger pipeline name
623 notice_type : int
624 GCN notice type integer
626 Returns
627 -------
628 skymap : array
629 Sky map array
631 """
632 # Dictionary definitions for core_weight, core, and tail values
633 # for different notice types.
634 # Flight notice: Values from first row of Table 7
635 # Ground notice: Values from first row of Table 3
636 # Final notice: Values from second row of Table 3
637 fermi_params = {
638 gcn.NoticeType.FERMI_GBM_FLT_POS: {"core_weight": 0.897,
639 "core_width": 7.52,
640 "tail_width": 55.6},
641 gcn.NoticeType.FERMI_GBM_GND_POS: {"core_weight": 0.804,
642 "core_width": 3.72,
643 "tail_width": 13.7},
644 gcn.NoticeType.FERMI_GBM_FIN_POS: {"core_weight": 0.900,
645 "core_width": 3.71,
646 "tail_width": 14.3},
647 }
649 # Correct 90% containment to 1-sigma for Swift
650 if pipeline == 'Swift':
651 error /= np.sqrt(-2 * np.log1p(-.9))
652 # Set minimum error radius so function does not return void
653 # FIXME: Lower this when fixes are made to ligo.skymap to fix nans when
654 # the error radius is too low
655 error = max(error, .08)
656 # This function adaptively refines the grid based on the given gaussian
657 # pdf to create the multi-ordered skymap.
658 if pipeline == 'Fermi' and notice_type is not None:
659 # Correct for Fermi systematics based on recommendations from GBM team
660 # Convolve with both a narrow core and wide tail Gaussian with error
661 # radius determined by the scales respectively, each comprising a
662 # fraction determined by the weights respectively.
663 if notice_type not in fermi_params:
664 raise AssertionError('Provide a supported Fermi notice type')
665 core_weight = fermi_params[notice_type]["core_weight"]
666 # Note that tail weight = 1 - core_weight
667 core_width = fermi_params[notice_type]["core_width"]
668 tail_width = fermi_params[notice_type]["tail_width"]
670 # Integrate the fermi_error_model using bayestar_adaptive_grid
671 skymap = bayestar_adaptive_grid(fermi_error_model, ra, dec, error,
672 core_width, tail_width, core_weight,
673 rounds=8)
674 else:
675 # Use generic cone method for Swift, INTEGRAL, etc.
676 skymap = bayestar_adaptive_grid(from_cone, ra, dec, error, rounds=8)
678 return skymap
681def write_to_fits(skymap, event, notice_type, notice_date):
682 """Write external sky map FITS file, populating the
683 header with relevant info.
685 Parameters
686 ----------
687 skymap : array
688 Sky map array
689 event : dict
690 Dictionary of external event
691 notice_type : int
692 GCN notice type integer
693 notice_date : str
694 External event trigger time in ISO format
696 Returns
697 -------
698 skymap_fits : str
699 Bytes string of sky map
701 """
703 if notice_type is None:
704 msgtype = event['pipeline'] + '_LVK_TARGETED_SEARCH'
705 else:
706 msgtype = NOTICE_TYPE_DICT[str(notice_type)]
708 gcn_id = event['extra_attributes']['GRB']['trigger_id']
709 with NamedTemporaryFile(suffix='.fits') as f:
710 fits.write_sky_map(f.name, skymap,
711 objid=gcn_id,
712 url=event['links']['self'],
713 instruments=event['pipeline'],
714 gps_time=event['gpstime'],
715 msgtype=msgtype,
716 msgdate=notice_date,
717 creator='gwcelery',
718 origin='LIGO-VIRGO-KAGRA',
719 vcs_version=_version.get_versions()['version'],
720 history='file only for internal use')
721 with open(f.name, 'rb') as file:
722 return file.read()
725@app.task(shared=False)
726def create_upload_external_skymap(event, notice_type, notice_date):
727 """Create and upload external sky map using
728 RA, dec, and error radius information.
730 Parameters
731 ----------
732 event : dict
733 Dictionary of external event
734 notice_type : int
735 GCN notice type integer
736 notice_date : str
737 External event trigger time in ISO format
739 """
740 graceid = event['graceid']
741 skymap_filename = event['pipeline'].lower() + '_skymap.multiorder.fits'
743 ra = event['extra_attributes']['GRB']['ra']
744 dec = event['extra_attributes']['GRB']['dec']
745 error = event['extra_attributes']['GRB']['error_radius']
746 pipeline = event['pipeline']
748 if not (ra or dec or error):
749 # Don't create sky map if notice only contains zeros, lacking info
750 return
751 skymap = create_external_skymap(ra, dec, error, pipeline, notice_type)
753 skymap_data = write_to_fits(skymap, event, notice_type, notice_date)
755 if notice_type is None:
756 extra_sentence = ' from {} via our joint targeted search'.format(
757 pipeline)
758 else:
759 msgtype = NOTICE_TYPE_DICT[str(notice_type)]
760 extra_sentence = ' from a {} type GCN notice'.format(msgtype)
762 message = (
763 'Mollweide projection of <a href="/api/events/{graceid}/files/'
764 '{filename}">{filename}</a>{extra_sentence}').format(
765 graceid=graceid, filename=skymap_filename,
766 extra_sentence=extra_sentence)
768 (
769 gracedb.upload.si(
770 skymap_data,
771 skymap_filename,
772 graceid,
773 'Sky map created from GCN RA, dec, and error{}.'.format(
774 extra_sentence),
775 ['sky_loc'])
776 |
777 skymaps.plot_allsky.si(skymap_data)
778 |
779 gracedb.upload.s(event['pipeline'].lower() + '_skymap.png',
780 graceid,
781 message,
782 ['sky_loc'])
783 |
784 gracedb.create_label.si('EXT_SKYMAP_READY', graceid)
785 ).delay()
788@app.task(shared=False)
789@closing_figures()
790def plot_overlap_integral(coinc_far_dict, superevent, ext_event,
791 var_label=r"\mathcal{I}_{\Omega}"):
792 """Plot and upload visualization of the sky map overlap integral computed
793 by ligo.search.overlap_integral.
795 Parameters
796 ----------
797 coinc_far_dict : dict
798 Dictionary containing coincidence false alarm rate results from
799 RAVEN
800 superevent : dict
801 Superevent dictionary
802 ext_event : dict
803 External event dictionary
804 var_label : str
805 The variable symbol used in plotting
807 """
808 if coinc_far_dict.get('skymap_overlap') is None:
809 return
810 if superevent['em_type'] != ext_event['graceid'] and \
811 'RAVEN_ALERT' in superevent['labels']:
812 return
814 superevent_id = superevent['superevent_id']
815 ext_id = ext_event['graceid']
817 log_overlap = np.log(coinc_far_dict['skymap_overlap'])
818 logI_string = np.format_float_positional(log_overlap, 1, trim='0',
819 sign=True)
820 # Create plot
821 fig, _ = plot_bayes_factor(
822 log_overlap, values=(1, 3, 5), xlim=7, var_label=var_label,
823 title=(r'Sky Map Overlap between %s and %s [$\ln\,%s = %s$]' %
824 (superevent_id, ext_id, var_label, logI_string)))
825 # Convert to bytes
826 outfile = io.BytesIO()
827 fig.savefig(outfile, format='png')
828 # Upload file
829 gracedb.upload.si(
830 outfile.getvalue(),
831 'overlap_integral.png',
832 superevent_id,
833 message='Sky map overlap integral between {0} and {1}'.format(
834 superevent_id, ext_id),
835 tags=['ext_coinc']
836 ).delay()