added Left and Right Circular Projection matrix

bin/pycbc_multi_inspiral

@@ -22,6 +22,7 @@ from collections import defaultdict
 import argparse
 from pycbc import vetoes, psd, waveform, events, strain, scheme, fft,\
     DYN_RANGE_FAC
+from pycbc.vetoes import sgchisq
 from pycbc.filter import MatchedFilterControl
 from pycbc.types import TimeSeries, zeros, float32, complex64
 from pycbc.types import MultiDetOptionAction
@@ -106,17 +107,20 @@ parser.add_argument("--null-step", type=float, default=20., help="""
                     according to the null_grad and null_min.""")
 parser.add_argument("--trigger-time", type=int, help="""
                     Time of the GRB, used to set the antenna patterns.""")
+parser.add_argument("--projection", default="standard",
+                    choices=["standard", "left", "right", "left+right"],
+                    help="""Choice of projection matrix. 'Left' and 'Right' 
+                    corresponds to face-away and face-on""")
 
 # Add options groups
 strain.insert_strain_option_group_multi_ifo(parser)
 strain.StrainSegments.insert_segment_option_group_multi_ifo(parser)
 psd.insert_psd_option_group_multi_ifo(parser)
 scheme.insert_processing_option_group(parser)
 fft.insert_fft_option_group(parser)
-from pycbc.vetoes.sgchisq import SingleDetSGChisq
 pycbc.opt.insert_optimization_option_group(parser)
 pycbc.inject.insert_injfilterrejector_option_group_multi_ifo(parser)
-SingleDetSGChisq.insert_option_group(parser)
+sgchisq.SingleDetSGChisq.insert_option_group(parser)
 opt = parser.parse_args()
 
 # Use the time in the middle of the segment to calculate the antenna patterns.
@@ -139,7 +143,6 @@ logging.basicConfig(format='%(asctime)s : %(message)s', level=log_level)
 inj_filter_rejector = pycbc.inject.InjFilterRejector.from_cli_multi_ifos(opt,opt.instruments)
 ctx = scheme.from_cli(opt)
 
-
 def network_chisq(chisq, chisq_dof, snr_dict):
     ifos = sorted(snr_dict.keys())
     chisq_per_dof = {}
@@ -161,7 +164,7 @@ def pycbc_reweight_snr(network_snr, network_chisq, a = 3, b = 1. / 6.):
                           trigger
             network_chisq: A chisq value for each trigger 
     """
-    denom = ((1 + network_chisq)**a) / 2
+    denom = (1 + (network_chisq)**a) / 2
     reweighted_snr = network_snr / denom**b
     return reweighted_snr
 
@@ -197,20 +200,33 @@ def get_weighted_antenna_patterns(Fp_dict, Fc_dict, sigma_dict):
     return wp, wc
 
 
-def get_projection_matrix(wp, wc):
+def get_projection_matrix(wp, wc, projection='standard'):
     """
     Output: projection_matrix: Projects the data onto the signal space
     Input:  wp,wc: The weighted antenna response fuctions to plus and cross
                    polarisations respectively
