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

Source Code for Module pylal.followup_missed

  1  # Copyright (C) 2006  Alexander Dietz 
  2  # 
  3  # This program is free software; you can redistribute it and/or modify it 
  4  # under the terms of the GNU General Public License as published by the 
  5  # Free Software Foundation; either version 2 of the License, or (at your 
  6  # option) any later version. 
  7  # 
  8  # This program is distributed in the hope that it will be useful, but 
  9  # WITHOUT ANY WARRANTY; without even the implied warranty of 
 10  # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General 
 11  # Public License for more details. 
 12  # 
 13  # You should have received a copy of the GNU General Public License along 
 14  # with this program; if not, write to the Free Software Foundation, Inc., 
 15  # 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA. 
 16   
 17  __author__ = "Darren Woods and Stephen Fairhurst <sfairhurs@gravity.phys.uwm.edu>" 
 18  __prog__ = "followup_missed.py" 
 19  __title__ = "Followup missed injections" 
 20   
 21  import os, sys, copy 
 22  from math import sqrt, pi 
 23   
 24  from pylab import rcParams, fill, figtext, figure, plot, axes, axis, xlabel, ylabel, title, close, grid, legend 
 25   
 26  from pylal import SnglInspiralUtils 
 27  from pylal import InspiralUtils 
 28  from pylal import SimInspiralUtils 
 29  from pylal import CoincInspiralUtils 
 30  from pylal import SearchSummaryUtils 
 31  from pylal import git_version 
 32  from glue import lal 
 33  from glue import markup 
 34  from glue import pipeline 
 35  from glue import segments 
 36  from glue import segmentsUtils 
 37  from glue.markup import oneliner as extra 
 38  from glue.ligolw import table 
 39  from glue.ligolw import lsctables 
 40  from glue.ligolw import utils 
 41  import numpy 
 42   
 43   
 44  ########################################################## 
