LALSimUtils

Top-level utility module for LALSimulation interfacing.

A collection of useful data analysis routines built from the SWIG wrappings of LAL and LALSimulation.

class RIFT.lalsimutils.ChooseWaveformParams(phiref=0.0, deltaT=0.000244140625, m1=1.9884098706980507e+31, m2=1.9884098706980507e+31, s1x=0.0, s1y=0.0, s1z=0.0, s2x=0.0, s2y=0.0, s2z=0.0, fmin=40.0, fref=0.0, dist=3.085677581491367e+22, incl=0.0, lambda1=0.0, lambda2=0.0, waveFlags=None, nonGRparams=None, ampO=0, phaseO=7, approx=84, theta=0.0, phi=0.0, psi=0.0, tref=0.0, radec=False, detector='H1', deltaF=None, fmax=0.0, taper=0, eccentricity=0.0, meanPerAno=0.0, E0=0.0, p_phi0=0.0, a6c=10000.0)

Bases: object

Class containing all the arguments needed for SimInspiralChooseTD/FDWaveform plus parameters theta, phi, psi, radec to go from h+, hx to h(t)

if radec==True: (theta,phi) = (DEC,RA) and strain will be computed using

XLALSimDetectorStrainREAL8TimeSeries

if radec==False: then strain will be computed using a simple routine

that assumes (theta,phi) are spherical coord. in a frame centered at the detector

FrequencyAtVelocity(v)

HardAlignedQ()

Test if L,S1,S2 all parallel to z

OrbitalAngularMomentumAtReference()

OrbitalAngularMomentumAtReferenceOverM2()

SoftAlignedQ()

Test if L,S1,S2 all parallel to one another

TotalAngularMomentumAtReference()

TotalAngularMomentumAtReferenceOverM2()

VelocityAtFrequency(f)

Xi()

assign_param(p, val)

Syntatic sugar to assign parameters to values. Provides ability to specify a binary by its values of

• mchirp, eta

• system frame parameters

VERY HELPFUL if you want to change just one parameter at a time (e.g., for Fisher )

copy()

Create a deep copy, so copy and original can be changed separately This does NOT work reliably (segfault side effects possible)

copy_lsctables_sim_inspiral(row)

Fill this ChooseWaveformParams with the fields of a row of an lsctables.SimInspiral table (i.e. SimInspiral table in the format as read from a file)

NB: SimInspiral table does not contain deltaT, deltaF, fref, fmax, lambda1, lambda2, waveFlags, nonGRparams, or detector fields, but ChooseWaveformParams does have these fields. This function will not alter these fields, so their values will be whatever values the instance previously had.

Adapted from code by Chris Pankow

copy_sim_inspiral(row)

Fill this ChooseWaveformParams with the fields of a row of a SWIG wrapped lalmetaio.SimInspiral table

NB: SimInspiral table does not contain deltaT, deltaF, fref, fmax, lambda1, lambda2, waveFlags, nonGRparams, or detector fields, but ChooseWaveformParams does have these fields. This function will not alter these fields, so their values will be whatever values the instance previously had.

create_sim_inspiral()

Create a sim_inspiral table from P. One element from it

extract_param(p)

assign_param Syntatic sugar to extract parameters to values. Necessary to make full use of assign_param on

• mchirp, eta

• system frame parameters

VERY HELPFUL if you want to change just one parameter at a time (e.g., for Fisher )

extract_system_frame(verbose=False)

Extract system frame angles. Returned in precisely the format needed for use in init_via_system_frame. P.init_via_system_frame(P.extract_system_frame()) should leave the state unchanged. CHECK DONE IN : LIGO-Development-UpdatedSpinningWaveforms/KISTI-MCMC/PrecessingFisher/Python/overlap_versus_systemframe PROBLEM: Polarization angle isn’t stable (oddly?)

init_via_system_frame(thetaJN=0.0, phiJL=0.0, theta1=0.0, theta2=0.0, phi12=0.0, chi1=0.0, chi2=0.0, psiJ=0.0)

Populate spin directions (and psi_L) using system frame parameters. Note this is NOT using eta Example: P.init_via_system_frame(thetaJN=0.1, phiJL=0.1, theta1=0.1, theta2=0.1, phi12=0.1, chi1=1., chi2=1., psiJ=0.)

manual_copy()

print_params(show_system_frame=False)

Print all key-value pairs belonging in the class instance

randomize(zero_spin_Q=False, aligned_spin_Q=False, volumetric_spin_prior_Q=False, default_inclination=None, default_phase=None, default_polarization=None, dMax=500.0, dMin=1.0, sMax=1)

scale_to_snr(new_SNR, psd, ifo_list, analyticPSD_Q=True, **kwargs)

• evaluates network SNR in the ifo list provided (assuming constant psd for all..may change)

• uses network SNR to rescale the distance of the source, so the SNR is now new_SNR

• returns current_SNR, for sanity

swap_components()

to_lal_dict()

to_lal_dict_extended(extra_args_dict=None)

to_lal_dict_extended:

Extended implementation, to address massive proliferation of options to control waveform version See https://git.ligo.org/lscsoft/bilby_pipe/-/merge_requests/371/diffs

RIFT.lalsimutils.ChooseWaveformParams_array_to_hdf5(P_list, fname='injections', deltaT=0.000244140625, fref=0.0, waveFlags=None, nonGRparams=None, detector='H1', deltaF=None, fMax=0.0)

HDF5 storage for parameters. Compare to lalinference HDF5 i/o https://github.com/lscsoft/lalsuite/blob/master/lalinference/python/lalinference/io/hdf5.py

TERRIBLE CODE: Assumes hardcoded order, does not embed metadata with field names SHOULD RESTRUCTURE to load into record-array like structure

RIFT.lalsimutils.ChooseWaveformParams_array_to_xml(P_list, fname='injections', minrow=None, maxrow=None, deltaT=0.000244140625, fref=0.0, waveFlags=None, nonGRparams=None, detector='H1', deltaF=None, fMax=0.0)

Standard XML storage for parameters. Note that lambda values are NOT stored in xml table entries — I have a special hack to do this

class RIFT.lalsimutils.ComplexIP(fLow=10.0, fMax=None, fNyq=2048.0, deltaF=0.125, psd=<built-in function SimNoisePSDaLIGOZeroDetHighPower>, analyticPSD_Q=True, inv_spec_trunc_Q=False, T_spec=0.0, waveform_is_psi4=False)

Bases: InnerProduct

Complex-valued inner product. self.ip(h1,h2) computes

fNyq

2 int h1(f) h2*(f) / Sn(f) df

-fNyq

And similarly for self.norm(h1)

N.B. DOES NOT assume h1, h2 are Hermitian - they should contain negative

and positive freqs. packed as

[ -N/2 * df, …, -df, 0, df, …, (N/2-1) * df ] DOES NOT maximize over time or phase

ip(h1, h2, include_epoch_differences=False)

Compute inner product between two COMPLEX16Frequency Series Accounts for time shfit

ip_matrix(listA, listB, chunk=262144)

Batched form of ip(): returns the matrix M[a,b] = self.ip(listA[a], listB[b]).

listA, listB are sequences of COMPLEX16FrequencySeries on the same 2-sided grid. The double loop over (a,b) is a single matrix product,

M = 2 df . conj(A) . (B . W)^T , A[a,f] = listA[a](f), B[b,f] = listB[b](f)

which turns O(Na.Nb) separate full-length reductions into one pass over the data plus a GEMM, and skips the frequency bins where the band weight is exactly zero. Accumulated in chunks of chunk bins so the working set stays bounded.

NOT bit-identical to the loop: the reduction order changes (pairwise np.sum over the full array vs. blocked GEMM accumulation over the band). Measured deviation is ~1e-15 relative to max|M| – see DESIGN_precompute_crossterm_batching.md. Callers that need the shipped rounding must keep using ip().

norm(h)

Compute norm of a COMPLEX16Frequency Series

class RIFT.lalsimutils.ComplexOverlap(fLow=10.0, fMax=None, fNyq=2048.0, deltaF=0.125, psd=<built-in function SimNoisePSDaLIGOZeroDetHighPower>, analyticPSD_Q=True, inv_spec_trunc_Q=False, T_spec=0.0, full_output=False, interpolate_max=False, waveform_is_psi4=False)

Bases: InnerProduct

Inner product maximized over time and polarization angle. This inner product does not assume Hermitianity and is therefore valid for waveforms that are complex in the TD, e.g. h+(t) + 1j hx(t). self.IP(h1,h2) computes:

fNyq

max 2 Abs int h1*(f,tc) h2(f) / Sn(f) df

tc -fNyq

h1, h2 must be COMPLEX16FrequencySeries defined in [-fNyq, fNyq-deltaF] At least one of which should be non-Hermitian for the maximization over phase to work properly.

If self.full_output==False: returns

The maximized overlap

If self.full_output==True: returns

The maximized overlap The entire COMPLEX16TimeSeries of overlaps for each possible time shift The index of the above time series at which the maximum occurs The phase rotation which maximizes the real-valued overlap

ip(h1, h2, **kwargs)

Compute inner product between two non-Hermitian COMPLEX16FrequencySeries

norm(h)

Compute norm of a non-Hermitian COMPLEX16FrequencySeries

wrap_times()

Return a vector of wrap-around time offsets, i.e. [ 0, dt, 2 dt, …, N dt, -(N-1) dt, -(N-1) dt, …, -2 dt, -dt ]

This is useful in conjunction with the ‘full_output’ option to plot the overlap vs timeshift. e.g. do:

IP = ComplexOverlap(full_output=True) t = IP.wrap_times() rho, ovlp, rhoIdx, rhoPhase = IP.ip(h1, h2) plot(t, abs(ovlp))

RIFT.lalsimutils.CreateCompatibleComplexIP(hlmf, **kwargs)

Creates complex IP (no maximization)

RIFT.lalsimutils.CreateCompatibleComplexOverlap(hlmf, **kwargs)

CreateCompatibleComplexOverlap: accepts dictionary or single instance of COMPLEX16FrequencySeries

RIFT.lalsimutils.DataFourier(ht)

RIFT.lalsimutils.DataFourierREAL8(ht)

RIFT.lalsimutils.DataInverseFourier(hf)

RIFT.lalsimutils.DataInverseFourierREAL8(hf)

RIFT.lalsimutils.DataRollBins(ht, nL)

RIFT.lalsimutils.DataRollTime(ht, DeltaT)

RIFT.lalsimutils.Fcross(theta, phi, psi)

Antenna pattern as a function of polar coordinates measured from directly overhead a right angle interferometer and polarization angle

RIFT.lalsimutils.Fplus(theta, phi, psi)

Antenna pattern as a function of polar coordinates measured from directly overhead a right angle interferometer and polarization angle

class RIFT.lalsimutils.HermitianComplexIP(fLow=10.0, fMax=None, fNyq=2048.0, deltaF=0.125, psd=<built-in function SimNoisePSDaLIGOZeroDetHighPower>, analyticPSD_Q=True, inv_spec_trunc_Q=False, T_spec=0.0, waveform_is_psi4=False)

Bases: InnerProduct

Complex-valued inner product. self.ip(h1,h2) computes

fNyq

4 int h1(f) h2*(f) / Sn(f) df

fLow

And similarly for self.norm(h1)

N.B. Assumes h1, h2 are Hermitian - i.e. they store only positive freqs.

with negative freqs. given by h(-f) = h*(f)

DOES NOT maximize over time or phase

ip(h1, h2)

Compute inner product between two COMPLEX16Frequency Series

norm(h)

Compute norm of a COMPLEX16Frequency Series

class RIFT.lalsimutils.InnerProduct(fLow=10.0, fMax=None, fNyq=2048.0, deltaF=0.125, psd=<built-in function SimNoisePSDaLIGOZeroDetHighPower>, analyticPSD_Q=True, inv_spec_trunc_Q=False, T_spec=0.0, waveform_is_psi4=False)

Bases: object

Base class for inner products

ip(h1, h2)

Compute inner product between two COMPLEX16Frequency Series

norm(h)

Compute norm of a COMPLEX16Frequency Series

RIFT.lalsimutils.Mceta(m1, m2)

Compute chirp mass and symmetric mass ratio from component masses

RIFT.lalsimutils.NINJA_data_to_hoff(fname, TDlen=0, scaleT=1.0, scaleH=1.0, Fp=1.0, Fc=0.0)

Function to read in data in the NINJA format, i.e. t_i h+(t_i) hx(t_i) and convert it to a COMPLEX16FrequencySeries holding the observed h(f) = FFT[ Fp*h+(t) + Fc*hx(t) ]

If TDlen == -1, do not zero-pad the TD waveform before FFTing If TDlen == 0 (default), zero-pad the TD waveform to the next power of 2 If TDlen == N, zero-pad the TD waveform to length N before FFTing

scaleT and scaleH can be used to rescale the time steps and strain resp. e.g. to get a waveform appropriate for a total mass M you would scale by scaleT = G M / c^3 scaleH = G M / (D c^2)

RIFT.lalsimutils.NINJA_data_to_hoft(fname, TDlen=-1, scaleT=1.0, scaleH=1.0, Fp=1.0, Fc=0.0)

Function to read in data in the NINJA format, i.e. t_i h+(t_i) hx(t_i) and convert it to a REAL8TimeSeries holding the observed h(t) = Fp*h+(t) + Fc*hx(t)

If TDlen == -1 (default), do not zero-pad the returned waveforms If TDlen == 0, zero-pad returned waveforms to the next power of 2 If TDlen == N, zero-pad returned waveforms to length N

scaleT and scaleH can be used to rescale the time steps and strain resp. e.g. to get a waveform appropriate for a total mass M you would scale by scaleT = G M / c^3 scaleH = G M / (D c^2)

RIFT.lalsimutils.NINJA_data_to_hp_hc(fname, TDlen=-1, scaleT=1.0, scaleH=1.0, deltaT=0)

Function to read in data in the NINJA format in file ‘fname’, i.e. t_i h+(t_i) hx(t_i) and convert it to two REAL8TimeSeries holding polarizations hp(t) and hc(t)

If TDlen == -1 (default), do not zero-pad the returned waveforms If TDlen == 0, zero-pad returned waveforms to the next power of 2 If TDlen == N, zero-pad returned waveforms to length N

scaleT and scaleH can be used to rescale the time steps and strain resp. e.g. if the input file provides time steps in t/M and the strain has mass and distance scaled out, to get a waveform appropriate for a total mass M and distance D you would scale by scaleT = G M / c^3 scaleH = G M / (D c^2)

Once time is properly scaled into seconds, you can interpolate the waveform to a different sample rate deltaT. e.g. if the file has time steps of t/M = 1 and you use scaleT to rescale for a 100 Msun binary, the time steps will be ~= 0.00049 s. If you provide the argument deltaT = 1./4096. ~= 0.00024 s the waveform will be interpolated and resampled at 4096 Hz.

If deltaT==0, then no interpolation will be done

NOTE: For improved accuracy, we convert h+ + i hx –> A e^(i Phi), interpolate and resample A and Phi, then convert back to h+, hx

RIFT.lalsimutils.NINJA_data_to_norm_hoff(fname, IP, TDlen=0, scaleT=1.0, scaleH=1.0, Fp=1.0, Fc=0.0)

Function to read in data in the NINJA format, i.e. t_i h+(t_i) hx(t_i) and convert it to a COMPLEX16FrequencySeries holding h(f) = FFT[ Fp*h+(t) + Fc*hx(t) ] normalized so (h|h)=1 for inner product IP

If TDlen == -1, do not zero-pad the TD waveform before FFTing If TDlen == 0 (default), zero-pad the TD waveform to the next power of 2 If TDlen == N, zero-pad the TD waveform to length N before FFTing

scaleT and scaleH can be used to rescale the time steps and strain resp. e.g. to get a waveform appropriate for a total mass M you would scale by scaleT = G M / c^3 scaleH = G M / (D c^2)

class RIFT.lalsimutils.Overlap(fLow=10.0, fMax=None, fNyq=2048.0, deltaF=0.125, psd=<built-in function SimNoisePSDaLIGOZeroDetHighPower>, analyticPSD_Q=True, inv_spec_trunc_Q=False, T_spec=0.0, full_output=False)

Bases: InnerProduct

Inner product maximized over time and phase. self.ip(h1,h2) computes:

fNyq

max 4 Abs int h1*(f,tc) h2(f) / Sn(f) df

tc fLow

h1, h2 must be COMPLEX16FrequencySeries defined in [0, fNyq] (with the negative frequencies implicitly given by Hermitianity)

If self.full_output==False: returns

The maximized (real-valued, > 0) overlap

If self.full_output==True: returns

The maximized overlap The entire COMPLEX16TimeSeries of overlaps for each possible time shift The index of the above time series at which the maximum occurs The phase rotation which maximizes the real-valued overlap

ip(h1, h2)

Compute inner product between two Hermitian COMPLEX16Frequency Series

norm(h)

Compute norm of a COMPLEX16Frequency Series

wrap_times()

Return a vector of wrap-around time offsets, i.e. [ 0, dt, 2 dt, …, N dt, -(N-1) dt, -(N-1) dt, …, -2 dt, -dt ]

This is useful in conjunction with the ‘full_output’ option to plot the overlap vs timeshift. e.g. do:

IP = Overlap(full_output=True) t = IP.wrap_times() rho, ovlp, rhoIdx, rhoPhase = IP.ip(h1, h2) plot(t, abs(ovlp))

RIFT.lalsimutils.RangeProtect(fn, val)

RangeProtect(fn) wraps fn, returning - np.inf in cases where an argument is not finite, and calling the function otherwise. Intended to be used with test_coord_output and coordinate conversion routines, to identify range errors, etc.

This function assumes fn acts on EACH ELEMENT INDIVIDUALLY, and does not reduce the dimension. Not to use.

RIFT.lalsimutils.RangeProtectReduce(fn, val)

RangeProtect(fn) wraps fn, returning - np.inf in cases where an argument is not finite, and calling the function otherwise. Intended to be used with test_coord_output and coordinate conversion routines, to identify range errors, etc.

This function assumes fn acts to REDUCE the data to one dimension (e.g., a list of points)

class RIFT.lalsimutils.RealIP(fLow=10.0, fMax=None, fNyq=2048.0, deltaF=0.125, psd=<built-in function SimNoisePSDaLIGOZeroDetHighPower>, analyticPSD_Q=True, inv_spec_trunc_Q=False, T_spec=0.0, waveform_is_psi4=False)

Bases: InnerProduct

Real-valued inner product. self.ip(h1,h2) computes

fNyq

4 Re int h1(f) h2*(f) / Sn(f) df

fLow

And similarly for self.norm(h1)

DOES NOT maximize over time or phase

ip(h1, h2)

Compute inner product between two COMPLEX16Frequency Series

norm(h)

Compute norm of a COMPLEX16Frequency Series

RIFT.lalsimutils.SphHarmFrequencySeries_to_dict(hlms, Lmax)

Convert a SphHarmFrequencySeries SWIG-wrapped linked list into a dictionary.

The keys are tuples of integers of the form (l,m) and a key-value pair will be created for each (l,m) with 2 <= l <= Lmax, |m| <= l for which lalsimulation.SphHarmFrequencySeriesGetMode(hlms, l, m) returns a non-null pointer.

RIFT.lalsimutils.SphHarmTimeSeries_to_dict(hlms, Lmax)

Convert a SphHarmTimeSeries SWIG-wrapped linked list into a dictionary.

The keys are tuples of integers of the form (l,m) and a key-value pair will be created for each (l,m) with 2 <= l <= Lmax, |m| <= l for which lalsimulation.SphHarmTimeSeriesGetMode(hlms, l, m) returns a non-null pointer.

RIFT.lalsimutils.VectorToFrame(vecRef)

Convert vector to a frame, by picking relatively arbitrary vectors perp to it. Used to convert J to a frame. spin_convention == radiation:

Uses the radiation frame conventions to do so., so the ‘z’ direction is orthogonal to the 2nd vector and ‘n’ (the line of sight) defines things

spin_convention == L:

Trickier, because we will want to define angles of L relative to J ! Pick a random alternate vector

RIFT.lalsimutils.Yagi13_fit_barlamdel(barlam2, ell)

Yagi 2013 fits for multipolar $bar{lambda}_ell$ = 2 k_ell/(C^{2ell+1} (2ell-1)!!)$ Eq.(10),(61); Tab.I; Fig.8 http://arxiv.org/abs/1311.0872

RIFT.lalsimutils.Yagi13_fitcoefs(ell)

Coefficients of Yagi 2013 fits for multipolar $bar{lambda}_ell = 2 k_ell/(C^{2ell+1} (2ell-1)!!)$ Tab.I (NS) http://arxiv.org/abs/1311.0872

RIFT.lalsimutils.barlamdel_to_kappal(q, barlaml, ell)

$kappa^{A,B}_ell(bar{lambda}_ell)$ Assume $q=M_A/M_B>=1$

RIFT.lalsimutils.check_FD_pending(code)

RIFT.lalsimutils.chieff_aligned_spins(m1, m2, s1z, s2z, xi)

s1z, s2z with xi = (m1 s1z + m2 s2z)/M and chiminus = (m2 s1z - m1 s2z)/M held from the inputs.

The inverse used by assign_param(‘chieff_aligned’); accepts arrays. Exchange-antisymmetric; for large mass ratio s1z -> xi, s2z -> -chiminus. Can give |s_iz| > 1.

RIFT.lalsimutils.complex_hoff(P, sgn=-1, fwdplan=None, **kwargs)

CURRENTLY ONLY WORKS WITH TD APPROXIMANTS

Generate a (non-Hermitian) FD waveform from ChooseWaveformParams P by creating a complex TD waveform of the form

h(t) = h+(t) + 1j sgn hx(t) where sgn = -1 (default) or 1

If P.deltaF==None, just pad up to next power of 2 If P.deltaF = 1/X, will generate a TD waveform, zero-pad to length X seconds

and then FFT. Will throw an error if waveform is longer than X seconds

If you do not provide a forward FFT plan, one will be created. If you are calling this function many times, it is best to create it once beforehand and pass it in, e.g.: fwdplan=lal.CreateForwardCOMPLEX16FFTPlan(TDlen,0)

Returns a COMPLEX16FrequencySeries object

RIFT.lalsimutils.complex_hoft(P, sgn=-1, **kwargs)

Generate a complex TD waveform from ChooseWaveformParams P Returns h(t) = h+(t) + 1j sgn hx(t) where sgn = -1 (default) or 1

Returns a COMPLEX16TimeSeries object

RIFT.lalsimutils.complex_hoft_IMRPv2(P_copy, sgn=-1)

Generate a complex TD waveform from ChooseWaveformParams P, using SimInspiralTD Returns h(t) = h+(t) + 1j sgn hx(t) where sgn = -1 (default) or 1

Returns a COMPLEX16TimeSeries object

Designed specifically to help interface with IMRPhenomPv2

RIFT.lalsimutils.complex_norm_hoff(P, IP, sgn=-1, fwdplan=None)

CURRENTLY ONLY WORKS WITH TD APPROXIMANTS

Generate a (non-Hermitian) FD waveform from ChooseWaveformParams P by creating a complex TD waveform of the form

h(t) = h+(t) + 1j sgn hx(t) where sgn = -1 (default) or 1

If P.deltaF==None, just pad up to next power of 2 If P.deltaF = 1/X, will generate a TD waveform, zero-pad to length X seconds

and then FFT. Will throw an error if waveform is longer than X seconds

If you do not provide a forward FFT plan, one will be created. If you are calling this function many times, it is best to create it once beforehand and pass it in, e.g.: fwdplan=lal.CreateForwardCOMPLEX16FFTPlan(TDlen,0)

Returns a COMPLEX16FrequencySeries object

RIFT.lalsimutils.conj_hlmoff(P, Lmax=2)

RIFT.lalsimutils.convert_waveform_coordinates(x_in, coord_names=['mc', 'eta'], low_level_coord_names=['m1', 'm2'], enforce_kerr=False, source_redshift=0, fref=None)

A wrapper for ChooseWaveformParams() ‘s coordinate tools (extract_param, assign_param) providing array-formatted coordinate changes. BE VERY CAREFUL, because coordinates may be defined inconsistently (e.g., holding different variables constant: M and eta, or mc and q). Note that if ChooseWaveformParam structuers are built ,the loops can be quite slow

Special cases:

• coordinates in x_out already in x_in are copied over directly

• xi==chi_eff, chiMinus, mu1,mu2 : transformed directly from mc, delta_mc, s1z,s2z coordinates, using fast vectorized transformations.

• source_redshift: if nonzero, convert m1 -> m1 (1+z)=m_z, as fit is done in the detector frame. We are assuming source-frame sampling

• rf_* (geometric4) coordinates: computed from the physical components; require fref.

RIFT.lalsimutils.convert_waveform_coordinates_with_eos(x_in, coord_names=['mc', 'eta'], low_level_coord_names=['m1', 'm2'], enforce_kerr=False, eos_class=None, no_matter1=False, no_matter2=False, source_redshift=0, backstop_novector=False)

A wrapper for ChooseWaveformParams() ‘s coordinate tools (extract_param, assign_param) providing array-formatted coordinate changes. BE VERY CAREFUL, because coordinates may be defined inconsistently (e.g., holding different variables constant: M and eta, or mc and q)

RIFT.lalsimutils.deltalamtilde_of_eta_lam1_lam2(eta, lam1, lam2)

This is the definition found in Les Wade’s paper. Les has factored out the quantity sqrt(1-4eta). It is different from Marc Favata’s paper. $deltatildeLambda(eta, Lambda_1, Lambda_2)$. Lambda_1 is assumed to correspond to the more massive (primary) star m_1. Lambda_2 is for the secondary star m_2.

RIFT.lalsimutils.eccentric_anomaly_from_mean(mean_anomaly, eccentricity)

RIFT.lalsimutils.enforce_swig_psd_fmin(raw_psd, fmin)

RIFT.lalsimutils.estimateDeltaF(P, LmaxEff=2)

Input: P Output:estimated duration (in s) based on Newtonian inspiral from P.fmin to infinite frequency

RIFT.lalsimutils.estimateWaveformDuration(P, LmaxEff=2)

Input: P Output:estimated duration (in s) based on Newtonian inspiral from P.fmin to infinite frequency

RIFT.lalsimutils.eta_crit(Mc, m2_min)

RIFT.lalsimutils.evaluate_fvals(lal_2sided_fseries, lal_convention=False)

Associates frequencies with a 2sided lal complex array. Compare with ‘self.longweights’ code Done by HAND in PrecessingOrbitModesOfFrequency Manually reverses convention re sign of omega used in lal!

lal_convention=False (default): RIFT’s own two-sided packing (arrays produced by

DataFourier/complex_hoff): even length, REVERSED, f[k] = deltaF*(npts/2 - k). WARNING: for odd-length arrays this assigns f on a grid offset by deltaF/2 – it is only correct for the even-length RIFT packing.

lal_convention=True: the packing of two-sided series returned directly by LAL

generators (e.g. SimInspiralChooseFDModes): ASCENDING [-fNyq, …, 0, …, +fNyq] with DC at index npts//2 (odd length TDlen+1; even length after truncation keeps DC at npts//2). f[k] = deltaF*(k - npts//2), exact for both parities of npts. Use this - not the default - on ChooseFDModes output: the default’s reversal + half-bin offset shifts any |f|-symmetric window by one bin between (l,m) and (l,-m) modes (which live at opposite signs of f), breaking the conjugate-pair identity h_{l,-m}(f) = (-1)^l conj(h_{lm}(-f)) at the percent level.

Notes:

1. XXXFrequencySeries have an f0 and a deltaF, and logically they should run from f0….f0+N*df

2. I will always use COMPLEX16FrequencySeries, which should run from -fNyq…fNyq-df

   with fNyq=fSample/2.=1/(2*deltaT)

3. In practice, I have a REVERSED frequency series associated…WHY? (finish doc review) Probably an overall sign difference in omega t vs -omega t in FFT definition? LAL documentation isn’t clear with the overall sign.

   lal: int dt exp( -i omega t) F(t) = tilde{F}(omega) [lal] me: int dt exp( i omega t) F(t) = tilde{F}(omega)

   following Poisson and Will, my old ‘DataFourier’ code, etc

4. The low-level code just rotates the original vector right by half the length – so the usual FFTW

   code data is mapped right.

Checks:

P=lalsimutils.ChooseWaveformParams() hf = lalsimutils.complex_hoff(P) lalsimutils.evaluate_fvals(hf) hf.f0

Notes:

‘’Physics convention<math> int dt exp(-i omega t) F(t) </math>

• my notes use -2 pi i f t convention

• my NR papers text uses int dt exp(-i omega t) convention

‘’Other convention (PN-convenient) <math> int dt exp(i omega t) F(t) </math>

• Prakash notes use opposite convention

• low-level mathematica code (DataFourier) uses int dt exp(i omega t) convention

• Poisson and Will uses opposite convention http://arxiv.org/pdf/gr-qc/9502040v1.pdf, as do all subsequent PN work

• Cutler and Flanagan use the opposite convention

RIFT.lalsimutils.evaluate_tvals(lal_tseries)

RIFT.lalsimutils.extend_swig_psd_series_to_sampling_requirements(raw_psd, dfRequired, fNyqRequired)

extend_psd_series_to_sampling_requirements: Takes a conventional 1-sided PSD and extends into a longer 1-sided PSD array by filling intervening samples. Currently only works by ‘interpolating’ between samples. Assumes fNyq does not change. Assumes the input and output are swig REAL8Timeseries objects

RIFT.lalsimutils.extract_JL_angles(P, return_inverse=True, phase_factor=1, opening_angle=None)

extract_JL_angles: finds Euler angles for transformation see eg C13, C14 in 2004.06503

Using directly-coded J,L expressions rather than calling the lalsuite routine, to make sure phases agree with what I intend.

RIFT.lalsimutils.findDeltaF(P, **kwargs)

Given ChooseWaveformParams P, generate the TD waveform, round the length to the next power of 2, and find the frequency bin size corresponding to this length. This is useful b/c deltaF is needed to define an inner product which is needed for norm_hoft and norm_hoff functions

RIFT.lalsimutils.frame_data_to_hoff(fname, channel, start=None, stop=None, TDlen=0, window_shape=0.0, verbose=True, **kwargs)

Function to read in data in the frame format and convert it to a COMPLEX16FrequencySeries holding h(f) = FFT[ h(t) ]

Before the FFT, applies a Tukey window with shape parameter ‘window_shape’. N.B. if window_shape=0, the window is the identity function

if window_shape=1, the window becomes a Hann window if 0<window_shape<1, the data will transition from zero to full

strength over that fraction of each end of the data segment.

If TDlen == -1, do not zero-pad the TD waveform before FFTing If TDlen == 0 (default), zero-pad the TD waveform to the next power of 2 If TDlen == N, zero-pad the TD waveform to length N before FFTing

RIFT.lalsimutils.frame_data_to_hoft(fname, channel, start=None, stop=None, window_shape=0.0, verbose=True, deltaT=None, deltaT_internal=None, upsample_method='cubic', use_gwpy=False, **kwargs)

Function to read in data in the frame format and convert it to a REAL8TimeSeries. fname is the path to a LIGO cache file.

Applies a Tukey window to the data with shape parameter ‘window_shape’. N.B. if window_shape=0, the window is the identity function

if window_shape=1, the window becomes a Hann window if 0<window_shape<1, the data will transition from zero to full

strength over that fraction of each end of the data segment.

Modified to rely on the lalframe read_timeseries function

https://github.com/lscsoft/lalsuite/blob/master/lalframe/python/lalframe/frread.py

RIFT.lalsimutils.frame_data_to_hoft_old(fname, channel, start=None, stop=None, window_shape=0.0, verbose=True)

Function to read in data in the frame format and convert it to a REAL8TimeSeries. fname is the path to a LIGO cache file.

Applies a Tukey window to the data with shape parameter ‘window_shape’. N.B. if window_shape=0, the window is the identity function

if window_shape=1, the window becomes a Hann window if 0<window_shape<1, the data will transition from zero to full

strength over that fraction of each end of the data segment.

RIFT.lalsimutils.frame_data_to_non_herm_hoff(fname, channel, start=None, stop=None, TDlen=0, window_shape=0.0, verbose=True, deltaT=None, **kwargs)

Function to read in data in the frame format and convert it to a COMPLEX16FrequencySeries h(f) = FFT[ h(t) ] Create complex FD data that does not assume Hermitianity - i.e. contains positive and negative freq. content

Before the FFT, applies a Tukey window with shape parameter ‘window_shape’. N.B. if window_shape=0, the window is the identity function

if window_shape=1, the window becomes a Hann window if 0<window_shape<1, the data will transition from zero to full

strength over that fraction of each end of the data segment.

If TDlen == -1, do not zero-pad the TD waveform before FFTing If TDlen == 0 (default), zero-pad the TD waveform to the next power of 2 If TDlen == N, zero-pad the TD waveform to length N before FFTing

RIFT.lalsimutils.get_intp_psd_series_from_xmldoc(fname, inst)

RIFT.lalsimutils.get_psd_series_from_xmldoc(fname, inst)

RIFT.lalsimutils.guess_mc_range(event_dict, force_mc_range=None)

RIFT.lalsimutils.guess_mc_range_mdc(P, force_mc_range=None)

RIFT.lalsimutils.hdf5_to_ChooseWaveformParams_array(fname='injections', deltaT=0.000244140625, fref=0.0, waveFlags=None, nonGRparams=None, detector='H1', deltaF=None, fMax=0.0)

HDF5 storage for parameters. Compare to lalinference HDF5 i/o https://github.com/lscsoft/lalsuite/blob/master/lalinference/python/lalinference/io/hdf5.py

TERRIBLE CODE: Assumes hardcoded order, does not embed metadata with field names SHOULD RESTRUCTURE to load into record-array like structure

RIFT.lalsimutils.hlmoff(P, Lmax=2)

Generate the FD h_lm -2-spin-weighted spherical harmonic modes of a GW with parameters P. Returns a SphHarmTimeSeries, a linked-list of modes with a COMPLEX16TimeSeries and indices l and m for each node.

The linked list will contain all modes with l <= Lmax and all values of m for these l.

RIFT.lalsimutils.hlmoft(P, Lmax=2, nr_polarization_convention=False, fixed_tapering=False, silent=True, fd_standoff_factor=0.964, no_condition=False, fd_L_frame=False, fd_centering_factor=0.9, fd_alignment_postevent_time=None, fd_end_taper_max_time=1.0, **kwargs)

Generate the TD h_lm -2-spin-weighted spherical harmonic modes of a GW with parameters P. Returns a SphHarmTimeSeries, a linked-list of modes with a COMPLEX16TimeSeries and indices l and m for each node.

The linked list will contain all modes with l <= Lmax and all values of m for these l.

nr_polarization_convention : apply a factor -1 to all modes provided by lalsuite

fd_standoff_factor: for ChooseFDWaveform, reduce starting frequency P.fmin by this factor internally when generating.

FD waveform calculation implicitly assumes fmin is in the inspiral regime (if not you are making bad life choices). Normal default value of 0.964 is based on changing inspiral duration by 10% (1.1^(3/8) ~ 1.036). Any constant factor will change the inspiral duration by a factor (1./fd_standoff_factor)^(8/3) or so.

WARNING: internal_hlm_generator will CHANGE this factor to 0.9, to better match complex_hoft. I’m keeping the default of this function here to match the nominal documentaiton of SimInspiralTD, while internal_hlm_generator has a different factor, chosen to match its actual output.

no_condition: disable conditioning: the ChooseFDModes high-pass/taper, the start tapers of the SEOB, SEOBv3 and IMRPhenomP branches, and the start taper applied to ChooseTDModes / TDModesFromPolarizations output at the end of this function. Used by the td-gs likelihood code path (no windows, no tapers) and for diagnostic plots on impact of/need for conditioning. Does not touch lalsuite’s own internal conditioning of FD-generated models. fd_L_frame: rotate hlmoft to L frame from J frame for return values from ChooseFDModes (XO4/XPHM). fd_end_taper_max_time: cap (seconds) on the ChooseFDModes end taper, which is otherwise 1% of the

buffer. With fd_alignment_postevent_time=2 an uncapped end taper covers the merger once the buffer exceeds ~185 s. None restores the uncapped 1%. Has no effect for buffers <= 100 s.

RIFT.lalsimutils.hlmoft_FromFD_dict(P, Lmax=2)

Uses Chris Pankow’s interface in lalsuite Do not redshift the source

RIFT.lalsimutils.hlmoft_IMRPv2_dict(P, sgn=-1)

Reconstruct hlm(t) from h(t,nhat) evaluated in 5 different nhat directions, specifically for IMRPhenomPv2 Using SimInspiralTD evaluated at nhat=zhat (zhat=Jhat), nhat=-zhat, and at three points in the equatorial plane. There is an analytic solution, but for code transparency I do the inverse numerically rather than code it in

ISSUES

• lots of instability in recovery depending on point set used (e.g., degeneracy in aligned spin case). Need to choose intelligently

• duration of hlm is longer than deltaF (because of SimInspiralTD) in general.

RIFT.lalsimutils.hlmoft_SEOB_dict(P, Lmax=2)