+            projection: standard, left or right
     """
-    denominator = np.dot(wp, wp) * np.dot(wc, wc) - np.dot(wp, wc)**2
-    projection_matrix = (np.dot(wc, wc)*np.outer(wp, wp) +
-                         np.dot(wp, wp)*np.outer(wc, wc) -
-                         np.dot(wp, wc)*(np.outer(wp, wc) +
-                         np.outer(wc, wp))) / denominator
-    return projection_matrix
-
+    if projection=='standard':
+        denominator = np.dot(wp, wp) * np.dot(wc, wc) - np.dot(wp, wc)**2
+        projection_matrix = (np.dot(wc, wc)*np.outer(wp, wp) +
+                             np.dot(wp, wp)*np.outer(wc, wc) -
+                             np.dot(wp, wc)*(np.outer(wp, wc) +
+                             np.outer(wc, wp))) / denominator
+    elif projection=='left':
+        projection_matrix = ((np.outer(wp, wp) + np.outer(wc, wc) +
+                             (np.outer(wp, wc) - np.outer(wc, wp)) *1j)
+                             / (np.dot(wp, wp) + np.dot(wc, wc)))
+    elif projection=='right':
+        projection_matrix = ((np.outer(wp, wp) + np.outer(wc, wc) +
+                             (np.outer(wc, wp) - np.outer(wp, wc)) *1j)
+                             / (np.dot(wp, wp) + np.dot(wc, wc)))
+    else:
+        raise ValueError('Unknown projection type: {}'.format(projection))
 
+    return projection_matrix
+
+
 def coherent_snr(snr_triggers, index, threshold, projection_matrix,
                 coinc_snr=[]):
     """
