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

Add mdc fetching code that was previously deleted, add lines which plot and...

Add mdc fetching code that was previously deleted, add lines which plot and calculate searched sky area for injections, remove pylal functions
parent 5af9ee2e
...@@ -7,15 +7,11 @@ import matplotlib.pyplot as plt ...@@ -7,15 +7,11 @@ import matplotlib.pyplot as plt
import glob import glob
import os import os
import sys import sys
from pylal.xlal.datatypes.ligotimegps import LIGOTimeGPS
#from pylal.xlal.date import XLALGreenwichSiderealTime
from lal import GreenwichSiderealTime from lal import GreenwichSiderealTime
from bayeswave_plot import find_trig as ft from bayeswave_plot import find_trig as ft
import healpy as hp import healpy as hp
from subprocess import call
# from scipy.optimize import curve_fit # from scipy.optimize import curve_fit
import acor #import acor
from pylal import bayespputils
import sky_area.sky_area_clustering as sky import sky_area.sky_area_clustering as sky
import getopt import getopt
from bayeswave_plot.readbwb import BwbParams from bayeswave_plot.readbwb import BwbParams
...@@ -164,7 +160,7 @@ def make_skyview(directory='.', mdc=None, NSIDE=128, inj=None, npost=5000): ...@@ -164,7 +160,7 @@ def make_skyview(directory='.', mdc=None, NSIDE=128, inj=None, npost=5000):
# -- Enter requested directory # -- Enter requested directory
os.chdir(directory) os.chdir(directory)
print "Entered", os.getcwd() print "Entered", os.getcwd()
jobname = glob.glob('*bayeswave.run')[0] #jobname = glob.glob('*bayeswave.run')[0]
try: try:
jobname = 'bayeswave.run' jobname = 'bayeswave.run'
...@@ -204,17 +200,11 @@ def make_skyview(directory='.', mdc=None, NSIDE=128, inj=None, npost=5000): ...@@ -204,17 +200,11 @@ def make_skyview(directory='.', mdc=None, NSIDE=128, inj=None, npost=5000):
# -------------------------------- # --------------------------------
# Create skymap summary statistics # Create skymap summary statistics
# -------------------------------- # --------------------------------
# -- Get run name
params = BwbParams()
jobname = params.jobname
print "got job name"
# -- Input skymap data # -- Input skymap data
# num, L, ralist, sin_dec, psi, e, dA, dphi, dt = np.loadtxt(filename, unpack=True) # num, L, ralist, sin_dec, psi, e, dA, dphi, dt = np.loadtxt(filename, unpack=True)
print "Extracting RA/DEC samples" print "Extracting RA/DEC samples"
filename = './chains/' + jobname + 'signal_params.dat.0' filename = './chains/' + 'signal_params.dat.0'
data = np.loadtxt(filename, unpack=True,usecols=(0,1,2)) data = np.loadtxt(filename, unpack=True,usecols=(0,1,2))
ralist = data[1] ralist = data[1]
sin_dec = data[2] sin_dec = data[2]
...@@ -240,9 +230,9 @@ def make_skyview(directory='.', mdc=None, NSIDE=128, inj=None, npost=5000): ...@@ -240,9 +230,9 @@ def make_skyview(directory='.', mdc=None, NSIDE=128, inj=None, npost=5000):
print "Finding sidereal time" print "Finding sidereal time"
print "Using GPS {0}".format(gps) print "Using GPS {0}".format(gps)
trigtime = float(gps) trigtime = float(gps)
lalgps = LIGOTimeGPS(trigtime) #lalgps = LIGOTimeGPS(trigtime)
#sidtime = XLALGreenwichSiderealTime(lalgps, 0) #sidtime = XLALGreenwichSiderealTime(lalgps, 0)
sidtime = GreenwichSiderealTime(lalgps, 0) sidtime = GreenwichSiderealTime(gps, 0)
sidtime = sidtime % (np.pi*2) sidtime = sidtime % (np.pi*2)
# -- Get the injection location # -- Get the injection location
...@@ -313,6 +303,16 @@ def make_skyview(directory='.', mdc=None, NSIDE=128, inj=None, npost=5000): ...@@ -313,6 +303,16 @@ def make_skyview(directory='.', mdc=None, NSIDE=128, inj=None, npost=5000):
print "Calculating searched area..." print "Calculating searched area..."
searcharea = (kde.searched_area( np.array([[injra, injdec]]) )*sq_deg)[0] searcharea = (kde.searched_area( np.array([[injra, injdec]]) )*sq_deg)[0]
print "Got searched area of {0}".format(searcharea) print "Got searched area of {0}".format(searcharea)
elif (inj is not None):
print "Calculating p value of command line injection"
injcontour = kde.p_values( np.array([inj['ra'], inj['dec']]) )[0]
print "Got contour of {0}".format(injcontour)
print "Calculating searched area..."
searcharea = (kde.searched_area( np.array([inj['ra'], inj['dec']]) )*sq_deg)[0]
print "Got searched area of {0}".format(searcharea)
else: else:
injcontour = 0 injcontour = 0
searcharea = 0 searcharea = 0
...@@ -373,6 +373,14 @@ if __name__ == "__main__": ...@@ -373,6 +373,14 @@ if __name__ == "__main__":
if ra is not None and dec is not None: if ra is not None and dec is not None:
injpos ={'ra': ra, 'dec':dec, 'id':0} injpos ={'ra': ra, 'dec':dec, 'id':0}
if mdc is not None:
print 'Reading mdc log {0}'.format(mdc)
try:
mdc = ft.Mdc(arg)
except:
print 'Warning! Failed to read mdc'
mdc = None
make_skyview(args[0], mdc, NSIDE, injpos) make_skyview(args[0], mdc, NSIDE, injpos)
......
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