45 -class FollowupMissed:
46 """ 47 This defines a class for followup missed injections or follow up any 48 injections in general. This class essentially creates time-series of triggers 49 around the time of a missed injection for the various steps of the pipeline. 50 51 Usage: 52 53 # first need to initialize the class 54 followup = followup_missed.FollowupMissed( cache, opts) 55 56 # later, one can do a followup of a missed injection 'inj', which returns the filename 57 # of the html file containing the tables and pictures created 58 followuphtml = followup.followup( inj, ifo ) 59 """ 60 61 62 # -----------------------------------------------------
63 - def __init__(self, cache, opts):
64 """ 65 Initialize this class and sets up all the cache files. 66 @param cache: The cache of all files 67 @param opts: The 'opts' structure from the main code 68 """ 69 rcParams.update({'text.usetex': False}) 70 71 # Check all the required options 72 option_list = ['verbose','followup_exttrig','output_path',\ 73 'followup_time_window','followup_tag','prefix',\ 74 'suffix','figure_resolution'] 75 for option in option_list: 76 if not hasattr(opts, option): 77 raise "Error: The following parameter is required in the "\ 78 "opts structure: ", option 79 80 81 # defining some colors and the stages 82 self.colors = {'H1':'r','H2':'b','L1':'g','V1':'m','G1':'c'} 83 self.stageLabels = ['INSPIRAL_FIRST', 'THINCA_FIRST',\ 84 'INSPIRAL_SECOND', 'THINCA_SECOND'] 85 self.orderLabels = copy.deepcopy(self.stageLabels) 86 self.orderLabels.extend( [ 'THINCA_SECOND_CAT_1','THINCA_SECOND_CAT_2', \ 87 'THINCA_SECOND_CAT_3','THINCA_SECOND_CAT_4'] ) 88 89 # set arguments from the options 90 self.opts = opts 91 self.verbose = opts.verbose 92 self.exttrig = opts.followup_exttrig 93 self.output_path = opts.output_path 94 self.time_window = opts.followup_time_window 95 96 # default value for the injection time-window. This value 97 # will be taken later from a processParams table (if possible) 98 self.injection_window = 0.050 99 100 # initialize a list of images created 101 self.fnameList = [] 102 103 # counter for the followups 104 self.number = 0 105 106 # for the estimated distances 107 self.flow = opts.followup_flow 108 self.spectrum = createSpectrum( self.flow, \ 109 sampleRate = 4096, \ 110 nPoints = 1048576) 111 112 if self.verbose: 113 print "\nStarting initializing the Followup class..." 114 115 # setting the caches 116 if opts.followup_tag: 117 self.cache = cache.sieve(description = opts.followup_tag) 118 else: 119 self.cache = cache 120 121 # retrieving the caches for the different stages 122 self.triggerCache = {} 123 for stage in self.stageLabels: 124 pattern = stage 125 self.triggerCache[stage] = self.cache.sieve(description=pattern) 126 if self.opts.verbose: 127 print "%d files found for stage %s" % (len(self.triggerCache[stage]), stage) 128 129 # sieve the injections 130 self.injectionCache = self.cache.sieve(description = "INJECTION").\ 131 sieve(ifos='HL') 132 133 # generate a dictionary based on the event-ID 134 self.injections=dict() 135 for file, entry in zip(self.injectionCache.pfnlist(), self.injectionCache): 136 injection_id = self.get_injection_id(cache_entry = entry) 137 self.injections[injection_id] = SimInspiralUtils.\ 138 ReadSimInspiralFromFiles( [file], verbose=False ) 139 if self.verbose: 140 print "parsing of cache files for the Followup class done..." 141 142 143 144 # get the process params table from one of the COIRE files 145 coire_file = self.cache.sieve(description = "FOUND").checkfilesexist()[0].pfnlist()[0] 146 try: 147 doc = SearchSummaryUtils.ReadTablesFromFiles([coire_file],[ lsctables.ProcessParamsTable]) 148 process_params = lsctables.ProcessParamsTable.get_table(doc) 149 except IOError: 150 sys.stderr.write("ERROR (IOError) while reading process_params table from file %s. "\ 151 "Does this file exist and does it contain a search_summary table?\n" %(coire_file)) 152 raise 153 except: 154 raise "Error while reading process_params table from file: ", coire_file 155 156 # and retrieve the time window from this file 157 found_flag = False 158 for tab in process_params: 159 if tab.param=='--injection-window': 160 found_flag = True 161 self.injection_window = float(tab.value)/1000.0 162 if not found_flag: 163 sys.stderr.write("WARNING: No entry '--injection-window' found in file %s"\ 164 "Value used is %.1f ms. If incorrect, please change file at %s\n" %\ 165 (coire_file, 1000.0*self.injection_window, __file__)) 166 if self.verbose: 167 print "Injection-window set to %.0f ms" % (1000*self.injection_window) 168 169 # read the veto files 170 self.readVetoFiles() 171 172 if self.verbose: 173 print "Initializing the Followup class done..."
174 175 176 # -----------------------------------------------------
177 - def get_injection_id(self, filename=None, url=None, cache_entry=None):
178 """ 179 Extracting the injection ID from the filename, using 180 the mechanism as used in lalapps_path2cache. You can specify the filename 181 itself, the url or the cache entry. You must not specify more than one input! 182 The injection-ID is calculated in the following way (exttrig only): 183 184 The code expects the INSPIRAL and THINCA files in the following scheme (example): 185 PREFIX-TAGPART_injections32_77-GPS-DURATION.xml 186 The number of the injection run is extracted (INJRUN) as well as the following 187 number (INJNUMBER). The injection ID is then calculated as: 188 INJID = 100000*INJRUN + INJNUMBER 189 so for this example the injectionj ID is 3200077. 190 191 @param file: filename from which the injection ID is extracted 192 @param url: url from which the injection ID is extracted 193 @param cache_entry: cache entry from which the injection ID is extracted 194 """ 195 196 # Check that only one input is given 197 if cache_entry: 198 if filename and url: 199 raise "Error in function get_injection_id: Only one input should be "\ 200 "specified. Now 'cache_entry' and another variable is specified. Check the code." 201 202 if cache_entry is None: 203 if filename and url: 204 raise "Error in function get_injection_id: Only one input should be "\ 205 "specified. Now 'filename' and 'url' is specified. Check the code." 206 207 if filename: 208 path, filename = os.path.split(filename.strip()) 209 url = "file://localhost%s" % os.path.abspath(os.path.join(path, filename)) 210 211 try: 212 cache_entry = lal.CacheEntry.from_T050017(url) 213 except ValueError, e: 214 raise "Error while extracting injection ID from file ", filename 215 216 # split the expression into different sub-pieces 217 pieces = cache_entry.description.split('_') 218 if self.exttrig: 219 220 # its easy for the exttrig case 221 injection_id = pieces[-2]+'_'+pieces[-1] 222 else: 223 224 # but need to check for the appearance of the CAT suffix else 225 index = 0 226 for ind, piece in enumerate(pieces): 227 if 'CAT' in piece: 228 index = ind 229 injection_id = pieces[index-1] 230 231 return injection_id
232 233 # -----------------------------------------------------
234 - def readVetoFiles( self ):
235 """ 236 Reads the veto segments given by veto-files (if any) 237 """ 238 self.vetodict = dict() 239 240 # loop over the IFO names 241 for ifoName in self.colors: 242 243 self.vetodict[ifoName]=None 244 # create the attribute name and check if it exists 245 attributeName = 'followup_vetofile_'+ifoName.lower() 246 if hasattr( self.opts, attributeName): 247 248 # get the filename 249 filename = getattr( self.opts, attributeName ) 250 251 if filename: 252 #read the segment lists 253 self.vetodict[ifoName] = segmentsUtils.fromsegwizard(open(filename))
254 255 # -----------------------------------------------------
256 - def reset( self ):
257 """ 258 Resets the counting number for the time-series plots generated. 259 """ 260 self.number=0
261 262 # -----------------------------------------------------
263 - def setTag( self, tag ):
264 """ 265 Just sets the tag, called from plotinspmissed (why needed?) 266 @param tag: tag for this followup 267 """ 268 269 self.tag = tag
270 271 # -----------------------------------------------------
272 - def print_inj( self, inj, injID ):
273 """ 274 Print some useful informations to the screen. 275 @param inj: the current missed injection 276 @param injID: the injection ID (used for exttrig only) 277 """ 278 279 if self.exttrig: 280 print "\nInjection details for injection %d with injID %s: " %\ 281 (self.number, injID) 282 else: 283 print "\nInjection details for injection %d:" % (self.number) 284 285 print "m1: %.2f m2:%.2f | end_time: %d.%d | "\ 286 "distance: %.2f eff_dist_h: %.2f eff_dist_l: %.2f" % \ 287 ( inj.mass1, inj.mass2, inj.geocent_end_time, inj.geocent_end_time_ns,\ 288 inj.distance, inj.eff_dist_h, inj.eff_dist_l )
289 290 # ----------------------------------------------------
291 - def savePlot( self, stage ):
292 """ 293 Saves the plots and store them in a seperate fnameList. 294 @param stage: the stage this plot belongs to (e.g. INSPIRAL, THINCA,...) 295 """ 296 fname = 'Images/'+self.opts.prefix + "_"+self.tag+"_map-"+\ 297 stage+"-"+str(self.number) +self.opts.suffix+'.png' 298 fname_thumb = InspiralUtils.\ 299 savefig_pylal( filename = self.output_path+fname,\ 300 doThumb = True, 301 dpi_thumb = self.opts.figure_resolution) 302 303 self.fnameList.append( fname ) 304 return fname
305 306 307 # -----------------------------------------------------
308 - def findInjection( self, missedInj ):
309 """ 310 Find the injection-ID corresponding to this particular 311 missed injection. 312 @param missedInj: the missed injection 313 """ 314 315 # injID: the ID (a number) 316 # groupInj: a list of SimInspiral tables 317 for injID, groupInj in self.injections.iteritems(): 318 for inj in groupInj: 319 if missedInj.geocent_end_time==inj.geocent_end_time and \ 320 missedInj.geocent_end_time_ns==inj.geocent_end_time_ns: 321 return injID 322 323 self.print_inj( missedInj, None) 324 raise "No injection ID found for the above particular missed Injection " 325 326 return None
327 328 # -----------------------------------------------------
329 - def getTimeTrigger( self, trig ):
330 """ 331 This is a helper function to return a GPS time as one float number 332 @param trig: a sngl_inspiral table entry 333 """ 334 335 return float(trig.end_time) + float(trig.end_time_ns) * 1.0e-9
336 337 338 # -----------------------------------------------------
339 - def getTimeSim(self, sim, ifo=None ):
340 """ 341 This is a helper function to return a GPS time as one float number 342 for a certain IFO. If no IFO is specified the injected geocentric 343 time is returned. 344 @param sim: a sim_inspiral table entry 345 @param ifo: the IFO for which we want the sim time 346 """ 347 348 time=0 349 nano=0 350 351 if not ifo: 352 time = sim.geocent_end_time 353 nano = sim.geocent_end_time_ns 354 if ifo: 355 time = sim.get(ifo[0].lower()+'_end_time' ) 356 nano = sim.get(ifo[0].lower()+'_end_time_ns' ) 357 358 return float(time) + float(nano) * 1.0e-9
359 360 # -----------------------------------------------------
361 - def highlightVeto( self, timeInjection, segLarge, ifoName, ylims ):
362 """ 363 Finds the intersection of the drawn triggers with the veto segments 364 for this IFO 365 """ 366 367 if not self.vetodict[ifoName]: 368 return 369 370 # find intersecting veto segments 371 segVeto = self.vetodict[ifoName] & segments.segmentlist( [segLarge] ) 372 373 # draw all intersecting segments 374 for seg1 in segVeto: 375 376 # after a tim-shift 377 seg=seg1.shift( -timeInjection) 378 vetox = [seg[0], seg[1], seg[1], seg[0], seg[0]] 379 vetoy = [ylims[0], ylims[0], ylims[1], ylims[1], ylims[0]] 380 fill ( vetox, vetoy, 'y', alpha=0.2)
381 382 383 # -----------------------------------------------------
384 - def isThereVeto( self, timeTrigger, ifoName ):
385 """ 386 This function checks if at the time 'timeTrigger' the 387 IFO 'ifoName' is vetoed. 388 @param timeTrigger: The time to be investigated 389 @param ifoName: The name of the IFO to be investigated 390 """ 391 392 if self.vetodict[ifoName] is None: 393 flag = False 394 else: 395 flag = timeTrigger in self.vetodict[ifoName] 396 397 return flag
398 399 # -----------------------------------------------------
400 - def getExpectedSNR(self, triggerFiles, inj, number ):
401 """ 402 Investigate template bank and returns exepcted horizon distance 403 @param triggerFiles: List of files containing the inspiral triggers 404 @param inj: The current injection 405 @param ifo: The current IFO 406 @param number: The consecutive number for this inspiral followup 407 """ 408 409 # read the inspiral file(s) 410 if self.verbose: print "Processing TMPLTBANK triggers from files ",\ 411 triggerFiles 412 inspiralSumm, massInspiralSumm = InspiralUtils.\ 413 readHorizonDistanceFromSummValueTable(\ 414 triggerFiles, self.verbose) 415 416 # get different informations on the missed injection 417 injMass = [inj.mass1, inj.mass2] 418 timeInjection = self.getTimeSim( inj ) 419 totalMass = inj.mass1 + inj.mass2 420 eta = inj.mass1 * inj.mass2 / totalMass / totalMass 421 422 # Given the masses (and therefore eta), the expected horizon distance(SNR) has to be scaled by this factor 423 factor = sqrt(4 * eta) 424 425 output = {} 426 # for each ifo that has been found 427 for ifo in massInspiralSumm.keys(): 428 # loop over the summ value table as long as horizon is not set 429 output[ifo] = [] 430 horizon = 0 431 for massNum in range(len(massInspiralSumm[ifo])): 432 # looking at the masses and horizon distance (and threshold) 433 if horizon > 0: 434 break 435 for this, horizon in zip(massInspiralSumm[ifo][massNum].\ 436 getColumnByName('comment'),\ 437 massInspiralSumm[ifo][massNum].getColumnByName('value').asarray()): 438 masses = this.split('_') 439 threshold = float(masses[2]) 440 readTotalMass = float(masses[0])+float(masses[1]) 441 # check if the current total Mass is greater tan the requested total Mass 442 # if so, then keep this horizon and break 443 if readTotalMass > totalMass: 444 startTimeSec = \ 445 massInspiralSumm[ifo][massNum].getColumnByName('start_time').asarray() 446 #sanity check 447 if len(startTimeSec)>1: 448 print >> sys.stderr, 'Warning in fct. expectedHorizonDistance: '\ 449 'More than 1 file found at particular GPS time. Using the first file.' 450 if startTimeSec[0] > inj.geocent_end_time: 451 text= """the start time of the template bank must be less than 452 the end_time of the injection. We found startTime of the bank to be %s and 453 the geocent end time of the injection to be %s""",startTimeSec[0], inj.geocent_end_time 454 raise ValueError, text 455 output[ifo] = horizon * factor * threshold / float(getattr(inj, 'eff_dist_'+ifo[0].lower() )) 456 break 457 #'otherwise, reset horizon distance 458 else: horizon = 0 459 460 return output
461 462 # -----------------------------------------------------
463 - def putText( self, text ):
464 """ 465 Puts some text into an otherwise empty plot. 466 @param text: text to put in the empty plot 467 """ 468 newText = '' 469 for i in range( int(len(text)/60.0)+1): 470 newText+=text[60*i:60*i+60]+'\n' 471 figtext(0.15,0.15, newText)
472 473 # -----------------------------------------------------
474 - def investigateTimeseries(self, triggerFiles, inj, ifoName, stage, number ):
475 """ 476 Investigate inspiral triggers and create a time-series 477 of the SNRs around the injected time 478 @param triggerFiles: List of files containing the inspiral triggers 479 @param inj: the current missed injection 480 @param ifoName: the IFO for which the plot is made 481 @param stage: the name of the stage (FIRST, SECOND) 482 @param number: the consecutive number for this inspiral followup 483 """ 484 485 # read the inspiral file(s) 486 if self.verbose: 487 print "Processing INSPIRAL triggers from files ", triggerFiles 488 489 snglTriggers = SnglInspiralUtils.ReadSnglInspiralFromFiles( \ 490 triggerFiles , verbose=False) 491 492 # create a figure and initialize some lists 493 fig=figure() 494 foundSet = set() 495 loudest_details = {} 496 noTriggersFound = True 497 498 if snglTriggers is None: 499 # put message on the plot instead 500 self.putText( 'No sngl_inspiral triggers in %s' % str(triggerFiles)) 501 502 else: 503 # selection segment 504 timeInjection = self.getTimeSim( inj ) 505 segSmall = segments.segment( timeInjection-self.injection_window, \ 506 timeInjection+self.injection_window ) 507 segLarge = segments.segment( timeInjection-self.time_window, \ 508 timeInjection+self.time_window ) 509 510 # create coincidences for THINCA stage 511 coincTriggers = None 512 if 'THINCA' in stage: 513 coincTriggers = CoincInspiralUtils.coincInspiralTable( snglTriggers, \ 514 CoincInspiralUtils.coincStatistic("snr") ) 515 selectedCoincs = coincTriggers.vetoed( segSmall ) 516 517 # loop over the IFOs (although this is a plot for IFO 'ifoName') 518 for ifo in self.colors.keys(): 519 520 # get the singles for this ifo 521 snglInspiral = snglTriggers.ifocut(ifo) 522 523 # select a range of triggers 524 selectedLarge = snglInspiral.vetoed( segLarge ) 525 timeLarge = [ self.getTimeTrigger( sel )-timeInjection \ 526 for sel in selectedLarge ] 527 528 selectedSmall = snglInspiral.vetoed( segSmall ) 529 timeSmall = [ self.getTimeTrigger( sel )-timeInjection \ 530 for sel in selectedSmall ] 531 532 # use the set of selected coincident triggers in the THINCA stages 533 if coincTriggers: 534 selectedSmall = selectedCoincs.cluster(2* self.injection_window).getsngls(ifo) 535 timeSmall = [ self.getTimeTrigger( sel )-timeInjection \ 536 for sel in selectedSmall ] 537 538 # skip if no triggers in the large time window 539 if len(timeLarge)==0: 540 continue 541 noTriggersFound = False 542 543 # add IFO to this set; the injection is found for this IFO and stage 544 if len(timeSmall)>0: 545 foundSet.add(ifo) 546 547 # record details of the loudest trigger 548 loudest_details[ifo] = {} 549 loudest = selectedSmall[selectedSmall.get_column('snr').argmax()] 550 loudest_details[ifo]["snr"] = loudest.snr 551 loudest_details[ifo]["mchirp"] = loudest.mchirp 552 loudest_details[ifo]["eff_dist"] = loudest.eff_distance 553 loudest_details[ifo]["chisq"] = loudest.chisq 554 loudest_details[ifo]["timeTrigger"] = self.getTimeTrigger( loudest ) 555 556 timeTrigger = self.getTimeTrigger( loudest ) 557 vetoSegs = self.vetodict[ifoName] 558 559 # plot the triggers 560 plot( timeLarge, selectedLarge.get_column('snr'),\ 561 self.colors[ifo]+'o', label="_nolegend_") 562 plot( timeSmall, selectedSmall.get_column('snr'), \ 563 self.colors[ifo]+'s', label=ifo) 564 565 # draw the injection times and other stuff 566 if noTriggersFound: 567 self.putText( 'No triggers/coincidences found within time window') 568 569 ylims=axes().get_ylim() 570 plot( [0,0], ylims, 'g--', label="_nolegend_") 571 plot( [-self.injection_window, -self.injection_window], ylims, 'c:',\ 572 label="_nolegend_") 573 plot( [+self.injection_window, +self.injection_window], ylims, 'c:',\ 574 label="_nolegend_") 575 576 self.highlightVeto( timeInjection, segLarge, ifoName, ylims ) 577 578 # save the plot 579 grid(True) 580 legend() 581 582 ylims=axes().get_ylim() 583 axis([-self.time_window, +self.time_window, ylims[0], ylims[1]]) 584 xlabel('time [s]') 585 ylabel('SNR') 586 title(stage+'_'+str(self.number)) 587 fname = self.savePlot( stage ) 588 close(fig) 589 590 result = {'filename':fname, 'foundset':foundSet, 'loudest_details':loudest_details} 591 return result
592 593 594 # -----------------------------------------------------
595 - def select_category(self, trigger_files, category):
596 """ 597 Return a trigger list that contains only files for the choosen category. 598 @param triggerList : a list of file names 599 @param category: a category number 600 @return: a sub list of filename corresponding to the category requested 601 """ 602 603 # there are two different labels used to denote 604 # the categories. THIS NEEDS TO BE UNIFIED 605 cat1 = 'CAT_'+str(category) 606 cat2 = 'CATEGORY_'+str(category) 607 608 if category==1: 609 # Category 1 files might not be labelled with a 'CAT_1' suffix. 610 # So, for now, all files NOT containing the expression 611 # 'CAT' in the filename are supposed to be CAT_1 files 612 new_list = [file for file in trigger_files \ 613 if 'CAT' not in file or cat1 in file or cat2 in file] 614 615 else: 616 cat = 'CAT_'+str(category) 617 new_list = [file for file in trigger_files if cat1 in file\ 618 or cat2 in file] 619 620 return new_list
621 622 623 # -----------------------------------------------------
624 - def followup(self, inj, selectIFO, description = None):
625 """ 626 Do the followup procedure for the missed injection 'inj' 627 and create the several time-series for INSPIRAL and THINCA. 628 The return value is the name of the created html file. 629 @param inj: sim_inspiral table of the injection that needs to be 630 followed up 631 @param selectIFO: the IFO that is investigated 632 @param description: Can be used to sieve further this pattern 633 from the description field. 634 """ 635 636 def fill_table(page, contents ): 637 """ 638 Making life easier... 639 """ 640 page.add('<tr>') 641 for content in contents: 642 page.add('<td>') 643 page.add( str(content) ) 644 page.add('</td>') 645 page.add('</tr>')
646 647 648 # get the ID corresponding to this injection 649 injection_id = self.findInjection( inj ) 650 651 # increase internal number: 652 self.number+=1 653 654 ## create the web-page and add a table 655 page = markup.page() 656 page.h1("Followup missed injection #"+str(self.number)+" in "+selectIFO ) 657 page.hr() 658 page.add('<table border="3" ><tr><td>') 659 page.add('<table border="2" >') 660 fill_table( page, ['<b>parameter','<b>value'] ) 661 fill_table( page, ['Number', self.number] ) 662 fill_table( page, ['inj ID', injection_id] ) 663 fill_table( page, ['mass1', '%.2f'% inj.mass1] ) 664 fill_table( page, ['mass2', '%.2f'%inj.mass2] ) 665 fill_table( page, ['mtotal', '%.2f' % (inj.mass1+inj.mass2)] ) 666 fill_table( page, ['mchirp', '%.2f' % (inj.mchirp)] ) 667 fill_table( page, ['end_time', inj.geocent_end_time] ) 668 fill_table( page, ['end_time_ns', inj.geocent_end_time_ns] ) 669 fill_table( page, ['distance', '%.1f' % inj.distance] ) 670 fill_table( page, ['eff_dist_h','%.1f' % inj.eff_dist_h] ) 671 fill_table( page, ['eff_dist_l','%.1f' % inj.eff_dist_l] ) 672 fill_table( page, ['eff_dist_v','%.1f' % inj.eff_dist_v] ) 673 fill_table( page, ['eff_dist_g','%.1f' % inj.eff_dist_g] ) 674 fill_table( page, ['playground','%s' % pipeline.s2play(inj.geocent_end_time)] ) 675 page.add('</table></td>') 676 677 # print infos to screen if required 678 if self.opts.verbose: 679 self.print_inj( inj, injection_id) 680 681 # sieve the cache for the required INSPIRAL and THINCA files 682 invest_dict = {} 683 for stage, cache in self.triggerCache.iteritems(): 684 685 trig_cache = lal.Cache() 686 for c in cache: 687 688 # check the time and the injection ID 689 if inj.geocent_end_time in c.segment: 690 if self.get_injection_id(url = c.url) == injection_id: 691 trig_cache.append( c ) 692 693 # create a filelist 694 file_list = trig_cache.sieve(description = description).pfnlist() 695 696 # check if the pfnlist is empty. ` 697 if len(file_list)==0: 698 print >>sys.stderr, "Error: No files found for stage %s in the "\ 699 "cache for ID %s and time %d; probably mismatch of a "\ 700 "pattern in the options. " % \ 701 ( stage, injection_id, inj.geocent_end_time) 702 continue 703 704 # if the stage if THINCA_SECOND... 705 if 'THINCA_SECOND' in stage: 706 707 # ... need to loop over the four categories 708 for cat in [1,2,3,4]: 709 710 select_list=self.select_category( file_list, cat) 711 if len(select_list)==0: 712 print "WARNING: No THINCA_SECOND files found for category ", cat 713 continue 714 715 modstage = stage+'_CAT_' + str(cat) 716 invest_dict[modstage] = self.investigateTimeseries( select_list, inj, selectIFO, modstage, self.number ) 717 718 #sys.exit(0) 719 else: 720 invest_dict[stage]=self.investigateTimeseries( file_list, inj, selectIFO, stage, self.number) 721 722 723 724 ## print out the result for this particular injection 725 page.add('<td><table border="2" >') 726 fill_table( page, ['<b>step','<b>F/M', '<b>Rec. SNR', '<b>Rec. mchirp', \ 727 '<b>Rec. eff_dist', '<b>Rec. chisq', '<b>Veto ON/OFF'] ) 728 729 # loop over the stages and create the table with 730 # the various data in it (when available) 731 for stage in self.orderLabels: 732 if stage in invest_dict: 733 result = invest_dict[stage] 734 735 # Fill in the details of the loudest found coinc. 736 #found_ifo='' 737 #if "INSPIRAL" in stage or "THINCA" in stage: 738 found_ifo='' 739 loudest_snr='' 740 loudest_mchirp='' 741 loudest_eff_dist='' 742 loudest_chisq='' 743 veto_onoff='' 744 745 # add all the IFO's for this coincident 746 for ifo in result['foundset']: 747 found_ifo += ifo+' ' 748 749 # Parameters of the loudest trigger, taken from the 750 # 'loudest-details' dictionary, created in 'investigateTimeseries' 751 loudest_snr += ifo + ': ' + str(result['loudest_details'][ifo]['snr'])+'<br>' 752 loudest_mchirp += ifo + ': ' + str(result['loudest_details'][ifo]['mchirp'])+'<br>' 753 loudest_eff_dist += ifo + ': ' + str(result['loudest_details'][ifo]['eff_dist'])+'<br>' 754 loudest_chisq += ifo + ': ' + str(result['loudest_details'][ifo]['chisq'])+'<br>' 755 756 # Check whether some of the ifo times is vetoed 757 timeTrigger = float(result['loudest_details'][ifo]['timeTrigger']) 758 if (self.vetodict[ifo]): 759 veto = self.isThereVeto (timeTrigger, ifo) 760 veto_txt = 'OFF' 761 if veto: 762 veto_txt = 'ON' 763 veto_onoff+=ifo+': '+veto_txt+'<br>' 764 else: 765 veto_onoff+=ifo+': No info<br>' 766 767 # Fill the table whether something is found or not 768 if len(result['foundset'])>0: 769 fill_table( page, [ stage, 'FOUND in '+found_ifo, 'loudest<br>'+loudest_snr, \ 770 'loudest<br>'+loudest_mchirp, 'loudest<br>'+loudest_eff_dist,\ 771 'loudest<br>'+loudest_chisq, veto_onoff]) 772 else: 773 fill_table( page, [ stage, '<font color="red">MISSED']) 774 775 page.add('</table>') 776 page.add('</td></tr></table><br><br>') 777 778 779 ## add the pictures to the webpage 780 for stage in self.orderLabels: 781 if stage in invest_dict: 782 result = invest_dict[stage] 783 784 ##if stage!="TMPLTBANK": 785 if True: 786 fname = result['filename'] 787 page.a(extra.img(src=[fname], width=400, \ 788 alt=fname, border="2"), title=fname, href=[ fname ]) 789 790 # add version information 791 page.add('<hr>Page created with %s Version %s' % \ 792 (__prog__, git_version.verbose_msg)) 793 794 # and write the html file 795 htmlfilename = self.opts.prefix + "_"+selectIFO+"_followup_"+str(self.number) +\ 796 self.opts.suffix+'.html' 797 file = open(self.opts.output_path+htmlfilename,'w') 798 file.write(page(False)) 799 file.close() 800 801 # store html file in fnameList 802 self.fnameList.append(htmlfilename) 803 804 # supply the output 805 return htmlfilename
806 807 # -----------------------------------------------------
808 - def estimatedDistance( self, mass1, mass2, distTarget):
809 """ 810 Calculates the approx. effective distance (in Mpc) for a 811 binary of two objects with masses 'mass1' and 'mass2', 812 when the efective distance for a 10/10 binary is 'distTarget'. 813 This distance is rescaled given the two masses and the spectrum. 814 Example 815 816 estimatedDistance(10.0, 10.0, 100.0) 817 should return exactly 100.0 again. 818 819 @param mass1: The mass of the first component (Solar masses) 820 @param mass2: The mass of the second component (Solar masses) 821 @param distTarget: Effective distance for a 10/10 binary 822 """ 823 824 snrActual = 8.0 825 distActual = computeCandleDistance( 10.0, 10.0, \ 826 self.flow, self.spectrum, \ 827 snrActual) 828 829 # rescale the SNR threshold 830 snrTarget = distActual / distTarget*snrActual 831 distance = computeCandleDistance( mass1, mass2, \ 832 self.flow, self.spectrum,\ 833 snrTarget) 834 return distance
835 836 837 # -----------------------------------------------------
838 -def createSpectrum( flow, sampleRate = 4096, nPoints = 1048576):
839 """ 840 Computes the approximate spectrum of LIGO-I. 841 @param flow: The lower cutoff frequency 842 @param sampleRate: The sample rate used (default: 4096) 843 @param nPoints The number of points (default: 1048576) 844 """ 845 846 def calcPSD( f ): 847 f0=150.0 848 a = pow(4.49*f/f0, -56.0) 849 b = 0.16*pow(f/f0, -4.52) 850 c = 0.52 851 d = 0.32*pow(f/f0,2) 852 return (a+b+c+d)*9.0e-46
853 854 chanDeltaT = 1.0/float(sampleRate) 855 deltaF = 1.0/( float(nPoints)*chanDeltaT ) 856 857 spectrum = [] 858 for k in range( nPoints/2+1): 859 f = k*deltaF 860 if f<flow: 861 spectrum.append(0) 862 else: 863 spectrum.append( calcPSD(f) ) 864 865 return spectrum 866 867 # -----------------------------------------------------
868 -def computeCandleDistance( candleM1, candleM2, fLow, spectrum = None, \ 869 snr=8.0, sampleRate = 4096, nPoints = 1048576):
870 """ 871 Computes the candle distance as computed in inspiralutils.c. 872 This code has been tested on 21 Feb 2008 and the derivation between 873 these calculated values with the LAL code is less than 2 Percent. 874 @param candleM1: The mass of the first component (Solar masses) 875 @param candleM2: The mass of the second component (Solar masses) 876 @param flow: The lower cutoff frequency 877 @param spectrum: The PSD that has been calculated earlier 878 @param snr: The SNR at which the distance is calculated (default: 8) 879 @param sampleRate: The sample rate used (default: 4096) 880 @param nPoints The number of points (default: 1048576) 881 """ 882 # TestCode: testEstimatedDistance.png 883 884 885 LAL_MTSUN_SI = 4.92549095e-6 886 chanDeltaT = 1.0/float(sampleRate) 887 negativeSevenOverThree = -7.0/3.0 888 totalMass = candleM1 + candleM2 889 mu = candleM1 * candleM2 / totalMass 890 distNorm = 9.5708317e-20 # = 2.0 * LAL_MRSUN_SI / (1.0e6 * LAL_PC_SI ) 891 a = sqrt( (5.0 * mu) / 96.0 ) * \ 892 pow( totalMass / ( pi*pi ), 1.0/3.0 ) *\ 893 pow( LAL_MTSUN_SI / chanDeltaT, -1.0/6.0 ) 894 sigmaSq = 4.0 * ( chanDeltaT / float(nPoints) ) * \ 895 distNorm * distNorm * a * a 896 897 fmax = 1.0 / (6.0 * sqrt(6.0) * pi * totalMass * LAL_MTSUN_SI) 898 deltaF = 1.0/( float(nPoints)*chanDeltaT ) 899 900 cut = int( fLow/deltaF ) 901 kmax = min( int( fmax/deltaF ), len(spectrum) ) 902 903 sigmaSqSum = 0 904 for k in range( cut, kmax): 905 sigmaSqSum += pow( float(k)/float(nPoints), negativeSevenOverThree ) /\ 906 spectrum[k] 907 908 sigmaSq *= sigmaSqSum 909 distance = sqrt( sigmaSq ) / snr 910 911 return distance
912