Merge branch 'snr-optimizer-dag-scripts' into 'main'

bayestar scripts

See merge request !50

scripts/snr_optimizer/calculate_fits_and_searched_area.py

+#!/usr/bin/env python3
+
+import copy
+import http.client
+from optparse import OptionParser
+import os
+import sys
+import tempfile
+import numpy as np
+import pickle
+
+from astropy.coordinates import SkyCoord
+from astropy.table import Table
+from astropy.io.registry.base import IORegistryError
+import astropy_healpix as ah
+from astropy import units as u
+from astropy.units.quantity import Quantity as Quantity
+
+from lal import LIGOTimeGPS, GreenwichMeanSiderealTime
+
+from ligo.gracedb.rest import GraceDb
+from ligo.gracedb.exceptions import HTTPError
+
+from ligo.lw import ligolw
+from ligo.lw import lsctables
+from ligo.lw import utils as ligolw_utils
+
+from ligo.skymap.bayestar import localize
+from ligo.skymap.io import events as LIGOSkymapEvents
+from ligo.skymap.io import read_sky_map
+from ligo.skymap.postprocess import crossmatch
+
+class LIGOLWContentHandler(ligolw.LIGOLWContentHandler):
+	pass
+lsctables.use_in(LIGOLWContentHandler)
+
+def parse_command_line():
+	parser = OptionParser()
+	parser.add_option('--input-coinc', help = 'Coinc filenames to process')
+	#parser.add_option('--input-dir', help = 'Coinc filenames to process')
+	parser.add_option('--output-skymap-dir', metavar = 'dir', default = 'fits_files', help = 'Directory to write out skymap fits files to.')
+	#parser.add_option('--output-data-dir', metavar = 'dir', default = 'data', help = 'Directory to write out searched are and searched prob data to.')
+	#parser.add_option('--injfile', metavar = 'url', default = 'injections.xml.gz', help = 'Location of inj file')
+	parser.add_option('--f_low', metavar = 'int', default = 15, help = 'f_low to use')
+	opts, args = parser.parse_args()
+
+	return opts, args
+
+## set up
+opts, args = parse_command_line()
+
+#def skymap_90_area(skymap):
+#	skymap.sort('PROBDENSITY', reverse=True)
+#	level, ipix = ah.uniq_to_level_ipix(skymap['UNIQ'])
+#	pixel_area = ah.nside_to_pixel_area(ah.level_to_nside(level))
+#	prob = pixel_area * skymap['PROBDENSITY']
+#	cumprob = np.cumsum(prob)
+#	#i = cumprob.searchsorted(0.9)
+#	i = cumprob.searchsorted(Quantity('0.9sr'))
+#	area_90 = pixel_area[:i].sum()
+#	return area_90.to_value(u.deg**2)
+
+#def skymap_searched_prob(skymap, this_inj):
+#	ra, dec = this_inj.ra_dec
+#	loc = SkyCoord(ra=ra, dec=dec, unit="rad")
+#	p = crossmatch(skymap, loc).searched_prob
+#	deg2 = crossmatch(skymap, loc).searched_area
+#	return p, deg2
+
+#inj_xmldoc = ligolw_utils.load_filename(opts.injfile, contenthandler = LIGOLWContentHandler)
+#inj_simtable = lsctables.SimInspiralTable.get_table(inj_xmldoc)
+
+event_source = LIGOSkymapEvents.ligolw.open(f"{opts.input_coinc}", psd_file=f"{opts.input_coinc}", coinc_def=None)
+for event_id, event in event_source.items():
+	skymap = localize(event, f_low = opts.f_low)
+	break # only one event per file
+skymap.write(f"{opts.output_skymap_dir}/{opts.input_coinc.split('/')[-1]}", format="fits", overwrite=True)
+
+#end_time = list(event_source.values())[0].singles[0].time
+#this_inj = None
+#for row in inj_simtable:
+#	this_end_time = row.geocent_end_time
+#	if this_end_time - 1 < end_time and this_end_time + 1 > end_time:
+#		this_inj = row
+#		break
+#if this_inj is None:
+#	raise ValueError(f'cannot find inj for {opts.input_coinc}')
+
+#searched_prob, searched_area = skymap_searched_prob(skymap, this_inj)
+#skymap_area = skymap_90_area(skymap)
+
+#with open(f"{opts.output_data_dir}/{opts.input_coinc}", 'wb') as f:
+#	pickle.dump((searched_prob, searched_area, skymap_area), f)

scripts/snr_optimizer/calculate_fits_and_searched_area_wrapper.sh

+#!/bin/bash
+
+source /cvmfs/oasis.opensciencegrid.org/ligo/sw/conda/etc/profile.d/conda.sh
+conda activate igwn
+./calculate_fits_and_searched_area.py "$@"

scripts/snr_optimizer/skymap.sub

+universe = vanilla
+executable = calculate_fits_and_searched_area_wrapper.sh
+arguments = --input-coinc $(input_coinc)
+periodic_release = (HoldReasonCode == 5) || ((CurrentTime - EnteredCurrentStatus > 180) && (HoldReasonCode != 34))
+request_cpus = 1
+request_memory = 6000
+request_disk = 3GB
+want_graceful_removal = True
+kill_sig = 15
+accounting_group = ligo.dev.o4.cbc.em.gstlalonline
+accounting_group_user = prathamesh.joshi
+environment = "LAL_DATA_PATH=/cvmfs/oasis.opensciencegrid.org/ligo/sw/pycbc/lalsuite-extra/current/share/lalsimulation"
+transfer_executable = False
+getenv = False
+requirements = has_avx2 && (CpuFamily =?= 6) && (HAS_SINGULARITY=?=True)
+success_exit_code = 0
+output = logs/$(nodename)-$(cluster)-$(process).out
+error = logs/$(nodename)-$(cluster)-$(process).err
+notification = never
+MY.MemoryUsage = ( 4000 ) * 2 / 3
+
+queue