@@ -269,7 +285,7 @@ def coincident_snr(snr_dict, index, threshold, time_delay_idx):
 
 
 def null_snr(rho_coh, rho_coinc, null_min=5.25, null_grad=0.2, null_step=20.,
-             index={}, snrv={}):
+             index={}, snrv={}, nullcut=False):
     """
     Output: null: null snr for surviving triggers
             rho_coh: Coherent snr for surviving triggers
@@ -298,11 +314,14 @@ def null_snr(rho_coh, rho_coinc, null_min=5.25, null_grad=0.2, null_step=20.,
     keep2 = np.logical_and(null < (rho_coh * null_grad + null_min),
                           rho_coh > null_step)
     keep = np.logical_or(keep1, keep2)
-    index = index[keep]
-    rho_coh  = rho_coh[keep]
-    snrv = {ifo : snrv[ifo][keep] for ifo in snrv.keys()}
-    rho_coinc = rho_coinc[keep]
-    null = null[keep]
+    if nullcut==True:
+        index = index[keep]
+        rho_coh  = rho_coh[keep]
+        snrv = {ifo : snrv[ifo][keep] for ifo in snrv.keys()}
+        rho_coinc = rho_coinc[keep]
+        null = null[keep]
+    else:
+        logging.info("nullcut is not applied")
     return null, rho_coh, rho_coinc, index, snrv
 
 
@@ -386,6 +405,8 @@ with ctx:
 
     # Matched filter each ifo. Don't cluster here for a coherent search.
     # Clustering happens at the end of the template loop.
+    # FIXME: The single detector SNR threshold should not necessarily be
+    #        applied to every IFO (usually only 2 most sensitive in network)
     matched_filter = {ifo : MatchedFilterControl(
         opt.low_frequency_cutoff, None, opt.snr_threshold, tlen, delta_f,
         complex64, segments[ifo], template_mem, use_cluster=False,
@@ -442,6 +463,7 @@ with ctx:
         'coherent_snr'   : float32,
         'null_snr'      : float32,
         'nifo'          : int,
+        'my_network_chisq':float32,
         'reweighted_snr' : float32
                }
     network_out_vals = {
@@ -451,6 +473,7 @@ with ctx:
         'coherent_snr'   : None,
         'null_snr'      : None,
         'nifo'          : None,
+        'my_network_chisq' : None, 
         'reweighted_snr' : None
                }
     network_names = sorted(network_out_vals.keys())
@@ -588,26 +611,48 @@ with ctx:
                 logging.info("%s points above coincident SNR threshold" % \
                              str(len(coinc_idx)))
                 if len(coinc_idx) != 0:
-                    logging.info("Max coincident SNR = %s"
-                                % str(max(rho_coinc)))
+                    logging.info("Max coincident SNR = %.2f", max(rho_coinc))
+
             # If there is only one ifo, then coinc_triggers is just the
             # triggers from ifo
             elif len(coinc_idx) != 0 and nifo == 1:
                 coinc_triggers = {opt.instruments[0]:
-                                      snr[opt.instruments[0]][coinc_idx_det_frame[opt.instruments[0]]]}
+                        snr[opt.instruments[0]][\
+                                coinc_idx_det_frame[opt.instruments[0]]]}
             else:
                 coinc_triggers = {}
                 logging.info("No triggers above coincident SNR threshold")
+
             # If we have triggers above coinc threshold and more than 2 ifos
             # then calculate the coherent statistics
             if len(coinc_idx) != 0 and nifo > 2:
                 wp,wc=get_weighted_antenna_patterns(Fp,Fc,sigma)
-                projection_matrix = get_projection_matrix(wp,wc)
-                rho_coh, coinc_idx, coinc_triggers, rho_coinc = coherent_snr(
-                    coinc_triggers, coinc_idx, opt.coinc_threshold,
-                    projection_matrix, rho_coinc)
-                logging.info("%s points above coherent threshold"
-                            % str(len(rho_coh)))
+                if opt.projection=='left+right':
+                    project_l = get_projection_matrix(wp, wc,
+                                                      projection='left')
+                    rho_coh_l, coinc_idx_l, coinc_triggers_l, rho_coinc_l = \
+                            coherent_snr(coinc_triggers, coinc_idx,
+                                         opt.coinc_threshold, project_l,
+                                         rho_coinc)
+                    project_r = get_projection_matrix(wp, wc,
+                                                      projection='right')
+                    rho_coh_r, coinc_idx_r, coinc_triggers_r, rho_coinc_r = \
+                            coherent_snr(coinc_triggers, coinc_idx,
+                                         opt.coinc_threshold, project_r,
+                                         rho_coinc)
+                    max_idx = np.argmax([rho_coh_l, rho_coh_r], axis=0)
+                    rho_coh = np.where(max_idx==0, rho_coh_l, rho_coh_r)
+                    coinc_idx = np.where(max_idx==0, coinc_idx_l, coinc_idx_r)
+                    coinc_triggers = np.where(max_idx==0, coinc_triggers_l,
+                                              coinc_triggers_r)
+                    rho_coinc = np.where(max_idx==0, rho_coinc_l, rho_coinc_r)
+                else:
+                    project = get_projection_matrix(wp, wc, projection=opt.projection)
+                    rho_coh, coinc_idx, coinc_triggers, rho_coinc = \
+                            coherent_snr(coinc_triggers, coinc_idx,
+                                         opt.coinc_threshold, project,
+                                         rho_coinc)
+                logging.info("%d points above coherent threshold", len(rho_coh))
                 if len(coinc_idx) != 0:
                     logging.info("Max coherent SNR = %s" % str(max(rho_coh)))
                     #Find the null snr
@@ -664,6 +709,7 @@ with ctx:
                 # entries that are single values need to be repeated once
                 # per event.
                 num_events = len(reweighted_snr)
+                # the output will only be possible if len(networkchi2) == num_events
 
                 for ifo in opt.instruments:
                     ifo_out_vals['bank_chisq'], ifo_out_vals['bank_chisq_dof'] = \
@@ -693,6 +739,8 @@ with ctx:
                     network_out_vals['coherent_snr'] = \
                         abs(snr[opt.instruments[0]][coinc_idx_det_frame[opt.instruments[0]]])
                 network_out_vals['reweighted_snr'] = reweighted_snr
+                network_out_vals['my_network_chisq'] = np.real(network_chisq_dict)
+
                 network_out_vals['time_index'] = \
                     coinc_idx+stilde[ifo].cumulative_index
                 network_out_vals['nifo'] = [nifo] * num_events
@@ -702,7 +750,7 @@ with ctx:
                                 [network_out_vals[n] for n in network_names])
         event_mgr.cluster_template_network_events(
             "time_index",
-            "reweighted_snr",
+            "coherent_snr",
             cluster_window
             )
         event_mgr.finalize_template_events()