Commit f24f6503 authored by Sudarshan Ghonge's avatar Sudarshan Ghonge
Browse files

Fix megasky, readbwb. Fix bayeswave_pipe and bayeswave_pipe_utisl to setup megasky accordingly

parent b872b3f2
......@@ -1164,10 +1164,12 @@ class bayeswave_fpeakNode(bayeswave_postNode):
class megaskyJob(pipeline.CondorDAGJob,pipeline.AnalysisJob):
def __init__(self, cp, dax=False):
def __init__(self, cp, inj, dax=False):
# --- Common workflow configuration
if inj is not None:
self.add_opt('inj', inj)
# --- OSG-specifics unique to this job
# Files to include in transfer
......@@ -1190,7 +1192,12 @@ class megaskyNode(pipeline.CondorDAGNode, pipeline.AnalysisNode):
pipeline.CondorDAGNode.__init__(self, megasky_job)
# Set eventnum if injection
def set_injevent(self, eventnum):
print 'eventnum:', eventnum
self.add_var_opt('eventnum', eventnum)
# Set work dir
def set_outputDir(self, outputDir):
......@@ -61,13 +61,15 @@ class BwbParams:
# This is my (Tyson's) horrendously hacky way to find the injection parameters.
# Can someon who knows Python fix this?
injline = lines[26].split(' ')
print injline[2]
if injline[2]=='Signal':
print "Found Signal Injection data!"
self.ra = float(lines[27].split()[1])
sin_dec = float(lines[27].split()[2])
self.dec = np.arcsin(sin_dec)
# Update: Fixed. The script now picks up the injection sky location from the xml file
# Sudarshan Ghonge
#injline = lines[26].split(' ')
#print injline[2]
#if injline[2]=='Signal':
# print "Found Signal Injection data!"
# self.ra = float(lines[27].split()[1])
# sin_dec = float(lines[27].split()[2])
# self.dec = np.arcsin(sin_dec)
......@@ -724,7 +724,7 @@ if opts.fpeak_analysis:
bayeswave_fpeak_job = pipe_utils.bayeswave_fpeakJob(cp, cache_files,
injfile=injfile, nrdata=nrdata)
megasky_job = pipe_utils.megaskyJob(cp)
megasky_job = pipe_utils.megaskyJob(cp, injfile)
megaplot_job = pipe_utils.megaplotJob(cp)
if opts.submit_to_gracedb: submitToGraceDB_job = pipe_utils.submitToGraceDB(cp)
......@@ -917,6 +917,10 @@ for t,trigger in enumerate(trigger_list.triggers):
# --- Add options for mega-scripts
if injfile is not None:
print 'Injfile is not none'
print 'adding event ', trigger.injevent
......@@ -20,14 +20,52 @@ import sky_area.sky_area_clustering as sky
import getopt
from bayeswave_plot.readbwb import BwbParams
import math
from optparse import OptionParser
from glue.ligolw import ligolw
from glue.ligolw import lsctables
from glue.ligolw import table
from glue.ligolw import utils
class LIGOLWContentHandler(ligolw.LIGOLWContentHandler):
#import lalinference.fits as lf
import as lf
import lalinference.plot as lp
import lalinference.plot.cmap
#import lalinference.cmap
def parser():
Parser for input (command line and ini file)
# --- cmd line
parser = OptionParser()
parser.add_option("-M", "--mdc", type=str, default=None)
parser.add_option("-N", "--nside", type=int, default=128)
parser.add_option("-I", "--inj", type=str, default=None, help='Injection xml')
parser.add_option("-e", "--eventnum", type=int, default=None, help='Injection event id')
parser.add_option("-r", "--ra", type=float, default=None)
parser.add_option("-d", "--dec", type=float, default=None)
parser.add_option("--geo", default=False, action="store_true")
(opts,args) = parser.parse_args()
if len(args)==0:
print >> sys.stderr, "ERROR: require BW directory"
if not os.path.isdir(args[0]):
print >> sys.stderr, "ERROR: BW dir %s does not exist"%args[0]
return opts, args,
def gps_to_mjd(gps_time):
Convert a floating-point GPS time in seconds to a modified Julian day.
......@@ -110,11 +148,10 @@ def gps_to_iso8601(gps_time):
# Works with single output directory from
# BayesWaveBurst
# -------------------------------------------
def make_skyview(directory='.', mdc=None, NSIDE=128, ra=None, dec=None, results=None, npost=5000):
def make_skyview(directory='.', mdc=None, NSIDE=128, inj=None, npost=5000):
# -- Close any open figures
# ----------------
# Read S6 MDC log file
# -----------------
......@@ -173,7 +210,6 @@ def make_skyview(directory='.', mdc=None, NSIDE=128, ra=None, dec=None, results=
jobname = params.jobname
print "got job name"
print "params.ra {0}".format(params.ra)
# -- Input skymap data
# num, L, ralist, sin_dec, psi, e, dA, dphi, dt = np.loadtxt(filename, unpack=True)
......@@ -224,49 +260,16 @@ def make_skyview(directory='.', mdc=None, NSIDE=128, ra=None, dec=None, results=
print injtheta, injra
# -- Special handling for injections drawn from prior
if params.ra is not None:
injdec = params.dec
injra = params.ra
injtheta = np.pi/2 - injdec
mdc = True
print "Psych! Got injection values from"
# -- Make plots directory, if needed
plotsDir = './plots'
if not os.path.exists(plotsDir):
# Add markers (e.g., for injections or external triggers).
if (ra is not None and dec is not None):
# Convert the right ascension to either a reference angle from -pi to pi
# or a wrapped angle from 0 to 2 pi, depending on the version of Matplotlib.
#for rarad, decrad in [np.deg2rad([ra, dec])]:
(rarad, decrad) = np.deg2rad([ra, dec])
if geo:
dlon = -sidtime # + or -?
#rarad = lalinference.plot.reference_angle(rarad + dlon)
rarad = rarad + dlon
while rarad > np.pi:
rarad = rarad - 2*np.pi
while rarad < -np.pi:
rarad = rarad + 2*np.pi
#rarad = lalinference.plot.wrapped_angle(rarad)
while rarad > 2*np.pi:
rarad = rarad - 2*np.pi
while rarad < 0:
rarad = rarad + 2*np.pi
# Shift the declination to match hp.projscatter conventions
decrad = np.pi*0.5 - decrad
# -- Plot the skymap and injection location
skymap = kde.as_healpix(NSIDE, nest=True)
#fig = plt.figure(frameon=False)
fig = plt.figure(figsize=(8,6), frameon=False)
ax = plt.subplot(111, projection='astro hours mollweide')
......@@ -275,8 +278,8 @@ def make_skyview(directory='.', mdc=None, NSIDE=128, ra=None, dec=None, results=
# TODO: our .png skymaps don't work anymore.
if (ra is not None and dec is not None):
plt.plot(rarad, dec, 'kx', ms=30, mew=1)
if (inj is not None):
plt.plot(inj['ra'], inj['dec'], 'kx', ms=30, mew=1)
if mdc is not None:
plt.plot(injra, injdec, 'kx', ms=30, mew=1)
......@@ -314,11 +317,6 @@ def make_skyview(directory='.', mdc=None, NSIDE=128, ra=None, dec=None, results=
injcontour = 0
searcharea = 0
pixarea = hp.nside2pixarea(NSIDE, degrees=True)
if results is not None:
for name in ['injcontour', 'area50', 'area90', 'searcharea']:
if name not in results: results[name] = []
# -- Make an output file with key statistics
outfile = open('skystats.txt', 'w')
......@@ -339,57 +337,44 @@ if __name__ == "__main__":
# Working directory is current directory unless it is specified by the
# first argument of the function call
inargs = sys.argv[1:]
if len(sys.argv) >= 2:
if os.path.isdir(sys.argv[1]):
workdir = sys.argv[1]
inargs = sys.argv[2:]
opts, args = parser()
opts, args = getopt.getopt(inargs, "", ['directory=', 'mdc=', 'NSIDE=', 'ra=', 'dec=', 'geo'])
# -- Set default argument values
ra = opts.ra
dec = opts.dec
geo = opts.dec
mdc = opts.mdc
NSIDE = opts.nside
injpos = None
if opts.inj is not None:
if(opts.eventnum is None):
print >> sys.stderr, "Provide event num if giving injfile"
print "Loading xml"
xmldoc = utils.load_filename(
opts.inj, contenthandler=LIGOLWContentHandler)
print('Checking if using a sim_inspiral table...')
injs = table.get_table(
xmldoc, lsctables.SimInspiralTable.tableName)
inj = injs[opts.eventnum]
injpos = {
'ra': inj.longitude,
'dec': inj.latitude,
'id': inj.simulation_id}
print(' yes')
print('Checking if using a sim_burst table...')
injs = table.get_table(xmldoc, lsctables.SimBurstTable.tableName)
inj = injs[opts.eventnum]
injpos = {'ra': inj.ra, 'dec': inj.dec, 'id': inj.simulation_id}
print(' yes')
print "Got these arguments"
print opts
print "With these values:"
print args
for opt, arg in opts:
print opt
if opt=='--mdc':
# -- mdc argument should be the name of the MDC log
mdc = ft.Mdc(arg)
print "Reading MDC log {0}".format(arg)
mdc = ft.Mdc(arg)
print "WARNING! Failed to read MDC log file"
#mdc = None
if opt == '--directory':
directory = arg
if opt == '--NSIDE':
NSIDE = int(arg)
if opt == '--sim':
sim = 'True' == arg
if opt == '--ra':
ra = float(arg)
if opt == '--dec':
dec = float(arg)
if opt == '--geo':
geo = 'True'
if ra is not None and dec is not None:
injpos ={'ra': ra, 'dec':dec, 'id':0}
make_skyview(directory, mdc, NSIDE, ra, dec, results)
make_skyview(args[0], mdc, NSIDE, injpos)
# Move back to original dir
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment