Loading [MathJax]/extensions/TeX/AMSsymbols.js
LALInference 4.1.9.1-3a66518
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
make_injtimes.py
Go to the documentation of this file.
1# -*- coding: utf-8 -*-
2"""
3Michalis Agathos (2014)
4
5This script provides a set of good gps times for injecting a GW signal
6into real multiple detector data.
7The input consists of a set of science segment and veto segment ASCII files
8having the 4-column format:
9
10segid gpsstart gpsend duration
11
12The output consists of an ASCII file containing a list of gps injection times.
13If the timeslides functionality is activated, an additional ASCII file is
14output, containing timeslide information (in seconds) for each IFO, with
15respect to the nominal times provided in the first file.
16"""
17
18__author__ = "Michalis Agathos"
19__maintainer__ = "Michalis Agathos"
20__email__ = "michalis.agathos@ligo.org"
21__version__ = "1.0.0"
22__status__ = "Production"
23
24import matplotlib as mpl
25import argparse
26import os
27import sys
28import random as rd
29mpl.use('Agg')
30from pylab import *
31py_version = sys.version_info[:2]
32np_version = np.__version__
33
34colordict={'H1':'r', 'L1':'g', 'V1':'b', 'H1L1':'c', 'L1H1':'c', 'H1V1':'m', 'V1H1':'m', 'L1V1':'y', 'V1L1':'y', 'H1L1V1':'k'}
35
36'''Exclude segments that contain the following gpstimes'''
37excludetimes = [966951349, 967054240, 968021010]
38
39###########################
40#
41# CLASS DEFINITIONS
42#
43###########################
44
45
46class IFO:
47 '''An interferometer. Can also refer to multiple interferometers to host doubles, triples etc.'''
48 def __init__(self, name, segments, vetoes, minlen=1):
49 self._name = name
50 self._minlen = minlen
51
52 if type(segments) is str:
53 print("Reading segments for " + self._name + "")
54 self._allsegments = self.readSegments(segments)
55 elif type(segments) is ndarray:
56 self._allsegments = segmentData(list(segments))
57 elif isinstance(segments, segmentData):
58 self._allsegments = segments
59 else:
60 print("Cannot recognize segments!")
61 return -1
62 self._segments = self._allsegments.fits(self._minlen)
63
64 if type(vetoes) is str:
65 print("Reading veto segments for " + self._name + "")
66 self._vetoes = self.readSegments(vetoes)
67 elif type(vetoes) is ndarray:
68 self._vetoes = segmentData(list(vetoes))
69 elif isinstance(vetoes, segmentData):
70 self._vetoes = vetoes
71 else:
72 print("Cannot recognize veto segments!")
73 return -1
74
75 self.setUnvetoed()
76 print("Number of unvetoed segments that fit " + str(self._minlen) + ": " + str(self._unvetoed.length()))
77
78 def setUnvetoed(self):
79 '''This is getting a list of unvetoed segments that fit the minimum length'''
80 self._unvetoed = self._segments.getUnvetoed(self._vetoes).fits(self._minlen)
81 #self._unvetoed = self._segments.fits(self._minlen).getUnvetoed(self._vetoes).fits(self._minlen)
82
83 def readSegments(self, segfname):
84 print("Reading " + segfname)
85 segdata = genfromtxt(segfname, dtype=int)
86 return segmentData(list(segdata))
87
88 def printUnvetoedToFile(self, outname):
89 self._unvetoed.printToFile(outname)
90
91 def getTrigTimes(self, whereInj='middle', interval=None, rmargin=2, lmargin=2, n=None, outfile=None):
92 '''Returns a list of gps times on which injections can be made'''
93 trigtimes = []
94 if interval is None:
95 l = self._minlen
96 else:
97 l = interval
98 for seg in self._unvetoed._seglist:
99 t = int(seg[1])
100 if whereInj == 'often':
101 '''Inject every l seconds until you reach rmargin seconds from end of segment.'''
102 t += lmargin
103 while t + rmargin <= seg[2]:
104 if t >= seg[1] + lmargin:
105 trigtimes.append(t)
106 t += l
107 elif whereInj == 'middle':
108 '''Inject in the middle of minlen-long intervals.'''
109 while t + self._minlen <= seg[2]:
110 if (t + self._minlen//2 >= seg[1] + lmargin) and (t + self._minlen//2 + rmargin <= seg[2]):
111 trigtimes.append(t + self._minlen//2)
112 t += self._minlen
113 else:
114 print("Invalid whereInj argument. Output times will be empty.")
115
116 if n is not None:
117 trigtimes = trigtimes[:n]
118 if outfile is None:
119 return array(trigtimes)
120 else:
121 savetxt(outfile, array(trigtimes), fmt='%i')
122
123 def plotCumulativeDurations(self, outfile, maxdur=None):
124 print("Plotting segment lengths distribution to file " + outfile)
125 fig = figure()
126 ax = fig.add_subplot(111)
127 ax.set_xlabel("segment length")
128 ax.set_ylabel("Cumulative # segments")
129 sdur = array(self._segments._seglist)[:,-1]
130 udur = array(self._unvetoed._seglist)[:,-1]
131 if maxdur is not None:
132 sdur = sdur[where(sdur <= maxdur)[0]]
133 udur = udur[where(udur <= maxdur)[0]]
134 ax.set_ylim((0,max([len(udur),len(sdur)])))
135 ax.set_title("Segment length distribution for " + self._name)
136 ax.hist(sdur, bins=sort(sdur), cumulative=True, histtype='stepfilled', alpha=0.3, label="All segments (" + str(self._segments.length()) + ")")
137 ax.hist(udur, bins=sort(udur), cumulative=True, histtype='stepfilled', alpha=0.3, label="Unvetoed segments (" + str(self._unvetoed.length()) + ")")
138 ax.legend(loc=4)
139 fig.savefig(outfile)
140
141 def plotSegments(self, outfile=None, lenperline=200000, segcolorlist=['b','r','g'], title=None):
142 '''Plot segments in rows of lenperline seconds'''
143 start = min(self._allsegments._start, self._vetoes._start)
144 end = max(self._allsegments._end, self._vetoes._end)
145 n = ceil((end-start)/lenperline)
146 plot_params = {'figure.figsize': [fig_width,fig_height*2]}
147 rcParams.update(plot_params)
148 fig = figure()
149 for k, segs in zip([0,1,2],[self._segments, self._vetoes, self._unvetoed]):
150 s = array(segs._seglist)
151 t = sort(hstack((s[:,1],s[:,1],s[:,2],s[:,2])))
152 y = array(len(s)*[0,1,1,0])
153 for i in arange(n):
154 ax=fig.add_subplot(n,1,i+1)
155 ax.tick_params(axis='y', which='both', left='off', right='off', labelleft='off')
156 ax.tick_params(axis='x', which='both', top='on', bottom='off', labelbottom='off')
157 ax.set_ylim(0,2)
158 ax.xaxis.set_major_formatter(FormatStrFormatter('%0.0f'))
159 t0 = start + i*lenperline
160 t1 = t0 + lenperline
161 ax.set_ylabel(int(t0), verticalalignment='top', rotation=45)
162 ax.set_xlim(t0, t1)
163 ax.fill_between(t, 0*y + 2*(k//2), y + 2*(1-y)*(k//2), alpha=0.3, color=segcolorlist[k])
164 if title is not None:
165 fig.suptitle(title)
166 if outfile is None:
167 return fig
168 else:
169 print('Plotting segments to ', outfile)
170 fig.savefig(outfile)
171 plot_params = {'figure.figsize': [fig_width,fig_height]}
172 rcParams.update(plot_params)
173
174
175# def distToVeto(self, time):
176# '''Returns the distance to the closest veto segment in seconds'''
177# dist = self._unvetoed
178
179class segment:
180 '''A time segment in an IFO'''
181
182 def __init__(self, data, gpsstart=None, gpsend=None):
183 '''Creates a segment'''
184 if gpsstart is not None and gpsend is not None:
185 self._id = data[0]
186 self._start = gpsstart
187 self._end = gpsend
188 self._len = max(gpsend-gpsstart, 0)
189 elif gpsstart is None and gpsend is None:
190 '''Overloading to generate segment from an array (id, start, end, length)'''
191 if len(data) != 4:
192 print("Segment data doesn't have the correct length!")
193 return -1
194 self._id = data[0]
195 self._start = data[1]
196 self._end = data[2]
197 self._len = data[3]
198 else:
199 print("Wrong arguments to create segment!")
200 return -1
201# if self._len != (self._end - self._start):
202# print("Error in segment data: inconsistent length! " + str(self._len) + " " + str(self._end - self._start))
203
204 def intersectWithList(self, id0, other, sortedlist=True, startid=None):
205 '''Intersect segment with list of (non-overlapping) segments'''
206 newlist = []
207 id = id0
208 minid = startid
209 '''Cut list short from the left'''
210 for oseg in other[startid:]:
211 if segment(oseg)._start > self._end and sortedlist:
212 minid = max((0,segment(oseg)._id - 1))
213 break
214 newseg = self.intersectWith(segment(oseg))
215 if newseg._start < newseg._end:
216 newseg._id = id
217 id += 1
218 newlist.append(newseg.toArray())
219 if startid is None:
220 return newlist
221 else:
222 return (newlist, minid)
223
224 def hasTime(self, time):
225 '''Checks if time is in (open) segment'''
226 if (time < self._start or time > self._end):
227 dist = min(self._start - time, time - self._end)
228 else:
229 dist = 0
230 return dist
231
232 def intersectWith(self, other):
233 '''Intersects with another segment'''
234 newstart = max(self._start, other._start)
235 newend = min(self._end, other._end)
236 newseg = segment([self._id], newstart, newend)
237 return newseg
238
239 def toArray(self):
240 a = array([self._id, self._start, self._end, self._len])
241 return a
242
243class segmentData:
244 '''Data that holds segments for one or more IFOs.'''
245
246 def __init__(self, seglist, gpsstart=None, gpsend=None):
247 if gpsstart is None:
248 gpsstart = seglist[0][1]
249 if gpsend is None:
250 gpsend = seglist[-1][2]
251 sortedlist = list(array(seglist)[argsort(array(seglist)[:,1])])
252 self._seglist = sortedlist
253 self._start = gpsstart
254 self._end = gpsend
255
256 def length(self):
257 return len(self._seglist)
258
259 def fits(self, minlen):
260 '''Returns the list of segments that fit a minimum required length.'''
261 data = array(self._seglist)
262 fitted = data[where(data[:,3] >= minlen)[0]]
263 fitted = vstack( (arange(len(fitted)),fitted[:, 1:].T) ).T
264 return segmentData(list(fitted), self._start, self._end)
265
266 def intersectSegments(self, other):
267 '''Returns the intersection with another list of segments.'''
268 newdata = []
269 startfrom = 0
270 id0 = self._seglist[0][0]
271 for s in self._seglist:
272 exclude = False
273 seg = segment(s)
274 for et in excludetimes:
275 exclude = exclude or (seg.hasTime(et) == 0)
276 if not exclude:
277 id = id0 + len(newdata)
278 nd, minid= seg.intersectWithList(id, other._seglist, startid=startfrom)
279 # nd= seg.intersectWithList(id, other._seglist, startid=None)
280 newdata += nd
281 startfrom = minid
282 return segmentData(newdata)
283
284 def notSegments(self, start, end):
285 '''Get the complement of the (non-overlapping) segments in this data'''
286 notlist = []
287 c = 0
288 times = array(self._seglist)[:,[1,2]]
289 # print(shape(times))
290 times = times[argsort(times[:,1])]
291 t1 = start
292 for i in (arange(len(times))):
293 t2 = times[i,0]
294 if t2 > t1 and t1 >= start and t2 <= end:
295 newseg = array([c, t1, t2, t2-t1]) #FIXME: check c initial value
296 notlist.append(newseg)
297 c += 1
298 t1 = times[i,1]
299 if t1 < end:
300 newseg = array([c, t1, end, end-t1])
301 notlist.append(newseg)
302 return segmentData(notlist, start, end)
303
304 def getUnvetoed(self, vetodata):
305 '''Returns a segmentData object with unvetoed segments'''
306 noveto = vetodata.notSegments(self._start, self._end)
307 unvetoed = self.intersectSegments(noveto)
308 return unvetoed
309
310 def unionSegments(self, other):
311 '''Take the union of segment lists. This is used for combining veto segments of different detectors.'''
312 '''Here I use AUB = not(notA^notB)'''
313 start = min(self._start, other._start)
314 end = max(self._end, other._end)
315 united = self.notSegments(start, end).intersectSegments(other.notSegments(start, end)).notSegments(start, end)
316 return united
317
318 def hasTime(self, time):
319 '''Checks whether time is in the segments. Gives 0 if true and the distance of time to the segment list otherwise.'''
320 dist = max(abs(self._start - time), abs(self._end - time))
321 for s in self._seglist:
322 dist = min(dist, s.hasTime(time))
323 return dist
324
325 def timeIn(self, time):
326 '''If t is in a segment returns the distance to the edges, otherwise returns False'''
327 sl = array(self._seglist)
328 ds = time - sl[:,1]
329 de = time - sl[:,2]
330 sgn = ds*de #should be negative if segment contains t
331 ids = where(sgn < 0)[0]
332 if len(ids) > 1:
333 print("Time " + time + " contained in more than one segment!")
334 if len(ids) == 1:
335 return [ds[ids[0]], -de[ids[0]]]
336 else:
337 return False
338
339
340 def printToFile(self, outfile):
341 print("Printing segment list to file " + outfile)
342 savetxt(outfile, array(self._seglist), fmt='%i')
343
344 def plotSegments(self, outfile=None, lenperline=200000, segcolor='b', title=None):
345 '''Plot segments in rows of lenperline seconds'''
346 s = array(self._seglist)
347 t = sort(hstack((s[:,1],s[:,1],s[:,2],s[:,2])))
348 y = len(s)*[0,1,1,0]
349 n = ceil((t[-1]-t[0])/lenperline)
350 fig = figure()
351 for i in arange(n):
352 ax=fig.add_subplot(n,1,i+1)
353 ax.tick_params(axis='y', which='both', left='off', right='off', labelleft='off')
354 ax.tick_params(axis='x', which='both', top='on', bottom='off', labelbottom='off')
355 ax.set_ylim(0,2)
356 ax.xaxis.set_major_formatter(FormatStrFormatter('%0.0f'))
357 t0 = t[0] + i*lenperline
358 t1 = t0 + lenperline
359 ax.set_ylabel(int(t0), verticalalignment='top', rotation=45)
360 ax.set_xlim(t0, t1)
361 ax.fill_between(t, 0*array(y), y, alpha=0.3, color=segcolor)
362 if title is not None:
363 fig.suptitle(title)
364 if outfile is None:
365 return fig
366 else:
367 print('Plotting segments to ', outfile)
368 fig.savefig(outfile)
369
370
371 def plotCumulativeDurations(self, outfile, maxdur=None):
372 print("Plotting segment lengths distribution to file " + outfile)
373 fig = figure()
374 ax = fig.add_subplot(111)
375 ax.set_xlabel("segment length")
376 ax.set_ylabel("# segments")
377 dur = array(self._seglist)[:,-1]
378 if maxdur is not None:
379 dur = dur[where(dur <= maxdur)[0]]
380 ax.hist(dur, bins=sort(dur), cumulative=True, histtype='stepfilled', alpha=0.3)
381 fig.savefig(outfile)
382
383
384
385###########################
386#
387# OTHER FUNCTIONS
388#
389###########################
390
391
392def ensure_dir(f):
393 """
394 CREATE FOLDER IF IT DOES NOT EXIST
395 """
396 if not os.path.exists(f):
397 os.makedirs(f)
398
399
400def getDoubles(ifo1, ifo2, unvetoed=False):
401 '''
402 Combines 2 interferometer objects into a new one, taking intersection of segments and union of vetoes
403 Setting unvetoed=True will speed up the calculations by ignoring
404 small/vetoed segments. Use with caution! (not good for plotting seglen distr)
405 '''
406 name1 = ifo1._name
407 name2 = ifo2._name
408 seg1 = ifo1._allsegments
409 seg2 = ifo2._allsegments
410 vet1 = ifo1._vetoes
411 vet2 = ifo2._vetoes
412 unv1 = ifo1._unvetoed
413 unv2 = ifo2._unvetoed
414 minl1 = ifo1._minlen
415 minl2 = ifo2._minlen
416
417 name12 = name1 + name2
418 print('Combining ' + name1 + ' and ' + name2 + ' into ' + name12)
419 if unvetoed:
420 seg12 = unv1.intersectSegments(unv2) # only uses unvetoed long segment list for each IFO
421 else:
422 seg12 = seg1.intersectSegments(seg2) # uses full segment lists including small&vetoed segments
423 vet12 = vet1.unionSegments(vet2)
424 minl12 = max(minl1, minl2)
425 ifo12 = IFO(name12, seg12, vet12, minl12)
426# if ifo12
427 return ifo12
428
429
430def getTriples(ifo1, ifo2, ifo3, unvetoed=False):
431 '''Combines 3 IFO objects into one'''
432 ifo12 = getDoubles(ifo1, ifo2, unvetoed)
433 ifo123 = getDoubles(ifo12, ifo3, unvetoed)
434 return ifo123
435
436
437def generateTimeslides(tdict, n, ref=None, outfolder=None, verbose=False):
438 '''
439 Generate a list of injection times and timeslides,
440 given a dictionary of (single) injection times
441 '''
442 ifos = tdict.keys()
443 # nifo = len(ifos)
444 if ref is None:
445 ref = ifos[0]
446 injtimes = []
447 injdict = {}
448 slidedict = {}
449 for k in tdict.keys():
450 injdict[k] = []
451 slidedict[k] = []
452 for i in arange(n):
453 for ifo in ifos:
454 injdict[ifo].append(rd.sample(tdict[ifo], 1)[0])
455
456 injtimes = injdict[ref]
457 for ifo in ifos:
458 slidedict[ifo] = array(injdict[ifo]) - array(injdict[ref])
459
460 if outfolder is None:
461 return injtimes, slidedict
462 else:
463 label=''.join(ifos) + '_' + str(n)
464 injtimesfile = os.path.join(outfolder,'injtimes_'+label+'.dat')
465 slidefilename = os.path.join(outfolder,'timeslides_'+label+'.dat')
466 slidefile = open(slidefilename, 'w')
467 header = " ".join(ifos)
468 slidefile.write(header+'\n')
469 slidedata = array(slidedict.values()).T
470 if verbose:
471 print(injtimes)
472 print(slidedict)
473 savetxt(injtimesfile, array(injtimes), fmt='%i')
474 savetxt(slidefile, slidedata, delimiter=' ', fmt='%i')
475 slidefile.close()
476 return (injtimesfile, slidefilename)
477
478
479
480###########################
481#
482# MAIN FUNCTION
483#
484###########################
485
486if __name__ == "__main__":
487
488###########################################
489#
490# Parse arguments (argparser)
491#
492###########################################
493
494 #parser = OptionParser()
495 parser = argparse.ArgumentParser(description="Reads segment files and outputs injection times in veto-free segments.")
496
497
498 parser.add_argument('-i', '--ifos', nargs=3, type=str, metavar='IFO', dest='ifos', help='List of IFOs to be considered (currently takes 3)') # FIXME: Variable length???
499 parser.add_argument('-s', '--segfiles', nargs=3, type=str, metavar='FILE', dest='segf', help='path to the segment files for IFOS in the order entered above')
500 parser.add_argument('-v', '--vetofiles', nargs=3, type=str, metavar='FILE', dest='vetof', help='path to the veto files for IFOS in the order entered above')
501 parser.add_argument('-o', '--outfolder', type=str, metavar='FILE', dest='outf', help='path to the output folder', default='.')
502 parser.add_argument('-l', '--length', type=int, metavar='INT', dest='length', help='length of signal segments in seconds', default=45.0)
503 parser.add_argument('-P', '--psdlength', type=int, metavar='INT', dest='psdlength', help='minimum length for calculating PSD in seconds', default=1024)
504 parser.add_argument('-N', '--Ninj', type=int, dest="Ninj", help="# injections", default=1000)
505 parser.add_argument('-I', '--interval', type=int, metavar='INT', dest='interval', help='minimum space between injections in seconds (defaults to injecting only in the middle of segments)', default=None)
506 parser.add_argument('-S', '--singles', action="store_true", dest="singles", help="process singles (default for timeslides)", default=False)
507 parser.add_argument('-D', '--doubles', action="store_true", dest="doubles", help="process doubles (N/A for timeslides)", default=False)
508 parser.add_argument('-T', '--triples', action="store_true", dest="triples", help="process triples (N/A for timeslides)", default=False)
509 parser.add_argument('-t', '--timeslides', action="store_true", dest="timeslides", help="Enable timeslides. This automatically deactivates --triples and --doubles.", default=False)
510 parser.add_argument('-p', '--plot', action="store_true", dest="plotsegdist", help="plot cumulative segment length distribution", default=False)
511 parser.add_argument('-c', '--checkslides', action="store_true", dest="check", help="verbose check that all timeslid times are in unvetoed segments", default=False)
512 parser.add_argument('-u', '--dumpunvetoed', action="store_true", dest="dumpunvetoed", help="dump list of unvetoed segments to file and plot segments", default=False)
513
514 #TODO:
515 # * add option to dump to pickle (IFOlist, doubleIFOs, tripleIFO)
516 # * add option to read IFOs from pickle
517
518 args = parser.parse_args()
519
520 outfolder = str(args.outf)
521 segfiles = args.segf
522 vetofiles = args.vetof
523 ifos = args.ifos
524 nifos = len(ifos)
525 Ninj = args.Ninj
526 seglen = args.length
527 psdlen = args.psdlength
528 plotsegdist = args.plotsegdist
529 dumpunvetoed = args.dumpunvetoed
530 interval = args.interval
531 singles = args.singles
532 doubles = args.doubles
533 triples = args.triples
534 timeslides = args.timeslides
535 check = args.check
536 maxplot = 40000
537
538 # Ensure that directories exist
539 ensure_dir(outfolder)
540 if dumpunvetoed:
541 ensure_dir(os.path.join(outfolder, 'segments'))
542
543 # Choose where to inject (currently middle of segment or every interval)
544 if interval is None:
545 whereInj = 'middle'
546 else:
547 whereInj = 'often'
548
549 # Read segments data in rawlist and only keep long enough segments in fitlist for each of the IFOs
550 timesdict = {}
551 minlen = max(psdlen, seglen)
552
553 # Read in segment files and veto files and create IFO list
554 IFOlist = []
555 for i in arange(nifos):
556 IFOlist.append(IFO(ifos[i], segfiles[i], vetofiles[i], minlen))
557 # Generate triggers for each IFO (for timeslides only)
558 if timeslides:
559 if plotsegdist:
560 IFOlist[-1].plotCumulativeDurations(os.path.join(outfolder,'singleseg_' + IFOlist[-1]._name +'.png'), maxplot)
561 timesdict[ifos[i]] = IFOlist[-1].getTrigTimes(whereInj=whereInj, interval=interval, lmargin=seglen)
562
563 # Generate timeslide output
564 if timeslides:
565 ensure_dir(os.path.join(outfolder, 'timeslides'))
566 generateTimeslides(timesdict, Ninj, ref='H1', outfolder=os.path.join(outfolder,'timeslides'))
567 singles = True
568 doubles = False
569 triples = False
570
571 if singles:
572 for single in IFOlist:
573 if plotsegdist:
574 single.plotCumulativeDurations(os.path.join(outfolder, 'singleseg_' + single._name + '.png'), maxplot)
575 if not timeslides:
576 single.getTrigTimes(whereInj=whereInj, n=Ninj, interval=interval, lmargin=seglen, outfile=os.path.join(outfolder, 'injtimes_' + single._name + '_' + str(single._minlen) + '.dat'))
577 # Check unvetoedness and prints the distance of trigger times to edges of unvetoed segments
578 if check:
579 print('checking single ' + single._name)
580 for time in single.getTrigTimes(whereInj=whereInj, n=Ninj, interval=interval, lmargin=seglen):
581 dist = single._unvetoed.timeIn(time)
582 print(dist, dist > 0)
583 # Dump unvetoed segment list to file and plot segments to file
584 if dumpunvetoed:
585 single._unvetoed.printToFile(os.path.join(outfolder, 'segments', single._name + '_unvetoed_segs.dat'))
586 figseg = single._unvetoed.plotSegments(os.path.join(outfolder, 'segments', single._name + '_unvetoed_segs.pdf'), title=single._name + ' unvetoed segments', segcolor=colordict[single._name])
587
588
589 # Combine IFOs into doubles
590 if doubles:
591 doubleIFOs = []
592 doubleIFOs.append( getDoubles(IFOlist[0], IFOlist[1]) ) #, unvetoed=not plotsegdist) ) #FIXME: test it with plots!
593 doubleIFOs.append( getDoubles(IFOlist[1], IFOlist[2]) ) #, unvetoed=not plotsegdist) ) # speedup is ~3s for full triples set
594 doubleIFOs.append( getDoubles(IFOlist[0], IFOlist[2]) ) #, unvetoed=not plotsegdist) ) # not worth overall?
595
596 # Generate output for doubles
597 for double in doubleIFOs:
598 double.getTrigTimes(whereInj=whereInj, n=Ninj, interval=interval, lmargin=seglen, outfile=os.path.join(outfolder, 'injtimes_' + double._name + '_' + str(double._minlen) +'.dat'))
599 if plotsegdist:
600 double.plotCumulativeDurations(os.path.join(outfolder,'doubleseg_' + double._name +'.png'), maxplot)
601 # Check unvetoedness and prints the distance of trigger times to edges of unvetoed segments
602 if check:
603 print('checking double ' + double._name)
604 for time in double.getTrigTimes(whereInj=whereInj, n=Ninj, interval=interval, lmargin=seglen):
605 dist = double._unvetoed.timeIn(time)
606 print(dist, dist > 0)
607 # Dump unvetoed segment list to file and plot segments to file
608 if dumpunvetoed:
609 double._unvetoed.printToFile(os.path.join(outfolder, 'segments', double._name + '_unvetoed_segs.dat'))
610 figseg = double._unvetoed.plotSegments(os.path.join(outfolder, 'segments', double._name + '_unvetoed_segs.pdf'), title=double._name + ' unvetoed segments')
611
612
613 # Combine IFOs into triples
614 if triples:
615 # If doubles are combined half of the work is already done
616 if doubles:
617 tripleIFO = getDoubles(doubleIFOs[0], IFOlist[2]) #, unvetoed=not plotsegdist)
618 else:
619 tripleIFO = getTriples(IFOlist[0], IFOlist[1], IFOlist[2]) #, unvetoed=not plotsegdist)
620
621 # Generate output for doubles
622 tripleIFO.getTrigTimes(whereInj=whereInj, n=Ninj, interval=interval, lmargin=seglen, outfile=os.path.join(outfolder,'injtimes_' + tripleIFO._name + '_' + str(tripleIFO._minlen) +'.dat'))
623 if plotsegdist:
624 tripleIFO.plotCumulativeDurations(os.path.join(outfolder,'tripleseg_' + tripleIFO._name +'.png'), maxplot)
625 # Check unvetoedness and prints the distance of trigger times to edges of unvetoed segments
626 if check:
627 print('checking triple ' + tripleIFO._name)
628 for time in tripleIFO.getTrigTimes(whereInj=whereInj, n=Ninj, interval=interval, lmargin=seglen):
629 dist = tripleIFO._unvetoed.timeIn(time)
630 print(dist, dist > 0)
631 # Dump unvetoed segment list to file and plot segments to file
632 if dumpunvetoed:
633 tripleIFO._unvetoed.printToFile(os.path.join(outfolder, 'segments', tripleIFO._name + '_unvetoed_segs.dat'))
634 figseg = tripleIFO._unvetoed.plotSegments(os.path.join(outfolder, 'segments', tripleIFO._name + '_unvetoed_segs.pdf'), title=tripleIFO._name + ' unvetoed segments')
#define max(a, b)
def setUnvetoed(self)
This is getting a list of unvetoed segments that fit the minimum length.
def plotSegments(self, outfile=None, lenperline=200000, segcolorlist=['b', 'r', 'g'], title=None)
Plot segments in rows of lenperline seconds.
def getTrigTimes(self, whereInj='middle', interval=None, rmargin=2, lmargin=2, n=None, outfile=None)
Returns a list of gps times on which injections can be made.
def printUnvetoedToFile(self, outname)
def plotCumulativeDurations(self, outfile, maxdur=None)
def __init__(self, name, segments, vetoes, minlen=1)
Data that holds segments for one or more IFOs.
def plotSegments(self, outfile=None, lenperline=200000, segcolor='b', title=None)
Plot segments in rows of lenperline seconds.
def getUnvetoed(self, vetodata)
Returns a segmentData object with unvetoed segments.
def timeIn(self, time)
If t is in a segment returns the distance to the edges, otherwise returns False.
def plotCumulativeDurations(self, outfile, maxdur=None)
def __init__(self, seglist, gpsstart=None, gpsend=None)
def notSegments(self, start, end)
Get the complement of the (non-overlapping) segments in this data.
def hasTime(self, time)
Checks whether time is in the segments.
def unionSegments(self, other)
Take the union of segment lists.
def fits(self, minlen)
Returns the list of segments that fit a minimum required length.
def intersectSegments(self, other)
Returns the intersection with another list of segments.
def intersectWith(self, other)
Intersects with another segment.
def __init__(self, data, gpsstart=None, gpsend=None)
Creates a segment.
def hasTime(self, time)
Checks if time is in (open) segment.
def intersectWithList(self, id0, other, sortedlist=True, startid=None)
Intersect segment with list of (non-overlapping) segments.
def getTriples(ifo1, ifo2, ifo3, unvetoed=False)
Combines 3 IFO objects into one.
def getDoubles(ifo1, ifo2, unvetoed=False)
Combines 2 interferometer objects into a new one, taking intersection of segments and union of vetoes...
def ensure_dir(f)
CREATE FOLDER IF IT DOES NOT EXIST.
def generateTimeslides(tdict, n, ref=None, outfolder=None, verbose=False)
Generate a list of injection times and timeslides, given a dictionary of (single) injection times.