Hardware injection waveform generation 

Introduction 

This page describes how to generate waveforms and save them as single-column ASCII waveform files that can be used by awgstream to inject into the detector.

There are two executables that can be used to generate single-column ASCII files; they are pycbc_generate_hwinj and pycbc_generate_hwinj_from_xml. Both executables use the PyCBC injection module (pycbc.inject) to inject the coherent waveform into a time series of zeroes.

The executable pycbc_generate_hwinj generates a waveform using parameters from the command line. The user inputs parameters such as --mass1, --mass2, etc. on the command line. This executable is useful for generating a specific coherent waveform for hardware injections.

The executable pycbc_generate_hwinj_from_xml generates all the waveforms in a LIGOLW sim_inspiral table. The output of lalapps_inspinj (an executable for generating a population of injections) is a LIGOLW sim_inspiral table. This executable is useful if you want to generate a population of coherent waveforms for hardware injections.

Generate waveform from command line (`pycbc_generate_hwinj`)

Here is a usage example for generating a CBC waveform using detector data with pycbc_generate_hwinj.

Select a time for the injection 

First on the command line set a variable for the GPS geocentric end time of the coherent injection

GEOCENT_END_TIME=1124381661

Select data for PSD estimation 

We need to set what data we will use to estimate the PSD. This should be a 2048 second interval. So in this example we will set

GPS_START_TIME=1124380361
GPS_END_TIME=1124382409

We will want to use data from when the detector was in science mode time to estimate the PSD. To check if the detector was in science mode you can query the segment database, check the section How to query the segment database for an example command.

Since we are using detector data we will need to read the data from frame files. We will need to specify the frame type and channel name. The frame type is a way to identify what list of channels is in the frame file. So on the command line we set

FRAME_TYPE=H1_HOFT_C00
CHANNEL_NAME=H1:GDS-CALIB_STRAIN

We can check that frames exist for this time by querying the LDR server, check the section How to query the LDR server for an example command.

Do not assume that the frame type and channel name are the same as in this section. These values are correct for this example.

Run `pycbc_generate_hwinj`

The --instruments option specifies the IFOs to include in the injection. The executable pycbc_generate_hwinj will write the waveform for all detectors listed with the --instruments option. There are a number of options, including --instruments that take a space-seperated list, example usage of this option is

--instruments H1 L1

Now let’s specify the parameters of the waveform. A full list of the command line options is available with pycbc_generate_hwinj --help. Here we provide an example for generating a coherent 1.4-1.4 component mass binary using the EOBNRv2 approximant. Here is an example of this command line case

--approximant SEOBNRv2 --order pseudoFourPN --mass1 25.0 --mass2 25.0 --inclination 0.0 --polarization 0.0 --ra 0.0 --dec 0.0

The sample rate of the output ASCII file is given by the --sample-rate option. The sample rate for each IFO must match, example usage of this option is

--sample-rate H1:16384 L1:16384

We have the option to query the LDR server directly with pycbc_generate_hwinj. So we only need to specify a --frame-type option to get the data. Alternatively you an pass a space-seperated list with --frame-files or a LAL frame cache with --frame-cache. We also need to say what channel to use for the PSD estimation with the --channel-name option. Example usage of --frame-type would be

--frame-type H1:${FRAME_TYPE} L1:${FRAME_TYPE} --channel-name H1:${CHANNEL_NAME} L1:${CHANNEL_NAME}

We do not want to inject a step-like response into the detector, therefore we taper the waveform at the beginning. The EOBNRv2 has a ringdown at the end so we do not need to taper the end. Example usage of this option is

--taper TAPER_START

We specify the network SNR we want the coherent injection to have on the command line. The network SNR calculation includes all IFOs specified in the --insturments option, example usage of this option is

--network-snr 28

In calculating the network SNR the executable pycbc_generate_hwinj will generate a PSD and calculate an SNR for the waveform. The options --low-frequency-cutoff and --high-frequency-cutoff set the min and max frequency for the SNR calculation. The waveform used in the SNR calculation is also generated at this low-frequency cutoff, note the waveform is not written to disk with this low-frequency cutoff. Example usage of the PSD options is