Generate the TD h_lm -2-spin-weighted spherical harmonic modes of a GW with parameters P. Returns a dictionary of modes. Just for SEOBNRv2 SEOBNRv1, and EOBNRv2. Uses aligned-spin trick to get (2,2) and (2,-2) modes. A hack. Works for any aligned-spin time-domain waveform with only (2,pm 2) modes though.

RIFT.lalsimutils.hlmoft_SEOBv3_dict(P, Lmax=2)

Generate the TD h_lm -2-spin-weighted spherical harmonic modes of a GW with parameters P. Returns a dictionary of modes. A hack for SEOBNRv3, because there’s not a natural dictionary output

Will NOT perform the -1 factor correction

RIFT.lalsimutils.hlmoft_SpinTaylorManual_dict(P, Lmax=2)

Generate the TD h_lm -2-spin-weighted spherical harmonic modes of a GW with parameters P. Returns a dictionary of modes. Just for SpinTaylor, using manual interface provided by Riccardo and Jake.

RIFT.lalsimutils.hoff(P, Fp=None, Fc=None, fwdplan=None, **kwargs)

Generate a FD waveform from ChooseWaveformParams P. Will return a COMPLEX16FrequencySeries object.

If P.approx is a FD approximant, hoff_FD is called. This path calls SimInspiralChooseFDWaveform

fwdplan must be None for FD approximants.

If P.approx is a TD approximant, hoff_TD is called. This path calls ChooseTDWaveform and performs an FFT.

The TD waveform will be zero-padded so it’s Fourier transform has frequency bins of size P.deltaT. If P.deltaF == None, the TD waveform will be zero-padded to the next power of 2.

RIFT.lalsimutils.hoff_FD(P, Fp=None, Fc=None, **kwargs)

Generate a FD waveform for a FD approximant. Note that P.deltaF (which is None by default) must be set

RIFT.lalsimutils.hoff_TD(P, Fp=None, Fc=None, fwdplan=None, **kwargs)

Generate a FD waveform from ChooseWaveformParams P by creating a TD waveform, zero-padding and then Fourier transforming with FFTW3 forward FFT plan fwdplan

If P.deltaF==None, just pad up to next power of 2 If P.deltaF = 1/X, will generate a TD waveform, zero-pad to length X seconds

and then FFT. Will throw an error if waveform is longer than X seconds

If you do not provide a forward FFT plan, one will be created. If you are calling this function many times, you may to create it once beforehand and pass it in, e.g.: fwdplan=lal.CreateForwardREAL8FFTPlan(TDlen,0)

You may pass in antenna patterns Fp, Fc. If none are provided, they will be computed from the information in ChooseWaveformParams

Returns a COMPLEX16FrequencySeries object

RIFT.lalsimutils.hoft(P, Fp=None, Fc=None, **kwargs)

Generate a TD waveform from ChooseWaveformParams P You may pass in antenna patterns Fp, Fc. If none are provided, they will be computed from the information in ChooseWaveformParams.

Returns a REAL8TimeSeries object

RIFT.lalsimutils.hoft_from_hlm(hlms, P, return_complex=False, extra_phase_shift=0)

hoft_from_hlm(hlm,P,Lmax): return hoft like output given hlmoft input. Important if hoft is not accessible (e.g., not provided by lalsuite)

RIFT.lalsimutils.hoft_to_frame_data(fname, channel, hoft)

Function to write h(t) to a frame file in a specific channel. NOT YET DEBUGGED

RIFT.lalsimutils.lam1_lam2_of_pe_params(eta, lamt, dlamt)

lam1 is for the the primary mass m_1. lam2 is for the the secondary mass m_2. m_1 >= m2.

RIFT.lalsimutils.lamtilde_of_eta_lam1_lam2(eta, lam1, lam2)

$tildeLambda(eta, Lambda_1, Lambda_2)$. Lambda_1 is assumed to correspond to the more massive (primary) star m_1. Lambda_2 is for the secondary star m_2.

RIFT.lalsimutils.load_resample_and_clean_psd(psd_fname, det, deltaF, verbose=False)

RIFT.lalsimutils.lsu_StringFromPNOrder(order)

RIFT.lalsimutils.m1m2(Mc, eta)

Compute component masses from Mc, eta. Returns m1 >= m2

RIFT.lalsimutils.mass1(Mc, eta)

Compute larger component mass from Mc, eta

RIFT.lalsimutils.mass2(Mc, eta)

Compute smaller component mass from Mc, eta

RIFT.lalsimutils.mchirp(m1, m2)

Compute chirp mass from component masses

RIFT.lalsimutils.mean_anomaly_from_true(true_anomaly, eccentricity)

RIFT.lalsimutils.modes_to_k(modes)

RIFT.lalsimutils.nextPow2(length)

Find next power of 2 <= length

RIFT.lalsimutils.nhat(th, ph)

RIFT.lalsimutils.non_herm_hoff(P)

Generate a FD waveform with two-sided spectrum. i.e. not assuming the Hermitian symmetry applies

RIFT.lalsimutils.norm_hoff(P, IP, Fp=None, Fc=None, fwdplan=None)

Generate a normalized FD waveform from ChooseWaveformParams P. Will return a COMPLEX16FrequencySeries object.

If P.approx is a FD approximant, norm_hoff_FD is called. This path calls SimInspiralChooseFDWaveform

fwdplan must be None for FD approximants.

If P.approx is a TD approximant, norm_hoff_TD is called. This path calls ChooseTDWaveform and performs an FFT.

The TD waveform will be zero-padded so it’s Fourier transform has frequency bins of size P.deltaT. If P.deltaF == None, the TD waveform will be zero-padded to the next power of 2.

RIFT.lalsimutils.norm_hoff_FD(P, IP, Fp=None, Fc=None)

Generate a FD waveform for a FD approximant normalized according to IP. Note that P.deltaF (which is None by default) must be set. IP and the waveform generated from P must have the same deltaF and

the waveform must extend to at least the highest frequency of IP’s PSD.

RIFT.lalsimutils.norm_hoff_TD(P, IP, Fp=None, Fc=None, fwdplan=None)

Generate a waveform from ChooseWaveformParams P normalized according to inner product IP by creating a TD waveform, zero-padding and then Fourier transforming with FFTW3 forward FFT plan fwdplan. Returns a COMPLEX16FrequencySeries object.

