Package pylal :: Package dq :: Module dqTriggerUtils
[hide private]
[frames] | no frames]

Source Code for Module pylal.dq.dqTriggerUtils

   1  #!/usr/bin/env python 
   2   
   3  # ============================================================================= 
   4  # Preamble 
   5  # ============================================================================= 
   6   
   7  from __future__ import division 
   8  import sys,os,re,math,datetime,glob,copy 
   9  from socket import getfqdn 
  10   
  11  from glue.ligolw import ligolw,table,lsctables,utils 
  12  from glue.ligolw.utils import process as ligolw_process 
  13  from glue import segments 
  14   
  15  from glue.lal import Cache as LALCache 
  16  from glue.lal import CacheEntry as LALCacheEntry 
  17   
  18  from pylal import date 
  19  from pylal.xlal.datatypes.ligotimegps import LIGOTimeGPS 
  20   
  21  from glue import git_version 
  22   
  23  from scipy import special 
  24  import numpy 
  25   
  26  __author__  = "Duncan Macleod <duncan.macleod@astro.cf.ac.uk>" 
  27  __version__ = "git id %s" % git_version.id 
  28  __date__    = git_version.date 
  29   
  30  """ 
  31  This module provides a bank of useful functions for manipulating triggers and trigger files for data quality investigations. 
  32  """ 
  33   
  34  # global regular expressions 
  35  trigsep = re.compile('[\t\s,]+') 
  36  _cchar_regex = re.compile('[-#%<!()_\[\]-{}:;\'\"\ ]') 
  37   
  38  # ============================================================================= 
  39  # Define ETG options 
  40  # ============================================================================= 
  41   
  42  _trig_regex  = re.compile('(burst|inspiral|ring)', re.I) 
  43  _burst_regex = re.compile('(burst|omega|kleine|kw|cwb|hacr)', re.I) 
  44  _cbc_regex   = re.compile('(ihope|inspiral|cbc)', re.I) 
  45  _ring_regex  = re.compile('(ring)', re.I) 
  46   
47 -def SnglTriggerTable(etg, columns=None):
48 49 """ 50 Handy function to return the correct type of ligolw table for the given ETG. 51 """ 52 53 if _burst_regex.search(etg): 54 t = lsctables.New(lsctables.SnglBurstTable, columns=columns) 55 elif _cbc_regex.search(etg): 56 t = lsctables.New(lsctables.SnglInspiralTable, columns=columns) 57 elif _ring_regex.search(etg): 58 t = lsctables.New(lsctables.SnglRingdownTable, columns=columns) 59 else: 60 raise AttributeError("etg=%s not recognised by SnglTriggerTable." % etg) 61 62 # set columns 63 if columns: 64 if isinstance(columns[0], str): columns = map(str.lower, columns) 65 if isinstance(columns[0], unicode): columns = map(unicode.lower, columns) 66 for c in t.columnnames: 67 if c.lower() not in columns: 68 idx = t.columnnames.index(c) 69 t.columnnames.pop(idx) 70 t.columntypes.pop(idx) 71 72 return t
73
74 -def SnglTrigger(etg):
75 """ 76 Handy function to return the correct type of ligolw table row object for 77 the given ETG 78 """ 79 80 if _burst_regex.search(etg): 81 return lsctables.SnglBurst() 82 elif _cbc_regex.search(etg): 83 return lsctables.SnglInspiral() 84 elif _ring_regex.search(etg): 85 return lsctables.SnglRingdown() 86 else: 87 raise AttributeError("etg=%s not recognised by SnglTrigger." % etg)
88 89 # ============================================================================= 90 # Define get_time choice 91 # ============================================================================= 92
93 -def def_get_time(tableName, ifo=None):
94 95 """ 96 Define the get_time() function for given table 97 """ 98 99 get_time = None 100 101 if ifo: ifo = ifo[0] 102 103 # if given an injection table: 104 if re.match('sim', tableName): 105 106 if re.search('inspiral',tableName): 107 get_time = lambda row: row.get_end(site=ifo) 108 else: 109 get_time = lambda row: row.get_time_geocent() 110 111 112 # if given a sngl trigger table 113 elif re.match('(sngl|multi|coinc)', tableName): 114 115 if re.search('inspiral',tableName): 116 get_time = lambda row: row.get_end() 117 elif re.search('ringdown',tableName): 118 get_time = lambda row: row.get_start() 119 else: 120 get_time = lambda row: row.get_peak() 121 122 return get_time
123 124 # ============================================================================= 125 # Convert list from text file into Sngl{Burst,Inspiral} object 126 # ============================================================================= 127
128 -def trigger(data, etg, ifo=None, channel=None, loadcolumns=None):
129 130 """ 131 Reads the list object data and returns a Sngl{Burst,Inspiral} object 132 with the correct columns seeded, based on the given string object etg. 133 134 Arguments: 135 136 data : [ string | list ] 137 string or list object containing the data to be parsed. 138 139 etg : [ "ihope" | "kw" | "omega" | "omegadq" ] 140 Defines how the data list is understood and parsed, based on the 141 standard column layouts (comma-, space-, or tab-delimited) for each ETG. 142 "ihope" is read as per the .csv files from ihope_daily. "omegadq" 143 is read as per the .txt or .clst produced by the detchar scripts used to 144 plot omega online events. 145 146 """ 147 148 # set up trig object 149 trig = SnglTrigger(etg) 150 151 # set up load columns 152 if loadcolumns==None: 153 loadcolumns = trig.__slots__ 154 155 # if given string, split on space, tab or comma 156 if isinstance(data,str): 157 data = trigsep.split(data.rstrip()) 158 159 # ===== 160 # ihope 161 # ===== 162 if re.match('ihope',etg): 163 # comma separated values are: 164 # end_time,end_time_ns,ifo,snr,mass1,mass2,mtotal,eta,event_duration, 165 # template_duration,eff_distance,chisq,chisq_dof,bank_chisq,bank_chisq_dof 166 # cont_chisq,cont_chisq_dof 167 168 trig.end_time = int(data[0]) 169 trig.end_time_ns = int(data[1]) 170 trig.ifo = str(data[2]) 171 trig.snr = float(data[3]) 172 trig.mass1 = float(data[4]) 173 trig.mass2 = float(data[5]) 174 trig.mtotal = float(data[6]) 175 trig.eta = float(data[7]) 176 trig.event_duration = float(data[8]) 177 trig.template_duration = float(data[9]) 178 trig.eff_distance = float(data[10]) 179 trig.chisq = float(data[11]) 180 trig.chisq_dof = float(data[12]) 181 trig.bank_chisq = float(data[13]) 182 try: 183 trig.bank_chisq_dof = float(data[14]) 184 except: 185 trig.bank_chisq_dof = None 186 try: 187 trig.cont_chisq = float(data[15]) 188 trig.cont_chisq_dof = float(data[16]) 189 except ValueError: 190 trig.cont_chisq = None 191 trig.cont_chisq_dof = None 192 193 # ============ 194 # kleine-welle 195 # ============ 196 elif etg=='kw': 197 # space separated values are: 198 # peak_time.peak_time_ns start_time.start_time_ns stop_time.stop_time_ns 199 # central_freq,energy,amplitude,n_pix,significance,N 200 201 # peak time 202 peak = LIGOTimeGPS(data[2]) 203 trig.peak_time = peak.seconds 204 trig.peak_time_ns = peak.nanoseconds 205 # start time 206 start = LIGOTimeGPS(data[0]) 207 ms_start = start 208 trig.start_time = start.seconds 209 trig.start_time_ns = start.nanoseconds 210 trig.ms_start_time = ms_start.seconds 211 trig.ms_start_time_ns = ms_start.nanoseconds 212 # end time 213 stop = LIGOTimeGPS(data[1]) 214 ms_stop = stop 215 trig.stop_time = stop.seconds 216 trig.stop_time_ns = stop.nanoseconds 217 trig.ms_stop_time = ms_stop.seconds 218 trig.ms_stop_time_ns = ms_stop.nanoseconds 219 # duration 220 trig.duration = stop-start 221 trig.ms_duration = ms_stop-ms_start 222 # others 223 trig.central_freq = float(data[3]) 224 trig.peak_frequency = float(data[3]) 225 energy = float(data[4]) 226 trig.amplitude = float(data[5]) 227 n_pix = float(data[6]) 228 #significance = float(data[7]) 229 #N = float(data[8]) 230 trig.snr = math.sqrt(trig.amplitude-n_pix) 231 232 # ===== 233 # omega 234 # ===== 235 if etg=='omega' or etg=='wpipe': 236 # space separated values are: 237 # peak_time.peak_time_ns peak_frequency duration bandwidth amplitude 238 # cluster_size cluster_norm_energy cluster_number 239 240 # peak time 241 peak = LIGOTimeGPS(data[0]) 242 trig.peak_time = peak.seconds 243 trig.peak_time_ns = peak.nanoseconds 244 # central frequency 245 trig.central_freq = float(data[1]) 246 trig.peak_frequency = trig.central_freq 247 # duration 248 trig.duration = LIGOTimeGPS(float(data[2])) 249 halfduration = LIGOTimeGPS(float(data[2])/2) 250 # start time 251 start = peak-halfduration 252 ms_start = start 253 trig.start_time = start.seconds 254 trig.start_time_ns = start.nanoseconds 255 trig.ms_start_time = ms_start.seconds 256 trig.ms_start_time_ns = ms_start.nanoseconds 257 # end time 258 stop = peak+halfduration 259 ms_stop = stop 260 trig.stop_time = stop.seconds 261 trig.stop_time_ns = stop.nanoseconds 262 trig.ms_stop_time = ms_stop.seconds 263 trig.ms_stop_time_ns = ms_stop.nanoseconds 264 265 trig.ms_duration = ms_stop-ms_start 266 # bandwidth and flow,fhigh 267 trig.bandwidth = float(data[3]) 268 trig.flow = trig.peak_frequency - 0.5*trig.bandwidth 269 trig.fhigh = trig.peak_frequency + 0.5*trig.bandwidth 270 # energy 271 trig.amplitude = float(data[4]) 272 273 # cluster parameters 274 if len(data)>5: 275 trig.param_one_name = 'cluster_size' 276 trig.param_one_value = float(data[5]) 277 trig.param_two_name = 'cluster_norm_energy' 278 trig.param_two_value = float(data[6]) 279 trig.param_three_name = 'cluster_number' 280 trig.param_three_value = float(data[7]) 281 282 # SNR 283 trig.snr = math.sqrt(2*trig.amplitude) 284 285 # ===== 286 # omega 287 # ===== 288 if etg=='omegaspectrum': 289 # space separated values are: 290 # peak_time.peak_time_ns peak_frequency duration bandwidth amplitude 291 # cluster_size cluster_norm_energy cluster_number 292 293 # peak time 294 peak = LIGOTimeGPS(data[0]) 295 trig.peak_time = peak.seconds 296 trig.peak_time_ns = peak.nanoseconds 297 # central frequency 298 trig.central_freq = float(data[1]) 299 trig.peak_frequency = trig.central_freq 300 301 trig.amplitude = float(data[2]) 302 303 # SNR 304 trig.snr = math.sqrt(trig.amplitude) 305 306 # use chisq to stor PSD variance, 307 if len(data)>3: 308 trig.chisq = math.sqrt(math.sqrt(float(data[3])))/trig.snr 309 310 # ============== 311 # omegadq 312 # ============== 313 # follow Cadonati's clustering output 314 if etg=='omegadq': 315 # start time 316 start = LIGOTimeGPS(data[0]) 317 trig.start_time = start.seconds 318 trig.start_time_ns = start.nanoseconds 319 # end time 320 stop = LIGOTimeGPS(data[1]) 321 trig.stop_time = stop.seconds 322 trig.stop_time_ns = stop.nanoseconds 323 # peak time 324 peak = LIGOTimeGPS(data[2]) 325 trig.peak_time = peak.seconds 326 trig.peak_time_ns = peak.nanoseconds 327 # duration 328 trig.duration = stop-start 329 # bandwidth, and flow,fhigh,central_freq 330 trig.flow = float(data[3]) 331 trig.fhigh = float(data[4]) 332 trig.bandwidth = trig.fhigh - trig.flow 333 trig.central_freq = trig.flow + 0.5*trig.bandwidth 334 335 # MS params 336 ms_start = LIGOTimeGPS(data[6]) 337 trig.ms_start_time = ms_start.seconds 338 trig.ms_start_time_ns = ms_start.nanoseconds 339 ms_stop = LIGOTimeGPS(data[7]) 340 trig.ms_stop_time = ms_stop.seconds 341 trig.ms_stop_time_ns = ms_stop.nanoseconds 342 trig.ms_duration = ms_stop-ms_start 343 trig.ms_flow = float(data[8]) 344 trig.ms_fhigh = float(data[9]) 345 trig.ms_bandwidth = trig.ms_fhigh - trig.ms_flow 346 trig.peak_frequency = trig.ms_flow + 0.5*trig.ms_bandwidth 347 348 # SNR 349 trig.amplitude = float(data[11]) 350 trig.snr = math.sqrt(2*float(data[11])) 351 trig.ms_snr = math.sqrt(2*float(data[12])) 352 353 # ==== 354 # HACR 355 # ==== 356 357 # based on output of hacr web interface 358 if etg=='hacr': 359 # peak time 360 trig.peak_time = int(data[0]) 361 trig.peak_time_ns = int(float(data[1])*math.pow(10,9)) 362 trig.param_one_name = 'peak_time_offset' 363 trig.param_one_value = float(data[1]) 364 # duration 365 trig.duration = float(data[4]) 366 # start time 367 start = trig.get_peak()-trig.duration 368 trig.start_time = start.seconds 369 trig.start_time_ns = start.nanoseconds 370 # end time 371 stop = trig.get_peak()+trig.duration 372 trig.stop_time = stop.seconds 373 trig.stop_time_ns = stop.nanoseconds 374 # bandwidth, and flow,fhigh,central_freq 375 trig.central_freq = float(data[2]) 376 trig.peak_frequency = trig.central_freq 377 trig.bandwidth = float(data[3]) 378 trig.peak_frequency_error = trig.bandwidth/trig.central_freq 379 trig.flow = trig.central_freq - 0.5*trig.bandwidth 380 trig.fhigh = trig.central_freq + 0.5*trig.bandwidth 381 #snr 382 trig.snr = float(data[6]) 383 trig.ms_snr = float(data[6]) 384 385 # ms extras 386 ms_start = start 387 trig.ms_start_time = ms_start.seconds 388 trig.ms_start_time_ns = ms_start.nanoseconds 389 ms_stop = stop 390 trig.ms_stop_time = ms_stop.seconds 391 trig.ms_stop_time_ns = ms_stop.nanoseconds 392 trig.ms_duration = ms_stop-ms_start 393 trig.ms_fhigh = trig.fhigh 394 trig.ms_flow = trig.flow 395 trig.ms_bandwidth = trig.ms_fhigh-trig.ms_flow 396 397 # others 398 trig.param_two_name = 'numPixels' 399 trig.param_two_value = int(data[5]) 400 trig.param_three_name = 'totPower' 401 trig.param_three_value = float(data[7]) 402 403 #trig.process_id = int(float(data[10])) 404 405 # sundries 406 if ifo: 407 trig.ifo = ifo 408 if channel: 409 trig.channel = channel 410 411 return trig
412 413 # ============================================================================= 414 # Write triggers to file in etg standard form 415 # ============================================================================= 416
417 -def totrigxml(file, table, program=None, params=[]):
418 419 """ 420 Write the given lsctables compatible table object to the file object file 421 in xml format. 422 423 Arguments: 424 425 file: file object 426 file object describing output file 427 428 table: lsctable 429 glue.ligolw.lsctables compatible table 430 431 Keyword Arguments: 432 433 program: string 434 name of generating executable, defaults to self 435 params: list 436 list of (param,type,value) tuples to append to process_params:table 437 438 """ 439 440 xmldoc = ligolw.Document() 441 xmlligolw = ligolw.LIGO_LW() 442 xmldoc.appendChild(ligolw.LIGO_LW()) 443 xmldoc.childNodes[-1].appendChild(lsctables.New(lsctables.ProcessTable)) 444 xmldoc.childNodes[-1].appendChild(lsctables.New(lsctables.ProcessParamsTable)) 445 446 # append process to table 447 if not program: 448 program='pylal.dq.dqDataUtils.totrigxml' 449 450 process = ligolw_process.append_process(xmldoc, program=program,\ 451 version=__version__,\ 452 cvs_repository = 'lscsoft',\ 453 cvs_entry_time = __date__) 454 455 ligolw_process.append_process_params(xmldoc, process, params) 456 457 # append trig table to file 458 xmldoc.childNodes[-1].appendChild(table) 459 460 # write triggers to file object file 461 ligolw_process.set_process_end_time(process) 462 with utils.SignalsTrap(): 463 utils.write_fileobj(xmldoc, file, gz=file.name.endswith('gz'))
464 465 # ============================================================================= 466 # Write triggers to text file in etg standard form 467 # ============================================================================= 468
469 -def totrigfile(file,table,etg,header=True,columns=None):
470 471 """ 472 Write the lsctables.table to the given file object file, in the standard 473 format for the given etg. 474 475 If columns is given as a list, the standard format will be overwritten 476 and only the given columns will be processed. 477 """ 478 479 etg = etg.lower() 480 481 # set columns 482 if not columns: 483 if re.match('ihope',etg): 484 # comma separated values are: 485 columns = ['end_time','end_time_ns','ifo','snr','mass1','mass2','mtotal',\ 486 'eta','event_duration','template_duration','eff_distance',\ 487 'chisq','chisq_dof','bank_chisq','bank_chisq_dof',\ 488 'cont_chisq','cont_chisq_dof'] 489 490 elif re.match('omegadq',etg) or re.match('omega_dq',etg): 491 columns = ['start_time','stop_time','peak_time','flow','fhigh',\ 492 'cluster_length','ms_start_time','ms_stop_time','ms_flow',\ 493 'ms_fhigh','cluster_size','amplitude','amplitude'] 494 495 elif re.match('omegaspectrum',etg): 496 columns = ['peak_time','peak_frequency','amplitude','chisq'] 497 498 elif re.match('omega',etg) or re.match('wpipe',etg): 499 if len(table) and hasattr(table[0], 'param_one_value'): 500 columns = ['peak_time','peak_frequency','duration',\ 501 'bandwidth','amplitude',\ 502 'param_one_value','param_two_value','param_three_value'] 503 else: 504 columns = ['peak_time','peak_frequency','duration',\ 505 'bandwidth','amplitude'] 506 507 elif re.match('kw',etg.lower()): 508 columns = ['peak_time','start_time','stop_time','peak_frequency',\ 509 'energy','amplitude','num_pixels','significance','N'] 510 511 elif re.match('hacr',etg.lower()): 512 columns = ['peak_time','param_one_value','peak_frequency','bandwidth',\ 513 'duration','param_two_value','snr','param_three_value'] 514 515 # set delimiter 516 if etg=='ihope': 517 d = ',' 518 else: 519 d=' ' 520 521 # print header 522 if header: 523 cols = [] 524 for c in columns: 525 if re.search('param_[a-z]+_value', c) and len(table)>0\ 526 and hasattr(table[0], c.replace('value','name')): 527 cols.append(getattr(table[0], c.replace('value','name'))) 528 else: 529 cols.append(c) 530 file.write('%s\n' % d.join(['#']+cols)) 531 532 # work out columns 533 if len(table)>0: 534 t = table[0] 535 columnnames = [t for t in table[0].__slots__ if hasattr(table[0], t)] 536 537 # print triggers 538 for row in table: 539 line = [] 540 for col in columns: 541 if col not in columnnames: 542 line.append('-1') 543 continue 544 entry = '' 545 # if ihope, print column 546 if re.match('(ihope|hacr)',etg.lower()): 547 # HACR default is to have peak_time_ns in seconds, not ns 548 if re.match('hacr',etg.lower()) and col=='peak_time_ns': 549 entry = str(getattr(row, col)) 550 else: 551 entry = str(getattr(row, col)) 552 # if not ihope, check for time and print full GPS 553 else: 554 if col=='peak_time': 555 entry = str(row.get_peak()) 556 elif col=='start_time': 557 entry = str(row.get_start()) 558 elif col=='ms_start_time': 559 entry = str(row.get_ms_start()) 560 elif col=='stop_time': 561 entry = str(row.get_stop()) 562 elif col=='ms_stop_time': 563 entry = str(row.get_ms_stop()) 564 else: 565 entry = str(getattr(row, col)) 566 567 line.append(entry) 568 569 file.write('%s\n' % d.join(line))
570 571 # ============================================================================= 572 # Function to load triggers from xml 573 # ============================================================================= 574
575 -def fromtrigxml(file,tablename='sngl_inspiral:table',start=None,end=None,\ 576 columns=None):
577 578 """ 579 Reads a trigger table from the given table from the xml 580 file object file 581 582 Arguments: 583 584 file : file object 585 586 Keyword arguments: 587 588 start : [ float | int | LIGOTimeGPS ] 589 GPS start time of requested period 590 end : [ float | int | LIGOTimeGPS ] 591 GPS end time of requested period 592 tablename : string 593 name of requested trigger table in xml file, defaults to 594 'sngl_inspiral:table' 595 """ 596 597 # set times 598 if not start: 599 start=0 600 if not end: 601 end=9999999999 602 603 span = segments.segment(start,end) 604 605 # set columns 606 if columns!=None: 607 if re.search('sngl_burst', tablename): 608 lsctables.SnglBurstTable.loadcolumns = columns 609 elif re.search('sngl_inspiral', tablename): 610 lsctables.SnglInspiralTable.loadcolumns = columns 611 if re.search('multi_burst', tablename): 612 lsctables.MultiBurstTable.loadcolumns = columns 613 elif re.search('multi_inspiral', tablename): 614 lsctables.MultiInspiralTable.loadcolumns = columns 615 616 # set tablename 617 if not tablename.endswith(':table'): 618 tablename = ':'.join([tablename,'table']) 619 620 # crack open xml file 621 xmldoc,digest = utils.load_fileobj(file,gz=file.name.endswith('gz')) 622 alltriggers = table.get_table(xmldoc,tablename) 623 624 triggers = table.new_from_template(alltriggers) 625 append = triggers.append 626 627 get_time = def_get_time(triggers.tableName) 628 # parse triggers in time 629 for row in alltriggers: 630 if float(get_time(row)) in span: 631 append(row) 632 633 # sort table in time 634 triggers.sort(key=lambda trig: float(get_time(trig))) 635 636 # reset columns 637 if columns: 638 type(triggers).loadcolumns = None 639 640 return triggers
641 642 # ============================================================================= 643 # Load triggers from text file 644 # ============================================================================= 645
646 -def fromtrigfile(file,etg,start=None,end=None,ifo=None,channel=None,\ 647 tabletype=None, columns=None, virgo=False):
648 649 """ 650 Reads the file object file containing standard columns for the given etg and 651 returns either a corresponding lsctable. 652 653 Arguments: 654 655 file : file object 656 etg : [ "ihope" | "kw" | "omega" | "omegadq" ] 657 string defining how to read text file. 658 "ihope" is read as per the .csv files from ihope_daily. "omegadq" 659 is read as per the .txt or .clst produced by the detchar scripts used to 660 plot omega online events. 661 662 Keyword arguments: 663 664 start : [ float | int | LIGOTimeGPS ] 665 GPS start time of requested period 666 end : [ float | int | LIGOTimeGPS ] 667 GPS end time of requested period 668 tabletype : type 669 Specific ligolw table type for output. By default tables will be 670 SnglInspiralTable or SnglBurstTable type depending on ETG 671 """ 672 673 if re.search('omegaspectrum', etg, re.I): 674 return fromomegaspectrumfile(file, start=start, end=end, ifo=ifo,\ 675 channel=channel, columns=columns) 676 elif re.search('omegadq', etg, re.I): 677 return fromomegadqfile(file, start=start, end=end, ifo=ifo,\ 678 channel=channel,columns=columns) 679 elif re.search('omega', etg, re.I): 680 return fromomegafile(file, start=start, end=end, ifo=ifo, channel=channel,\ 681 columns=columns, virgo=virgo) 682 elif re.search('kw', etg, re.I): 683 return fromkwfile(file, start=start, end=end, ifo=ifo, channel=channel,\ 684 columns=columns) 685 elif re.search('hacr', etg, re.I): 686 return fromhacrfile(file, start=start, end=end, ifo=ifo, channel=channel,\ 687 columns=columns) 688 elif re.search('ihope', etg, re.I): 689 return fromihopefile(file, start=start, end=end, ifo=ifo, channel=channel,\ 690 columns=columns)
691 692 # ============================================================================= 693 # Load triggers from a cache 694 # ============================================================================= 695
696 -def fromLALCache(cache, etg, start=None, end=None, columns=None,\ 697 virgo=False, verbose=False, snr=False):
698 699 """ 700 Extract triggers froa given ETG from all files in a glue.lal.Cache object. 701 Returns a glue.ligolw.Table relevant to the given trigger generator etg. 702 """ 703 704 # set up counter 705 if verbose: 706 sys.stdout.write("Extracting %s triggers from %d files... "\ 707 % (etg, len(cache))) 708 sys.stdout.flush() 709 delete = '\b\b\b' 710 num = len(cache)/100 711 712 trigs = SnglTriggerTable(etg, columns=columns) 713 714 # load files 715 for i,e in enumerate(cache): 716 if re.search('(xml|xml.gz)\Z', e.path): 717 trigsTmp = fromtrigxml(open(e.path), tablename=trigs.tableName,\ 718 start=start, end=end, columns=columns) 719 else: 720 trigsTmp = fromtrigfile(open(e.path), etg=etg, start=start, end=end,\ 721 columns=columns, virgo=virgo) 722 # keep only triggers above SNR threshold if requested 723 if snr: 724 trigs.extend(t for t in trigsTmp if t.snr > snr) 725 else: 726 trigs.extend(trigsTmp) 727 728 # print verbose message 729 if verbose and len(cache)>1: 730 progress = int((i+1)/num) 731 sys.stdout.write('%s%.2d%%' % (delete, progress)) 732 sys.stdout.flush() 733 734 if verbose: sys.stdout.write("\n") 735 return trigs
736 737 # ============================================================================= 738 # Generate a daily ihope cache 739 # ============================================================================= 740
741 -def daily_ihope_cache(start,end,ifo,cluster=None,filetype='xml',cat=0):
742 743 """ 744 Generates glue.lal.Cache containing CacheEntires for all daily ihope 745 INSPIRAL files for given ifo and clustering between start and end time. 746 747 Arguments: 748 749 start : [ float | int | LIGOTimeGPS ] 750 GPS start time of requested period 751 end : [ float | int | LIGOTimeGPS ] 752 GPS end time of requested period 753 ifo : [ "H1" | "L1" | "V1" ] 754 IFO 755 756 Keyword arguments: 757 cluster : [ "unclustered" | "100ms" | "30ms" | "16s" ] 758 clustering time in human format 759 filetype : [ "xml" | "csv" ] 760 file format of desired cache contents 761 cat : [ 0 | 1 | 2 | 3 | 4 ] 762 int veto category of trigger files requested 763 """ 764 765 # daily path 766 ihope_daily_path = os.path.expanduser('~cbc/ihope_daily') 767 768 # set clustering tag 769 if cluster==None or cluster.upper()=='UNCLUSTERED': 770 cluster_tag='UNCLUSTERED' 771 elif cluster.upper()=='100MS': 772 cluster_tag='100MILLISEC_CLUSTERED' 773 elif cluster.upper()=='30MS': 774 cluster_tag='30MILLISEC_CLUSTERED' 775 elif cluster.upper()=='16S': 776 cluster_tag='16SEC_CLUSTERED' 777 778 # work out days 779 days = gps_day_list(start,end) 780 span = segments.segment(start,end) 781 cache = LALCache() 782 # loop over days gathering files 783 for day in days: 784 utc = datetime.datetime(*date.XLALGPSToUTC(day)[:6]) 785 day_path = os.path.join(ihope_daily_path,utc.strftime("%Y%m"), 786 utc.strftime("%Y%m%d")) 787 788 if filetype=='xml': 789 filenames = glob.glob(os.path.join(day_path, 790 ifo+'-INSPIRAL_'+cluster_tag+'*.xml.gz')) 791 792 for filename in filenames: 793 e = LALCacheEntry.from_T050017(filename) 794 if span.intersects(e.segment): cache.append(e) 795 796 elif filetype=='csv': 797 csvfile = os.path.join(day_path,ifo+'-'+str(cat)+'-INSPIRAL_'+\ 798 cluster_tag+'.csv') 799 if os.path.isfile(csvfile): 800 e = LALCacheEntry.from_T050017(csvfile) 801 if span.intersects(e.segment): cache.append(e) 802 803 cache.sort(key=lambda e: e.path) 804 805 return cache
806 807 # ============================================================================= 808 # Function to generate an omega online cache 809 # ============================================================================= 810
811 -def omega_online_cache(start, end, channel, mask='DOWNSELECT',\ 812 check_files_exist=False, **kwargs):
813 814 """ 815 Returns a glue.lal.Cache contatining CacheEntires for all omega online 816 trigger files between the given start and end time for the given ifo. 817 For S6 triggers are only available for each IFO on it's own site cluster. 818 819 Arguments: 820 821 start : [ float | int | LIGOTimeGPS ] 822 GPS start time of requested period 823 end : [ float | int | LIGOTimeGPS ] 824 GPS end time of requested period 825 channel name or ifo : [ "H1" | "L1" | "V1" | "G1:SEI_TCC_STS2x" | ... ] 826 IFO 827 """ 828 829 # format channel 830 if re.match('\w\d:', channel): 831 ifo, channel = channel.split(':', 1) 832 else: 833 ifo = channel 834 channel = None 835 836 cache = LALCache() 837 838 # verify host 839 host = getfqdn() 840 ifo_host = { 'G1':'atlas', 'H1':'ligo-wa', 'H2':'ligo-wa', 'L1':'ligo-la'} 841 if (not kwargs.has_key('directory') and not re.search(ifo_host[ifo],host)): 842 print >>sys.stderr, "Error: Omega online files are not available for "+\ 843 "IFO=%s on this host." % ifo 844 return cache 845 846 span = segments.segment(start,end) 847 if ifo == 'G1': 848 if channel: 849 kwargs.setdefault('directory', '/home/omega/online/%s_%s/segments' % (ifo, channel)) 850 else: 851 kwargs.setdefault('directory', '/home/omega/online/%s/segments' % ifo) 852 kwargs.setdefault('epoch', 0) 853 else: 854 kwargs.setdefault('directory',\ 855 '/home/omega/online/%s/archive/S6/segments' % ifo) 856 kwargs.setdefault('epoch', 931211808) 857 kwargs.setdefault('duration', 64) 858 kwargs.setdefault('overlap', 8) 859 860 # optimise 861 append = cache.append 862 splitext = os.path.splitext 863 isfile = os.path.isfile 864 intersects = span.intersects 865 segment = segments.segment 866 from_T050017 = LALCacheEntry.from_T050017 867 basedir = kwargs['directory'] 868 basetime = kwargs['epoch'] 869 triglength = kwargs['duration'] 870 overlap = kwargs['overlap'] 871 872 # get times 873 start_time = int(start-math.fmod(start-basetime,triglength-overlap)) 874 t = start_time 875 876 # loop over time segments constructing file paths and appending to the cache 877 while t<end: 878 if ifo == 'G1': 879 trigfile = '%s/%d/%d-%d/%s-OMEGA_TRIGGERS_%s-%.10d-%d.txt'\ 880 % (basedir, t/100000, t, t+triglength, ifo, mask, t, triglength) 881 else: 882 trigfile = '%s/%.10d-%10.d/%s-OMEGA_TRIGGERS_%s-%.10d-%d.txt'\ 883 % (basedir, t, t+triglength, ifo, mask, t, triglength) 884 if intersects(segment(t, t+triglength))\ 885 and (not check_files_exist or isfile(trigfile)): 886 append(from_T050017(trigfile)) 887 t+=triglength-overlap 888 889 cache.sort(key=lambda e: e.path) 890 891 return cache
892 893 # ============================================================================= 894 # Function to generate an omega spectrum online cache 895 # ============================================================================= 896
897 -def omega_spectrum_online_cache(start, end, ifo, **kwargs):
898 899 """ 900 Returns a glue.lal.Cache contatining CacheEntires for all omega online 901 trigger files between the given start and end time for the given ifo. 902 For S6 triggers are only available for each IFO on it's own site cluster. 903 904 Arguments: 905 906 start : [ float | int | LIGOTimeGPS ] 907 GPS start time of requested period 908 end : [ float | int | LIGOTimeGPS ] 909 GPS end time of requested period 910 ifo : [ "H1" | "L1" | "V1" ] 911 IFO 912 """ 913 914 return omega_online_cache(start, end, ifo, mask='SPECTRUM', **kwargs)
915 916 # ============================================================================= 917 # DetChar 'omegadq' cache 918 # ============================================================================= 919
920 -def omega_dq_cache(start,end,ifo):
921 922 """ 923 Arguments: 924 925 start : [ float | int | LIGOTimeGPS ] 926 GPS start time of requested period 927 end : [ float | int | LIGOTimeGPS ] 928 GPS end time of requested period 929 ifo : [ "H1" | "L1" | "V1" ] 930 IFO 931 """ 932 933 # verify host 934 host = getfqdn() 935 ifo_host = {'H1':'ligo-wa','H2':'ligo-wa','L1':'ligo-la'} 936 if not re.search(ifo_host[ifo],host): 937 print >>sys.stderr, "Error: OmegaClustered files are not available for "+\ 938 "IFO="+ifo+" on this host." 939 return LALCache() 940 941 cache = LALCache() 942 basedir = os.path.expanduser('~detchar/public_html/S6/glitch/Wdata') 943 944 # work out days 945 days = gps_day_list(start,end) 946 span = segments.segment(start,end) 947 triglength=86400 948 cache = LALCache() 949 950 for day in days: 951 dirstart = day 952 dirend = day+triglength 953 dirpath = os.path.join(basedir,'%s_%s' % (dirstart,dirend)) 954 trigfile = os.path.join(dirpath,'clusters.txt') 955 if os.path.isfile(trigfile): 956 957 e = LALCacheEntry(ifo,'OMEGADQ',segments.segment(dirstart,dirend),\ 958 os.path.realpath(trigfile)) 959 if span.intersects(e.segment): cache.append(e) 960 961 cache.sort(key=lambda e: e.path) 962 963 return cache
964 965 # ============================================================================== 966 # Class for KWCacheEntry 967 # ============================================================================== 968
969 -class KWCacheEntry(LALCacheEntry):
970 971 _regex = re.compile(r"\A\s*(?P<obs>\S+)\s+(?P<dsc>\S+)\s*\Z") 972
973 - def from_KWfilename(cls, url, coltype = LIGOTimeGPS):
974 """ 975 Parse a URL in the style of KW filenames into a FrameCacheEntry. 976 The KW file name format is, essentially, 977 978 /path/to/start_end/observatory_description.txt 979 980 981 """ 982 983 try: 984 head,tail = os.path.split(url) 985 observatory,description = re.split('_',os.path.splitext(tail)[0],\ 986 maxsplit=1) 987 observatory = observatory 988 start,end = [coltype(t) for t in os.path.basename(head).split('_')] 989 duration = end-start 990 991 segment = segments.segment(start,end) 992 993 except: 994 raise ValueError, "could not convert %s to KWCacheEntry" % repr(url) 995 996 return cls(observatory, description, segment, url)
997 998 from_KWfilename = classmethod(from_KWfilename) 999
1000 - def from_T050017(cls,url,coltype = LIGOTimeGPS):
1001 """ 1002 Redirects to from_KWfilename for KWCacheEntry objects due to KW not 1003 following T50017-00 conventions. 1004 """ 1005 return KWCacheEntry.from_KWfilename(url,coltype=coltype)
1006 1007 from_T050017 = classmethod(from_T050017)
1008 1009 # ============================================================================= 1010 # Function to generate a KW DARM_ERR cache 1011 # ============================================================================= 1012
1013 -def kw_cache(start, end, channel='H1:LSC-DARM_ERR', frequency=None):
1014 1015 """ 1016 Returns a list of KW trigger files between the given start and end 1017 time for the given ifo. For S6 triggers are only available for each IFO on 1018 it's own site cluster. 1019 1020 Arguments: 1021 1022 start : [ float | int | LIGOTimeGPS ] 1023 GPS start time of requested period 1024 end : [ float | int | LIGOTimeGPS ] 1025 GPS end time of requested period 1026 ifo : [ "H1" | "L1" | "V1" ] 1027 IFO 1028 """ 1029 1030 # format channel 1031 if re.match('\w\d:', channel): 1032 ifo, channel = channel.split(':', 1) 1033 #_cchar_regex.sub('_', channel) 1034 else: 1035 raise ValueError("Please give channel in the form \"IFO:CHANNEL-NAME\"") 1036 1037 # verify host 1038 host = getfqdn() 1039 ifo_host = {'H0':'(ligo-wa|ligo\.)', 'H1':'(ligo-wa|ligo\.)',\ 1040 'H2':'ligo-wa', 'L0':'(ligo-la|ligo\.)', 'L1':'(ligo-la|ligo\.)',\ 1041 'V1':'ligo\.'} 1042 if not re.search(ifo_host[ifo],host): 1043 print >>sys.stderr, "Warning: KW files are not available for "+\ 1044 "IFO="+ifo+" on this host." 1045 return LALCache() 1046 1047 cache = LALCache() 1048 if ifo == 'V1': 1049 basedir = os.path.expanduser('~mabizoua/public_html/KW') 1050 elif re.search('ligo\.', host): 1051 basedir = os.path.expanduser('~lindy/public_html/triggers/s6-merged') 1052 else: 1053 basedir = os.path.expanduser('~lindy/public_html/triggers/s6') 1054 1055 # times are numbere from a given start, which for S6 is: 1056 if ifo=='V1': 1057 base = LIGOTimeGPS(938736015) 1058 else: 1059 base = LIGOTimeGPS(938736000) 1060 triglength = 86400 1061 1062 # construct span 1063 span = segments.segment(start, end) 1064 1065 # get fbin 1066 if not frequency: 1067 f = '*_*' 1068 elif frequency.lower() == 'low': 1069 f = '32_2048' 1070 elif frequency.lower() == 'high': 1071 f = '1024_4096' 1072 elif not isinstance(frequency ,str) and len(frequency)==2: 1073 f = '%d_%d' % tuple(frequency) 1074 else: 1075 f = '*_*' 1076 1077 # get times 1078 start_time = int(start-math.fmod(start-base,triglength)) 1079 t = start_time 1080 1081 # loop over time segments constructing file paths and appending to the cache 1082 while t<end: 1083 dirstart = str(t) 1084 dirend = str(t+triglength) 1085 dirpath = os.path.join(basedir,dirstart+'_'+dirend) 1086 trigfile = '%s/%s_%s_%s.trg' % (dirpath, ifo, channel, f) 1087 if '*' in trigfile: 1088 trigfiles = glob.glob(trigfile) 1089 else: 1090 trigfiles = [trigfile] 1091 for trigfile in trigfiles: 1092 if os.path.isfile(trigfile): 1093 e = KWCacheEntry.from_KWfilename(trigfile) 1094 if span.intersects(e.segment): 1095 cache.append(e) 1096 1097 t+=triglength 1098 1099 cache.sort(key=lambda e: e.path) 1100 1101 return cache
1102 1103 # ============================================================================== 1104 # Function to construct list of days from start and end times 1105 # ============================================================================== 1106
1107 -def gps_day_list(start,end):
1108 1109 """ 1110 This script will construct a list of LIGOTimeGPS days encompassing the start and end GPS times. 1111 """ 1112 1113 start=LIGOTimeGPS(start) 1114 1115 start_d = date.XLALUTCToGPS(datetime.datetime(*date.XLALGPSToUTC(start)[:6])\ 1116 .replace(hour=0,minute=0,second=0)\ 1117 .timetuple()) 1118 days = [] 1119 day = start_d 1120 while day<=end: 1121 days.append(day) 1122 day+=86400 1123 1124 return days
1125 1126 # ============================================================================= 1127 # Function to cluster triggers 1128 # ============================================================================= 1129
1130 -def cluster(triggers,params=[('time',1)],rank='snr'):
1131 1132 """ 1133 Cluster the lsctable triggers in the each of the pairs (column,width), 1134 using the rank column. 1135 1136 Arguments: 1137 1138 triggers: glue.ligowl.Table 1139 Table containing trigger columns for clustering 1140 1141 params : list 1142 List object containing (column,width) pairs. Clustering is nested in 1143 order 1144 1145 rank: string 1146 Column by which to rank clusters 1147 test 1148 """ 1149 1150 outtrigs = table.new_from_template(triggers) 1151 1152 j = 0 1153 1154 cols = [p[0] for p in params] 1155 coldata = dict((p, get_column(triggers, p)) for p in cols+[rank]) 1156 # need bandwidth and duration for all burst triggers 1157 if _burst_regex.search(triggers.tableName): 1158 coldata['stop_time'] = get_column(triggers, 'stop_time') +\ 1159 get_column(triggers, 'stop_time_ns')*1e-9 1160 coldata['start_time'] = get_column(triggers, 'start_time') +\ 1161 get_column(triggers, 'start_time_ns')*1e-9 1162 coldata['flow'] = get_column(triggers, 'flow') 1163 coldata['fhigh'] = get_column(triggers, 'fhigh') 1164 # need time for time clustering 1165 elif 'time' in cols: 1166 coldata['stop_time'] = coldata['time'] 1167 coldata['start_time'] = coldata['time'] 1168 1169 for key in coldata.keys(): 1170 coldata[key] = coldata[key].astype(float) 1171 1172 # for each parameter break the clusters generated using the previous 1173 # parameter into smaller clusters by sorting triggers and clustering 1174 # when all parameters have been used, pick the loudest in each cluster 1175 1176 clusters = [range(len(triggers))] 1177 1178 while j < len(params): 1179 1180 col,width = params[j] 1181 1182 newclusters = [] 1183 1184 for k,subcluster in enumerate(clusters): 1185 1186 # sort triggers incluster parameter 1187 subcluster.sort(key=lambda i: coldata[col][i]) 1188 subsubcluster = [] 1189 1190 for i in subcluster: 1191 1192 # get value of param 1193 if col=='time': 1194 valueStop = coldata['stop_time'][i] 1195 valueStart = coldata['start_time'][i] 1196 elif col=='peak_frequency': 1197 valueStop = coldata['fhigh'][i] 1198 valueStart = coldata['flow'][i] 1199 else: 1200 valueStop = coldata[col][i] 1201 valueStart = valueStop 1202 1203 # if subcluster is empty, simply add the first trigger 1204 if not subsubcluster: 1205 subsubcluster = [i] 1206 prevStop = valueStop 1207 prevStart = valueStart 1208 continue 1209 1210 # if current trig is inside width, append to cluster 1211 if (valueStart-prevStop)<width: 1212 subsubcluster.append(i) 1213 1214 # if not the subcluster is complete, append it to list and start again 1215 else: 1216 newclusters.append(subsubcluster) 1217 subsubcluster=[i] 1218 1219 prevStart = valueStart 1220 prevStop = valueStop 1221 1222 # append final subsubcluster 1223 newclusters.append(subsubcluster) 1224 1225 clusters = copy.deepcopy(newclusters) 1226 j += 1 1227 1228 # process clusters 1229 for cluster in clusters: 1230 if len(cluster)==1: 1231 outtrigs.append(copy.deepcopy(triggers[cluster[0]])) 1232 elif len(cluster) > 1: 1233 carray = numpy.asarray(cluster) 1234 cluster.sort(key=lambda i: coldata[rank][i], reverse=True) 1235 t = copy.deepcopy(triggers[cluster[0]]) 1236 # reset burst params for a clustered event 1237 if _burst_regex.search(triggers.tableName): 1238 # record most significant trigger 1239 t.ms_start_time = t.start_time 1240 t.ms_start_time_ns = t.start_time_ns 1241 t.ms_stop_time = t.stop_time 1242 t.ms_stop_time_ns = t.stop_time_ns 1243 t.ms_duration = t.duration 1244 t.ms_bandwidth = t.bandwidth 1245 t.ms_flow = t.flow 1246 t.ms_fhigh = t.fhigh 1247 t.ms_snr = t.snr 1248 # record cluster 1249 start = LIGOTimeGPS(min(coldata['start_time'][carray])) 1250 t.start_time = start.seconds 1251 t.start_time_ns = start.nanoseconds 1252 stop = LIGOTimeGPS(max(coldata['stop_time'][carray])) 1253 t.stop_time = stop.seconds 1254 t.stop_time_ns = stop.nanoseconds 1255 t.duration = float(t.get_stop()-t.get_start()) 1256 t.flow = min(coldata['flow'][carray]) 1257 t.fhigh = max(coldata['fhigh'][carray]) 1258 t.bandwidth = t.fhigh-t.flow 1259 t.tfvolume = t.bandwidth * t.duration 1260 outtrigs.append(t) 1261 1262 # resort trigs in first parameter 1263 outtrigs.sort(key=lambda t: get(t, cols[0])) 1264 1265 return outtrigs
1266 1267 # ============================================================================= 1268 # Compute trigger auto-correlation 1269 # ============================================================================= 1270
1271 -def autocorr(triggers,column='time',timeStep=0.02,timeRange=60):
1272 1273 """ 1274 Compute autocorrelation of lsctable triggers in the each of the pairs 1275 (column,width), using the rank column. 1276 1277 Arguments: 1278 1279 triggers: glue.ligowl.Table 1280 Table containing trigger columns for clustering 1281 1282 column: 1283 On which trigger column to auto-correlate, almost always time 1284 1285 timeStep: 1286 Step (bin width) for the autocorrelation 1287 1288 timeRange: 1289 Longest time to consider for autocorrelation 1290 1291 """ 1292 1293 # time sort triggers before proceeding 1294 get_time = def_get_time(triggers.tableName) 1295 triggers.sort(key=lambda trig: get_time(trig)) 1296 1297 previousTimes = [] 1298 histEdges = numpy.arange(timeStep,timeRange,timeStep); 1299 delayHist = numpy.zeros(int(math.ceil(timeRange/timeStep))) 1300 1301 for trig in triggers: 1302 curTime = trig.peak_time + 1e-9*trig.peak_time_ns 1303 # remove previous times which are beyond the considered timeRange 1304 while len(previousTimes) > 0 and curTime - previousTimes[0] > timeRange: 1305 previousTimes.pop() 1306 for t in previousTimes: 1307 pos = int(math.floor((curTime - t)/timeStep)) 1308 if pos < len(delayHist): 1309 delayHist[pos] += 1 1310 previousTimes.append(curTime) 1311 1312 delayHistFFT = numpy.abs(numpy.fft.fft(delayHist)) 1313 freqBins = numpy.fft.fftfreq(len(delayHist), d=timeStep) 1314 1315 return delayHistFFT, freqBins, delayHist, histEdges
1316 1317 # ============================================================================= 1318 # Get coincidences between two tables 1319 # ============================================================================= 1320
1321 -def get_coincs(table1, table2, dt=1, returnsegs=False, timeshift=0):
1322 1323 """ 1324 Returns the table of those entries in table1 whose time is within +-dt of 1325 an entry in table2. 1326 """ 1327 1328 t1 = get_column(table1, 'time') 1329 t2 = get_column(table2, 'time') 1330 1331 coincsegs = segments.segmentlist(segments.segment(t-dt+timeshift, t+dt+timeshift) for t in t2)\ 1332 .coalesce() 1333 coincsegs.sort() 1334 coinctrigs = table.new_from_template(table1) 1335 coinctrigs.extend(t for i,t in enumerate(table1) if t1[i] in coincsegs) 1336 1337 if returnsegs: 1338 return coinctrigs,coincsegs 1339 else: 1340 return coinctrigs
1341 1342 # ============================================================================= 1343 # Get number of coincidences between two tables 1344 # ============================================================================= 1345
1346 -def get_number_coincs(table1, table2, dt=1, timeshift=0, tabletype='trigger'):
1347 1348 """ 1349 Returns the numbers of entries in table1 whose time is within +-dt of 1350 an entry in table2. 1351 """ 1352 1353 # t1 = get_column(table1, 'time') 1354 # t2 = get_column(table2, 'time') 1355 1356 # coincsegs = segments.segmentlist(segments.segment(t-dt+timeshift, t+dt+timeshift) for t in t2)\ 1357 # .coalesce() 1358 # coincsegs.sort() 1359 # coinctrigs = table.new_from_template(table1) 1360 # coinctrigs.extend(t for i,t in enumerate(table1) if t1[i] in coincsegs) 1361 1362 # return len(coinctrigs) 1363 1364 if tabletype == 'trigger': 1365 time1 = get_column(table1, 'time') 1366 time2 = get_column(table2, 'time') 1367 elif tabletype == 'time': 1368 time1 = table1 1369 time2 = table2 1370 else: 1371 raise ValueError("Unrecognized table type for coincidence number: %s" % tabletype) 1372 1373 time1.sort() 1374 time2.sort() 1375 1376 i2 = 0 1377 ncoinc = 0 1378 for i1,t1 in enumerate(time1): 1379 while i2 < len(time2) and t1 > time2[i2]+timeshift+dt: 1380 i2 = i2 + 1 1381 if i2 >= len(time2): 1382 break 1383 if t1 >= time2[i2]+timeshift-dt: 1384 ncoinc = ncoinc + 1 1385 1386 return ncoinc
1387 1388 # ============================================================================== 1389 # Calculate poisson significance of time coincidences 1390 # ============================================================================== 1391
1392 -def coinc_significance_times(gwtrigtime, auxtrigtime, window=1, livetime=None):
1393 1394 1395 # get livetime 1396 if not livetime: 1397 start = min(gwtrigtime) 1398 end = max(gwtrigtime) 1399 livetime = end-start 1400 1401 # calculate probability of a GW trigger falling within the window 1402 gwprob = len(gwtrigtime) * 2.0 * float(window) / float(livetime) 1403 1404 # calculate mean of Poisson distribution 1405 mu = gwprob * len(auxtrigtime) 1406 1407 # get coincidences 1408 ncoinc = get_number_coincs(gwtrigtime, auxtrigtime, dt=window, tabletype='time') 1409 1410 g = special.gammainc(ncoinc, mu) 1411 1412 # if no coincidences, set significance to zero 1413 if ncoinc<1: 1414 significance = 0 1415 # if significance would blow up, use other formula (ref. hveto_significance.m) 1416 elif g == 0: 1417 significance = -ncoinc * math.log10(mu) + \ 1418 mu * math.log10(math.exp(1)) +\ 1419 special.gammaln(ncoinc + 1) / math.log(10) 1420 # otherwise use the standard formula 1421 else: 1422 significance = -math.log(g, 10) 1423 1424 return significance
1425 1426 1427 # ============================================================================== 1428 # Calculate poisson significance of trigger coincidences 1429 # ============================================================================== 1430
1431 -def coinc_significance(gwtriggers, auxtriggers, window=1, livetime=None,\ 1432 returnsegs=False):
1433 1434 get_time = def_get_time(gwtriggers.tableName) 1435 aux_get_time = def_get_time(auxtriggers.tableName) 1436 1437 # get livetime 1438 if not livetime: 1439 start = min([get_time(t) for t in gwtriggers]) 1440 end = max([get_time(t) for t in gwtriggers]) 1441 livetime = end-start 1442 1443 # calculate probability of a GW trigger falling within the window 1444 gwprob = len(gwtriggers) * 2.0 * float(window) / float(livetime) 1445 1446 # calculate mean of Poisson distribution 1447 mu = gwprob * len(auxtriggers) 1448 1449 # get coincidences 1450 if returnsegs: 1451 coinctriggers,coincsegs = get_coincs(gwtriggers, auxtriggers, dt=window,\ 1452 returnsegs=True) 1453 ncoinc = len(coinctriggers) 1454 else: 1455 ncoinc = get_number_coincs(gwtriggers, auxtriggers, dt=window) 1456 1457 g = special.gammainc(ncoinc, mu) 1458 1459 # if no coincidences, set significance to zero 1460 if ncoinc<1: 1461 significance = 0 1462 # if significance would blow up, use other formula (ref. hveto_significance.m) 1463 elif g == 0: 1464 significance = -ncoinc * math.log10(mu) + \ 1465 mu * math.log10(math.exp(1)) +\ 1466 special.gammaln(ncoinc + 1) / math.log(10) 1467 # otherwise use the standard formula 1468 else: 1469 significance = -math.log(g, 10) 1470 1471 if returnsegs: 1472 return significance,coincsegs 1473 else: 1474 return significance
1475 1476 # ============================================================================= 1477 # Extract column from generic table 1478 # ============================================================================= 1479
1480 -def get_column(lsctable, column):
1481 1482 """ 1483 Extract column from the given glue.ligolw.table lsctable as numpy array. 1484 Tries to use a 'get_col() function if available, otherwise uses 1485 getColumnByName(), treating 'time' as a special case for the known 1486 Burst/Inspiral/Ringdown tables. 1487 """ 1488 1489 # format column 1490 column = str(column).lower() 1491 obj_type = str(type(lsctable)) 1492 1493 # if there's a 'get_' function, use it 1494 if hasattr(lsctable, 'get_%s' % column): 1495 return numpy.asarray(getattr(lsctable, 'get_%s' % column)()) 1496 1497 # treat 'time' as a special case 1498 elif column == 'time' and _trig_regex.search(obj_type): 1499 if _burst_regex.search(obj_type): 1500 tcol = 'peak_time' 1501 if _cbc_regex.search(obj_type): 1502 tcol = 'end_time' 1503 if _ring_regex.search(obj_type): 1504 tcol = 'start_time' 1505 return numpy.asarray(lsctable.getColumnByName(tcol)) + \ 1506 numpy.asarray(lsctable.getColumnByName('%s_ns' % tcol))*10**-9 1507 1508 return numpy.asarray(lsctable.getColumnByName(column))
1509
1510 -def get(self, parameter):
1511 1512 """ 1513 Extract parameter from given ligolw table row object. 1514 """ 1515 1516 # format 1517 parameter = parameter.lower() 1518 obj_type = str(type(self)) 1519 1520 # if there's a 'get_' function, use it 1521 if hasattr(self, 'get_%s' % parameter): 1522 return getattr(self, 'get_%s' % parameter)() 1523 1524 # treat 'time' as a special case 1525 elif parameter == 'time' and _trig_regex.search(obj_type): 1526 if _burst_regex.search(obj_type): 1527 tcol = 'peak_time' 1528 elif _cbc_regex.search(obj_type): 1529 tcol = 'end_time' 1530 elif _ring_regex.search(obj_type): 1531 tcol = 'start_time' 1532 return getattr(self, tcol)+getattr(self, '%s_ns' % tcol)*10**-9 1533 1534 else: 1535 return getattr(self, parameter)
1536 1537 # ============================================================================= 1538 # Veto triggers from a ligolw table 1539 # ============================================================================= 1540
1541 -def veto(self, seglist, inverse=False):
1542 1543 """ 1544 Returns a ligolw table of those triggers outwith the given seglist. 1545 If inverse=True is given, the opposite is returned, i.e. those triggers 1546 within the given seglist. 1547 """ 1548 1549 get_time = def_get_time(self.tableName) 1550 1551 keep = table.new_from_template(self) 1552 if inverse: 1553 keep.extend(t for t in self if float(get_time(t)) in seglist) 1554 else: 1555 keep.extend(t for t in self if float(get_time(t)) not in seglist) 1556 1557 return keep
1558
1559 -def vetoed(self, seglist):
1560 1561 """ 1562 Returns the opposite of veto, i.e. those triggers that lie within the given 1563 seglist. 1564 """ 1565 1566 return veto(self, seglist, inverse=True)
1567 1568 # ============================================================================= 1569 # read triggers from file 1570 # ============================================================================= 1571
1572 -def fromomegafile(fname, start=None, end=None, ifo=None, channel=None,\ 1573 columns=None, virgo=False):
1574 1575 """ 1576 Load triggers from an Omega format text file into a SnglBurstTable object. 1577 Use start and end to restrict the returned triggers, and give ifo and 1578 channel to fill those columns in the table. 1579 1580 If columns is given as a list, only those columns in the table will be 1581 filled. This is advisable to speed up future operations on this table. 1582 1583 Arguments : 1584 1585 fname : file or str 1586 file object or filename path to read with numpy.loadtext 1587 1588 Keyword arguments : 1589 1590 start : float 1591 minimum peak time for returned triggers 1592 end : float 1593 maximum peak time for returned triggers 1594 ifo : str 1595 name of IFO to fill in table 1596 channel : str 1597 name of channel to fill in table 1598 columns : iterable 1599 list of columnnames to populate in table 1600 """ 1601 1602 # set columns 1603 if columns==None: columns = lsctables.SnglBurst.__slots__ 1604 if start or end: 1605 if not start: 1606 start = 0 1607 if not end: 1608 end = numpy.inf 1609 span = segments.segment(start, end) 1610 if 'peak_time' not in columns: columns.append('peak_time') 1611 if 'peak_time_ns' not in columns: columns.append('peak_time_ns') 1612 check_time = True 1613 else: 1614 check_time = False 1615 1616 if 'snr' in columns and not 'amplitude' in columns: 1617 columns.append('amplitude') 1618 1619 # generate table 1620 out = SnglTriggerTable('omega', columns=columns) 1621 1622 # force filename not file object 1623 if hasattr(fname, 'readline'): 1624 fh = fname 1625 else: 1626 fh = open(fname, 'r') 1627 1628 dat = loadtxt(fh) 1629 1630 if not hasattr(fname, 'readline'): 1631 fh.close() 1632 1633 if numpy.shape(dat) == (0,): 1634 return out 1635 1636 if virgo: 1637 start, stop, peak, freq, bandwidth, cln, cle, snr = dat 1638 duration = stop-start 1639 amplitude = snr**2/2 1640 omega_clusters = False 1641 if len(dat)==11: 1642 peak, freq, duration, bandwidth, amplitude, cls, cle, cln, av_freq, av_bandwidth, err_freq = dat 1643 omega_clusters = True 1644 else: 1645 if len(dat)==8: 1646 peak, freq, duration, bandwidth, amplitude, cls, cle, cln = dat 1647 omega_clusters = True 1648 elif len(dat)==5: 1649 peak, freq, duration, bandwidth, amplitude = dat 1650 omega_clusters = False 1651 else: 1652 raise ValueError("Wrong number of columns in omega format file. "\ 1653 "Cannot read.") 1654 av_freq = freq 1655 av_bandwidth = bandwidth 1656 err_freq = av_bandwidth/av_freq 1657 1658 numtrigs = len(peak) 1659 attr_map = dict() 1660 1661 if 'start_time' in columns or 'start_time_ns' in columns: 1662 start = map(LIGOTimeGPS, peak - duration/2) 1663 attr_map['start_time'], attr_map['start_time_ns'] =\ 1664 zip(*[(s.seconds, s.nanoseconds) for s in start]) 1665 if 'stop_time' in columns or 'stop_time_ns' in columns: 1666 stop = map(LIGOTimeGPS, peak + duration/2) 1667 attr_map['stop_time'], attr_map['stop_time_ns'] =\ 1668 zip(*[(s.seconds, s.nanoseconds) for s in stop]) 1669 if 'peak_time' in columns or 'peak_time_ns' in columns: 1670 peak = map(LIGOTimeGPS, peak) 1671 attr_map['peak_time'], attr_map['peak_time_ns'] =\ 1672 zip(*[(s.seconds, s.nanoseconds) for s in peak]) 1673 1674 if 'ms_start_time' in columns or 'ms_start_time_ns' in columns: 1675 ms_start = map(LIGOTimeGPS, peak-duration/2) 1676 attr_map['ms_start_time'], attr_map['ms_start_time_ns'] =\ 1677 zip(*[(s.seconds, s.nanoseconds) for s in ms_start]) 1678 if 'ms_stop_time' in columns or 'ms_stop_time_ns' in columns: 1679 ms_stop = map(LIGOTimeGPS, peak+duration/2) 1680 attr_map['ms_stop_time'], attr_map['ms_stop_time_ns'] =\ 1681 zip(*[(s.seconds, s.nanoseconds) for s in ms_stop]) 1682 1683 if 'central_freq' in columns: attr_map['central_freq'] = freq 1684 if 'peak_frequency' in columns: attr_map['peak_frequency'] = av_freq 1685 if 'peak_frequency_error' in columns: attr_map['peak_frequency_error'] = err_freq 1686 if 'bandwidth' in columns: attr_map['bandwidth'] = av_bandwidth 1687 if 'ms_bandwidth' in columns: attr_map['ms_bandwidth'] = bandwidth 1688 if 'flow' in columns: attr_map['flow'] = freq-bandwidth/2 1689 if 'fhigh' in columns: attr_map['fhigh'] = freq+bandwidth/2 1690 if 'ms_flow' in columns: attr_map['ms_flow'] = freq-bandwidth/2 1691 if 'ms_fhigh' in columns: attr_map['ms_fhigh'] = freq+bandwidth/2 1692 1693 if 'duration' in columns: attr_map['duration'] = duration 1694 if 'ms_duration' in columns: attr_map['ms_duration'] = duration 1695 if 'snr' in columns: attr_map['snr'] = (2*amplitude)**(1/2) 1696 if 'ms_snr' in columns: attr_map['ms_snr'] = (2*amplitude)**(1/2) 1697 1698 if 'cluster_size' in columns or 'param_one_value' in columns: 1699 attr_map['param_one_name'] = ['cluster_size'] * numtrigs 1700 if omega_clusters: 1701 attr_map['param_one_value'] = cls 1702 else: 1703 attr_map['param_one_value'] = [numpy.NaN] * numtrigs 1704 if 'cluster_norm_energy' in columns or 'param_two_value' in columns: 1705 attr_map['param_two_name'] = ['cluster_norm_energy'] * numtrigs 1706 if omega_clusters: 1707 attr_map['param_two_value'] = cls 1708 else: 1709 attr_map['param_two_value'] = [numpy.NaN] * numtrigs 1710 if 'cluster_size' in columns or 'param_three_value' in columns: 1711 attr_map['param_three_name'] = ['cluster_number'] * numtrigs 1712 if omega_clusters: 1713 attr_map['param_three_value'] = cls 1714 else: 1715 attr_map['param_three_value'] = [numpy.NaN] * numtrigs 1716 1717 cols = attr_map.keys() 1718 append = out.append 1719 for i in range(numtrigs): 1720 t = lsctables.SnglBurst() 1721 for c in cols: setattr(t, c, attr_map[c][i]) 1722 if not check_time or (check_time and float(t.get_peak()) in span): 1723 append(t) 1724 1725 return out
1726
1727 -def fromkwfile(fname, start=None, end=None, ifo=None, channel=None,\ 1728 columns=None):
1729 1730 """ 1731 Load triggers from a KW format text file into a SnglBurstTable object. 1732 Use start and end to restrict the returned triggers, and give ifo and 1733 channel to fill those columns in the table. 1734 1735 If columns is given as a list, only those columns in the table will be 1736 filled. This is advisable to speed up future operations on this table. 1737 1738 Arguments : 1739 1740 fname : file or str 1741 file object or filename path to read with numpy.loadtext 1742 1743 Keyword arguments : 1744 1745 start : float 1746 minimum peak time for returned triggers 1747 end : float 1748 maximum peak time for returned triggers 1749 ifo : str 1750 name of IFO to fill in table 1751 channel : str 1752 name of channel to fill in table 1753 columns : iterable 1754 list of columnnames to populate in table 1755 """ 1756 1757 # set columns 1758 if columns==None: columns = lsctables.SnglBurst.__slots__ 1759 if start or end: 1760 if not start: 1761 start = 0 1762 if not end: 1763 end = numpy.inf 1764 span = segments.segment(start, end) 1765 if 'peak_time' not in columns: columns.append('peak_time') 1766 if 'peak_time_ns' not in columns: columns.append('peak_time_ns') 1767 check_time = True 1768 else: 1769 check_time = False 1770 1771 # generate table 1772 out = SnglTriggerTable('kw', columns=columns) 1773 1774 # force filename not file object 1775 if hasattr(fname, 'readline'): 1776 fh = fname 1777 else: 1778 fh = open(fname, 'r') 1779 1780 # load data from file 1781 dat = loadtxt(fh, usecols=[0,1,2,3,4,5,6,7]) 1782 1783 # close file if we opened it 1784 if not hasattr(fname, 'readline'): 1785 fh.close() 1786 1787 if numpy.shape(dat) == (0,): 1788 return out 1789 1790 if len(dat)==8: 1791 st, stop, peak, freq, energy, amplitude, n_pix, sig = dat 1792 else: 1793 raise ValueError("Wrong number of columns in KW format file. "\ 1794 "Cannot read.") 1795 1796 numtrigs = len(peak) 1797 1798 attr_map = dict() 1799 1800 if 'duration' in columns: attr_map['duration'] = stop-st 1801 if 'ms_duration' in columns: attr_map['ms_duration'] = stop-st 1802 1803 if 'start_time' in columns or 'start_time_ns' in columns: 1804 start = map(LIGOTimeGPS, st) 1805 attr_map['start_time'], attr_map['start_time_ns'] =\ 1806 zip(*[(s.seconds, s.nanoseconds) for s in start]) 1807 if 'stop_time' in columns or 'stop_time_ns' in columns: 1808 stop = map(LIGOTimeGPS, stop) 1809 attr_map['stop_time'], attr_map['stop_time_ns'] =\ 1810 zip(*[(s.seconds, s.nanoseconds) for s in stop]) 1811 if 'peak_time' in columns or 'peak_time_ns' in columns: 1812 peak = map(LIGOTimeGPS, peak) 1813 attr_map['peak_time'], attr_map['peak_time_ns'] =\ 1814 zip(*[(s.seconds, s.nanoseconds) for s in peak]) 1815 1816 if 'ms_start_time' in columns or 'ms_start_time_ns' in columns: 1817 ms_start = map(LIGOTimeGPS, st) 1818 attr_map['ms_start_time'], attr_map['ms_start_time_ns'] =\ 1819 zip(*[(s.seconds, s.nanoseconds) for s in ms_start]) 1820 if 'ms_stop_time' in columns or 'ms_stop_time_ns' in columns: 1821 ms_stop = map(LIGOTimeGPS, stop) 1822 attr_map['ms_stop_time'], attr_map['ms_stop_time_ns'] =\ 1823 zip(*[(s.seconds, s.nanoseconds) for s in ms_stop]) 1824 1825 if 'central_freq' in columns: attr_map['central_freq'] = freq 1826 if 'peak_frequency' in columns: attr_map['peak_frequency'] = freq 1827 if 'bandwidth' in columns: attr_map['bandwidth'] = numpy.zeros(len(freq)) 1828 if 'ms_bandwidth' in columns: attr_map['ms_bandwidth'] = attr_map['bandwidth'] 1829 if 'flow' in columns: attr_map['flow'] = freq 1830 if 'fhigh' in columns: attr_map['fhigh'] = freq 1831 if 'ms_flow' in columns: attr_map['ms_flow'] = freq 1832 if 'ms_fhigh' in columns: attr_map['ms_fhigh'] = freq 1833 1834 if 'duration' in columns: attr_map['duration'] = stop-st 1835 if 'ms_duration' in columns: attr_map['ms_duration'] = stop-st 1836 if 'snr' in columns: attr_map['snr'] = (amplitude-n_pix)**(1/2) 1837 if 'ms_snr' in columns: attr_map['ms_snr'] = (amplitude-n_pix)**(1/2) 1838 if 'confidence' in columns: attr_map['confidence'] = sig 1839 1840 if 'n_pix' in columns or 'param_one_value' in columns: 1841 attr_map['param_one_name'] = ['n_pix'] * numtrigs 1842 attr_map['param_one_value'] = n_pix 1843 """ 1844 if 'significance' in columns or 'param_two_value' in columns: 1845 attr_map['param_two_name'] = ['significance'] * numtrigs 1846 attr_map['param_two_value'] = sig 1847 """ 1848 1849 cols = attr_map.keys() 1850 append = out.append 1851 for i in range(numtrigs): 1852 t = lsctables.SnglBurst() 1853 for c in cols: setattr(t, c, attr_map[c][i]) 1854 if ifo!=None: 1855 t.ifo = ifo 1856 if channel!=None: 1857 t.channel = channel 1858 if not check_time or (check_time and t.get_peak() in span): 1859 append(t) 1860 1861 return out
1862
1863 -def fromomegaspectrumfile(fname, start=None, end=None, ifo=None, channel=None,\ 1864 columns=None):
1865 1866 """ 1867 Load triggers from an OmegaSpectrum format text file into a SnglBurstTable 1868 object. 1869 Use start and end to restrict the returned triggers, and give ifo and 1870 channel to fill those columns in the table. 1871 1872 If columns is given as a list, only those columns in the table will be 1873 filled. This is advisable to speed up future operations on this table. 1874 1875 Arguments : 1876 1877 fname : file or str 1878 file object or filename path to read with numpy.loadtext 1879 1880 Keyword arguments : 1881 1882 start : float 1883 minimum peak time for returned triggers 1884 end : float 1885 maximum peak time for returned triggers 1886 ifo : str 1887 name of IFO to fill in table 1888 channel : str 1889 name of channel to fill in table 1890 columns : iterable 1891 list of columnnames to populate in table 1892 """ 1893 1894 # set columns 1895 if columns==None: columns = lsctables.SnglBurst.__slots__ 1896 if start or end: 1897 if not start: 1898 start = 0 1899 if not end: 1900 end = numpy.inf 1901 span = segments.segment(start, end) 1902 if 'peak_time' not in columns: columns.append('peak_time') 1903 if 'peak_time_ns' not in columns: columns.append('peak_time_ns') 1904 check_time = True 1905 else: 1906 check_time = False 1907 1908 # generate table 1909 out = SnglTriggerTable('omegaspectrum', columns=columns) 1910 1911 # force filename not file object 1912 if hasattr(fname, 'readline'): 1913 fh = fname 1914 else: 1915 fh = open(fname, 'r') 1916 1917 # load data from file 1918 dat = loadtxt(fh) 1919 1920 # close file if we opened it 1921 if not hasattr(fname, 'readline'): 1922 fh.close() 1923 1924 if numpy.shape(dat) == (0,): 1925 return out 1926 1927 if len(dat)==3: 1928 peak, freq, amplitude = dat 1929 elif len(dat)==4: 1930 peak, freq, amplitude, chisq = dat 1931 else: 1932 raise ValueError("Wrong number of columns in omega spectrum format file. "\ 1933 "Cannot read.") 1934 numtrigs = len(peak) 1935 1936 attr_map = dict() 1937 1938 if 'peak_time' in columns or 'peak_time_ns' in columns: 1939 peak = map(LIGOTimeGPS, peak) 1940 attr_map['peak_time'], attr_map['peak_time_ns'] =\ 1941 zip(*[(s.seconds, s.nanoseconds) for s in peak]) 1942 1943 if 'central_freq' in columns: attr_map['central_freq'] = freq 1944 if 'peak_frequency' in columns: attr_map['peak_frequency'] = freq 1945 if 'amplitude' in columns: attr_map['amplitude'] = amplitude 1946 if 'snr' in columns: attr_map['snr'] = amplitude**(1/2) 1947 if 'chisq' in columns: attr_map['chisq'] = chisq**(0.25)/attr_map['snr'] 1948 1949 cols = attr_map.keys() 1950 append = out.append 1951 for i in range(numtrigs): 1952 t = lsctables.SnglBurst() 1953 for c in cols: setattr(t, c, attr_map[c][i]) 1954 if ifo!=None: 1955 t.ifo = ifo 1956 if channel!=None: 1957 t.channel = channel 1958 if not check_time or (check_time and t.get_peak() in span): 1959 append(t) 1960 1961 return out
1962
1963 -def fromomegadqfile(fname, start=None, end=None, ifo=None, channel=None,\ 1964 columns=None):
1965 1966 """ 1967 Load triggers from an OmegaDQ format text file into a SnglBurstTable object. 1968 Use start and end to restrict the returned triggers, and give ifo and 1969 channel to fill those columns in the table. 1970 1971 If columns is given as a list, only those columns in the table will be 1972 filled. This is advisable to speed up future operations on this table. 1973 1974 Arguments : 1975 1976 fname : file or str 1977 file object or filename path to read with numpy.loadtext 1978 1979 Keyword arguments : 1980 1981 start : float 1982 minimum peak time for returned triggers 1983 end : float 1984 maximum peak time for returned triggers 1985 ifo : str 1986 name of IFO to fill in table 1987 channel : str 1988 name of channel to fill in table 1989 columns : iterable 1990 list of columnnames to populate in table 1991 """ 1992 1993 # set columns 1994 if columns==None: columns = lsctables.SnglBurst.__slots__ 1995 if start or end: 1996 if not start: 1997 start = 0 1998 if not end: 1999 end = numpy.inf 2000 span = segments.segment(start, end) 2001 if 'peak_time' not in columns: columns.append('peak_time') 2002 if 'peak_time_ns' not in columns: columns.append('peak_time_ns') 2003 check_time = True 2004 else: 2005 check_time = False 2006 2007 # generate table 2008 out = SnglTriggerTable('omegadq', columns=columns) 2009 2010 # force filename not file object 2011 if hasattr(fname, 'readline'): 2012 fh = fname 2013 else: 2014 fh = open(fname, 'r') 2015 2016 # load data from file 2017 dat = loadtxt(fh) 2018 2019 # close file if we opened it 2020 if not hasattr(fname, 'readline'): 2021 fh.close() 2022 2023 if numpy.shape(dat) == (0,): 2024 return out 2025 2026 if len(dat)==13: 2027 st, stop, peak, flow, fhigh, nev, ms_start, ms_stop, ms_flow, ms_fhigh,\ 2028 cls, cle, ms_cle = dat 2029 else: 2030 raise ValueError("Wrong number of columns in OmegaDQ format file. "\ 2031 "Cannot read.") 2032 numtrigs = len(peak) 2033 2034 attr_map = dict() 2035 2036 if 'duration' in columns: attr_map['duration'] = stop-st 2037 if 'ms_duration' in columns: attr_map['ms_duration'] = ms_stop-ms_start 2038 2039 if 'start_time' in columns or 'start_time_ns' in columns: 2040 start = map(LIGOTimeGPS, st) 2041 attr_map['start_time'], attr_map['start_time_ns'] =\ 2042 zip(*[(s.seconds, s.nanoseconds) for s in start]) 2043 if 'stop_time' in columns or 'stop_time_ns' in columns: 2044 stop = map(LIGOTimeGPS, stop) 2045 attr_map['stop_time'], attr_map['stop_time_ns'] =\ 2046 zip(*[(s.seconds, s.nanoseconds) for s in stop]) 2047 if 'peak_time' in columns or 'peak_time_ns' in columns: 2048 peak = map(LIGOTimeGPS, peak) 2049 attr_map['peak_time'], attr_map['peak_time_ns'] =\ 2050 zip(*[(s.seconds, s.nanoseconds) for s in peak]) 2051 2052 if 'ms_start_time' in columns or 'ms_start_time_ns' in columns: 2053 ms_start = map(LIGOTimeGPS, ms_start) 2054 attr_map['ms_start_time'], attr_map['ms_start_time_ns'] =\ 2055 zip(*[(s.seconds, s.nanoseconds) for s in ms_start]) 2056 if 'ms_stop_time' in columns or 'ms_stop_time_ns' in columns: 2057 ms_stop = map(LIGOTimeGPS, ms_stop) 2058 attr_map['ms_stop_time'], attr_map['ms_stop_time_ns'] =\ 2059 zip(*[(s.seconds, s.nanoseconds) for s in ms_stop]) 2060 2061 if 'flow' in columns: attr_map['flow'] = flow 2062 if 'fhigh' in columns: attr_map['fhigh'] = fhigh 2063 if 'bandwidth' in columns: attr_map['bandwidth'] = fhigh-flow 2064 if 'ms_flow' in columns: attr_map['ms_flow'] = flow 2065 if 'ms_fhigh' in columns: attr_map['ms_fhigh'] = fhigh 2066 if 'ms_bandwidth' in columns: attr_map['ms_bandwidth'] = ms_fhigh-ms_flow 2067 2068 if 'central_freq' in columns: attr_map['central_freq'] = (flow+fhigh)/2 2069 if 'peak_frequency' in columns: attr_map['peak_frequency'] = (flow+fhigh)/2 2070 2071 if 'snr' in columns: attr_map['snr'] = cle**(1/2) 2072 if 'ms_snr' in columns: attr_map['ms_snr'] = ms_cle**(1/2) 2073 2074 if 'cluster_size' in columns or 'param_one_value' in columns: 2075 attr_map['param_one_name'] = ['cluster_size'] * numtrigs 2076 attr_map['param_one_value'] = cls 2077 if 'cluster_number' in columns or 'param_two_value' in columns: 2078 attr_map['param_two_name'] = ['cluster_number'] * numtrigs 2079 attr_map['param_two_value'] = nev 2080 2081 cols = attr_map.keys() 2082 append = out.append 2083 for i in range(numtrigs): 2084 t = lsctables.SnglBurst() 2085 for c in cols: setattr(t, c, attr_map[c][i]) 2086 if ifo!=None: 2087 t.ifo = ifo 2088 if channel!=None: 2089 t.channel = channel 2090 if not check_time or (check_time and t.get_peak() in span): 2091 append(t) 2092 2093 return out
2094
2095 -def fromhacrfile(fname, start=None, end=None, ifo=None, channel=None,\ 2096 columns=None):
2097 2098 """ 2099 Load triggers from a HACR format text file into a SnglBurstTable object. 2100 Use start and end to restrict the returned triggers, and give ifo and 2101 channel to fill those columns in the table. 2102 2103 If columns is given as a list, only those columns in the table will be 2104 filled. This is advisable to speed up future operations on this table. 2105 2106 Arguments : 2107 2108 fname : file or str 2109 file object or filename path to read with numpy.loadtext 2110 2111 Keyword arguments : 2112 2113 start : float 2114 minimum peak time for returned triggers 2115 end : float 2116 maximum peak time for returned triggers 2117 ifo : str 2118 name of IFO to fill in table 2119 channel : str 2120 name of channel to fill in table 2121 columns : iterable 2122 list of columnnames to populate in table 2123 """ 2124 2125 # set columns 2126 if columns==None: columns = lsctables.SnglBurst.__slots__ 2127 if start or end: 2128 if not start: 2129 start = 0 2130 if not end: 2131 end = numpy.inf 2132 span = segments.segment(start, end) 2133 if 'peak_time' not in columns: columns.append('peak_time') 2134 if 'peak_time_ns' not in columns: columns.append('peak_time_ns') 2135 check_time = True 2136 else: 2137 check_time = False 2138 2139 # generate table 2140 out = SnglTriggerTable('hacr', columns=columns) 2141 2142 # force filename not file object 2143 if hasattr(fname, 'readline'): 2144 fh = fname 2145 else: 2146 fh = open(fname, 'r') 2147 2148 # load data from file 2149 dat = loadtxt(fh) 2150 2151 # close file if we opened it 2152 if not hasattr(fname, 'readline'): 2153 fh.close() 2154 2155 if numpy.shape(dat) == (0,): 2156 return out 2157 elif numpy.shape(dat) == (8,): 2158 peak_time, peak_time_offset, freq, bandwidth, duration, n_pix, snr,\ 2159 totPower = map(lambda n: numpy.asarray([n]), dat) 2160 elif len(dat)==8: 2161 peak_time, peak_time_offset, freq, bandwidth, duration, n_pix, snr,\ 2162 totPower = dat 2163 else: 2164 raise ValueError("Wrong number of columns in HACR format file. "\ 2165 "Cannot read.") 2166 2167 numtrigs = len(peak_time) 2168 2169 if columns==None: columns = lsctables.SnglBurst.__slots__ 2170 if start or end: 2171 if not start: 2172 start = 0 2173 if not end: 2174 end = numpy.inf 2175 span = segments.segment(start, end) 2176 if 'peak_time' not in columns: columns.append('peak_time') 2177 if 'peak_time_ns' not in columns: columns.append('peak_time_ns') 2178 check_time = True 2179 else: 2180 check_time = False 2181 2182 attr_map = dict() 2183 2184 peak = peak_time+peak_time_offset 2185 if 'start_time' in columns or 'start_time_ns' in columns: 2186 start = map(LIGOTimeGPS, peak-duration/2) 2187 attr_map['start_time'], attr_map['start_time_ns'] =\ 2188 zip(*[(s.seconds, s.nanoseconds) for s in start]) 2189 if 'stop_time' in columns or 'stop_time_ns' in columns: 2190 stop = map(LIGOTimeGPS, peak+duration/2) 2191 attr_map['stop_time'], attr_map['stop_time_ns'] =\ 2192 zip(*[(s.seconds, s.nanoseconds) for s in stop]) 2193 if 'peak_time' in columns or 'peak_time_ns' in columns: 2194 peak = map(LIGOTimeGPS, peak) 2195 attr_map['peak_time'], attr_map['peak_time_ns'] =\ 2196 zip(*[(s.seconds, s.nanoseconds) for s in peak]) 2197 2198 if 'ms_start_time' in columns or 'ms_start_time_ns' in columns: 2199 ms_start = map(LIGOTimeGPS, peak-duration/2) 2200 attr_map['ms_start_time'], attr_map['ms_start_time_ns'] =\ 2201 zip(*[(s.seconds, s.nanoseconds) for s in ms_start]) 2202 if 'ms_stop_time' in columns or 'ms_stop_time_ns' in columns: 2203 ms_stop = map(LIGOTimeGPS, peak+duration/2) 2204 attr_map['ms_stop_time'], attr_map['ms_stop_time_ns'] =\ 2205 zip(*[(s.seconds, s.nanoseconds) for s in ms_stop]) 2206 2207 if 'duration' in columns: attr_map['duration'] = duration 2208 if 'ms_duration' in columns: attr_map['ms_duration'] = duration 2209 if 'central_freq' in columns: attr_map['central_freq'] = freq 2210 if 'peak_frequency' in columns: attr_map['peak_frequency'] = freq 2211 if 'peak_frequency_error' in columns: attr_map['peak_frequency_error'] = bandwidth/freq 2212 2213 if 'flow' in columns: attr_map['flow'] = freq-bandwidth/2 2214 if 'fhigh' in columns: attr_map['fhigh'] = freq+bandwidth/2 2215 if 'bandwidth' in columns: attr_map['bandwidth'] = bandwidth 2216 if 'ms_flow' in columns: attr_map['ms_flow'] = freq-bandwidth/2 2217 if 'ms_fhigh' in columns: attr_map['ms_fhigh'] = freq+bandwidth/2 2218 if 'ms_bandwidth' in columns: attr_map['ms_bandwidth'] = bandwidth 2219 2220 if 'snr' in columns: attr_map['snr'] = snr 2221 if 'ms_snr' in columns: attr_map['ms_snr'] = snr 2222 2223 if 'peak_time_offset' in columns or 'param_one_value' in columns\ 2224 or 'peak_time_ns' in columns: 2225 attr_map['param_one_name'] = ['peak_time_offset'] * numtrigs 2226 attr_map['param_one_value'] = peak_time_offset 2227 if 'numPixels' in columns or 'param_two_value' in columns: 2228 attr_map['param_two_name'] = ['numPixels'] * numtrigs 2229 attr_map['param_two_value'] = n_pix 2230 if 'totPower' in columns or 'param_three_value' in columns: 2231 attr_map['param_three_name'] = ['totPower'] * numtrigs 2232 attr_map['param_three_value'] = totPower 2233 2234 cols = attr_map.keys() 2235 append = out.append 2236 for i in range(numtrigs): 2237 t = lsctables.SnglBurst() 2238 for c in cols: setattr(t, c, attr_map[c][i]) 2239 if ifo!=None: 2240 t.ifo = ifo 2241 if channel!=None: 2242 t.channel = channel 2243 if not check_time or (check_time and t.get_peak() in span): 2244 append(t) 2245 2246 return out
2247
2248 -def fromihopefile(fname, start=None, end=None, ifo=None, channel=None,\ 2249 columns=None):
2250 2251 """ 2252 Load triggers from an iHope format CSV file into a SnglInspiralTable object. 2253 Use start and end to restrict the returned triggers, and give ifo and 2254 channel to fill those columns in the table. 2255 2256 If columns is given as a list, only those columns in the table will be 2257 filled. This is advisable to speed up future operations on this table. 2258 2259 Arguments : 2260 2261 fname : file or str 2262 file object or filename path to read with numpy.loadtext 2263 2264 Keyword arguments : 2265 2266 start : float 2267 minimum peak time for returned triggers 2268 end : float 2269 maximum peak time for returned triggers 2270 ifo : str 2271 name of IFO to fill in table 2272 channel : str 2273 name of channel to fill in table 2274 columns : iterable 2275 list of columnnames to populate in table 2276 """ 2277 2278 # get columns 2279 def_cols = dict(enumerate(['end_time','end_time_ns','ifo','snr','mass1',\ 2280 'mass2', 'mtotal','eta','event_duration',\ 2281 'template_duration','eff_distance','chisq',\ 2282 'chisq_dof','bank_chisq','bank_chisq_dof',\ 2283 'cont_chisq','cont_chisq_dof'])) 2284 2285 if columns==None: columns = lsctables.SnglInspiral.__slots__ 2286 if start or end: 2287 if not start: 2288 start = 0 2289 if not end: 2290 end = numpy.inf 2291 span = segments.segment(start, end) 2292 if 'end_time' not in columns: columns.append('end_time') 2293 if 'end_time_ns' not in columns: columns.append('end_time_ns') 2294 check_time = True 2295 else: 2296 check_time = False 2297 2298 usecols = [t for t in def_cols if def_cols[t] in columns] 2299 2300 # force filename not file object 2301 if hasattr(fname, 'readline'): 2302 fh = fname 2303 else: 2304 fh = open(fname, 'r') 2305 2306 # load data from file 2307 dat = loadtxt(fh, usecols) 2308 2309 # close file if we opened it 2310 if not hasattr(fname, 'readlin'): 2311 fh.close() 2312 2313 if usecols: 2314 numtrigs = len(dat[0]) 2315 else: 2316 numtrigs = 0 2317 2318 # generate table 2319 out = SnglTriggerTable('ihope', columns=columns) 2320 2321 cols = numpy.arange(len(dat)) 2322 append = out.append 2323 for i in range(numtrigs): 2324 t = lsctables.SnglInspiral() 2325 for c in cols: 2326 setattr(t, def_cols[usecols[c]],\ 2327 int(data[c][i]) if re.search('time', def_cols[usecols[c]])\ 2328 else data[c][i]) 2329 if ifo!=None: 2330 t.ifo = ifo 2331 if channel!=None: 2332 t.channel = channel 2333 if not check_time or (check_time and t.get_peak() in span): 2334 append(t) 2335 2336 return out
2337 2338 # ============================================================================== 2339 # Time shift trigger table 2340 # ============================================================================== 2341
2342 -def time_shift(lsctable, dt=1):
2343 2344 """ 2345 Time shift lsctable by time dt. 2346 """ 2347 2348 out = table.new_from_template(lsctable) 2349 get_time = def_get_time(lsctable.tableName) 2350 2351 for t in lsctable: 2352 t2 = copy.deepcopy(t) 2353 if re.search('inspiral', lsctable.tableName, re.I): 2354 t2.set_end(t2.get_end()+dt) 2355 elif re.search('burst', lsctable.tableName, re.I): 2356 t2.set_peak(t2.get_peak()+dt) 2357 elif re.search('ringdown', lsctable.tableName, re.I): 2358 t2.set_start(t2.get_start()+dt) 2359 out.append(t2) 2360 2361 return out
2362 2363 # ============================================================================= 2364 # Read file 2365 # ============================================================================= 2366
2367 -def loadtxt(fh, usecols=None):
2368 2369 """ 2370 Stripped down version of numpy.loadtxt to work with empty files. 2371 """ 2372 2373 _comment = re.compile('[#%]') 2374 _delim = re.compile('[\t\,\s]+') 2375 output = [] 2376 nVals = 0 2377 for i,line in enumerate(fh): 2378 if _comment.match(line): continue 2379 vals = _delim.split(line.rstrip()) 2380 if nVals>0 and len(vals) != nVals: 2381 print "Warning, line %d of file %s was skipped, uncorrect column number" % (i, fh) 2382 continue 2383 nVals = len(vals) 2384 if usecols is not None: 2385 output.append(tuple(map(float, [vals[j] for j in usecols]))) 2386 else: 2387 output.append(tuple(map(float, vals))) 2388 return numpy.squeeze(numpy.array(output, float)).T
2389