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