If P.deltaF==None, just pad up to next power of 2 If P.deltaF = 1/X, will generate a TD waveform, zero-pad to length X seconds

and then FFT. Will throw an error if waveform is longer than X seconds

If you do not provide a forward FFT plan, one will be created. If you are calling this function many times, you may to create it once beforehand and pass it in, e.g.: fwdplan=lal.CreateForwardREAL8FFTPlan(TDlen,0)

You may pass in antenna patterns Fp, Fc. If none are provided, they will be computed from the information in ChooseWaveformParams.

N.B. IP and the waveform generated from P must have the same deltaF and

the waveform must extend to at least the highest frequency of IP’s PSD.

RIFT.lalsimutils.polar_angles_in_frame(frm, vec)

Take a vector in the default frame. Evaluate the polar angles of that unit vector in a new (orthonormal) frame ‘frm’. Not the fastest solution. Not naturally vectorizable.

RIFT.lalsimutils.polar_angles_in_frame_alt(frm, theta, phi, xpy=<module 'numpy' from '/opt/hostedtoolcache/Python/3.11.18/x64/lib/python3.11/site-packages/numpy/__init__.py'>)

Take polar angles in the default frame. Evaluate the polar angles of that unit vector in a new (orthonormal) frame ‘frmInverse’. Probably easier to vectorize

RIFT.lalsimutils.psd_windowing_factor(window_shape, TDlen)

[Implementing <w^2>, due to https://dcc.ligo.org/LIGO-T1900249] If a PSD is calculated using repeated instances of windowed data, the PSD will be biased downward (i.e., we have less noise than we expect) That’s fine and self-consistent if we will be analyzing data which has the same windowing applied. But if we’re not applying the same window shape to the data we analyze, we need to correct for the difference in noise actually present to the noise used to create the PSD.

Ths routine creates <w^2> for w a window shape. This could be computed much more efficiently (it is not memory-efficient), but this code is cleaner

RIFT.lalsimutils.pylal_psd_to_swig_psd(raw_pylal_psd)

Why do I do a conversion? I am having trouble returning modified PSDs

RIFT.lalsimutils.regularize_swig_psd_series_near_nyquist(raw_psd, DeltaFToZero)

regularize_psd_series_near_nyquist Near nyquist, some PSD bins are insane. As a temporary measure, I am going to use this routine to explicitly zero out those frequency bins, using a 30 Hz window near nyquist

RIFT.lalsimutils.resample_psd_series(psd, df=None, fmin=None, fmax=None)

RIFT.lalsimutils.rotate_hlm_static(hlm, alphaA, betaA, gammaA, extra_polarization=None)

Rotate output hlm by these Euler angles Note this is a full system rotation, and therefore should be associated with rotation of the physical quantities (e.g., spin vectors, J,L)

extra_polarization is applied at the end, as a multiplicative factor, if nonzero

RIFT.lalsimutils.rotation_matrix(axis, theta)

RIFT.lalsimutils.safe_int(mystr)

RIFT.lalsimutils.sanitize_eta(eta, tol=1e-10, exception='error')

If ‘eta’ is slightly outside the physically allowed range for symmetric mass ratio, push it back in. If ‘eta’ is further outside the physically allowed range, throw an error or return a special value. Explicitly:

• If ‘eta’ is in [tol, 0.25], return eta.

• If ‘eta’ is in [0, tol], return tol.

• If ‘eta’ in is (0.25, 0.25+tol], return 0.25

• If ‘eta’ < 0 OR eta > 0.25+tol,

  • if exception==’error’ raise a ValueError

  • if exception is anything else, return exception

RIFT.lalsimutils.singleIFOSNR(data, psd, fNyq, fmin=None, fmax=None)

Calculate single IFO SNR using inner product class.

RIFT.lalsimutils.std_and_conj_hlmoff(P, Lmax=2, **kwargs)

RIFT.lalsimutils.stringGPSNice(tgps)

Return a string with nice formatting displaying the value of a LIGOTimeGPS

RIFT.lalsimutils.symRatio(m1, m2)

Compute symmetric mass ratio from component masses

RIFT.lalsimutils.symmetry_sign_exchange(coord_names)

RIFT.lalsimutils.test_coord_output(x_out)

Checks if any of the x_out are -np.inf. Returns a boolean array [ True, False, False, …] with True if the corresponding coordinate row is ok, false otherwise

RIFT.lalsimutils.tidal_lambda_from_tilde(mass1, mass2, lam_til, dlam_til)

Determine physical lambda parameters from effective parameters.

RIFT.lalsimutils.tidal_lambda_tilde(mass1, mass2, lambda1, lambda2)

‘Effective’ lambda parameters. Lackey et al https://arxiv.org/pdf/1402.5156.pdf, Eq. (5,6).

RIFT.lalsimutils.true_anomaly_from_eccentric(eccentric_anomaly, eccentricity)

RIFT.lalsimutils.unit_frame()

RIFT.lalsimutils.unwind_phase(phase, thresh=5.0)

Unwind an array of values of a periodic variable so that it does not jump discontinuously when it hits the periodic boundary, but changes smoothly outside the periodic range.

Note: ‘thresh’, which determines if a discontinuous jump occurs, should be somewhat less than the periodic interval. Empirically, 5 is usually a safe value of thresh for a variable with period 2 pi.

Fast method: take element-by-element differences, use mod 2 pi, and then add

RIFT.lalsimutils.vecCross(v1, v2)

RIFT.lalsimutils.vecDot(v1, v2)

RIFT.lalsimutils.vecUnit(v)

RIFT.lalsimutils.xml_to_ChooseWaveformParams_array(fname, minrow=None, maxrow=None, deltaT=0.000244140625, fref=0.0, lambda1=0.0, lambda2=0.0, waveFlags=None, nonGRparams=None, detector='H1', deltaF=None, fmax=0.0)

Function to read an xml file ‘fname’ containing a SimInspiralTable, convert rows of the SimInspiralTable into ChooseWaveformParams instances and return an array of the ChooseWaveformParam instances

Can optionally give ‘minrow’ and ‘maxrow’ to convert only rows in the range (starting from zero) [minrow, maxrow). If these arguments are not given, this function will convert the whole SimInspiral table.

The rest of the optional arguments are the fields in ChooseWaveformParams that are not present in SimInspiral tables. Any of these arguments not given values will use the standard default values of ChooseWaveformParams.