diff --git a/pycbc/events/stat.py b/pycbc/events/stat.py index 2b10f0b1fd2..a8e52c66d49 100644 --- a/pycbc/events/stat.py +++ b/pycbc/events/stat.py @@ -760,8 +760,10 @@ def coinc_lim_for_thresh( if not self.has_hist: self.get_hist() - fixed_stat = [b["snglstat"] for a, b in sngls_list if a != limifo][0] - s1 = thresh**2.0 - fixed_stat**2.0 + fixed_stat_sq = sum( + [b["snglstat"] ** 2 for a, b in sngls_list if a != limifo] + ) + s1 = thresh ** 2.0 - fixed_stat_sq # Assume best case scenario and use maximum signal rate s1 -= 2.0 * self.hist_max s1[s1 < 0] = 0 @@ -838,6 +840,9 @@ def __init__(self, sngl_ranking, files=None, ifos=None, **kwargs): self.kwargs.get("minimum_statistic_cutoff", -30.0) ) + # This will be used to keep track of the template number being used + self.curr_tnum = None + # Go through the keywords and add class information as needed: if self.kwargs["sensitive_volume"]: # Add network sensitivity beckmark @@ -850,7 +855,6 @@ def __init__(self, sngl_ranking, files=None, ifos=None, **kwargs): axis=0, ) self.benchmark_logvol = 3.0 * numpy.log(hl_net_med_sigma) - self.curr_tnum = None if self.kwargs["dq"]: # Reweight the noise rate by the dq reweighting factor @@ -959,11 +963,6 @@ def setup_segments(self, key): def find_dq_noise_rate(self, trigs, dq_state): """Get dq values for a specific ifo and dq states""" - try: - tnum = trigs.template_num - except AttributeError: - tnum = trigs["template_id"] - try: ifo = trigs.ifo except AttributeError: @@ -976,10 +975,10 @@ def find_dq_noise_rate(self, trigs, dq_state): if ifo in self.dq_rates_by_state: for i, st in enumerate(dq_state): - if isinstance(tnum, numpy.ndarray): - bin_name = self.dq_bin_by_tid[ifo][tnum[i]] + if isinstance(self.curr_tnum, numpy.ndarray): + bin_name = self.dq_bin_by_tid[ifo][self.curr_tnum[i]] else: - bin_name = self.dq_bin_by_tid[ifo][tnum] + bin_name = self.dq_bin_by_tid[ifo][self.curr_tnum] dq_val[i] = self.dq_rates_by_state[ifo][bin_name][st] return dq_val @@ -1100,8 +1099,8 @@ def find_fits(self, trigs): # fits_by_tid is a dictionary of dictionaries of arrays # indexed by ifo / coefficient name / template_id - alphai = self.fits_by_tid[ifo]["smoothed_fit_coeff"][tnum] - ratei = self.fits_by_tid[ifo]["smoothed_rate_above_thresh"][tnum] + alphai = self.fits_by_tid[ifo]["smoothed_fit_coeff"][self.curr_tnum] + ratei = self.fits_by_tid[ifo]["smoothed_rate_above_thresh"][self.curr_tnum] thresh = self.fits_by_tid[ifo]["thresh"] return alphai, ratei, thresh @@ -1162,6 +1161,15 @@ def lognoiserate(self, trigs): lognoisel: numpy.array Array of log noise rate density for each input trigger. """ + # What is the template number currently being used? + try: + # exists if accessed via coinc_findtrigs, this is a number + self.curr_tnum = trigs.template_num + except AttributeError: + # exists for SingleDetTriggers & pycbc_live get_coinc, + # this is a numpy array + self.curr_tnum = trigs["template_id"] + alphai, ratei, thresh = self.find_fits(trigs) sngl_stat = self.get_sngl_ranking(trigs) lognoisel = ( @@ -1237,14 +1245,6 @@ def single(self, trigs): if self.kwargs["sensitive_volume"]: # populate fields to allow sensitive volume factor calculation singles["sigmasq"] = trigs["sigmasq"][:] - try: - # exists if accessed via coinc_findtrigs - self.curr_tnum = trigs.template_num - except AttributeError: - # exists for SingleDetTriggers - self.curr_tnum = trigs["template_id"] - # Should only be one ifo fit file provided - assert len(self.ifos) == 1 # Store benchmark log volume as single-ifo information since # the ranking methods do not have access to template id singles["benchmark_logvol"] = self.benchmark_logvol[self.curr_tnum] @@ -1256,9 +1256,14 @@ def single(self, trigs): if self.kwargs["chirp_mass"]: from pycbc.conversions import mchirp_from_mass1_mass2 - self.curr_mchirp = mchirp_from_mass1_mass2( - trigs.param["mass1"], trigs.param["mass2"] - ) + try: + mass1 = trigs.param['mass1'] + mass2 = trigs.param['mass2'] + except AttributeError: + mass1 = trigs['mass1'] + mass2 = trigs['mass2'] + self.curr_mchirp = mchirp_from_mass1_mass2(mass1, mass2) + return numpy.array(singles, ndmin=1) def sensitive_volume_factor(self, sngls): @@ -1323,7 +1328,11 @@ def logsignalrate_shared(self, sngls_info): if self.kwargs["chirp_mass"]: # chirp mass reweighting - mchirp = min(self.curr_mchirp, self.mcm) + if isinstance(self.curr_mchirp, numpy.ndarray): + mchirp = numpy.minimum(self.curr_mchirp, self.mcm) + else: + # curr_mchirp will be a number + mchirp = min(self.curr_mchirp, self.mcm) sr_factor += numpy.log((mchirp / 20.0) ** (11.0 / 3.0)) if self.kwargs["kde"]: @@ -1400,8 +1409,11 @@ def rank_stat_coinc( # Basic noise rate: sum of noise rates multiplied by the # window they can form coincidences in ln_noise_rate = coinc_rate.combination_noise_lograte( - sngl_dict, kwargs["time_addition"] + sngl_dict, + kwargs["time_addition"], + dets=kwargs.get('dets', None), ) + ln_noise_rate -= self.benchmark_lograte # Basic option is not to have any signal-based assumptions, @@ -1414,7 +1426,9 @@ def rank_stat_coinc( # normalized # Extent of time-difference space occupied noise_twindow = coinc_rate.multiifo_noise_coincident_area( - self.hist_ifos, kwargs["time_addition"] + self.hist_ifos, + kwargs["time_addition"], + dets=kwargs.get('dets', None), ) # Volume is the allowed time difference window, multiplied by 2pi # for each phase difference dimension and by allowed range of SNR