--low-frequency-cutoff 40.0 --high-frequency-cutoff 1000.0 --psd-estimation median --psd-segment-length 16 --psd-segment-stride 8 --pad-data 8

The additional PSD options dictate how the PSD will be calculated, ie. how many seconds per FFT and how much overlap. The --pad-data option is how much data to disgard at the edges of our time series used in PSD estimation to avoid data corruption.

The --waveform-low-frequency-cutoff option is the frequency that pycbc_generate_hwinj will begin generating the waveform that is written to file.

Here is a full example command for generating an injection in only H1

pycbc_generate_hwinj --high-frequency-cutoff 1000.0 --geocentric-end-time ${GEOCENT_END_TIME} --gps-start-time H1:${GPS_START_TIME} --gps-end-time H1:${GPS_END_TIME} --frame-type H1:${FRAME_TYPE} --channel-name H1:${CHANNEL_NAME} --approximant SEOBNRv2 --order pseudoFourPN --mass1 25.0 --mass2 25.0 --inclination 0.0 --polarization 0.0 --ra 0.0 --dec 0.0 --taper TAPER_START --network-snr 28 --waveform-low-frequency-cutoff 10.0 --low-frequency-cutoff 40.0 --sample-rate 16384 --pad-data 8 --strain-high-pass 30.0 --psd-estimation median --psd-segment-length 16 --psd-segment-stride 8 --instruments H1

This will generate a single-column ASCII files that contains the h(t) time series for each detector and a LIGOLW XML file with the waveform parameters. The output filenames are not specified on the command line, they are determined internally by pycbc_generate_hwinj. In this example the ASCII file with the waveform will be named hwinjcbc_${START}_H1.txt where ${START} is the start time stamp of the time series. The LIGOLW XML file will be named hwinjcbc_${START}.xml.gz.

The LIGOLW XML file contains a process_params table that saves the command line that was used to generate the waveform for future reference. It also includes a sim_inspiral table and a sngl_inspiral table. The sim_inspiral entry allows us to use the parameters of the waveform as a software injection in the PyCBC matched filtering executable pycbc_inspiral. The sngl_inspiral entry allows us to use the parameters of the waveform as the filter in pycbc_inspiral.

The user should inspect the waveforms. For a waveform plotting executable see section Plot ASCII waveform files with pycbc_plot_hwinj.

Generate waveform from `lalapps_inspinj` output (`pycbc_generate_hwinj_from_xml`)

Here is a usage case for generating a population of waveforms with lalapps_inspinj. This example generates an injection every Tuesday for three months.

Run `lalapps_inspinj`

Here we show an example on how to use lalapps_inspinj to generate a population of injections.

In this example we will select distributions for time, distance, inclination, mass, spin, and sky location. Below is a explaination of the command line options in this example. A full list of command line options can be found with lalapps_inspinj --help.

Our time distribution will be a fixed time step to perform an injection every Tuesday. We will allow the injections to be anytime during the day (86400 seconds) and have a minimum one week (604800 seconds) between injections. The command line options will be

--time-interval 86400 --time-step 604800 --gps-start-time 1126368017 --gps-end-time 1130371217

Our distance distribution will be uniformly distributed in volume. We set the minimum and maximum chirp distance in units of kpc. The command line options will be

--d-distr volume --min-distance 10000 --max-distance 40000

Our mass distribution will be uniform in total mass. We can select the minimum and maximum component masses in units of solar masses. The command line options will be

--m-distr totalMass --min-mass1 1.0 --max-mass1 2.0 --min-mass2 1.0 --max-mass2 2.0

We can select the minimum and maximum component spins. The command line options will be

--enable-spin --min-spin1 0.0 --max-spin1 0.04 --min-spin2 0.0 --max-spin2 0.04

Our inclination distribution will be uniform. The command line option will be

--i-distr uniform

Our source distribution will be random. The command line option will be

--l-distr random

