1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
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
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
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
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
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
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
97
98 self.injection_window = 0.050
99
100
101 self.fnameList = []
102
103
104 self.number = 0
105
106
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
116 if opts.followup_tag:
117 self.cache = cache.sieve(description = opts.followup_tag)
118 else:
119 self.cache = cache
120
121
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
130 self.injectionCache = self.cache.sieve(description = "INJECTION").\
131 sieve(ifos='HL')
132
133
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
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
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
170 self.readVetoFiles()
171
172 if self.verbose:
173 print "Initializing the Followup class done..."
174
175
176
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
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
217 pieces = cache_entry.description.split('_')
218 if self.exttrig:
219
220
221 injection_id = pieces[-2]+'_'+pieces[-1]
222 else:
223
224
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
235 """
236 Reads the veto segments given by veto-files (if any)
237 """
238 self.vetodict = dict()
239
240
241 for ifoName in self.colors:
242
243 self.vetodict[ifoName]=None
244
245 attributeName = 'followup_vetofile_'+ifoName.lower()
246 if hasattr( self.opts, attributeName):
247
248
249 filename = getattr( self.opts, attributeName )
250
251 if filename:
252
253 self.vetodict[ifoName] = segmentsUtils.fromsegwizard(open(filename))
254
255
257 """
258 Resets the counting number for the time-series plots generated.
259 """
260 self.number=0
261
262
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
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
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
309 """
310 Find the injection-ID corresponding to this particular
311 missed injection.
312 @param missedInj: the missed injection
313 """
314
315
316
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
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
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
371 segVeto = self.vetodict[ifoName] & segments.segmentlist( [segLarge] )
372
373
374 for seg1 in segVeto:
375
376
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
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
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
410 if self.verbose: print "Processing TMPLTBANK triggers from files ",\
411 triggerFiles
412 inspiralSumm, massInspiralSumm = InspiralUtils.\
413 readHorizonDistanceFromSummValueTable(\
414 triggerFiles, self.verbose)
415
416
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
423 factor = sqrt(4 * eta)
424
425 output = {}
426
427 for ifo in massInspiralSumm.keys():
428
429 output[ifo] = []
430 horizon = 0
431 for massNum in range(len(massInspiralSumm[ifo])):
432
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
442
443 if readTotalMass > totalMass:
444 startTimeSec = \
445 massInspiralSumm[ifo][massNum].getColumnByName('start_time').asarray()
446
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
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
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
486 if self.verbose:
487 print "Processing INSPIRAL triggers from files ", triggerFiles
488
489 snglTriggers = SnglInspiralUtils.ReadSnglInspiralFromFiles( \
490 triggerFiles , verbose=False)
491
492
493 fig=figure()
494 foundSet = set()
495 loudest_details = {}
496 noTriggersFound = True
497
498 if snglTriggers is None:
499
500 self.putText( 'No sngl_inspiral triggers in %s' % str(triggerFiles))
501
502 else:
503
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
511 coincTriggers = None
512 if 'THINCA' in stage:
513 coincTriggers = CoincInspiralUtils.coincInspiralTable( snglTriggers, \
514 CoincInspiralUtils.coincStatistic("snr") )
515 selectedCoincs = coincTriggers.vetoed( segSmall )
516
517
518 for ifo in self.colors.keys():
519
520
521 snglInspiral = snglTriggers.ifocut(ifo)
522
523
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
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
539 if len(timeLarge)==0:
540 continue
541 noTriggersFound = False
542
543
544 if len(timeSmall)>0:
545 foundSet.add(ifo)
546
547
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
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
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
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
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
604
605 cat1 = 'CAT_'+str(category)
606 cat2 = 'CATEGORY_'+str(category)
607
608 if category==1:
609
610
611
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
649 injection_id = self.findInjection( inj )
650
651
652 self.number+=1
653
654
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
678 if self.opts.verbose:
679 self.print_inj( inj, injection_id)
680
681
682 invest_dict = {}
683 for stage, cache in self.triggerCache.iteritems():
684
685 trig_cache = lal.Cache()
686 for c in cache:
687
688
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
694 file_list = trig_cache.sieve(description = description).pfnlist()
695
696
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
705 if 'THINCA_SECOND' in stage:
706
707
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
719 else:
720 invest_dict[stage]=self.investigateTimeseries( file_list, inj, selectIFO, stage, self.number)
721
722
723
724
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
730
731 for stage in self.orderLabels:
732 if stage in invest_dict:
733 result = invest_dict[stage]
734
735
736
737
738 found_ifo=''
739 loudest_snr=''
740 loudest_mchirp=''
741 loudest_eff_dist=''
742 loudest_chisq=''
743 veto_onoff=''
744
745
746 for ifo in result['foundset']:
747 found_ifo += ifo+' '
748
749
750
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
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
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
780 for stage in self.orderLabels:
781 if stage in invest_dict:
782 result = invest_dict[stage]
783
784
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
791 page.add('<hr>Page created with %s Version %s' % \
792 (__prog__, git_version.verbose_msg))
793
794
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
802 self.fnameList.append(htmlfilename)
803
804
805 return htmlfilename
806
807
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
830 snrTarget = distActual / distTarget*snrActual
831 distance = computeCandleDistance( mass1, mass2, \
832 self.flow, self.spectrum,\
833 snrTarget)
834 return distance
835
836
837
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
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
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