1 import sys
2 from glue import segments
3 from glue.ligolw import ligolw
4 from glue.ligolw import table
5 from glue.ligolw import lsctables
6 from glue.ligolw import utils
7 from pylal.tools import XLALCalculateEThincaParameter
8 from pylal.xlal import date
9 from pylal.xlal.datatypes.ligotimegps import LIGOTimeGPS
10 import glue.iterutils
11 import numpy
12 import cmath
13
14
15
16 ifos = ("G1", "H1", "H2", "L1", "T1", "V1")
17
18
19
20
21 -class ExtractCoincInspiralTableLIGOLWContentHandler(ligolw.PartialLIGOLWContentHandler):
22 """
23 LIGOLWContentHandler that will extract only the CoincInspiralTable from a document.
24 See glue.ligolw.LIGOLWContentHandler help for more info.
25 """
26 - def __init__(self,document):
27 def filterfunc(name,attrs):
28 return name==ligolw.Table.tagName and attrs.has_key('Name') and table.Table.TableName(attrs.get('Name'))==lsctables.CoincInspiralTable.tableName
29 ligolw.PartialLIGOLWContentHandler.__init__(self,document,filterfunc)
30
31
33 ifo_combos = []
34 for num_ifos in range(2, len(ifo_list) + 1):
35 ifo_combos.extend(list(glue.iterutils.choices(ifo_list, num_ifos)))
36
37 return ifo_combos
38
40 """
41 Return the e-thinca corresponding to the distance in parameter space between two inspiral triggers.
42
43 The average distance defined below is only an approximation to the true distance and is
44 valid whenever two triggers are nearby. The simplified version of
45 the e-thinca parameter is calculated based on that definition of the distance.
46
47 d_average=(1/2)[(Gamma(x1)_{ij}(x2-x1)^i(x2-x1)^j)^(1/2) + (Gamma(x2)_{ij}(x2-x1)^i(x2-x1)^j)^(1/2)]
48 then simple_ethinca= d_average^2/4
49
50 @param trigger1: is a single inspiral triggers.
51 @param trigger2: is a single inspiral triggers.
52 """
53 dt = 0
54 dtau0 = trigger2.tau0-trigger1.tau0
55 dtau3 = trigger2.tau3-trigger1.tau3
56
57 dist1 = dt * (dt * trigger1.Gamma0 + dtau0 * trigger1.Gamma1 + dtau3 * trigger1.Gamma2) + \
58 dtau0 * (dt * trigger1.Gamma1 + dtau0 * trigger1.Gamma3 + dtau3 * trigger1.Gamma4) + \
59 dtau3 * (dt * trigger1.Gamma2 + dtau0 * trigger1.Gamma4 + dtau3 * trigger1.Gamma5)
60
61 dist2 = dt * (dt * trigger2.Gamma0 + dtau0 * trigger2.Gamma1 + dtau3 * trigger2.Gamma2) + \
62 dtau0 * (dt * trigger2.Gamma1 + dtau0 * trigger2.Gamma3 + dtau3 * trigger2.Gamma4) + \
63 dtau3 * (dt * trigger2.Gamma2 + dtau0 * trigger2.Gamma4 + dtau3 * trigger2.Gamma5)
64
65 average_distance = 0.5 * (numpy.sqrt(dist1) + numpy.sqrt(dist2))
66
67 simple_ethinca = (average_distance * average_distance) / 4.0
68
69 return simple_ethinca
70
72 """
73 read in the Sngl and SimInspiralTables from a list of files
74 if Sngls are found, construct coincs, add injections (if any)
75 also return Sims (if any)
76 @param fileList: list of input files
77 @param statistic: instance of coincStatistic, to use in creating coincs
78 """
79 if not fileList:
80 return coincInspiralTable(), None
81
82 if not (isinstance(statistic,coincStatistic)):
83 raise TypeError, "invalid statistic, must be coincStatistic"
84
85 sims = None
86 coincs = None
87
88 lsctables.use_in(ExtractCoincInspiralTableLIGOLWContentHandler)
89 for thisFile in fileList:
90 doc = utils.load_filename(thisFile, gz = (thisFile or "stdin").endswith(".gz"), contenthandler=ExtractCoincInspiralTableLIGOLWContentHandler)
91
92 try:
93 simInspiralTable = lsctables.SimInspiralTable.get_table(doc)
94 if sims: sims.extend(simInspiralTable)
95 else: sims = simInspiralTable
96 except: simInspiralTable = None
97
98
99 try: snglInspiralTable = lsctables.SnglInspiralTable.get_table(doc)
100 except: snglInspiralTable = None
101 if snglInspiralTable:
102 coincFromFile = coincInspiralTable(snglInspiralTable,statistic)
103 if simInspiralTable:
104 coincFromFile.add_sim_inspirals(simInspiralTable)
105 if coincs: coincs.extend(coincFromFile)
106 else: coincs = coincFromFile
107
108 doc.unlink()
109
110 return coincs, sims
111
112
113
115 """
116 This class specifies the statistic to be used when dealing with coincident events.
117 It also contains parameter for such cases as the BBH bitten-L statistics.
118 """
119
120 __slots__ = ["name","a","b","rsq","bl","eff_snr_denom_fac","new_snr_index"]
121
122 - def __init__(self, name, a=0, b=0, denom_fac=250.0, chisq_index=6.0):
123 self.name=name
124 self.a=a
125 self.b=b
126 self.rsq=0
127 self.bl=0
128 self.eff_snr_denom_fac = denom_fac
129 self.new_snr_index = chisq_index
130
133
135 blx=self.a*snr-self.b
136 if bl==0:
137 return blx
138 else:
139 return min(bl, blx)
140
141
142
144 """
145 Table to hold coincident inspiral triggers. Coincidences are reconstructed
146 by making use of the event_id contained in the sngl_inspiral table.
147 The coinc is a dictionary with entries: event_id, numifos, stat, and
148 each available IFO (G1, H1, etc.).
149 The stat is set by default to the snrsq: the sum of the squares of the snrs
150 of the individual triggers.
151 """
153 __slots__ = ["event_id", "numifos", "stat", "likelihood",
154 "sim", "rsq", "bl", "fap"] + list(ifos)
155
156 - def __init__(self, event_id, numifos = 0, stat = 0, likelihood = 0):
157 self.event_id = event_id
158 self.numifos = numifos
159 self.stat = stat
160 self.likelihood = likelihood
161 self.rsq=0
162 self.bl=0
163 self.fap = 0.0
164
166
167
168
169
170
171 assert not hasattr(self, trig.ifo), "Trying to add %s trigger to a"\
172 " coincidence for the second time. Coincidence so far:\n%s"\
173 "\n\nTrigger:\n%s" % (trig.ifo, dict([(x, getattr(self, x)) for x in \
174 self.__slots__ if hasattr(self, x)]), trig.event_id)
175
176 self.numifos +=1
177 if statistic.name == 'effective_snr':
178 self.stat = (self.stat**2 + trig.get_effective_snr(statistic.eff_snr_denom_fac)**2)**(1./2)
179 elif statistic.name == 'new_snr':
180 self.stat = (self.stat**2 + trig.get_new_snr(statistic.new_snr_index)**2)**0.5
181 elif 'bitten_l' in statistic.name:
182 snr=trig.snr
183 self.rsq= (self.rsq**2 + snr**2)**(1./2)
184 self.bl=statistic.get_bittenl( self.bl, snr )
185 self.stat=min( self.bl, self.rsq )
186 if statistic.name == 'bitten_lsq' and self.numifos >2:
187 self.stat = self.rsq
188 elif 'far' == statistic.name:
189 self.stat = trig.get_far()
190 elif 'ifar' == statistic.name:
191 self.stat = trig.get_ifar()
192 elif 'lvS5stat' == statistic.name:
193 self.stat = trig.get_lvS5stat()
194 else:
195 self.stat = (self.stat**2 + getattr(trig,statistic.name)**2)**(1./2)
196
197
198 setattr(self,trig.ifo,trig)
199
201 setattr(self,"sim",sim)
202
205 ifos = property(fget=_get_ifo_set)
206
208 return lsctables.ifos_from_instrument_set(self.ifos)
209 ifostring = property(fget=_get_ifo_string)
210
212 ifo_string = ""
213 ifolist_in_coinc = []
214 for ifo in ifos:
215 if hasattr(self,ifo):
216 ifo_string = ifo_string + ifo
217 ifolist_in_coinc.append(ifo)
218
219 return ifo_string, ifolist_in_coinc
220
225 slide_num = property(fget=_get_slide_num)
226
228 for ifo in ifos:
229 if hasattr(self,ifo):
230 return getattr(self, ifo).end_time+getattr(self, ifo).end_time_ns*1.0e-9
231 raise ValueError, "This coincident trigger does not contain any "\
232 "single trigger. This should never happen."
233
235 """
236 Return a dictionary of the GPS times associated with each
237 trigger in the coincidence. The time is stored in LIGOTimeGPS format.
238 """
239 gpstimes={}
240 for ifo in ifos:
241 if hasattr(self,ifo):
242 gpstimes[ifo]= LIGOTimeGPS(getattr(self, ifo).end_time, getattr(self, ifo).end_time_ns)
243 return gpstimes
244
246 """
247 Return an iterator over the triggers in this coinc.
248 """
249 for ifo in ifos:
250 if hasattr(self, ifo):
251 yield getattr(self, ifo)
252
253 - def __init__(self, inspTriggers = None, stat = None):
254 """
255 @param inspTriggers: a metaDataTable containing inspiral triggers
256 from which to construct coincidences
257 @param stat: an instance of coincStatistic
258 """
259 self.stat = stat
260 self.sngl_table = inspTriggers
261 self.sim_table = None
262 self.rows = []
263 if inspTriggers is None:
264 return
265
266
267
268 row_dict = {}
269 unique_id_list = []
270 for trig in inspTriggers:
271 event_id = trig.event_id
272 if event_id not in row_dict:
273 unique_id_list.append(event_id)
274 row_dict[event_id] = self.row(event_id)
275 row_dict[event_id].add_trig(trig, stat)
276
277
278 pruned_rows = [row_dict[k] for k in unique_id_list \
279 if row_dict[k].numifos > 1]
280
281 self.rows = pruned_rows
282
284 return len(self.rows)
285
288
291
293 """
294 Retrieve the value in this column in row i.
295 """
296 return self.rows[i]
297
299 return numpy.array([c.stat for c in self.rows], dtype=float)
300
301 - def sort(self, descending = True):
302 """
303 Sort the list based on stat value
304 default is to descending
305 """
306 stat_list = [ (coinc.stat, coinc) for coinc in self.rows ]
307 stat_list.sort()
308 if descending:
309 stat_list.reverse()
310 self.rows = [coinc for (stat,coinc) in stat_list]
311
313 """
314 Calculates false alarm probability for each coinc using stats array.
315 @param stats: array of loudest statistics forund in each of the time slides
316 """
317
318 for coinc in self:
319 if not use_likelihood:
320 index = numpy.searchsorted(stats, coinc.stat)
321 if index == 0:
322 coinc.fap = 100.0
323 elif index == len(stats):
324 coinc.fap = 0.0
325 else:
326 coinc.fap = 100.0 * float((len(stats) - index)) /float(len(stats))
327
329 """
330 Return the triggers with a specific slide number.
331 @param slide_num: the slide number to recover (contained in the event_id)
332 """
333 slide_coincs = coincInspiralTable(stat=self.stat)
334 slide_coincs.sngl_table = self.sngl_table
335 slide_coincs.extend([c for c in self.rows if c.slide_num == slide_num])
336 return slide_coincs
337
339 """
340 Return the coincs which have triggers from the ifos in ifolist.
341 @param ifolist: a list of ifos
342 """
343 selected_coincs = coincInspiralTable(stat=self.stat)
344 selected_coincs.sngl_table = self.sngl_table
345 for coinc in self:
346 keep_trig = True
347 for ifo in ifolist:
348 if hasattr(coinc,ifo) == False:
349 keep_trig = False
350 break
351
352 if keep_trig == True:
353 selected_coincs.append(coinc)
354
355 return selected_coincs
356
358 """
359 Return the coincs which are from ifos.
360 @param ifolist: a list of ifos
361 """
362 coincs = self.coincinclude(ifolist)
363 selected_coincs = coincInspiralTable(stat=self.stat)
364 selected_coincs.sngl_table = self.sngl_table
365 for coinc in coincs:
366 if coinc.numifos == len(ifolist):
367 selected_coincs.append(coinc)
368
369 return selected_coincs
370
372 """
373 Return the coincs which are NOT from the coinc type made from combining all
374 the ifos in the ifolist.
375 @param ifolist: a list of ifos
376 """
377 ifolist.sort()
378 selected_coincs = coincInspiralTable(stat=self.stat)
379 selected_coincs.sngl_table = self.sngl_table
380 for coinc in self:
381 thiscoinctype = coinc.get_ifos()[1]
382 thiscoinctype.sort()
383 if not (ifolist == thiscoinctype):
384 selected_coincs.append(coinc)
385
386 return selected_coincs
387
389 """
390 Return the sngls for a specific ifo.
391 @param ifo: ifo for which to retrieve the single inspirals.
392 """
393 from glue.ligolw import table
394 try: ifoTrigs = table.new_from_template(self.sngl_table)
395 except: ifoTrigs = lsctables.New(lsctables.SnglInspiralTable)
396 for coinc in self:
397 if hasattr(coinc,ifo):
398 ifoTrigs.append(getattr(coinc,ifo))
399
400 return ifoTrigs
401
403 """
404 Returns a list of coincident triggers that are vetoed
405 entirely by a segment of the segment list.
406 A coincident trigger is added to the list only
407 if all triggers lie within a segment. If any
408 trigger lies outside any segment, it is not added.
409 @param seglist: segment list used to veto coincidences
410 """
411 vetoed_coincs = coincInspiralTable(stat=self.stat)
412
413
414 for coinc in self:
415 flagVeto = True
416 for ifo in ifos:
417 if hasattr(coinc, ifo):
418
419 if getattr(coinc,ifo).get_end() not in seglist:
420 flagVeto = False
421 break
422
423
424
425 if flagVeto:
426 vetoed_coincs.append( coinc )
427
428 return vetoed_coincs
429
431 """
432 Return the clustered triggers, returning the one with the largest stat in
433 each fixed cluster_window, exactly as in lalapps_coire
434 (in package/tools/CoincInspiralUtils.c:XLALClusterCoincInspiralTable)
435
436 @param cluster_window: length of time over which to cluster (seconds)
437 """
438
439
440 stat_list = [ (coinc.get_time(), coinc) for coinc in self.rows ]
441 stat_list.sort()
442 self.rows = [coinc for (t,coinc) in stat_list]
443
444
445
446 thisCoinc = 0
447 nextCoinc = 1
448 while nextCoinc<len(self.rows):
449
450
451 thisTime = self.rows[thisCoinc].get_time()
452 nextTime = self.rows[nextCoinc].get_time()
453
454
455 if nextTime-thisTime<cluster_window:
456
457
458 thisStat = self.rows[thisCoinc].stat
459 nextStat = self.rows[nextCoinc].stat
460
461
462 if (nextStat>thisStat):
463 self.rows.pop(thisCoinc)
464 else:
465 self.rows.pop(nextCoinc)
466
467 else:
468
469
470 thisCoinc+=1
471 nextCoinc+=1
472
473
474 - def cluster(self, cluster_window, numSlides = None):
475 """
476 New clustering method working the same way as the clustering method
477 used in lalapps_coire
478 (in package/tools/CoincInspiralUtils.c:XLALClusterCoincInspiralTable)
479
480 @param cluster_window: length of time over which to cluster (seconds)
481 @param numSlides: number of slides if this are time-slide triggers
482 """
483
484 if not numSlides:
485
486
487 self.cluster_core(cluster_window)
488
489 else:
490
491 cluster = coincInspiralTable()
492
493
494 for slide in range(-numSlides, numSlides+1):
495
496
497 slideCoinc = self.getslide( slide )
498
499
500 slideCoinc.cluster_core( cluster_window )
501
502
503 cluster.extend( slideCoinc )
504
505
506 self.rows = cluster.rows
507
508
509 return self
510
511
513 """
514 FIXME: We should really store the sim coincidence info in the event_id
515 Method to add simulated injections to a list of coincs
516
517 @param sim_inspiral: a simInspiralTable
518 """
519 self.sim_table = sim_inspiral
520
521 if len(self) != len(sim_inspiral):
522 raise ValueError, "Number of injections doesn't match number of coincs"
523
524 for i in range(len(self)):
525 self[i].add_sim(sim_inspiral[i])
526
527
529 """
530 Add missed sim inspirals to the list of coincs, set the stat = -1
531 @param sim_inspiral: a simInspiralTable
532 """
533 for sim in sim_inspiral:
534 row = coincInspiralTable.row(-1,stat=-1)
535 row.add_sim(sim)
536 self.append(row)
537
539 """
540 Method to return the sim_inspiral table associated to the coincs.
541 If thresh is specified, only return sims from those coincs whose stat
542 exceeds thresh (or is under thresh if statistic == far).
543
544 @param statistic: the statistic to use
545 @param thresh: the threshold on the statistic
546 """
547 from glue.ligolw import table
548 try: simInspirals = table.new_from_template(self.sim_table)
549 except: simInspirals = lsctables.New(lsctables.SimInspiralTable)
550 for coinc in self:
551 if statistic == 'far':
552 if (hasattr(coinc,"sim")) and \
553 (((coinc.stat <= thresh) and (coinc.stat >= 0)) or (thresh < 0)):
554 simInspirals.append(coinc.sim)
555 else:
556 if (hasattr(coinc,"sim")) and (coinc.stat >= thresh):
557 simInspirals.append(coinc.sim)
558
559 return simInspirals
560
562 """
563 Return (triggers with stat < threshold,
564 triggers with stat == threshold,
565 triggers with stat > threshold).
566
567 The set of (stat == threshold) is of zero measure, but often, as when
568 doing something like the loudest event statistic, the threshold is taken
569 from a coinc in self.
570 """
571 stats = self.getstat()
572
573 lesser_coincs = coincInspiralTable(stat=self.stat)
574 lesser_coincs.extend([self[i] for i in (stats < threshold).nonzero()[0]])
575
576 equal_coincs = coincInspiralTable(stat=self.stat)
577 equal_coincs.extend([self[i] for i in (stats == threshold).nonzero()[0]])
578
579 greater_coincs = coincInspiralTable(stat=self.stat)
580 greater_coincs.extend([self[i] for i in (stats > threshold).nonzero()[0]])
581
582 return lesser_coincs, equal_coincs, greater_coincs
583
585 """
586 Return triggers with mLow <= mean total mass < mHigh
587 @param mLow: a float
588 @param mHigh: a float
589 """
590 triggers_in_mass_range = coincInspiralTable()
591 for coinc in self:
592 mass_numer = 0.0
593 mass_denom = float(coinc.numifos)
594 for ifo in ifos:
595 if hasattr(coinc,ifo):
596 mass_numer += getattr(coinc,ifo).mass1 + getattr(coinc,ifo).mass2
597 mean_total_mass = mass_numer / mass_denom
598 if (mean_total_mass >= mLow) and (mean_total_mass < mHigh):
599 triggers_in_mass_range.append(coinc)
600
601 return triggers_in_mass_range
602
604 """
605 Return triggers with mLow <= mean chirp mass < mHigh
606 @param mLow: a float
607 @param mHigh: a float
608 """
609 triggers_in_mass_range = coincInspiralTable()
610 for coinc in self:
611 mass_numer = 0.0
612 mass_denom = float(coinc.numifos)
613 for ifo in ifos:
614 if hasattr(coinc,ifo):
615 mass_numer += getattr(coinc,ifo).mchirp
616 mean_total_mass = mass_numer / mass_denom
617 if (mean_total_mass >= mLow) and (mean_total_mass < mHigh):
618 triggers_in_mass_range.append(coinc)
619
620 return triggers_in_mass_range
621
623 """
624 Return ethinca values for the coincidences
625 @param ifos: a list of the 2 ifos
626 """
627 ethinca = numpy.zeros(len(self), dtype=float)
628 for i,coinc in enumerate(self):
629 if hasattr(coinc, ifos[0]) and hasattr(coinc, ifos[1]):
630 try:
631 ethinca[i] = XLALCalculateEThincaParameter(getattr(coinc, ifos[0]),
632 getattr(coinc, ifos[1]))
633 except ValueError:
634 ethinca[i] = 2.0
635
636 return ethinca
637
639 """
640 Return simple ethinca values for the coincidences of the ifo pair ifos.
641
642 For coincs that do not have both ifos specified, return 0.0.
643 """
644 ethinca = numpy.zeros(len(self), dtype=float)
645 for i,coinc in enumerate(self):
646 if hasattr(coinc, ifos[0]) and hasattr(coinc, ifos[1]):
647 ethinca[i] = simpleEThinca(getattr(coinc, ifos[0]),
648 getattr(coinc, ifos[1]))
649 return ethinca
650
652 """
653 Return a new coincInspiralTable with triggers whose end_times lie within
654 segment; always use alphabetically first ifo's end_time.
655 """
656 triggers_within_segment = coincInspiralTable(stat=self.stat)
657
658 for trig in self:
659 end_time = getattr(trig, trig.get_ifos()[1][0]).end_time
660 if end_time in segment:
661 triggers_within_segment.append(trig)
662
663 return triggers_within_segment
664