We select to use the SpinTaylorT4 approximant and begin the waveforms at 10.0Hz. Here we taper the injection at the start and end of the injection. The command line options will be

--waveform SpinTaylorT4threePointFivePN --f-lower 10 --taper-injection startend --band-pass-injection

Now we can combine all the options above and run lalapps_inspinj as

lalapps_inspinj --time-interval 86400 --time-step 604800 --gps-start-time 1126368017 --gps-end-time 1130371217 --d-distr volume --min-distance 10000 --max-distance 40000 --m-distr totalMass --min-mass1 1.0 --max-mass1 2.0 --min-mass2 1.0 --max-mass2 2.0 --enable-spin --min-spin1 0.0 --max-spin1 0.04 --min-spin2 0.0 --max-spin2 0.04 --i-distr uniform --l-distr random --waveform SpinTaylorT4threePointFivePN --f-lower 10 --taper-injection startend --band-pass-injection

In this example lalapps_inspinj will write a LIGOLW XML file called HL-INJECTIONS_1-1126368017-4003200.xml that has a sim_inspiral table with the population of injections.

Run `pycbc_generate_hwinj_from_xml`

Running lalapps_inspinj has written a LIGOLW XML file with a sim_inspiral table. Now we can run pycbc_generate_hwinj_from_xml to write single-column ASCII waveform files for the population of injections.

There are just two command line options --injection-file (path to the LIGOLW XML file that lalapps_inspinj had written) and --sample-rate (the sample rate of the ASCII waveform files).

In this example we set the sample rate to 16384Hz so on the command line do

pycbc_generate_hwinj_from_xml --injection-file HL-INJECTIONS_1-1126368017-4003200.xml --sample-rate 16384

As this command runs it will generate a H1 and L1 ASCII waveform file for each row in the sim_inspiral table.

The ASCII waveform files will be named ${IFO}-HWINJ_CBC_SIMULATION_ID_${SIMID}-${START}-${DURATION}.txt where where ${SIMID} is the simulation_id number for the sim_inspiral row, ${START} is the GPS start time of the ASCII waveform file, and ${DURATION} is the duration of the file in seconds.

The user should inspect the waveforms. For a waveform plotting executable see section Plot ASCII waveform files with pycbc_plot_hwinj.

Checks for the hardware injection output 

Here are some follow-up checks the user can do.

Plot ASCII waveform files with `pycbc_plot_hwinj`

You can plot the ASCII waveform files with an X11 connection. It’s strongly recommended to use the X11 connection instead of saving a static image of the entire waveform. The X11 connection allows the user to zoom in and inspect the waveform more closely. A basic inspection would include checking the amplitude, the tapering, and the ringdown of the waveforms are reasonable. For the pycbc_generate_hwinj example above one would do

pycbc_plot_hwinj --input-file hwinjcbc_${START}_H1.txt --output-file ${OUTPUT_PATH}

where ${OUTPUT_PATH} is the path to the output plot.

If you are using ssh to log into a cluster, you can provide the -Y option to open an X11 connection. For example

gsissh -Y ldas-pcdev1.ligo.caltech.edu

Recover software injection with `pycbc_inspiral`

The executable pycbc_generate_hwinj will create an XML file with both a sim_inspiral and sngl_inspiral table. Therefore we can inject the exact waveform parameters and recover them with the exact template.

The analogous software injection command for the example above would be

TMPLTBANK_FILE=hwinjcbc_${START}.xml.gz
INSPIRAL_FILE=H1-INSPIRAL_PYCBC-${GPS_START_TIME}-$((${GPS_END_TIME}-${GPS_START_TIME})).hdf
pycbc_inspiral --segment-end-pad 64  --segment-length 256 --segment-start-pad 64 --psd-estimation median --psd-segment-length 16 --psd-segment-stride 8 --psd-inverse-length 16 --pad-data 8 --sample-rate 4096 --low-frequency-cutoff 40 --strain-high-pass 30 --filter-inj-only --processing-scheme cpu --cluster-method template --approximant SEOBNRv2 --order 8 --snr-threshold 5.5 --chisq-bins 16 --channel-name ${CHANNEL_NAME} --gps-start-time ${GPS_START_TIME} --gps-end-time ${GPS_END_TIME} --trig-start-time $(($GEOCENT_END_TIME - 2)) --trig-end-time $(($GEOCENT_END_TIME + 2)) --frame-type ${FRAME_TYPE} --injection-file ${TMPLTBANK_FILE}  --bank-file ${TMPLTBANK_FILE} --output ${INSPIRAL_FILE} --verbose

Where ${START} is the start of the injection. We kept the same PSD options (eg. --psd-segment-length, etc.), data, high-pass filter, and low-frequency-cutoff.

You can print out the recovered SNR and other parameters as follows

echo `python -c "import numpy;from pycbc.io.hdf import SingleDetTriggers; \
h1_triggers=SingleDetTriggers('${INSPIRAL_FILE}', 'H1'); \
imax=numpy.argmax(h1_triggers.snr); max_snr=h1_triggers.snr[imax]; \
time=h1_triggers.end_time[imax]; print(time, max_snr)"`

Recover ASCII file injection with `pycbc_inspiral`

There is an executable pycbc_insert_frame_hwinj that will read the single-column ASCII file and insert it into frame data. An example command is here

HWINJ_FILE=hwinjcbc_${START}_H1.txt
pycbc_insert_frame_hwinj --frame-type ${FRAME_TYPE} --channel-name H1:${CHANNEL_NAME} --gps-start-time $((${GPS_START_TIME} - 16)) --gps-end-time $((${GPS_END_TIME} + 16)) --pad-data 8 --strain-high-pass 30.0 --sample-rate 16384 --hwinj-file ${HWINJ_FILE} --hwinj-start-time ${START} --ifo H1 --output-file H1-HWINJ.gwf

Where ${START} is the start time of the injection.

Then you can run pycbc on the output frame file H1-HWINJ.gwf.

How to query the segment database 

Here is an example on how to check if the detector was in science mode for a GPS time interval. To do this we query the segment database. A command line tool to do check the pycbc_generate_hwinj example above is

ligolw_segment_query_dqsegdb --query-segments --segment-url https://dqsegdb5.phy.syr.edu --gps-start-time 1124380361 --gps-end-time 1124382409 --include-segments L1:DMT-ANALYSIS_READY:1 --output-file L1-SEGMENTS.xml
ligolw_segment_query_dqsegdb --query-segments --segment-url https://dqsegdb5.phy.syr.edu --gps-start-time 1124380361 --gps-end-time 1124382409 --include-segments H1:DMT-ANALYSIS_READY:1 --output-file H1-SEGMENTS.xml

This should write two XML files L1-SEGMENTS.xml and another H1-SEGMENTS.xml. You can check the segment table to see if the detector was in science mode for this time. A command line tool that helps is

ligolw_print --table segment --column start_time --column end_time L1-SEGMENTS.xml
ligolw_print --table segment --column start_time --column end_time H1-SEGMENTS.xml

The output should be 1124380361,1124382409 for both ligolw_print commands. This tells us that the detector was in science mode for the entire time since there is one segment that is the equal to the interval of --gps-start-time to --gps-end-time.

Do not assume that the segment databse URL and science-mode segment names are the same as in this section. These values are correct for this example.

How to query the LDR server 

Here is an example on how to check if frame files exist for a GPS time interval. To do this we query the LDR server. A command line tool to do check the pycbc_generate_hwinj example above is

gw_data_find --observatory L --type L1_RDS --gps-start-time 1124380361  --gps-end-time 1124382409 --url-type file --gaps
gw_data_find --observatory H --type H1_RDS --gps-start-time 1124380361  --gps-end-time 1124382409 --url-type file --gaps

A list of frame files will be print to your terminal if they are accessible. If Missing segments is printed, then you will not be able to access all the frame files.

Do not assume that the frame type and channel name are the same as in this section. These values are correct for this example.