diff --git a/bin/pycbc_multi_inspiral b/bin/pycbc_multi_inspiral index ca4a6d0598c..b5b11d3a44d 100755 --- a/bin/pycbc_multi_inspiral +++ b/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,6 +107,10 @@ 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) @@ -113,10 +118,9 @@ 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()