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