Commit 648ad75c authored by Salvatore Vitale's avatar Salvatore Vitale
Browse files

Merge branch 'lalmaster' into mergelal

parents 9e0f2e5e b657ecbb
......@@ -24,6 +24,10 @@
#
try:
from fpconst import NegInf
except ImportError:
NegInf = float("-inf")
import itertools
import math
import scipy.stats
......@@ -156,7 +160,7 @@ class LnLRDensity(snglcoinc.LnLRDensity):
except AttributeError:
self.mkinterps()
interps = self.interps
return sum(interps[param](*value) for param, value in params.items())
return sum(interps[param](*value) for param, value in params.items()) if params else NegInf
def __iadd__(self, other):
if type(self) != type(other) or set(self.densities) != set(other.densities):
......
......@@ -22,6 +22,9 @@
# Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
# MA 02110-1301, USA.
import matplotlib
matplotlib.use('Agg') #sets backend to not need to open windows
import numpy as np
import astropy.table as apt
import scipy.integrate as si
......
......@@ -5403,6 +5403,7 @@ def readCoincXML(xml_file, trignum):
from glue.ligolw import ligolw
from glue.ligolw import lsctables
from glue.ligolw import utils
coincXML = utils.load_filename(xml_file, contenthandler = lsctables.use_in(ligolw.LIGOLWContentHandler))
coinc = lsctables.CoincTable.get_table(coincXML)
coincMap = lsctables.CoincMapTable.get_table(coincXML)
......
......@@ -18,6 +18,7 @@ import random
from itertools import permutations
import shutil
import numpy as np
import math
# We use the GLUE pipeline utilities to construct classes for each
# type of job. Each class has inputs and outputs, which are used to
......@@ -138,71 +139,81 @@ def readLValert(threshold_snr=None,gid=None,flow=20.0,gracedb="gracedb",basepath
Based on Chris Pankow's script
"""
output=[]
from glue.ligolw import utils
from glue.ligolw import utils as ligolw_utils
from glue.ligolw import lsctables
from glue.ligolw import ligolw
class PSDContentHandler(ligolw.LIGOLWContentHandler):
class LIGOLWContentHandler(ligolw.LIGOLWContentHandler):
pass
lsctables.use_in(PSDContentHandler)
from glue.ligolw import param
from glue.ligolw import array
import subprocess
lsctables.use_in(LIGOLWContentHandler)
from lal import series as lalseries
from subprocess import Popen, PIPE
import lal
from lalsimulation import SimInspiralChirpTimeBound, GetApproximantFromString, IMRPhenomDGetPeakFreq
from ligo.gracedb.rest import GraceDb, HTTPError
try:
from gstlal import reference_psd
except ImportError:
reference_psd = None
cwd=os.getcwd()
os.chdir(basepath)
print "%s download %s coinc.xml"%(gracedb,gid)
subprocess.call([gracedb,"download", gid ,"coinc.xml"])
xmldoc=utils.load_filename("coinc.xml",contenthandler = PSDContentHandler)
coinctable = lsctables.CoincInspiralTable.get_table(xmldoc)
coinc_events = [event for event in coinctable]
sngltable = lsctables.SnglInspiralTable.get_table(xmldoc)
sngl_events = [event for event in sngltable]
ifos = coinctable[0].ifos.split(",")
trigSNR = coinctable[0].snr
print "Download %s coinc.xml" % gid
client = GraceDb()
xmldoc = ligolw_utils.load_fileobj(client.files(gid, "coinc.xml"), contenthandler = LIGOLWContentHandler)[0]
ligolw_utils.write_filename(xmldoc, "coinc.xml")
coinc_events = lsctables.CoincInspiralTable.get_table(xmldoc)
sngl_event_idx = dict((row.event_id, row) for row in lsctables.SnglInspiralTable.get_table(xmldoc))
ifos = sorted(coinc_events[0].instruments)
trigSNR = coinc_events[0].snr
# Parse PSD
srate_psdfile=16384
fhigh=None
psdfileobj = None
if downloadpsd:
print "gracedb download %s psd.xml.gz" % gid
subprocess.call([gracedb,"download", gid ,"psd.xml.gz"])
xmlpsd = lalseries.read_psd_xmldoc(utils.load_filename('psd.xml.gz',contenthandler = lalseries.PSDContentHandler))
if os.path.exists("psd.xml.gz"):
psdasciidic=get_xml_psds(os.path.realpath("./psd.xml.gz"),ifos,os.path.realpath('./PSDs'),end_time=None)
combine=np.loadtxt(psdasciidic[psdasciidic.keys()[0]])
print "Download %s psd.xml.gz" % gid
try:
psdfileobj = client.files(gid, "psd.xml.gz")
except HTTPError:
print "Failed to download %s psd.xml.gz. lalinference will estimate the psd itself." % gid
if psdfileobj is not None:
if reference_psd is not None:
xmlpsd = ligolw_utils.load_fileobj(psdfileobj, contenthandler = lalseries.PSDContentHandler)[0]
psd = lalseries.read_psd_xmldoc(xmlpsd)
ligolw_utils.write_filename(xmlpsd, "psd.xml.gz", gz = True)
else:
open("psd.xml.gz", "w").write(psdfileobj.read())
psdasciidic = get_xml_psds(os.path.realpath("./psd.xml.gz"),ifos,os.path.realpath('./PSDs'),end_time=None)
combine = np.loadtxt(psdasciidic[psdasciidic.keys()[0]])
srate_psdfile = pow(2.0, ceil( log(float(combine[-1][0]), 2) ) ) * 2
else:
print "Failed to gracedb download %s psd.xml.gz. lalinference will estimate the psd itself." % gid
coinc_map = lsctables.CoincMapTable.get_table(xmldoc)
for coinc in coinc_events:
these_sngls = [e for e in sngl_events if e.event_id in [c.event_id for c in coinc_map if c.coinc_event_id == coinc.coinc_event_id] ]
these_sngls = [sngl_event_idx[c.event_id] for c in coinc_map if c.coinc_event_id == coinc.coinc_event_id]
dur=[]
srate=[]
horizon_distance=[]
for e in these_sngls:
# Review: Replace this with a call to LALSimulation function at some point
if roq==False:
try:
p=Popen(["lalapps_chirplen","--flow",str(flow),"-m1",str(e.mass1),"-m2",str(e.mass2)],stdout=PIPE, stderr=PIPE, stdin=PIPE)
strlen = p.stdout.read()
dur.append(pow(2.0, ceil( log(max(8.0,float(strlen.splitlines()[2].split()[5]) + 2.0), 2) ) ) )
srate.append(pow(2.0, ceil( log(float(strlen.splitlines()[1].split()[5]), 2) ) ) * 2 )
except:
print "WARNING: lalapps_chirplen --flow",str(flow),"-m1",str(e.mass1),"-m2",str(e.mass2)," failed."
print "WARNING: make sure you have set manually seglen and srate in the .ini file."
snr = e.snr
eff_dist = e.eff_distance
chirplen = SimInspiralChirpTimeBound(flow, e.mass1 * lal.MSUN_SI, e.mass2 * lal.MSUN_SI, 0.0, 0.0)
fstop = IMRPhenomDGetPeakFreq(e.mass1, e.mass2, 0.0, 0.0)
dur.append(pow(2.0, ceil( log(max(8.0, chirplen + 2.0), 2) ) ) )
srate.append(pow(2.0, ceil( log(fstop, 2) ) ) * 2)
# determine horizon distance
if threshold_snr is not None:
if snr > threshold_snr:
horizon_distance.append(eff_dist * snr/threshold_snr)
if e.eff_distance is not None and not math.isnan(e.eff_distance):
if e.snr > threshold_snr:
horizon_distance.append(e.eff_distance * e.snr / threshold_snr)
else:
horizon_distance.append(2 * eff_dist)
horizon_distance.append(2 * e.eff_distance)
else:
if reference_psd is not None and psdfileobj is not None:
if not roq==False:
fstop = IMRPhenomDGetPeakFreq(e.mass1, e.mass2, 0.0, 0.0)
HorizonDistanceObj = reference_psd.HorizonDistance(f_min = flow, f_max = fstop, delta_f = 1.0 / 32.0, m1 = e.mass1, m2 = e.mass2)
horizon_distance.append(HorizonDistanceObj(psd[e.ifo], snr = threshold_snr)[0])
if srate:
if max(srate)<srate_psdfile:
srate = max(srate)
else:
srate = srate_psdfile
if downloadpsd:
if psdfileobj is not None:
fhigh = srate_psdfile/2.0 * 0.95 # Because of the drop-off near Nyquist of the PSD from gstlal
else:
srate = None
......@@ -435,13 +446,13 @@ def get_xml_psds(psdxml,ifos,outpath,end_time=None):
(end_time): trigtime for this event. Will be used a part of the PSD file name
"""
lal=1
from glue.ligolw import utils
from glue.ligolw import utils as ligolw_utils
try:
#from pylal import series
from lal import series as series
from lal import series as lalseries
lal=0
except ImportError:
print "ERROR, cannot import pylal.series in bppu/get_xml_psds()\n"
print "ERROR, cannot import lal.series in bppu/get_xml_psds()\n"
exit(1)
out={}
......@@ -470,7 +481,7 @@ def get_xml_psds(psdxml,ifos,outpath,end_time=None):
if not os.path.isfile(psdxml):
print "ERROR: impossible to open the psd file %s. Exiting...\n"%psdxml
sys.exit(1)
xmlpsd = series.read_psd_xmldoc(utils.load_filename(psdxml,contenthandler = series.PSDContentHandler))
xmlpsd = lalseries.read_psd_xmldoc(ligolw_utils.load_filename(psdxml,contenthandler = lalseries.PSDContentHandler))
# Check the psd file contains all the IFOs we want to analize
for ifo in ifos:
if not ifo in [i.encode('ascii') for i in xmlpsd.keys()]:
......@@ -512,24 +523,22 @@ def get_xml_psds(psdxml,ifos,outpath,end_time=None):
def get_trigger_chirpmass(gid=None,gracedb="gracedb"):
from glue.ligolw import lsctables
from glue.ligolw import ligolw
from glue.ligolw import utils
class PSDContentHandler(ligolw.LIGOLWContentHandler):
from glue.ligolw import utils as ligolw_utils
class LIGOLWContentHandler(ligolw.LIGOLWContentHandler):
pass
lsctables.use_in(PSDContentHandler)
import subprocess
lsctables.use_in(LIGOLWContentHandler)
from ligo.gracedb.rest import GraceDb
cwd=os.getcwd()
subprocess.call([gracedb,"download", gid ,"coinc.xml"])
xmldoc=utils.load_filename("coinc.xml",contenthandler = PSDContentHandler)
coinctable = lsctables.CoincInspiralTable.get_table(xmldoc)
coinc_events = [event for event in coinctable]
sngltable = lsctables.SnglInspiralTable.get_table(xmldoc)
sngl_events = [event for event in sngltable]
client = GraceDb()
xmldoc = ligolw_utils.load_fileobj(client.files(gid, "coinc.xml"), contenthandler = LIGOLWContentHandler)[0]
ligolw_utils.write_filename(xmldoc, "coinc.xml")
coinc_events = lsctables.CoincInspiralTable.get_table(xmldoc)
sngl_event_idx = dict((row.event_id, row) for row in lsctables.SnglInspiralTable.get_table(xmldoc))
coinc_map = lsctables.CoincMapTable.get_table(xmldoc)
mass1 = []
mass2 = []
for coinc in coinc_events:
these_sngls = [e for e in sngl_events if e.event_id in [c.event_id for c in coinc_map if c.coinc_event_id == coinc.coinc_event_id] ]
these_sngls = [sngl_event_idx[c.event_id] for c in coinc_map if c.coinc_event_id == coinc.coinc_event_id]
for e in these_sngls:
mass1.append(e.mass1)
mass2.append(e.mass2)
......@@ -998,13 +1007,13 @@ class LALInferencePipelineDAG(pipeline.CondorDAG):
if self.config.has_option('input','burst-injection-file'):
#from pylal import SimBurstUtils
from glue.ligolw import lsctables
from glue.ligolw import utils
from glue.ligolw import utils as ligolw_utils
from glue.ligolw import ligolw
injfile=self.config.get('input','burst-injection-file')
class PSDContentHandler(ligolw.LIGOLWContentHandler):
class LIGOLWContentHandler(ligolw.LIGOLWContentHandler):
pass
lsctables.use_in(PSDContentHandler)
injTable=lsctables.SimBurstTable.get_table(utils.load_filename(injfile,contenthandler = PSDContentHandler))
lsctables.use_in(LIGOLWContentHandler)
injTable=lsctables.SimBurstTable.get_table(ligolw_utils.load_filename(injfile,contenthandler = LIGOLWContentHandler))
events=[Event(SimBurst=inj) for inj in injTable]
self.add_pfn_cache([create_pfn_tuple(self.config.get('input','burst-injection-file'))])
# SnglInspiral Table
......
......@@ -40,7 +40,11 @@ parser.add_argument('--analytic-csv-dir', type=str,
parser.add_argument('--pptest', action='store_true',
default=False,
help='Runs a P-P analysis.')
help='Runs a P-P analysis. Must specify a single engine.')
parser.add_argument('--pp-approximant', action='store', type=str,
default='IMRPhenomPv2pseudoFourPN,SEOBNRv4_ROMpseudoFourPN',
help='Approximant(s) used for the P-P analysis.')
parser.add_argument('--bbh-injection', type=str, nargs='?',
default=False,
......@@ -93,7 +97,7 @@ lalinf_prefix=''
try:
lalinf_prefix=os.environ['LALINFERENCE_PREFIX']
except KeyError:
print 'LALINFERENCE_PREFIX variable not defined, could not find LALInference installation.'
print('LALINFERENCE_PREFIX variable not defined, could not find LALInference installation.')
sys.exit()
def init_ini_file(file=args.ini_file):
......@@ -278,7 +282,6 @@ if args.analytic_tests:
os.chdir(args.output+'/' + test_func + '/')
shutil.copy(args.bbh_injection,args.output+'/'+test_func+'/')
analytic_ini_file=os.path.join(args.output,test_func,'analytic.ini')
cpanalytic=set_analytic_test(init_ini_file(), test_func)
......@@ -299,17 +302,20 @@ if args.analytic_tests:
############################################################
def set_pptest(cp):
def set_pptest(cp, engine, approximant):
cp.set('paths','webdir',web_outputdir+'/pptest/webdir/')
cp.set('paths','webdir', os.path.join(web_outputdir, 'pptest', engine, approximant, 'webdir/'))
cp.set('ppanalysis','webdir', os.path.join(web_outputdir, 'PPcheck', engine, approximant + '/'))
cp.set('lalinference','fake-cache',"{'H1':'LALSimAdLIGO','L1':'LALSimAdLIGO','V1':'LALSimAdVirgo'}")
cp.set('analysis','dataseed','1234')
cp.set('analysis', 'engine', engine)
cp.remove_option('engine','margphi')
cp.set('engine','margtime','')
cp.set('engine','amporder','-1')
cp.set('engine','fref','0')
cp.set('engine','distance-max','2000')
cp.set('engine', 'approx', approximant)
cp.set('resultspage','deltaLogP','7')
......@@ -317,24 +323,29 @@ def set_pptest(cp):
if args.pptest:
os.makedirs(args.output+'/pptest/')
os.chdir(args.output+'/pptest/')
# The PP test needs a single engine to be specified
for engine in args.engine.split(','):
for approximant in args.pp_approximant.split(','):
os.makedirs(os.path.join(args.output, 'pptest', engine, approximant + '/'))
os.chdir(os.path.join(args.output, 'pptest', engine, approximant + '/'))
pptest_ini_file=os.path.join(args.output,'pptest','pptest.ini')
pptest_ini_file=os.path.join(args.output, 'pptest', engine, approximant, 'pptest.ini')
cppptest=set_pptest(init_ini_file())
with open(pptest_ini_file,'w') as cpfile:
cppptest.write(cpfile)
cppptest=set_pptest(init_ini_file(), engine=engine,
approximant=approximant)
with open(pptest_ini_file,'w') as cpfile:
cppptest.write(cpfile)
lalinferenceargs = [ 'lalinference_pp_pipe'
, '-r'
, './run'
, '-N'
, '100'
, pptest_ini_file ]
lalinferenceargs = [ 'lalinference_pp_pipe'
, '-r'
, './run'
, '-N'
, '100'
, pptest_ini_file ]
subprocess.call(lalinferenceargs)
subprocess.call(lalinferenceargs)
if args.condor_submit:
condor_submit_dag = ['condor_submit_dag', os.path.join(args.output, 'pptest', engine, approximant, 'run/priortest.dag')]
subprocess.call(condor_submit_dag)
if args.condor_submit:
condor_submit_dag = ['condor_submit_dag',args.output+'/pptest/run/priortest.dag']
subprocess.call(condor_submit_dag)
Supports Markdown
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