[pycbc live] Allowing the use of template fits in the pycbc live ranking statistic

bin/pycbc_live

@@ -968,6 +968,9 @@ parser.add_argument('--gracedb-labels', metavar='LABEL', nargs='+',
 parser.add_argument('--ifar-upload-threshold', type=float, required=True,
                     help='Inverse-FAR threshold for uploading candidate '
                          'triggers to GraceDB, in years.')
+parser.add_argument('--coinc-threshold', type=float,
+                    help="The time window allowed between coincident triggers",
+                    default=0.002)
 parser.add_argument('--file-prefix', default='Live')
 
 parser.add_argument('--round-start-time', type=int, metavar='X',

pycbc/events/coinc.py

@@ -27,6 +27,7 @@
 
 import numpy, logging, pycbc.pnutils, pycbc.conversions, copy, lal
 from pycbc.detector import Detector, ppdets
+from pycbc.conversions import mchirp_from_mass1_mass2
 from .eventmgr_cython import coincbuffer_expireelements
 from .eventmgr_cython import coincbuffer_numgreater
 from .eventmgr_cython import timecoincidence_constructidxs
@@ -877,6 +878,7 @@ def __init__(self, num_templates, analysis_block, background_statistic,
 
         self.timeslide_interval = timeslide_interval
         self.return_background = return_background
+        self.coinc_threshold = coinc_threshold
 
         self.ifos = ifos
         if len(self.ifos) != 2:
@@ -971,6 +973,7 @@ def from_cli(cls, args, num_templates, analysis_chunk, ifos):
                    ifar_limit=args.background_ifar_limit,
                    timeslide_interval=args.timeslide_interval,
                    ifos=ifos,
+                   coinc_threshold=args.coinc_threshold,
                    **kwargs)
 
     @staticmethod
@@ -1097,9 +1100,11 @@ def _add_singles_to_buffer(self, results, ifos):
 
             if len(trigs['snr'] > 0):
                 trigsc = copy.copy(trigs)
+                trigsc['ifo'] = ifo
                 trigsc['chisq'] = trigs['chisq'] * trigs['chisq_dof']
                 trigsc['chisq_dof'] = (trigs['chisq_dof'] + 2) / 2
                 single_stat = self.stat_calculator.single(trigsc)
+                del trigsc['ifo']
             else:
                 single_stat = numpy.array([], ndmin=1,
                               dtype=self.stat_calculator.single_dtype)
@@ -1170,6 +1175,9 @@ def _find_coincs(self, results, valid_ifos):
                 trig_stat = trigs['stat'][i]
                 trig_time = trigs['end_time'][i]
                 template = trigs['template_id'][i]
+                mass1 = trigs['mass1'][i]
+                mass2 = trigs['mass2'][i]
+                mchirp = mchirp_from_mass1_mass2(mass1, mass2)
 
                 # Get current shift_ifo triggers in the same template
                 times = self.singles[shift_ifo].data(template)['end_time']
@@ -1206,12 +1214,18 @@ def _find_coincs(self, results, valid_ifos):
                 # ranking statistic values.
                 sngls_list = [[fixed_ifo, self.trig_stat_memory[:len(i1)]],
                               [shift_ifo, stats[i1]]]
-                c = self.stat_calculator.rank_stat_coinc(
-                    sngls_list,
-                    slide,
-                    self.timeslide_interval,
-                    shift_vec
-                )
+
+                if len(i1):
+                    c = self.stat_calculator.rank_stat_coinc(
+                        sngls_list,
+                        slide,
+                        self.timeslide_interval,
+                        shift_vec,
+                        time_addition=self.coinc_threshold,
+                        mchirp=mchirp
+                    )
+                else:
+                    c = numpy.array([])
 
                 # Store data about new triggers: slide index, stat value and
                 # times.

pycbc/events/stat.py

@@ -837,13 +837,16 @@ def find_fits(self, trigs):
             The thresh fit value(s)
         """
         try:
-            tnum = trigs.template_num  # exists if accessed via coinc_findtrigs
+            tnum = trigs.template_num
+        except AttributeError:
+            tnum = trigs['template_id']
+
+        try:
             ifo = trigs.ifo
         except AttributeError:
-            tnum = trigs['template_id']  # exists for SingleDetTriggers
-            assert len(self.ifos) == 1
-            # Should be exactly one ifo provided
-            ifo = self.ifos[0]
+            ifo = trigs['ifo']
+            assert ifo in self.ifos
+
         # 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]
@@ -1473,10 +1476,9 @@ def single(self, trigs):
             # exists if accessed via coinc_findtrigs
             self.curr_tnum = trigs.template_num
         except AttributeError:
-            # exists for SingleDetTriggers
+            # exists for SingleDetTriggers & pycbc_live get_coinc
             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 coinc
         # method does not have access to template id
         singles['benchmark_logvol'] = self.benchmark_logvol[self.curr_tnum]
@@ -1555,6 +1557,12 @@ def rank_stat_coinc(self, s, slide, step, to_shift,
 
         # First get signal PDF logr_s
         stat = {ifo: st for ifo, st in s}
+
+        try:
+            self.curr_mchirp = kwargs['mchirp']
+        except KeyError as e:
+            print(e)
+
         logr_s = self.logsignalrate(stat, slide * step, to_shift)
 
         # Find total volume of phase-time-amplitude space occupied by noise
@@ -1779,8 +1787,16 @@ def single(self, trigs):
             The array of single detector values
         """
         from pycbc.conversions import mchirp_from_mass1_mass2
-        self.curr_mchirp = mchirp_from_mass1_mass2(trigs.param['mass1'],
-                                                   trigs.param['mass2'])
+        try:
+            self.curr_mchirp = mchirp_from_mass1_mass2(trigs.param['mass1'],
+                                                       trigs.param['mass2'])
+        except AttributeError:
+            try:
+                self.curr_mchirp = mchirp_from_mass1_mass2(trigs['mass1'],
+                                                           trigs['mass2'])
+            except KeyError as e:
+                print(e)
+
         if self.mcm is not None:
             # Careful - input might be a str, so cast to float
             self.curr_mchirp = min(self.curr_mchirp, float(self.mcm))

test/test_live_coinc_compare.py

@@ -73,7 +73,8 @@ def setUp(self, *args):
             statistic_keywords=None,
             timeslide_interval=0.1,
             background_ifar_limit=100,
-            store_background=True
+            store_background=True,
+            coinc_threshold=0.002
         )
 
         # number of templates in the bank