Skip to content
Snippets Groups Projects
Commit 011b6a9a authored by Chad Hanna's avatar Chad Hanna
Browse files

bin/gstlal_inspiral_create_dt_dphi_snr_ratio_pdfs: program to figure out what...

bin/gstlal_inspiral_create_dt_dphi_snr_ratio_pdfs: program to figure out what probability you have of forming coincs with certain parameters
parent 1229493c
No related branches found
No related tags found
No related merge requests found
......@@ -16,72 +16,7 @@
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
import sys
import numpy
import numpy.random
import lal
import lal.series
import argparse
from glue.ligolw import utils as ligolw_utils
from gstlal.stats import inspiral_extrinsics
#import cProfile
import itertools
parser = argparse.ArgumentParser(description='Generate PDFs of extrinsic parameters')
parser.add_argument('--horizons', action = "append", help='add horizon distances in the form of e.g, H1=200.0')
parser.add_argument('--min-instruments', type=int, default = 1, help = 'min instruments to consider')
parser.add_argument('--snr-thresh', type=float, default = 4.0, help = 'set the snr minimum to define found')
parser.add_argument('--num-samples', type=int, default = 1000000, help = 'set the number of samples to compute')
parser.add_argument('--num-samples-random-time-phase', type=int, default = 10000, help = 'set the number of samples to compute')
parser.add_argument('--output-file', default = 'extparams.h5', help = 'set the output hdf5 file. Default extparams.h5')
args = parser.parse_args()
horizons = dict((x, float(y)) for (x,y) in [h.split("=") for h in args.horizons])
RS = inspiral_extrinsics.RandomSource(horizons)
DB = {}
IE = inspiral_extrinsics.InspiralExtrinsicParameters(horizons, args.snr_thresh, args.min_instruments)
for hist in IE.histograms.values():
print hist
i = 0
last = dict((k,0) for k in IE.histograms)
num = dict((k,0.) for k in IE.histograms)
#pr = cProfile.Profile()
#pr.enable()
while i < args.num_samples:
time, phase, snr, dist, prob = RS()
detected_ifos = tuple(sorted([ifo for ifo in snr if snr[ifo] >= args.snr_thresh]))
if any(abs(time[choice[0]] - time[choice[1]]) > IE.max_dt for choice in itertools.combinations(detected_ifos, 2)):
continue
# make sure at least one detector found it
if len(detected_ifos) >= args.min_instruments:
i += 1
else:
continue
num[detected_ifos]+=1.
missed_ifos = tuple(sorted([ifo for ifo in snr if snr[ifo] < args.snr_thresh]))
for ifo in missed_ifos:
del time[ifo]
del phase[ifo]
del snr[ifo]
try:
IE.insert(time, phase, snr, prob)
except ValueError as e:
print time, phase, snr, " was rejected"
if not(i % 10000):
print >> sys.stderr, "\nSIMULATION %d \n " % i
print "%20s | %20s | %20s | %20s" % ("IFO combos", "number of cells", "number of points", "fraction change")
print "".join(["="] * 88)
for combo in IE.histograms:
print "%20s | %20s | %20s | %20s" % (str(combo), IE.histograms[combo].num_cells, IE.histograms[combo].num_points, (IE.histograms[combo].num_cells - last[combo]) / num[combo])
last[combo] = IE.histograms[combo].num_cells
num[combo] = 0.
#pr.disable()
#pr.print_stats(sort='time')
IE.to_hdf5(args.output_file)
IE = inspiral_extrinsics.InspiralExtrinsics()
IE.to_hdf5("inspiral_dtdphi_pdf.h5")
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment