The code version for the submission of soichio kuwahara master thesis

gstlal-burst/bin/gstlal_cherenkov_burst

@@ -116,7 +116,7 @@ def parse_command_line():
 	parser.add_option("--reference-psd", metavar = "filename", help = "Use fixed PSDs from this file (one for each detector) instead of measuring them from the data.  Optional.")
 	parser.add_option("--output", metavar = "filename", help = "Name of output xml file.")
 	parser.add_option("--rankingstat-output", metavar = "filename", help = "Name of output xml file to record rankingstat object.")
-	parser.add_option("--threshold", metavar = "snr_threshold", type = "float", default = 4.0, help = "SNR threshold (default = 4.0).")
+	parser.add_option("--threshold", metavar = "snr_threshold", type = "float", default = 3.0, help = "SNR threshold (default = 3.0).")
 	parser.add_option("--cluster-events", metavar = "cluster_events", type = "float", default = 0.1, help = "Cluster events with input timescale in seconds (default = 0.1).")
 	parser.add_option("--min-instruments", metavar = "min_instruments", type = "int", default = 2, help = "Minimum number of instruments (default is 2).")
 	parser.add_option("--snr-dump", metavar = "start,stop", help = "Turn on SNR time series dumping to a text file within the GPS time range [start, stop).")
@@ -133,7 +133,7 @@ def parse_command_line():
 	required_options = set(("template_bank", "time_slide_file", "output", "rankingstat_output"))
 	missing_options = set(option for option in required_options if getattr(options, option) is None)
 	if missing_options:
-		raise ValueError("missing required option(s) %s" % ", ".join("--%s" % option.subst("_", "-") for option in (required_options - missing_options)))
+		raise ValueError("missing required option(s) %s" % ", ".join("--%s" % option.replace("_", "-") for option in (required_options - missing_options)))
 
 	options.gps_start_time = lal.LIGOTimeGPS(options.gps_start_time)
 	options.gps_end_time = lal.LIGOTimeGPS(options.gps_end_time)
@@ -145,6 +145,7 @@ def parse_command_line():
 
 	if options.hoft_dump is not None:
 		options.hoft_dump = segments.segment(*map(lal.LIGOTimeGPS, options.hoft_dump.split(",")))
+
 	return options, process_params, datasourceinfo, filenames
 
 
@@ -258,7 +259,7 @@ class PipelineHandler(simplehandler.Handler):
 		# make templates, whiten, put into firbank
 		for i, row in enumerate(self.template_bank):
 			# obtain cherenkov radiation-like waveform. Since this wave form is linealy polarized, only the plus mode will be considered.
-			template_t[i], _ = lalsimulation.SimBurstCherenkovRadiation(row.central_freq, row.bandwidth, row.amplitude, 1.0 / sample_rate)
+			template_t[i], _ = lalsimulation.SimBurstCherenkovRadiation(row.bandwidth, 1.0, 1.0 / sample_rate)
 
 			# Resize the template as its deltaF matches to the PSD.
 			assert template_t[i].data.length <= fir_length, "data length is too huge comparted to the zero pad."
@@ -504,7 +505,7 @@ for instrument in all_inst:
 
 
 	with tempfile.NamedTemporaryFile(suffix = ".xml.gz") as fp:
-		template_bank.getColumnByName("ifo")[:] = [instrument for i in range(len(template_bank))]
+		template_bank.getColumnByName("ifo")[:] = [instrument] * len(template_bank)
 		xmldoc = ligolw.Document()
 		xmldoc.appendChild(ligolw.LIGO_LW())
 		xmldoc.childNodes[-1].appendChild(template_bank)

gstlal-burst/bin/gstlal_cherenkov_inj

@@ -66,9 +66,12 @@ def parse_command_line():
 		description = "GstLAL-based cherenkov burst injection pipeline."
 	)
 	parser.add_option("--dEoverdA", type = "float", help = "Set dEoverdA. Non input of this generates uniform random amplitude.")
-	parser.add_option("--gps-geocent-time", metavar = "s", help = "Set the start time of the tiling in GPS seconds (required).")
 	parser.add_option("--deltaT", metavar = "s", default = 0, type = "float", help = "Set the time interval of injections. Default value is 0 for injecting only one.")
-	parser.add_option("--num", type = "int", default = 1, help = "Set the number of injections. Default value is 1.")
+	parser.add_option("--start", type = "float", help = "start gps time of the injection. Non input will generate the injection in whole segements.")
+	parser.add_option("--stop", type = "float", help = "end gps time of the injection. Non input will generate the injection in whole segements.")
+	parser.add_option("--min-instruments", metavar = "num", default = 1, type = "int", help = "When constructing the segment list, require at least this many instruments to be on (default = 1).  See also --segments-file and --segments-name.")
+	parser.add_option("--segments-file", metavar = "filename", help = "Load the segment list for the injections from this file (required).  See also --segments-name and --min-instruments.")
+	parser.add_option("--segments-name", metavar = "name", default = "datasegments", help = "Use the segments with this name in the segments file (default = \"datasegments\")).  See also --segments-file and --min-instruments")
 	parser.add_option("--time-slide-file", metavar = "filename", help = "Associate injections with the first time slide ID in this XML file (required).")
 	parser.add_option("--output", metavar = "filename", help = "Set the name of the output file (default = stdout).")
 	parser.add_option("--verbose", action = "store_true", help = "Be verbose.")
@@ -78,14 +81,11 @@ def parse_command_line():
 	process_params = dict(options.__dict__)
 
 	# check for params
-	required_options = set(("output", "time_slide_file", "gps_geocent_time"))
+	required_options = set(("min_instruments", "segments_file", "segments_name", "output", "time_slide_file"))
 	missing_options = set(option for option in required_options if getattr(options, option) is None)
 	if missing_options:
 		raise ValueError("missing required option(s) %s" % ", ".join("--%s" % option.replace("_", "-") for option in missing_options))
 
-	# type-cast
-	options.gps_geocent_time = lal.LIGOTimeGPS(options.gps_geocent_time)
-
 	return options, process_params, filenames
 
 
@@ -124,8 +124,34 @@ process = ligolw_process.register_to_xmldoc(xmldoc, "cherenkov_burst_inj", proce
 
 time_slide_table = lsctables.TimeSlideTable.get_table(xmldoc)
 time_slide_id = time_slide_table[0].time_slide_id
+offset_vector = time_slide_table.as_dict()[time_slide_id]
 if options.verbose:
-	print("associating injections with time slide (%d) %s" % (time_slide_id, time_slide_table.as_dict()[time_slide_id]), file = sys.stderr)
+	print("associating injections with time slide (%d) %s" % (time_slide_id, offset_vector), file = sys.stderr)
+
+#
+# construct the segment list for injection times
+#
+
+
+segmentlists = ligolw_segments.segmenttable_get_by_name(ligolw_utils.load_filename(options.segments_file, contenthandler = LIGOLWContentHandler, verbose = options.verbose), options.segments_name)
+segmentlists.coalesce()
+
+if set(segmentlists) != set(offset_vector):
+	raise ValueError("time slide had instruments %s but segment list has instruments %$.  must be the same." % (", ".join(sorted(offset_vector)), ", ".join(sorted(segmentlists))))
+
+# injections are done so that they are coincident when the data has the
+# given offset vector applied.  so we apply the offset vector to the
+# segment lists, determine when the number of on instruments meets the
+# --min-instruments criterion, and then *without restoring the offsets* use
+# that segment list for the injections.  when, later, the injection code
+# performs these injections the injection times will have the opposite time
+# shift applied which will place them back into the available segments in
+# the respective detectors.
+
+segmentlists.offsets.update(offset_vector)
+segmentlist = segments_utils.vote(segmentlists.values(), options.min_instruments)
+if not abs(segmentlist):
+	raise ValueError("the --min-instruments criterion cannot be satisfied")
 
 
 #
@@ -148,7 +174,19 @@ sim_burst_tbl = xmldoc.childNodes[-1].appendChild(lsctables.New(lsctables.SimBur
 # populate the sim_burst table with injections
 #
 
-for i in range(options.num):
+if options.start is None and options.stop is None:
+	start, stop = segmentlist.extent()
+else:
+	start = options.start
+	stop = options.stop
+for i in range(math.ceil(float(start) / options.deltaT), math.ceil(float(stop) / options.deltaT)):
+	# we use floating point arithmetic to compute the times because we
+	# don't really care when they occur, not to the nanosecond.
+	time_geocent = lal.LIGOTimeGPS(i * options.deltaT)
+	# skip injections outside the desired segments
+	if time_geocent not in segmentlist:
+		continue
+	# add an injection at the desired time
 	sim_burst_tbl.append(sim_burst_tbl.RowType(
 		# metadata
 		process_id = process.process_id,
@@ -157,10 +195,10 @@ for i in range(options.num):
 		waveform = "Cherenkov",
 
 		# waveform parameters
-		time_geocent = options.gps_geocent_time + i * options.deltaT,
+		time_geocent = time_geocent,
 		frequency = math.nan,
 		# bandwidth parameter is actually not used. Hence, the length of Enterprise in StarTrek.is inputted as a fixed parameter.
-		bandwidth = 642.,
+		bandwidth = 641.,
 		egw_over_rsquared = math.nan,
 
 		# sky location and polarization axis orientation
@@ -171,7 +209,7 @@ for i in range(options.num):
 		# unnecessary columns
 		waveform_number = 0.,
 		# Create the uniform and random amplitude in log10.
-		amplitude = pow(10, random.uniform(-1, 4.5)) if options.dEoverdA is None else options.dEoverdA,
+		amplitude = pow(10, random.uniform(-1., 5.)) if options.dEoverdA is None else options.dEoverdA,
 		duration = math.nan,
 		q = math.nan,
 		hrss = math.nan,

gstlal-burst/bin/gstlal_cherenkov_plot_rankingstat

@@ -70,12 +70,24 @@ def parse_command_line():
 		description = "GstLAL-based cherenkov burst ranking statistic plotting tool."
 	)
 	parser.add_option("--ranking-stat-cache", metavar = "filename", help = "Also load the ranking statistic likelihood ratio data files listed in this LAL cache.  See lalapps_path2cache for information on how to produce a LAL cache file.")
+	parser.add_option("--instrument", type = "string", help = "choose the instrument for the pdf calculation.")
+	parser.add_option("--width", type = "float", help = "width of the gaussian distribution in rankingstat in log scale.")
+	parser.add_option("--slope", type = "float", help = "slope of the pdf in the region where gaussian noise is dominant.")
+	parser.add_option("--threshold", type = "float", help = "threshold of the SNR where the dominance changes.")
+	parser.add_option("--matching", type = "float", help = "The chi^2 / rho value when the unmatch between template and injections become dominant.")
 	parser.add_option("--output", metavar = "string", help = "Name of output png.")
 	parser.add_option("-v", "--verbose", action = "store_true", help = "Be verbose.")
+	parser.add_option("--is-calculate", action = "store_true", help = "Whether plot calculated PDF for adjustment.")
 	options, urls = parser.parse_args()
 
 	if options.ranking_stat_cache is not None:
 		urls += [CacheEntry(line).url for line in open(options.ranking_stat_cache)]
+	# check for input parameters
+	required_options = set(("instrument", "width", "slope", "threshold", "matching"))
+	if options.is_calculate:
+		missing_options = set(option for option in required_options if getattr(options, option) is None)
+		if missing_options:
+			raise ValueError("missing required option(s) %s" % ", ".join("--%s" % option.replace("_", "-") for option in missing_options))
 
 	# save for the process_params table
 	options.options_dict = dict(options.__dict__)
@@ -101,6 +113,26 @@ options, urls = parse_command_line()
 
 rankingstat = cherenkov_rankingstat.marginalize_rankingstat_urls(urls, verbose = options.verbose)
 
+#
+# create artficial PDF
+#
+
+def gaussian(x, mu, sig):
+	return numpy.exp(-(x - mu)**2. / (2. * sig**2.))
+
+def calculated_pdf(x, y, width, slope, snr, matching):
+	pdf = numpy.zeros((len(x), len(y)))
+
+	for i in range(len(x)):
+		if x[i] >= snr:
+			# Gaussian in log scale
+			pdf[i][:] = gaussian(numpy.log10(y), matching, width) + (-4. * math.log10(x[i]))
+
+		else:
+			# Gaussian with decreasing width
+			pdf[i][:] = gaussian(numpy.log10(y), matching + slope * math.log10(x[i] / snr), 0.05 + (width - 0.05) * math.log10(x[i]) / math.log10(snr)) + (-4. * math.log10(x[i]))
+	return pdf
+
 
 #
 # plot denominator snr, chisq PDFs
@@ -110,18 +142,22 @@ rankingstat = cherenkov_rankingstat.marginalize_rankingstat_urls(urls, verbose =
 for instrument, pdf in rankingstat.denominator.lnpdfs.items():
 	x,y = pdf.bins.centres()
 	lnpdf = pdf.at_centres()
+	if options.is_calculate:
+		calpdf = calculated_pdf(x, y, options.width, options.slope, options.threshold, options.matching)
 
 	fig = figure.Figure()
 	FigureCanvas(fig)
 	axes = fig.gca()
 
 	norm = matplotlib.colors.Normalize(vmin = -30., vmax = lnpdf.max())
-	mesh = axes.pcolormesh(x, y, lnpdf.T, norm = None, cmap = "afmhot", shading = "gouraud")
+	mesh = axes.pcolormesh(x, y, lnpdf.T, norm = norm, cmap = "afmhot", shading = "gouraud")
+	if options.is_calculate:
+		axes.pcolormesh(x, y, calpdf.T, norm = norm, cmap = "winter", shading = "gouraud")
 	axes.contour(x, y, lnpdf.T, colors = "k", linestyles = "-", linewidths = .5, alpha = .3)
 	axes.loglog()
 	axes.grid(which = "both")
-	axes.set_xlim((4, 500))
-	axes.set_ylim((.00001, 5))
+	axes.set_xlim((3., 500.))
+	axes.set_ylim((.000001, 5.))
 	fig.colorbar(mesh, ax = axes)
 	axes.set_xlabel("$\mathrm{SNR}$")
 	axes.set_ylabel("$\chi^{2}/ \mathrm{SNR}^{2}$")
@@ -144,8 +180,8 @@ for instrument, pdf in rankingstat.denominator.lnpdfs.items():
 	axes.contour(x, y, lnpdf.T, norm = norm, colors = "k", linestyles = "-", linewidths = .5, alpha = .3)
 	axes.loglog()
 	axes.grid(which = "both")
-	axes.set_xlim((4, 500))
-	axes.set_ylim((.00001, 5))
+	axes.set_xlim((3., 500.))
+	axes.set_ylim((.000001, 5))
 	colorbar = fig.colorbar(mesh, ax = axes)
 	colorbar.set_clim(-30, lnpdf.max())
 	axes.set_xlabel("$\mathrm{SNR}$")

gstlal-burst/bin/gstlal_cherenkov_plot_summary

@@ -47,12 +47,14 @@ matplotlib.rcParams.update({
 from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas
 import numpy
 from scipy import interpolate
+from scipy import optimize
 import sys
 
 from optparse import OptionParser
 
 from lal import rate
 from lalburst import binjfind
+import lalsimulation
 from lal.utils import CacheEntry
 from gstlal import cherenkov
 from gstlal import stats as gstlalstats
@@ -75,8 +77,9 @@ from ligo.lw import lsctables
 
 def parse_command_line():
 	parser = OptionParser()
-	parser.add_option("--ranking-stat-pdf", metavar = "filename", help = "Set the name of the xml file rankingstatpdf.")
 	parser.add_option("--candidates-cache", metavar = "filename", help = "Also load the candidates from files listed in this LAL cache.  See lalapps_path2cache for information on how to produce a LAL cache file.")
+	parser.add_option("--power", type = "float", help = "Power of the source. It is used for making the contour graph of true rate.")
+	parser.add_option("--ranking-stat-pdf", metavar = "filename", help = "Set the name of the xml file rankingstatpdf.")
 	parser.add_option("--verbose", "-v", action = "store_true", help = "Be verbose.")
 
 	options, urls = parser.parse_args()
@@ -147,7 +150,7 @@ class FAPFAR(object):
 		# the zero-lag LR distribution
 		rate_normalization_lr_threshold, = rankingstatpdf.zl_lr_lnpdf.argmax()
 		# FIXME:  override with a manual cut until an automatic algorithm can be found
-		rate_normalization_lr_threshold = -28.
+		rate_normalization_lr_threshold = -19.
 
 		# record trials factor, with safety checks
 		counts = rankingstatpdf.zl_lr_lnpdf.count
@@ -224,7 +227,7 @@ def rankingstat_pdfs_plot(rankingstatpdf):
 	line1, = axes.semilogy(rankingstatpdf.noise_lr_lnpdf.bins[0].centres(), numpy.exp(rankingstatpdf.noise_lr_lnpdf.at_centres()), linewidth = 1)
 	rankingstatpdf.zl_lr_lnpdf.normalize()
 	axes.semilogy(rankingstatpdf.zl_lr_lnpdf.bins[0].centres(), numpy.exp(rankingstatpdf.zl_lr_lnpdf.at_centres()), color = 'red')
-	axes.set_xlim(-35., -15)
+	axes.set_xlim(-27., 5)
 	axes.set_ylim(1e-4, 1e1) 
 	axes.set_xlabel("$\ln\mathcal{L}$")
 	axes.set_ylabel("$P$")
@@ -297,20 +300,19 @@ def create_rate_vs_lnL_plot(axes, stats, stats_label, fapfar):
 
 
 def RateVsThreshold(fapfar, zerolag_stats, background_stats):
-	zerolag_fig, axes = create_plot("$\ln \mathcal{L}$ Threshold", "Number of Events $\geq \ln \mathcal{L}$")
+	zerolag_fig, axes = create_plot(r"$\ln \mathcal{L}$ Threshold", "Number of Events $\geq \ln \mathcal{L}$")
 	zerolag_stats = numpy.array(sorted(zerolag_stats, reverse = True))
 	create_rate_vs_lnL_plot(axes, zerolag_stats, r"Observed", fapfar)
-	axes.set_title("Event Count vs.\ Ranking Statistic Threshold")
+	axes.set_title(r"Event Count vs.\ Ranking Statistic Threshold")
 
-	background_fig, axes = create_plot("$\ln \mathcal{L}$ Threshold", "Number of Events $\geq \ln \mathcal{L}$")
+	background_fig, axes = create_plot(r"$\ln \mathcal{L}$ Threshold", "Number of Events $\geq \ln \mathcal{L}$")
 	background_stats = numpy.array(sorted(background_stats, reverse = True))
 	create_rate_vs_lnL_plot(axes, background_stats, r"Observed (time shifted)", fapfar)
-	axes.set_title("Event Count vs.\ Ranking Statistic Threshold (Closed Box)")
+	axes.set_title(r"Event Count vs.\ Ranking Statistic Threshold (Closed Box)")
 	detection_threshold = zerolag_stats[0]
 	return zerolag_fig, background_fig, detection_threshold
 
 
-
 #
 # =============================================================================
 #
@@ -320,6 +322,81 @@ def RateVsThreshold(fapfar, zerolag_stats, background_stats):
 #
 
 
+def slope(x, y):
+	"""
+	From the x and y arrays, compute the slope at the x co-ordinates
+	using a first-order finite difference approximation.
+	"""
+	slope = numpy.zeros((len(x),), dtype = "double")
+	slope[0] = (y[1] - y[0]) / (x[1] - x[0])
+	for i in range(1, len(x) - 1):
+		slope[i] = (y[i + 1] - y[i - 1]) / (x[i + 1] - x[i - 1])
+	slope[-1] = (y[-1] - y[-2]) / (x[-1] - x[-2])
+	return slope
+
+
+def upper_err(y, yerr, deltax):
+	z = y + yerr
+	deltax = int(deltax)
+	upper = numpy.zeros((len(yerr),), dtype = "double")
+	for i in range(len(yerr)):
+		upper[i] = max(z[max(i - deltax, 0) : min(i + deltax, len(z))])
+	return upper - y
+
+
+def lower_err(y, yerr, deltax):
+	z = y - yerr
+	deltax = int(deltax)
+	lower = numpy.zeros((len(yerr),), dtype = "double")
+	for i in range(len(yerr)):
+		lower[i] = min(z[max(i - deltax, 0) : min(i + deltax, len(z))])
+	return y - lower
+
+
+def render_from_bins(axes, efficiency_numerator, efficiency_denominator, cal_uncertainty, filter_width, colour = "k", erroralpha = 0.3, linestyle = "-"):
+	# extract array of x co-ordinates, and the factor by which x
+	# increases from one sample to the next.
+	(x,) = efficiency_denominator.centres()
+	x_factor_per_sample = efficiency_denominator.bins[0].delta
+
+	# compute the efficiency, the slope (units = efficiency per
+	# sample), the y uncertainty (units = efficiency) due to binomial
+	# counting fluctuations, and the x uncertainty (units = samples)
+	# due to the width of the smoothing filter.
+	eff = efficiency_numerator.array / efficiency_denominator.array
+	dydx = slope(numpy.arange(len(x), dtype = "double"), eff)
+	yerr = numpy.sqrt(eff * (1. - eff) / efficiency_denominator.array)
+	xerr = numpy.array([filter_width / 2.] * len(yerr))
+
+	# compute the net y err (units = efficiency) by (i) multiplying the
+	# x err by the slope, (ii) dividing the calibration uncertainty
+	# (units = percent) by the fractional change in x per sample and
+	# multiplying by the slope, (iii) adding the two in quadradure with
+	# the y err.
+	net_yerr = numpy.sqrt((xerr * dydx)**2. + yerr**2. + (cal_uncertainty / x_factor_per_sample * dydx)**2.)
+
+	# compute net xerr (units = percent) by dividing yerr by slope and
+	# then multiplying by the fractional change in x per sample.
+	net_xerr = net_yerr / dydx * x_factor_per_sample
+
+	# plot the efficiency curve and uncertainty region
+	#patch = patches.Polygon(zip(numpy.concatenate((x, x[::-1])), numpy.concatenate((eff + upper_err(eff, yerr, filter_width / 2.), (eff - lower_err(eff, yerr, filter_width / 2.))[::-1]))), edgecolor = colour, facecolor = colour, alpha = erroralpha)
+	#axes.add_patch(patch)
+	line, = axes.plot(x, eff, colour + linestyle)
+
+	# compute 50% point and its uncertainty
+	A50 = None #optimize.bisect(interpolate.interp1d(x, eff - 0.5), x[0], x[-1], xtol = 1e-40)
+	A50_err = None #interpolate.interp1d(x, net_xerr)(A50)
+
+	# print some analysis FIXME:  this calculation needs attention
+	num_injections = efficiency_denominator.array.sum()
+	num_samples = len(efficiency_denominator.array)
+	#print("Bins were %g samples wide, ideal would have been %g" % (filter_width, (num_samples / num_injections / interpolate.interp1d(x, dydx)(A50)**2.0)**(1.0/3.0)), file=sys.stderr)
+	#print("Average number of injections in each bin = %g" % efficiency_denominator.array.mean(), file=sys.stderr)
+
+	return line, A50, A50_err
+
+
 def create_efficiency_plot(axes, all_injections, found_injections, detection_threshold):
 	# formats
 	axes.semilogx()
@@ -327,43 +404,53 @@ def create_efficiency_plot(axes, all_injections, found_injections, detection_thr
 	
 	# set desired yticks
 	axes.set_yticks((0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0))
-	axes.set_yticklabels(("\(0\)", "\(0.1\)", "\(0.2\)", "\(0.3\)", "\(0.4\)", "\(0.5\)", "\(0.6\)", "\(0.7\)", "\(0.8\)","\(0.9\)", "\(1.0\)"))
+	axes.set_yticklabels((r"\(0\)", r"\(0.1\)", r"\(0.2\)", r"\(0.3\)", r"\(0.4\)", r"\(0.5\)", r"\(0.6\)", r"\(0.7\)", r"\(0.8\)", r"\(0.9\)", r"\(1.0\)"))
 
 	axes.xaxis.grid(True, which = "major, minor")
 	axes.yaxis.grid(True, which = "major, minor")
 
 	# the denominator and numerator binning for the plot
-	efficiency_numerator = rate.BinnedArray(rate.NDBins((rate.LogarithmicBins(min(sim.egw_over_rsquared for sim in all_injections), max(sim.egw_over_rsquared for sim in all_injections), 400),)))
-	efficiency_denominator = rate.BinnedArray(rate.NDBins((rate.LogarithmicBins(min(sim.egw_over_rsquared for sim in all_injections), max(sim.egw_over_rsquared for sim in all_injections), 400),)))
+	efficiency_numerator = rate.BinnedArray(rate.NDBins((rate.LogarithmicBins(min(sim.amplitude for sim in all_injections), max(sim.amplitude for sim in all_injections), 400),)))
+	efficiency_denominator = rate.BinnedArray(efficiency_numerator.bins)
 
 	# add into the vins
 	for sim in found_injections:
-		efficiency_numerator[sim.egw_over_rsquared,] += 1
+		efficiency_numerator[sim.amplitude,] += 1
 	for sim in all_injections:
-		efficiency_denominator[sim.egw_over_rsquared,] += 1
+		efficiency_denominator[sim.amplitude,] += 1
 
 	# adjust window function normalization
-	#filter_width = 16.7
-	#windowfunc = rate.gaussian_window(filter_width)
-	#windowfunc /= windowfunc[len(windowfunc) / 2 + 1]
-	#rate.filter_array(efficiency_numerator.array, windowfunc)
-	#rate.filter_array(efficiency_denominator.array, windowfunc)
+	filter_width = 10
+	windowfunc = rate.gaussian_window(filter_width)
+	windowfunc /= windowfunc[int((len(windowfunc) + 1) / 2)]
+	rate.filter_array(efficiency_numerator.array, windowfunc)
+	rate.filter_array(efficiency_denominator.array, windowfunc)
 
 	# regularize: adjust unused bins so that the efficiency is 0 avoid Nan
 	assert(efficiency_numerator.array <= efficiency_denominator.array).all()
 	efficiency_denominator.array[(efficiency_numerator.array == 0) & (efficiency_denominator.array == 0)] = 1
 	
-	(x,) = efficiency_denominator.centres()
-	eff = efficiency_numerator.array / efficiency_denominator.array
-	print(eff, file = sys.stderr)
-	line, = axes.plot(x, eff)
+	# FIXME: for output checking only
+	eff = rate.BinnedArray(efficiency_numerator.bins, array = efficiency_numerator.array / efficiency_denominator.array)
+	#print(eff, file = sys.stderr)
 
-def Efficiency(all_injections, found_injections, detection_threshold):
-	fig, axes = create_plot("Injection Amplitude (\(\mathrm{J/s^2}\))", "Detection Efficiency")
-	create_efficiency_plot(axes, all_injections, found_injections, detection_threshold)
-	axes.set_title("Detection Efficiency vs.\ Amplitude")
-	return fig
+	line1, A50, A50_error = render_from_bins(axes, efficiency_numerator, efficiency_denominator, 0.2, filter_width)
+
+	# add a legend to the axis
+	axes.legend((line1,), (r"\noindent Injections with $\log \Lambda > %.2f$" % detection_threshold,), loc = "lower right")
 
+	# adjust limits
+	axes.set_ylim([0.0, 1.0])
+
+	# return eff for true rate contour
+	return eff
+
+
+def Efficiency(all_injections, found_injections, detection_threshold):
+	fig, axes = create_plot(r"Injection Amplitude ($\mathrm{J/m^2}$)", "Detection Efficiency")
+	eff = create_efficiency_plot(axes, all_injections, found_injections, detection_threshold)
+	axes.set_title(r"Detection Efficiency vs.\ Amplitude")
+	return fig, eff
 
 
 #
@@ -374,6 +461,39 @@ def Efficiency(all_injections, found_injections, detection_threshold):
 # =============================================================================
 #
 
+
+def create_truerate_contour(axes, power, eff):
+	# set beta and impact parameter axis
+	rate = numpy.zeros((300,300))
+	beta = 10.**numpy.linspace(0.01, 2., 300)
+	r = 10.**numpy.linspace(5., 9., 300)
+
+	print(eff, file = sys.stderr)
+	# calculate
+	for i in range(len(beta)):
+		for j in range(len(r)):
+			dE_dA = lalsimulation.SimBurstCherenkov_dE_dA(power, beta[i], r[j])
+			dE_dA = max(dE_dA, eff.bins[0].min)
+			dE_dA = min(dE_dA, eff.bins[0].max)
+			rate[i][j] = 3.890 / eff[dE_dA,]
+
+	# set color maps
+	norm = matplotlib.colors.LogNorm(vmin = 1., vmax = 10000.)
+	cmap = matplotlib.cm.viridis
+	cmap.set_bad("k")
+	mesh = axes.pcolormesh(beta, r, rate.T, norm = norm, cmap = cmap, shading = "gouraud")
+	contour = axes.contour(beta, r, rate.T, norm = norm, colors = "black", linestyle = "-", linewidth = .5, alpha = 1., levels = numpy.logspace(0., 3., 7, base = 10.))
+	fmt = matplotlib.ticker.LogFormatterMathtext()
+	axes.clabel(contour, contour.levels, colors = "black", fmt = fmt)
+	axes.loglog()
+	axes.grid(which = "both")
+
+def Truerate(power, eff):
+	fig, axes = create_plot(r"$\beta$", "Impact parameter(metres)")
+	create_truerate_contour(axes, power, eff)
+	axes.set_title(r"90\%% Confidence Rate Upper Limit for NCC-1701-D Enterprise ($P=%s$ W)" % gstlalplotutil.latexnumber("%.4g" % power))
+	return fig
+
 #
 # =============================================================================
 #
@@ -463,6 +583,7 @@ for filename, fig in [("rate_vs_threshold_open", rate_vs_thresh_zl), ("rate_vs_t
 
 all_injections = []
 found_injections = []
+is_injection = False
 
 for n, url in enumerate(urls, 1):
 	if options.verbose:
@@ -471,6 +592,7 @@ for n, url in enumerate(urls, 1):
 
 	if not is_injection_document(xmldoc):
 		continue
+	is_injection = True
 	coinc_def_id = lsctables.CoincDefTable.get_table(xmldoc).get_coinc_def_id(search = cherenkov.CherenkovBBCoincDef.search, search_coinc_type = cherenkov.CherenkovBBCoincDef.search_coinc_type, create_new = False)
 	lnL_index = dict((coinc.coinc_event_id, coinc.likelihood) for coinc in lsctables.CoincTable.get_table(xmldoc) if coinc.coinc_def_id == coinc_def_id)
 
@@ -492,9 +614,18 @@ for n, url in enumerate(urls, 1):
 		if sim_lnL_index.get(sim.simulation_id, float("-inf")) >= threshold:
 			found_injections.append(sim)
 
-efficiency = Efficiency(all_injections, found_injections, threshold)
+if is_injection:
+	efficiencyfig, eff = Efficiency(all_injections, found_injections, threshold)
 
-for ext in (".png", ".pdf"):
-	if options.verbose:
-		print("writing Efficiency%s with threshold of %s", (ext,threshold), file = sys.stderr)
-	efficiency.savefig("Efficiency%s" % ext)
+	for ext in (".png", ".pdf"):
+		if options.verbose:
+			print("writing Efficiency%s with threshold of %s" % (ext,threshold), file = sys.stderr)
+		efficiencyfig.savefig("Efficiency%s" % ext)
+
+	if options.power is not None:
+		rate = Truerate(options.power, eff)
+
+		for ext in (".png", ".pdf"):
+			if options.verbose:
+				print("writing True rate contour graph", file = sys.stderr)
+			rate.savefig("Truerate%s" % ext)

gstlal-burst/bin/gstlal_cherenkov_zl_rank_pdfs

@@ -60,6 +60,7 @@ def parse_command_line():
 	parser.add_option("--candidates-cache", metavar = "filename", help = "Also load the candidates from files listed in this LAL cache.  See lalapps_path2cache for information on how to produce a LAL cache file.")
 	parser.add_option("-v", "--verbose", action = "store_true", help = "Be verbose.")
 	parser.add_option("--ranking-stat-pdf", metavar = "filename", help = "Write zero-lag ranking statistic PDF to this file.  Must contain exactly one Cherenkov burst ranking statistic PDF object.")
+	parser.add_option("--is-timeshifted", action = "store_true", help = "Whether the open(zerolag) or closed(time shifted) box.")
 	options, urls = parser.parse_args()
 
 	paramdict = options.__dict__.copy()
@@ -119,7 +120,10 @@ for n, url in enumerate(urls, 1):
 		print("%d/%d: " % (n, len(urls)), end = "", file = sys.stderr)
 	xmldoc = ligolw_utils.load_url(url, contenthandler = cherenkov_rankingstat.LIGOLWContentHandler, verbose = options.verbose)
 	coinc_def_id = lsctables.CoincDefTable.get_table(xmldoc).get_coinc_def_id(search = cherenkov.CherenkovBBCoincDef.search, search_coinc_type = cherenkov.CherenkovBBCoincDef.search_coinc_type, create_new = False)
-	zl_time_slide_ids = frozenset(time_slide_id for time_slide_id, offsetvector in lsctables.TimeSlideTable.get_table(xmldoc).as_dict().items() if not any(offsetvector.values()))
+	if options.is_timeshifted:
+		zl_time_slide_ids = frozenset(time_slide_id for time_slide_id, offsetvector in lsctables.TimeSlideTable.get_table(xmldoc).as_dict().items() if any(offsetvector.values()))
+	else:
+		zl_time_slide_ids = frozenset(time_slide_id for time_slide_id, offsetvector in lsctables.TimeSlideTable.get_table(xmldoc).as_dict().items() if not any(offsetvector.values()))
 
 	for coinc in lsctables.CoincTable.get_table(xmldoc):
 		if coinc.coinc_def_id == coinc_def_id and coinc.time_slide_id in zl_time_slide_ids:

gstlal-burst/python/cherenkov/rankingstat.py

@@ -53,15 +53,17 @@ from gstlal.stats import trigger_rate
 # ============================================================================
 #
 
+
 #
 # Numerator
 #
 
 class LnSignalDensity(snglcoinc.LnLRDensity):
-	ln_chisq_const = {"H1": -1.5, "L1": -1.}
-	threshold = {"H1": 80., "L1": 50.}
-	width = {"H1": 2.3, "L1": 2.3}
-	ln_matching = {"H1": -7.26443022292087, "L1": -5.80914299031403}
+	ln_chisq_const = {"H1": -1.9, "L1": -2.0}
+	threshold = {"H1": 120., "L1": 100.}
+	width = {"H1": 1.0, "L1": 1.5}
+	ln_matching = {"H1": -9.5, "L1": -9.5}
+
 
 	def __init__(self, instruments, delta_t, min_instruments):
 		self.instruments = instruments