1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 """Utilities for the coherent PTF inspiral analysis.
18 """
19
20 import os
21 import sys
22 import numpy
23 import glob
24 import math
25 import re
26
27 from pylal import grbsummary, antenna, InspiralUtils, SimInspiralUtils
28 from lal import PI as LAL_PI
29 from lal import MTSUN_SI as LAL_MTSUN_SI
30
31 from glue import segmentsUtils
32 from glue.ligolw import lsctables, table
33 from glue.ligolw.utils import process as ligolw_process
34
35
37 """Returns the chisq value needed to weight snr into new_snr
38 """
39 chisqnorm = (snr/new_snr)**q
40 if chisqnorm <= 1:
41 return 1E-20
42 return chisq_dof * (2*chisqnorm - 1)**(n/q)
43
44
45 -def get_bestnr( trig, q=4.0, n=3.0, null_thresh=(4.25,6), snr_threshold=6.,\
46 sngl_snr_threshold=4., chisq_threshold = None,\
47 null_grad_thresh=20., null_grad_val = 1./5.):
48 """
49 Calculate BestNR (coh_PTF detection statistic) through signal based vetoes:
50
51 The signal based vetoes are as follows:
52 * Coherent SNR < 6
53 * Bank chi-squared reduced (new) SNR < 6
54 * Auto veto reduced (new) SNR < 6
55 * Single-detector SNR (from two most sensitive IFOs) < 4
56 * Null SNR (CoincSNR^2 - CohSNR^2)^(1/2) < nullthresh
57 Returns BestNR as float
58 """
59
60 snr = trig.snr
61 if not chisq_threshold:
62 chisq_threshold = snr_threshold
63
64
65 if (snr < snr_threshold) \
66 or (trig.get_new_snr(index=q, nhigh=n, column='bank_chisq')\
67 < chisq_threshold) \
68 or (trig.get_new_snr(index=q, nhigh=n, column='cont_chisq')\
69 < chisq_threshold):
70 return 0
71
72
73 ifos = map(str,trig.get_ifos())
74
75
76 sens = {}
77 fPlus,fCross = get_det_response(numpy.degrees(trig.ra),\
78 numpy.degrees(trig.dec),\
79 trig.get_end())
80 for ifo in ifos:
81 if ifo.lower()[0] == 'h' :
82 i = ifo.lower()
83 else:
84 i = ifo[0].lower()
85 sens[ifo] = getattr(trig, 'sigmasq_%s' % i.lower()) * \
86 sum(numpy.array([fPlus[ifo], fCross[ifo]])**2)
87 ifos.sort(key=lambda ifo: sens[ifo], reverse=True)
88 for i in xrange(0,2):
89 if ifos[i].lower()[0] == 'h' :
90 i = ifos[i].lower()
91 else:
92 i = ifos[i][0].lower()
93 if getattr(trig, 'snr_%s' % i) < sngl_snr_threshold:
94 return 0
95
96
97 bestNR = trig.get_bestnr(index=1, nhigh=n, \
98 null_snr_threshold=null_thresh[0], \
99 null_grad_thresh=null_grad_thresh, null_grad_val = null_grad_val)
100
101
102
103 if trig.chisq == 0:
104
105 print >> sys.stderr,\
106 "Chisq not calculated for trigger with end time and snr:"
107 print >> sys.stderr, trig.get_end(),trig.snr
108 raise ValueError("Chisq has not been calculated for trigger.")
109
110 return bestNR
111
112 -def calculate_contours(q=4.0, n=3.0, null_thresh=6., null_grad_snr=20,\
113 new_snr_thresh=6.0, new_snrs=[5.5,6,6.5,7,8,9,10,11],\
114 null_grad_val = 0.2, chisq_dof = 60,\
115 bank_chisq_dof = 40, cont_chisq_dof = 160):
116 """Generate the plot contours for chisq variable plots
117 """
118
119 num_vals = len(new_snrs)
120 colors = ["k-" if snr == new_snr_thresh else
121 "y-" if snr == int(snr) else
122 "y--" for snr in new_snrs]
123
124
125 snr_low_vals = numpy.arange(4,30,0.1)
126 snr_high_vals = numpy.arange(30,500,1)
127 snr_vals = numpy.asarray(list(snr_low_vals) + list(snr_high_vals))
128
129
130 bank_conts = numpy.zeros([len(new_snrs),len(snr_vals)],
131 dtype=numpy.float64)
132 auto_conts = numpy.zeros([len(new_snrs),len(snr_vals)],
133 dtype=numpy.float64)
134 chi_conts = numpy.zeros([len(new_snrs),len(snr_vals)],
135 dtype=numpy.float64)
136 null_cont = []
137
138
139 for j,snr in enumerate(snr_vals):
140 for i,new_snr in enumerate(new_snrs):
141 bank_conts[i][j] = new_snr_chisq(snr, new_snr, bank_chisq_dof, q, n)
142 auto_conts[i][j] = new_snr_chisq(snr, new_snr, cont_chisq_dof, q, n)
143 chi_conts[i][j] = new_snr_chisq(snr, new_snr, chisq_dof, q, n)
144
145 if snr > null_grad_snr:
146 null_cont.append(null_thresh + (snr-null_grad_snr)*null_grad_val)
147 else:
148 null_cont.append(null_thresh)
149 null_cont = numpy.asarray(null_cont)
150
151 return bank_conts,auto_conts,chi_conts,null_cont,snr_vals,colors
152
153
155 for i in range(len(contours)):
156 plot_vals_x = []
157 plot_vals_y = []
158 for j in range(len(snr_vals)):
159 if contours[i][j] > 1E-15:
160
161
162 if not plot_vals_x:
163 plot_vals_x.append(snr_vals[j])
164 plot_vals_y.append(0.1)
165 plot_vals_x.append(snr_vals[j])
166 plot_vals_y.append(contours[i][j])
167 axis.plot(plot_vals_x,plot_vals_y,colors[i])
168
169
171 times = {}
172 for name,fileName in\
173 zip(["buffer", "off", "on"],\
174 ["bufferSeg.txt","offSourceSeg.txt","onSourceSeg.txt"]):
175
176 segs = segmentsUtils.fromsegwizard(open(os.path.join(segdir,fileName),
177 'r'))
178 if len(segs)>1:
179 raise AttributeError, 'More than one segment, an error has occured.'
180 times[name] = segs[0]
181 return times
182
183
185 import pylab
186 pylab.rcParams.update({
187 "text.usetex": True,
188 "text.verticalalignment": "center",
189
190
191
192 "figure.figsize": [8.0, 6.0],
193 "font.size": 20,
194 "axes.titlesize": 16,
195 "axes.labelsize": 24,
196 "xtick.labelsize": 18,
197 "ytick.labelsize": 18,
198 "legend.fontsize": 18,
199 "font.family": "serif",
200 "font.weight": "bold",
201 })
202
203
205 """Return detector response for complete set of IFOs for given sky
206 location and time. Inclination and polarization are unused so are
207 arbitrarily set to 0.
208 """
209 f_plus = {}
210 f_cross = {}
211 inclination = 0
212 polarization = 0
213 for ifo in ['G1','H1','H2','L1','T1','V1']:
214 f_plus[ifo],f_cross[ifo],_,_ = antenna.response(trigTime, ra, dec,\
215 inclination,
216 polarization, 'degree',
217 ifo)
218 return f_plus,f_cross
219
220
222 """FIXME
223 """
224 if re.search('SimInspiral', str(self)):
225 ra = numpy.degrees(self.longitude)
226 dec = numpy.degrees(self.latitude)
227 t = self.get_time_geocent()
228 else:
229 ra = numpy.degrees(self.ra)
230 dec = numpy.degrees(self.dec)
231 t = self.get_end()
232
233 fplus, fcross = get_det_response(ra, dec, t)
234 return dict((ifo, fplus[ifo]**2 + fcross[ifo]**2) for ifo in fplus.keys())
235
236
238 """Construct and append process and process_params tables to
239 ligolw.Document object xmldoc, using the given sys.argv variable
240 args and other parameters.
241 """
242 xmldoc.childNodes[-1].appendChild(lsctables.New(lsctables.ProcessTable))
243 xmldoc.childNodes[-1].appendChild(
244 lsctables.New(lsctables.ProcessParamsTable))
245
246
247 progName = args[0]
248 process = ligolw_process.append_process(xmldoc, program=progName,
249 version=version, cvs_repository='lscsoft',
250 cvs_entry_time=date)
251 params = []
252 for i in range(len(args)):
253 p = args[i]
254 if not p.startswith('-'):
255 continue
256 v = ''
257 if i < len(sys.argv)-1:
258 v = sys.argv[i+1]
259 params.append( map( unicode, (p,'string',v) ) )
260
261 ligolw_process.append_process_params(xmldoc,process,params)
262
263 return xmldoc
264
265
267 files = glob.glob(os.path.join(log_dir,"*err"))
268
269 badInjs = []
270
271 for file in files:
272 if os.stat(file)[6] != 0:
273 fp = open(file,"r")
274 conts = fp.read()
275 if conts.find('terminated') != -1:
276 conts=conts.split('\n')
277 for line in conts:
278 line = line.split(' ')
279 line = [entry.replace(',','') for entry in line if entry]
280 if 'terminated' in line:
281 injDict = {}
282 injDict['mass1'] = float(line[6])
283 injDict['mass2'] = float(line[8])
284 injDict['spin1x'] = float(line[10])
285 injDict['spin1y'] = float(line[12])
286 injDict['spin1z'] = float(line[14])
287 injDict['spin2x'] = float(line[16])
288 injDict['spin2y'] = float(line[18])
289 injDict['spin2z'] = float(line[20])
290 if not injDict in badInjs:
291 badInjs.append(injDict)
292 return badInjs
293
294
296 new_sims = []
297 for sim in sims:
298 for badInj in badInjs:
299 if ((abs(sim.mass1-badInj['mass1']) < 0.001) and
300 (abs(sim.mass2-badInj['mass2']) < 0.001) and
301 (abs(sim.spin1x-badInj['spin1x']) < 0.001) and
302 (abs(sim.spin1y-badInj['spin1y']) < 0.001) and
303 (abs(sim.spin1z-badInj['spin1z']) < 0.001) and
304 (abs(sim.spin2x-badInj['spin2x']) < 0.001) and
305 (abs(sim.spin2y-badInj['spin2y']) < 0.001) and
306 (abs(sim.spin2z-badInj['spin2z']) < 0.001)):
307 print("Removing injection:", sim.mass1, sim.mass2, sim.spin1x,
308 sim.spin1y, sim.spin1z, sim.spin2x, sim.spin2y,
309 sim.spin2z)
310 break
311 else:
312 new_sims.append(sim)
313
314 return new_sims
315
316
318
319 angmomfac = self.mass1 * self.mass2 * \
320 numpy.power(LAL_PI * LAL_MTSUN_SI * (self.mass1 + self.mass2) *
321 self.f_lower, -1.0/3.0)
322 m1sq = self.mass1 * self.mass1
323 m2sq = self.mass2 * self.mass2
324
325
326 L = numpy.zeros(3)
327 L[0] = angmomfac * numpy.sin(self.inclination)
328 L[1] = 0
329 L[2] = angmomfac * numpy.cos(self.inclination)
330
331
332 S = numpy.zeros(3)
333 S[0] = m1sq * self.spin1x + m2sq * self.spin2x
334 S[1] = m1sq * self.spin1y + m2sq * self.spin2y
335 S[2] = m1sq * self.spin1z + m2sq * self.spin2z
336
337
338 J = L + S
339
340 theta = math.atan2(math.sqrt(J[0]*J[0] + J[1]*J[1]),J[2])
341 if theta > math.pi/2.:
342 theta = math.pi - theta
343
344 if theta < 0 or theta > math.pi/2.:
345 raise Error("Theta is too big or too small")
346
347 return theta
348
349
351 """Veto events in a MultiInspiralTable based on their (coherent) SNR
352 value.
353
354 @param mi_table
355 a MultiInspiralTable from which to veto events
356 @param snr
357 the value of coherent SNR on which to threshold
358 @param return_index
359 boolean to return the index array of non-vetoed elements rather
360 than a new table containing the elements themselves
361
362 @returns
363 a new MultiInspiralTable with those events not vetoed OR
364 the indices of the original mi_table not vetoed if return_index=True
365 """
366 mi_snr = numpy.asarray(mi_table.get_column("snr"))
367 keep = mi_snr >= snr
368 if return_index:
369 return keep
370 else:
371 out = table.new_from_template(mi_table)
372 out.extend(numpy.asarray(mi_table)[keep])
373 return out
374
375
377 """Veto events in a MultiInspiralTable based on their \f$\chi^2\f$
378 re-weighted coherent SNR.
379
380 @param mi_table
381 a MultiInspiralTable from which to veto events
382 @param snr
383 the value of coherent new SNR on which to threshold
384 @param chisq_index
385 the index \f$\iota\f$ used in the newSNR calculation:
386 \f[\rho_{\mbox{new}} =
387 \frac{\rho}{\left[\frac{1}{2}
388 \left(1 + \left(\frac{\chi^2}{n_\mbox{dof}}\right)^{\iota/3}
389 \right)\right]^{1/\iota}}
390 \f]
391 @param return_index
392 boolean to return the index array of non-vetoed elements rather
393 than a new table containing the elements themselves
394 """
395 new_snr = numpy.asarray(mi_table.get_new_snr(column="chisq"))
396 keep = new_snr >= snr
397 if return_index:
398 return keep
399 else:
400 out = table.new_from_template(mi_table)
401 out.extend(numpy.asarray(mi_table)[keep])
402 return out
403
404
405 -def apply_bank_veto(mi_table, snr=6.0, chisq_index=4.0, return_index=False):
406 """Veto events in a MultiInspiralTable based on their bank chisq-
407 weighted (new) coherent SNR.
408
409 @param mi_table
410 a MultiInspiralTable from which to veto events
411 @param snr
412 the value of coherent new SNR on which to threshold
413 @param chisq_index
414 the index \f$\iota\f$ used in the newSNR calculation:
415 \f[\rho_{\mbox{new}} =
416 \frac{\rho}{\left[\frac{1}{2}
417 \left(1 + \left(\frac{\chi^2}{n_\mbox{dof}}\right)^{\iota/3}
418 \right)\right]^{1/\iota}}
419 \f]
420 @param return_index
421 boolean to return the index array of non-vetoed elements rather
422 than a new table containing the elements themselves
423
424 @returns
425 a new MultiInspiralTable with those events not vetoed OR
426 the indices of the original mi_table not vetoed if return_index=True
427 """
428 bank_new_snr = numpy.asarray(mi_table.get_new_snr(column="bank_chisq"))
429 keep = bank_new_snr >= snr
430 if return_index:
431 return keep
432 else:
433 out = table.new_from_template(mi_table)
434 out.extend(numpy.asarray(mi_table)[keep])
435 return out
436
437
438 -def apply_auto_veto(mi_table, snr=6.0, chisq_index=4.0, return_index=False):
439 """Veto events in a MultiInspiralTable based on their auto chisq-
440 weighted (new) coherent SNR.
441
442 @param mi_table
443 a MultiInspiralTable from which to veto events
444 @param snr
445 the value of coherent new SNR on which to threshold
446 @param chisq_index
447 the index \f$\iota\f$ used in the newSNR calculation:
448 \f[\rho_{\mbox{new}} =
449 \frac{\rho}{\left[\frac{1}{2}
450 \left(1 + \left(\frac{\chi^2}{n_\mbox{dof}}\right)^{\iota/3}
451 \right)\right]^{1/\iota}}
452 \f]
453 @param return_index
454 boolean to return the index array of non-vetoed elements rather
455 than a new table containing the elements themselves
456
457 @returns
458 a new MultiInspiralTable with those events not vetoed OR
459 the indices of the original mi_table not vetoed if return_index=True
460 """
461 cont_new_snr = numpy.asarray(mi_table.get_new_snr(column="cont_chisq"))
462 keep = cont_new_snr >= snr
463 if return_index:
464 return keep
465 else:
466 out = table.new_from_template(mi_table)
467 out.extend(numpy.asarray(mi_table)[keep])
468 return out
469
470
472 """Veto events in a MultiInspiralTable based on their single-detector
473 snr in the most sensitive detectors.
474
475 @param mi_table
476 a MultiInspiralTable from which to veto events
477 @param snrs
478 an X-element list of single-detector SNRs on which to threshold
479 for the X most sensitive detectors (in sensitivity order)
480 @param return_index
481 boolean to return the index array of non-vetoed elements rather
482 than a new table containing the elements themselves
483
484 @returns
485 a new MultiInspiralTable with those events not vetoed OR
486 the indices of the original mi_table not vetoed if return_index=True
487 """
488 if len(mi_table) == 0:
489 if return_index:
490 return numpy.zeros(0).astype(bool)
491 else:
492 return mi_table
493
494 ifos = lsctables.instrument_set_from_ifos(mi_table[0].ifos)
495 if "H1" in ifos and "H2" in ifos:
496 ifos.remove("H2")
497 mi_sngl_snr = numpy.asarray([numpy.asarray(mi_table.get_sngl_snr(ifo))
498 for ifo in ifos])
499 if len(snrs) > len(ifos):
500 raise ValueError("%s single-detector thresholds given, but only %d "
501 "detectors found." % (len(snrs), len(ifos)))
502
503 snr_array = numpy.zeros(len(ifos))
504 snr_array[:len(snrs)] = snrs
505 snr_array.sort()
506
507
508 mi_sngl_snr.sort(axis=0)
509
510
511 keep = (mi_sngl_snr.T > snr_array).all(axis=1)
512 if return_index:
513 return keep
514 else:
515 out = table.new_from_template(mi_table)
516 out.extend(numpy.asarray(mi_table)[keep])
517 return out
518
519
522 """Veto events in a MultiInspiralTable based on their null SNR.
523
524 @param mi_table
525 a MultiInspiralTable from which to veto events
526 @param null_snr
527 the value of null SNR on which to threshold
528 @param snr
529 the value of coherent SNR on above which to grade the null SNR
530 threshold
531 @param grad
532 the rate at which to increase the null SNR threshold above snr
533 @param return_index
534 boolean to return the index array of non-vetoed elements rather
535 than a new table containing the elements themselves
536
537 @returns
538 a new MultiInspiralTable with those events not vetoed OR
539 the indices of the original mi_table not vetoed if return_index=True
540 """
541 mi_snr = mi_table.get_column("snr")
542 mi_null_snr = mi_table.get_null_snr()
543
544 null_thresh = numpy.ones(len(mi_table)) * null_snr
545 grade = mi_snr >= snr
546 null_thresh[grade] += (mi_snr[grade] - snr)*5.0
547
548 keep = mi_null_snr < null_thresh
549 if return_index:
550 return keep
551 else:
552 out = table.new_from_template(mi_table)
553 out.extend(numpy.asarray(mi_table)[keep])
554 return out
555
556 -def veto(self, seglist, time_slide_table=None):
557 """Return a MultiInspiralTable with those row from self not lying
558 inside (i.e. not vetoed by) any elements of seglist.
559
560 If time_slide_table is not None, any time shifts will be undone and each
561 detector checked individually
562 """
563 seglist = type(seglist)([type(seg)(*map(float, seg)) for seg in seglist])
564 keep = table.new_from_template(self)
565 if time_slide_table:
566 slides = time_slide_table.as_dict()
567 for id_,vector in slides.iteritems():
568 idx = str(id_).split(":")[-1]
569 slides["multi_inspiral:time_slide_id:%s" % idx] = vector
570 del slides[id_]
571 for row in self:
572 ifos = row.get_ifos()
573 for i,ifo in enumerate(ifos):
574 ifo_time = float(row.get_end() +
575 slides[str(row.time_slide_id)][ifo])
576 if ifo_time in seglist:
577 i = -1
578 break
579 if i != -1:
580 keep.append(row)
581 else:
582 for row in self:
583 time = float(row.get_end())
584 if time in seglist:
585 continue
586 else:
587 keep.append(row)
588 return keep
589
590
591 -def vetoed(self, seglist, time_slide_table=None):
592 """Return a MultiInspiralTable with those row from self lying
593 inside (i.e. vetoed by) any elements of seglist.
594
595 If time_slide_table is not None, any time shifts will be undone and each
596 detector checked individually
597 """
598 seglist = type(seglist)(map(float, seglist))
599 vetoed = table.new_from_template(self)
600 if time_slide_table:
601 slides = time_slide_table.as_dict()
602 for id_,vector in slides.iteritems():
603 idx = str(id_).split(":")[-1]
604 slides["multi_inspiral:time_slide_id:%s" % idx] = vector
605 del slides[id_]
606 slides = get_slide_vectors(time_slide_table)
607 for row in self:
608 ifos = row.get_ifos()
609 for i,ifo in enumerate(ifos):
610 ifo_time = (float(row.get_end()) -
611 slides[str(row.time_slide_id)][ifo])
612 if ifo_time in seglist:
613 vetoed.append(row)
614 break
615 else:
616 for row in self:
617 time = float(row.get_end())
618 if time in seglist:
619 vetoed.append(row)
620 return vetoed
621