Time of coalescence incorrect for time-domain waveforms

This is an issue I encountered when running with parallel bilby but I think the underlying cause is in Bilby itself. If one uses a time-domain model for analysis, the coalescence time (i.e. geocentric end time) can be off by almost 100 ms. I have tested this by injecting an SEOBNRv4 signal and recovering with SEOBNRv4 and SEOBNRv4_ROM. While the ROM result is close to the true value, the time-domain is way off. With @roberto.cotesta we have looked at the difference in handling the time-domain vs frequency-domain waveforms. One thing that we think is missing is that XLALSimInspiralFD does not take into account the fact that the waveform does not start at 0, but rather at a negative time such that the peak of the waveform occurs at 0. Thus, we think there is an additional rotation in the Fourier domain that should be applied. I tried a preliminary version of the code and obtained the following results for the SEOBNRv4 injection:

As one can see, the geocentric end time posterior is now much closer to the injected value than before. However, it's still not centred on the true value. Thus I think there might be something else missing.

Could the core developers take a look?

added High priority label

@serguei.ossokine is this related to #443 (closed)? I.e. is it the same issue, or different?

I'm currently struggling to understand if this is a waveforms interface issue or a bilby issue. I personally can't see a "bug" right now in bilby with how we use lalsimulation. I'm going to email the waveform and PE groups to try and get some perspective.

Hi Greg, apologies for not being clearer. I think this is more of a waveform interface issue. This issue is related but not the same as #443 (closed) . Both issue stem from using XLALSimInspiralFD. In the present case, I think this is what causes the difference in the time of coalescence for time-domain waveforms. Now lalinference calls XLALChooseTDWaveform directly and then does all the conditioning, padding, FFTing and time-shifting internally. I think it would be good to take a good look at what they do, but unfortunately I didn't have enough time to parse it.

For what it's worth, my 2¢ are that all waveforms (FD and TD) called via XLALSimInspiralFD should produce similar outputs (meaning they have a common reference time). Given that the FD waveforms do have a common reference time, what is the rationale for not having one with TD waveforms? My concern with considering this an interface-issue is that it falls on, e.g., PE-code developers to figure out what convention particular TD waveform families have, which can easily lead to errors creeping in. Because of this, I think we should consider modifying XLALSimInspiralFD so that all of the TD waveforms are aligned in time in a way that's consistent with the FD waveforms.

I'm pretty sure that XLALSimInspiralFD does not assume that waveforms start at 0. In fact, it is the convention that 0 is roughly the coalescence / peak time. Is the problem actually within the ROM? Also, please check to make sure that you use the .epoch of the returned timeseries as this can be shifted in order to compensate with phase rotation in the frequency domain that is needed to keep the waveform from wrapping in time.

Interesting, thanks @jolien-creighton ! It would be good if someone more familiar with the bilby likelihood could confirm if .epoch is used to shift the signals, because there's potentially a v. easy fix for the TD-waveform likelihoods here.

@rory-smith the wrapping of all lalsim interfaces is done in the utils: https://git.ligo.org/lscsoft/bilby/blob/master/bilby/gw/utils.py#L793

I don't see anywhere that epoch is used

Thanks @gregory.ashton I would expect that epoch would be used in the likelihood itself, as it's supposed to be used to shift the waveform to whatever epoch happens to be in the time series (I recall something similar is done in LI). Could you point me to the relevant line(s) in the likelihood function where the time shifting is done?

I believe the time-shifting happens here: https://git.ligo.org/lscsoft/bilby/blob/master/bilby/gw/detector/interferometer.py#L316

Thanks, @gregory.ashton I think the issue is on exactly line 316. I'd bet that

dt = dt_geocent + time_shift

Just needs to be changed to something like

dt = dt_geocent + time_shift ± output_of_XLALSimInspiralFD.epoch

The ± sign being there because I can never remember which way the shift is applied. Do you know a convenient way of adding output_of_XLALSimInspiralFD.epoch to the interferometer class?

I wrote a simple script below that computes the epoch for a waveform computed using SimInspiralFD and SimInspiralChooseFDWaveform. I found that the the epoch is indeed different between the two waveform interfaces. Moreover, the definition of epoch is consistent for all FD waveforms generated using SimInspiralChooseFDWaveform (but not SimInspiralFD). For reference, the standard definition in SimInspiralChooseFDWaveform is epoch=-T where T is the data/signal duration. In light of this, there's a time offset that needs to be applied to all waveforms from SimInspiralFD which is just dt = T - abs(waveform.epoch) By applying this offset, the definition of epoch should be the same between SimInspiralFD and SimInspiralChooseFDWaveform

import lal
import lalsimulation as lalsim

longAscNodes = 0
eccentricity = 0 # Assume circular orbit
meanPerAno = 0

WFdict = lal.CreateDict()

lalsim.SimInspiralWaveformParamsInsertFrameAxis(WFdict, 2)
lalsim.SimInspiralWaveformParamsInsertPNSpinOrder(WFdict, -1)
lalsim.SimInspiralWaveformParamsInsertPNTidalOrder(WFdict, -1)
lalsim.SimInspiralWaveformParamsInsertPNPhaseOrder(WFdict, -1)
lalsim.SimInspiralWaveformParamsInsertPNAmplitudeOrder(WFdict, -1)
lalsim.SimInspiralWaveformParamsInsertTidalLambda1(WFdict, 0.0)
lalsim.SimInspiralWaveformParamsInsertTidalLambda2(WFdict, 0.0)
    
s1 = [0.0, 0, 0.0]
s2 = [0.0, 0.0, 0.0]
dist = 500 # Distance in Mpc
iota = np.pi*0.4  
phi_c = 0. # Overall phase
T = 4.
deltaF = 1./T 
f_high = 4096. 
Fs = 2*f_high
deltaT = 1./Fs
f_ref = 40. # Reference frequency at which spins are calculated
f_low = 20.

approx = lalsim.IMRPhenomD # Approximation scheme in time domain

m1=30
m2=30

hplus_INSPIRAL_FD, hcross_INSPIRAL_FD = lalsim.SimInspiralFD(m1 * lal.MSUN_SI, m2 * lal.MSUN_SI,\
                                                          s1[0], s1[1], s1[2],\
                                                          s2[0], s2[1], s2[2],\
                                                          dist * 1e6 * lal.PC_SI, iota, phi_c,\
                                                          longAscNodes, eccentricity, meanPerAno,
                                                          deltaF, f_low, f_high, f_ref,\
                                                          WFdict, approx)





hplus_CHOOSE_FD, hcross_CHOOSE_FD = lalsim.SimInspiralChooseFDWaveform(m1 * lal.MSUN_SI, m2 * lal.MSUN_SI,\
                                                          0, 0., 0,\
                                                          0, 0, 0,\
                                                          dist * 1e6 * lal.PC_SI, iota, phi_c,\
                                                          longAscNodes, eccentricity, meanPerAno,\
                                                          deltaF, f_low, f_high, f_ref, WFdict, approx)
print(hplus_CHOOSE_FD.epoch)
print(hplus_INSPIRAL_FD.epoch)

For info: There are some time-domain waveforms where "t=0" does not at all correspond to merger time. Specifically, the surrogate models NRsur7dq2, NRsur7dq4 set t=0 at the start of the waveforms, some 30-40 GW cycles before merger.

Incidentally, the correct way to do this IMO is to use the LALSimulation routines XLALSimDetectorStrainREAL8TimeSeries to generate the strain on a detector and XLALSimAddInjectionREAL8TimeSeries to add the strain to the data (I guess for PE purposes, a routine that subtracted the strain from the data would be in order...). These are designed to work with the other routines in LALSimulation. The first one additionally takes into account Doppler shifting due to the Earth's rotation (which affects long duration signals) as well as the non-long-wavelength-limit effects on the response pattern at high frequencies (not too important for CBCs in aLIGO but will become more important especially for IFOs with longer arms). The second one deals with the subsample interpolation.

See also the stand-alone programs lalsim-inspiral, lalsim-detector-strain, lalsim-inject, and specifically the example given in the latter.

Hi @serguei.ossokine @gregory.ashton I implemented my patch and performed a run on with SEOBNRv4 on GW150914. The tc recovery using my patch are consistent with the IMRPhenomPv2 results. In contrast, the master branch of the code shows a significant bias the tc recovery:

The full set of comparisons plots are here on the summarypage : https://ldas-jobs.ligo.caltech.edu/~rory.smith/pbilby_review/150914/pb_li_comparison_with_SEOBNRv4_2/html/Comparison_geocent_time.html

@rory-smith is the point that pbilby_seobnrv4_fix is closer to the IMR results? Is it not concerning that it still seems to be shifted by a small fraction?

Not really. If you look at a typical event, the spread on tc from different waveform families is quite broad, e.g., https://ldas-jobs.ligo.caltech.edu/~mpitkin/projects/O3/S190521g/offline/html/Comparison_geocent_time.html

@rory-smith @gregory.ashton I think the most convincing thing would be to run SEOBNRv4 with lalinference on this event. Then we can know for sure. @rory-smith could you set up such a run?

@serguei.ossokine : @roberto.cotesta has a nice demonstration using one of the NR surrogates which shows that the LI and pbilby tc posteriors now agree.

Hi all,

I've used the NRHybSurr3dq8 both with LI MCMC and with pbilby and I get a very consistent tc. This is with Rory's patch.

I think that we can close this issue if this patch is merged to master.

The PE run is here. Unfortunately the priors for the other parameters are a bit different so we cannot really compare the two runs for the other parameters.

@gregory.ashton can we close this issue and merge the commits? It is important for the exceptional events and we have already checked that is working properly

@roberto.cotesta this was fixed in !736 (merged) and has been merged to master. Closing.

closed

@gregory.ashton thanks!

mentioned in commit b3cb81f4

mentioned in issue #551 (closed)