Commit 19a6703c authored by Francesco Pannarale's avatar Francesco Pannarale
Browse files

Let megasky.py place all plots in the plots/ folder and have it mark the...

Let megasky.py place all plots in the plots/ folder and have it mark the injected sky location if ra and dec are passed on as an explicit argument to the script

git-svn-id: https://svn.ligo.caltech.edu/svn/bayeswave/trunk@341 c56465c9-8126-4a4f-9d7d-ac845eff4865
parent 3195f223
......@@ -22,7 +22,7 @@ from readbwb import BwbParams
# Works with single output directory from
# BayesWaveBurst
# -------------------------------------------
def make_skyview(directory='.', mdc=None, NSIDE=256, results=None):
def make_skyview(directory='.', mdc=None, NSIDE=256, ra=None, dec=None, results=None):
# -- Close any open figures
plt.close('all')
......@@ -33,7 +33,6 @@ def make_skyview(directory='.', mdc=None, NSIDE=256, results=None):
if mdc is None:
print "No MDC log provided."
# -- Save PWD for future reference
topdir = os.getcwd()
......@@ -145,19 +144,44 @@ def make_skyview(directory='.', mdc=None, NSIDE=256, results=None):
print injtheta, injra
# -- Make plots directory, if needed
if not os.path.exists('./plots'):
os.makedirs('./plots')
plotsDir = './plots'
if not os.path.exists(plotsDir):
os.makedirs(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
else:
#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)
hp.mollview(skymap, title='Skymap from KDE clusters', nest=True)
if (ra is not None and dec is not None):
hp.projscatter(decrad, rarad, s=200, c='yellow', marker=(5,1), alpha=0.5)
if mdc is not None:
hp.projscatter(injtheta, injra, s=200, c='yellow', marker=(5,1), alpha=0.5)
plt.savefig('plots/skymap.png')
plt.savefig(plotsDir+'/skymap.png')
plt.close()
hp.write_map("skymap_{0}.fits".format(gps), skymap, nest=True)
# -- Calculate the 50 and 90% credible intervals
sq_deg = (180.0/np.pi)**2
print "Calculating sky area ..."
......@@ -175,7 +199,7 @@ def make_skyview(directory='.', mdc=None, NSIDE=256, results=None):
print "Got searched area of {0}".format(searcharea)
else:
injcontour = 0
searcharea = 0
searcharea = 0
pixarea = hp.nside2pixarea(NSIDE, degrees=True)
......@@ -189,14 +213,14 @@ def make_skyview(directory='.', mdc=None, NSIDE=256, results=None):
plt.figure()
plt.plot(acf,'o', markersize=2)
plt.title('RA autocorrelation', fontsize=14)
plt.savefig('ra_acf.png')
plt.savefig(plotsDir+'/ra_acf.png')
plt.close()
# -- Try scatter plot of ra
plt.figure()
plt.plot(ralist, 'o', markersize=2)
plt.title('RA chain: {0} samples'.format(ralist.size), fontsize=14)
plt.savefig('ra_chain.png')
plt.savefig(plotsDir+'/ra_chain.png')
# -- Make an output file with key statistics
outfile = open('skystats.txt', 'w')
......@@ -218,9 +242,9 @@ def make_skyview(directory='.', mdc=None, NSIDE=256, results=None):
skyview.write('The searched area is {0:.2f} sq. degrees.<br/>\n'.format(searcharea))
skyview.write('The injection was found at the {0:.2f} percent area contour.<br/>\n'.format(injcontour * 100.0))
skyview.write('<h2>Plots</h2>\n')
skyview.write('<img src="plots/skymap.png" width=600>\n')
skyview.write('<img src="ra_chain.png" width=600>\n')
skyview.write('<img src="ra_acf.png" width=600>\n')
skyview.write('<img src="'+plotsDir+'/skymap.png" width=600>\n')
skyview.write('<img src="'+plotsDir+'/ra_chain.png" width=600>\n')
skyview.write('<img src="'+plotsDir+'/ra_acf.png" width=600>\n')
skyview.write('</body></html>')
skyview.close()
......@@ -240,12 +264,15 @@ def make_skyview(directory='.', mdc=None, NSIDE=256, results=None):
# -- Write main script for command line running
if __name__ == "__main__":
opts, args = getopt.getopt(sys.argv[1:], "", ['directory=', 'mdc=', 'NSIDE='])
opts, args = getopt.getopt(sys.argv[1:], "", ['directory=', 'mdc=', 'NSIDE=', 'ra=', 'dec=', 'geo'])
# -- Set default argument values
directory='.'
mdc=None
NSIDE=256
ra=None
dec=None
geo=False
results=None
print "Got these arguments"
......@@ -270,6 +297,12 @@ if __name__ == "__main__":
NSIDE = int(arg)
if opt == '--sim':
sim = 'True' == arg
make_skyview(directory, mdc, NSIDE, results)
if opt == '--ra':
ra = float(arg)
if opt == '--dec':
dec = float(arg)
if opt == '--geo':
geo = 'True'
make_skyview(directory, mdc, NSIDE, ra, dec, results)
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