Package pylal :: Module extended_background_utils
[hide private]
[frames] | no frames]

Source Code for Module pylal.extended_background_utils

  1  #!/usr/bin/env python 
  2  # read in the summary files and make coincidences etc. 
  3  from glue.ligolw import ligolw, lsctables, utils 
  4  from pylal.tools import XLALCalculateEThincaParameter 
  5  from glue import segments, lal, git_version 
  6  from glue.ligolw.utils import process 
  7  import copy 
  8  import lal 
  9  import numpy 
 10  import math 
 11  import glob 
 12  import optparse 
 13  import sys 
 14   
15 -def max_dt (row, eMatch):
16 # Calculate the timing uncertainty + earth travel time 17 # We use this as a conservative (but inexpensive to calculate) 18 # ethinca test. 19 a11 = row.Gamma0 / eMatch; 20 a12 = row.Gamma1 / eMatch; 21 a13 = row.Gamma2 / eMatch; 22 a22 = row.Gamma3 / eMatch; 23 a23 = row.Gamma4 / eMatch; 24 a33 = row.Gamma5 / eMatch; 25 26 x = (a23 * a23 - a22 * a33) * a22; 27 denom = (a12*a23 - a22*a13) * (a12*a23 - a22*a13) - (a23*a23 - a22*a33) * (a12*a12 - a22*a11); 28 29 return math.sqrt( x / denom ) + 2. * lal.REARTH_SI / lal.C_SI;
30
31 -def find_slide_coincs(htrigs, ltrigs, min_mchirp, max_mchirp, 32 ethinca, slide, threshold):
33 """ Calculate coincs from single detector triggers. We window based on 34 mchirp, apply a new snr threshold, and only keep triggers that lie on a 35 slide boundary. 36 """ 37 38 coinc_sngls = lsctables.New(lsctables.SnglInspiralTable) 39 num_trig = 0 40 t_snrsq = float(threshold) ** 2 41 42 # Collect the mchirp, new snr, and end time of the second detector's 43 # triggers into numpy arrays for faster math operations later. 44 l_ends = [] 45 lt_mc = [] 46 lt_sn = [] 47 lt_en = [] 48 for l in ltrigs: 49 lt_mc.append(l.mchirp) 50 lt_sn.append(l.get_new_snr()**2) 51 lt_en.append(float(l.get_end())) 52 lt_mc = numpy.array(lt_mc) 53 lt_sn = numpy.array(lt_sn) 54 lt_en = numpy.array(lt_en) 55 56 # The list of index location of "good" triggers in the second detector 57 # good triggers are ones that might possibly still form coincs with the 58 # first detector's triggers. Initially, this is all of them. 59 loc = numpy.arange(0, len(ltrigs)) 60 61 # Sort the first detector's triggers by newsnr ascending 62 stats = [] 63 htrigs.sort(lambda a, b: cmp(a.get_new_snr(), b.get_new_snr()), reverse=True) 64 65 # Iterate over each of the first detector's triggers and calculate what 66 # coincs it will generate with the second detector's triggers 67 for h in htrigs: 68 h_end = float(h.get_end()) 69 h_mchirp = h.mchirp 70 h_max_td = max_dt(h, ethinca) 71 h_snrsq = h.get_new_snr() ** 2 72 r_snrsq = t_snrsq - h_snrsq 73 lind = 0 74 75 # Grab only the information about triggers that could still form 76 # coincs 77 lt_snl = lt_sn[loc] 78 lt_mcl = lt_mc[loc] 79 lt_enl = lt_en[loc] 80 81 # Determine the time difference between these triggers and the 82 # considered "h" trigger 83 num_slides = numpy.round((h_end - lt_enl)/slide) 84 lt_en_tmp = lt_enl + num_slides * slide 85 td = abs(h_end - lt_en_tmp) 86 87 # Get location of coincs that are above the coinc new snr threshold 88 # Triggers that fall below this threshold are never considered again 89 goods = (lt_snl >= r_snrsq) 90 91 # Get the location of coincs that pass the rough time coincidence 92 # and are within the mchirp bin. We calculate ethinca only for these. 93 good = goods & (td <=h_max_td) & (lt_mcl >= 2*min_mchirp - h_mchirp) & (lt_mcl <= 2*max_mchirp - h_mchirp) 94 95 # For each remaining trigger we calculate the ethinca test and add 96 # to the list of coincs if they pass 97 for lind, ns in zip(loc[good], num_slides[good]): 98 99 # don't allow zerolag or offsets too close to it! 100 if abs(slide * ns) < 0.1: 101 continue 102 103 # Get the single inspiral trigger and apply the time offset 104 # so that we can perform coincidence 105 l = ltrigs[lind] 106 l_end_old = l.get_end() 107 l_end_tmp = l_end_old + float( slide * ns ) 108 109 l.set_end(l_end_tmp) 110 epar = XLALCalculateEThincaParameter(h, l) 111 112 # If we pass coincidence we add the first and second detector 113 # trigger pair to the output list 114 if epar < ethinca: 115 # This should never happen 116 if abs(slide*ns) < 0.1: 117 print 'I USED ZEROLAG!!' 118 exit(1) 119 hcopy = copy.deepcopy(h) 120 l.set_end(h.get_end()) 121 lcopy = copy.deepcopy(l) 122 hcopy.event_id = lsctables.SnglInspiralID(num_trig) 123 lcopy.event_id = lsctables.SnglInspiralID(num_trig) 124 num_trig += 1 125 coinc_sngls.append(hcopy) 126 coinc_sngls.append(lcopy) 127 l.set_end(l_end_old) 128 129 # Finally set the list of second triggers to consider to the reduced 130 # set that are still above the coinc new snr threshold 131 loc = loc[goods] 132 133 # No more triggers to consider 134 if len(loc) == 0: 135 break 136 return( coinc_sngls )
137
138 -def parse_mchirp_bins(param_range):
139 slots = param_range.split(";") 140 bins = [] 141 for slot in slots: 142 minb, maxb = slot[1:-1].split(',') 143 bins.append((float(minb), float(maxb))) 144 return bins
145
146 -def get_mchirp_bin(mchirp, mchirp_bins):
147 for mchirp_low, mchirp_high in mchirp_bins: 148 if mchirp >= mchirp_low and mchirp <= mchirp_high: 149 return mchirp_low, mchirp_high 150 return None, None
151
152 -def get_veto_coincs(coincs, veto_start, veto_end):
153 # Create a segment to veto the time given 154 veto_seg = segments.segmentlist() 155 veto_seg.append(segments.segment(veto_start, veto_end)) 156 157 # Remove the events around that time 158 veto_coincs = coincs.veto(veto_seg) 159 return veto_coincs
160
161 -def get_background_livetime(lt1, lt2, slide_step, veto_window=0):
162 """ 163 Returns the expected total time-slid coincident livetime for two detectors 164 having single-ifo livetimes lt1, lt2, when doing all possible slides with 165 time offsets being multiples of slide_step. 166 Specifying veto_window > 0 removes a given duration from each of the 167 single-ifo livetimes. 168 169 Derivation: Consider an analysis of length T with two ifos having duty 170 factors d1, d2. Then the expected coinc livetime per slide is T*d1*d2. 171 The maximum possible number of slides is T/s where s is the step. The 172 total time is the product which can be expressed as 173 (T*d1)(T*d2)/s = lt1*lt2/s. 174 """ 175 lt1 = float(lt1) 176 lt2 = float(lt2) 177 slide_step = float(slide_step) 178 veto_window = float(veto_window) 179 return (lt1 - veto_window*2) * (lt2 - veto_window*2) / slide_step
180