Add a feature to fit noise trigger rates with respect to eccentricity

bin/all_sky_search/pycbc_fit_sngls_binned

@@ -184,21 +184,20 @@ for veto_file, veto_segment_name in zip(args.veto_file, args.veto_segment_name):
 
 ### Functions for doing the pruning (removal of trigs at loudest times)
 
-def get_pars(args, tag, m1, m2, s1z, s2z):
+def get_pars(args, tag, bank, tid):
     # here used for both pruning and binning params
     paramarg = getattr(args, tag+'_param')
     try:
-        # will fail if m1 is a float rather than a sequence
-        logging.info('Getting %s values for %i triggers' % (paramarg, len(m1)))
+        # will fail if mass1 is not in the bank
+        logging.info('Getting %s values for %i triggers' % (paramarg, len(bank['mass1'])))
     except:
         pass
-    return triggers.get_param(paramarg, args, m1, m2, s1z, s2z)
+    return triggers.get_param(paramarg, args, bank, tid)
 
 if args.prune_param:
     logging.info('Getting min and max param values')
-    prpars = get_pars(args, 'prune',
-                      templatef['mass1'][:], templatef['mass2'][:],
-                      templatef['spin1z'][:], templatef['spin2z'][:])
+    templateid = np.arange(len(templatef['mass1'])) # get the full range
+    prpars = get_pars(args, 'prune', templatef, templateid)
     minprpar = min(prpars)
     maxprpar = max(prpars)
     del prpars
@@ -232,8 +231,7 @@ if args.prune_param:
         lstat = statpruneall[loudest]
         ltid = tidpruneall[loudest]
         ltime = timepruneall[loudest]
-        m1, m2, s1z, s2z = triggers.get_mass_spin(templatef, ltid)
-        lbin = trstats.which_bin(get_pars(args, 'prune', m1, m2, s1z, s2z),
+        lbin = trstats.which_bin(get_pars(args, 'prune', templatef, ltid),
                                  minprpar, maxprpar,
                                  args.prune_bins, log=args.log_prune_param)
         # is the bin where the loudest trigger lives full already?
@@ -275,8 +273,7 @@ if args.prune_param:
 if trig_dur:
     binpars = tdur + args.min_duration
 else:
-    m1, m2, s1z, s2z = triggers.get_mass_spin(templatef, tid)
-    binpars = get_pars(args, 'bin', m1, m2, s1z, s2z)
+    binpars = get_pars(args, 'bin', templatef, tid)
 logging.info("Parameter range of triggers: %f - %f" %
                                                   (min(binpars), max(binpars)))
 

bin/all_sky_search/pycbc_fit_sngls_by_template

@@ -261,9 +261,8 @@ total_time = abs(all_segments)
 # do pruning (removal of trigs at N loudest times defined over param bins)
 if args.prune_param:
     logging.info('Getting min and max param values')
-    pars = triggers.get_param(args.prune_param, args,
-                     templatef['mass1'][:], templatef['mass2'][:],
-                     templatef['spin1z'][:], templatef['spin2z'][:])
+    templateid = np.arange(len(templatef['mass1'])) # get the full range
+    pars = triggers.get_param(args.prune_param, args, templatef, templateid)
     minpar = min(pars)
     maxpar = max(pars)
     del pars
@@ -297,9 +296,8 @@ if args.prune_param:
         lstat = statpruneall[loudest]
         ltid = tidpruneall[loudest]
         ltime = timepruneall[loudest]
-        m1, m2, s1z, s2z = triggers.get_mass_spin(templatef, ltid)
         lbin = trstats.which_bin(triggers.get_param(args.prune_param, args,
-                                                    m1, m2, s1z, s2z),
+                                                    templatef, ltid),
                                  minpar, maxpar, args.prune_bins,
                                                       log=args.log_prune_param)
         # is the bin where the loudest trigger lives full already?

bin/all_sky_search/pycbc_fit_sngls_over_multiparam

@@ -193,6 +193,7 @@ parser.add_argument("--fit-param", nargs='+',
                     help="Parameter(s) over which to regress the background "
                          "fit coefficients. Required. Either read from "
                          "template fit file or choose from mchirp, mtotal, "
+                         "q, eccentricity, "
                          "chi_eff, eta, tau_0, tau_3, template_duration, "
                          "a frequency cutoff in pnutils or a frequency function"
                          "in LALSimulation. To regress the background over "
@@ -388,13 +389,12 @@ if args.smoothing_method == 'n_closest' and n_required > nabove.sum():
 
 logging.info('Calculating template parameter values')
 bank = HFile(args.bank_file, 'r')
-m1, m2, s1z, s2z = triggers.get_mass_spin(bank, tid)
 
 parvals = []
 parnames = []
 
 for param, slog in zip(args.fit_param, args.log_param):
-    data = triggers.get_param(param, args, m1, m2, s1z, s2z)
+    data = triggers.get_param(param, args, bank, tid)
     if slog in ['false', 'False', 'FALSE']:
         logging.info('Using param: %s', param)
         parvals.append(data)
@@ -502,6 +502,14 @@ else:
             **kwarg_dict
         )
 
+non_positive_mask = smoothed_vals[:, 0] <= 0
+if numpy.any(non_positive_mask):
+    n_bad = numpy.sum(non_positive_mask)
+    raise ValueError(
+        f"{n_bad} template(s) have zero triggers above threshold. "
+        "This will cause problems for the fit coefficient values, and is likely "
+        "due to too small smoothing width or too large fitting threshold. "
+    )
 logging.info("Writing output")
 outfile = HFile(args.output, 'w')
 outfile['template_id'] = tid

pycbc/events/triggers.py

@@ -104,32 +104,40 @@ def get_mass_spin(bank, tid):
     return m1, m2, s1z, s2z
 
 
-def get_param(par, args, m1, m2, s1z, s2z):
+def get_param(par, args, bank, tid):
     """
-    Helper function
+    Helper function to extract parameters from bank and calculate
+    derived parameters
 
     Parameters
     ----------
     par : string
         Name of parameter to calculate
     args : Namespace object returned from ArgumentParser instance
         Calling code command line options, used for f_lower value
-    m1 : float or array of floats
-        First binary component mass (etc.)
+    bank : h5py File object
+        Bank parameter file
+    tid : integer or array of int
+        Indices of the entries to be returned
 
     Returns
     -------
     parvals : float or array of floats
         Calculated parameter values
     """
+    m1, m2, s1z, s2z = get_mass_spin(bank, tid)
     if par == 'mchirp':
         parvals = conversions.mchirp_from_mass1_mass2(m1, m2)
     elif par == 'mtotal':
         parvals = m1 + m2
     elif par == 'eta':
         parvals = conversions.eta_from_mass1_mass2(m1, m2)
+    elif par == 'q':
+        parvals = conversions.q_from_mass1_mass2(m1, m2)
     elif par in ['chi_eff', 'effective_spin']:
         parvals = conversions.chi_eff(m1, m2, s1z, s2z)
+    elif par == 'eccentricity':
+        parvals = bank['eccentricity'][:][tid]
     elif par == 'template_duration':
         # default to SEOBNRv4 duration function
         if not hasattr(args, 'approximant') or args.approximant is None: