Reproducing an Event Analysis

Reproducing a gravitational-wave event analysis means recovering a declared result from specified public inputs with documented numerical checks. Begin here with an executable calculation: correlate one supplied waveform against the GW150914 strain in each LIGO detector, recover its timing peak, and investigate the noise estimate. The calculation is small enough to inspect line by line. A later section develops the more demanding protocol for GW170817 parameter inference, which is a proposed project rather than a completed result on this page.

Helpful background. Gravitational-Wave Astronomy supplies noise-weighted inner products, PSD estimation, detector response, and Bayesian likelihoods. Relativistic Cosmology supplies the distance and redshift frame conventions used in released parameters. Black-Hole Astrophysics supplies waveform regimes and physically informed source priors.

A single-template GW150914 calculation

A matched filter asks, “At what time does this waveform best line up with the data, after giving less weight to noisy frequencies?” The GW150914 release, strain-data and tutorial sections supplies both detector records and a numerical-relativity waveform. We use its V2, 4096 Hz, 32-second H1 and L1 files, beginning at GPS 1126259446. V2 matters: it incorporates a calibration revision and corrects a roughly millisecond timing offset in the older 4096 Hz files GWOSC, technical details, “Version 2 files”.

This waveform is SXS:BBH:0305, rescaled to detector-frame total mass $74.6\,M_\odot$. Its amplitude and phase were adjusted using the event data. Consequently, this exercise teaches the filtering calculation and recovers an event feature; it does not independently discover GW150914, estimate its masses, or reproduce its false-alarm probability. Maximizing over a time interval chosen around a known event is also not a blind search.

The calculation below was executed on 2026-09-23 using this historical release. The volume’s observational evidence cutoff remains 2026-08-18; the new run does not revise the catalog or the published interpretation of the event.

Download these four files into one directory:

• Python calculation and independent checks.
• Input URLs, hashes, provenance, and analysis settings.
• Pinned Python dependencies.
• Reference numerical results.

Use a Python 3.9–3.12 interpreter with the pinned dependencies. The reference run used macOS 27.0 on arm64, CPython 3.9.6, NumPy 2.0.2, SciPy 1.13.1, and h5py 3.14.0. The Python file includes its MIT license; the input specification records the separate data attribution and terms. For example, in a terminal in that directory:

python3 -m venv .venv
.venv/bin/python -m pip install -r gw150914-matched-filter-requirements.txt
.venv/bin/python gw150914-matched-filter.py --self-test
.venv/bin/python gw150914-matched-filter.py \
  --compare gw150914-matched-filter-results.json --output result.json

The script downloads the three inputs from GWOSC into a temporary-directory cache, checks their SHA-256 hashes and the strain files’ released MD5 checksums, and writes the small result file you requested. For an existing download, use --cache-dir /path/to/data --offline. Without --output, JSON is printed to the terminal. There are no random draws. The reference results include hashes of the code and input specification, so changing the calculation requires a new reference rather than silently passing the old comparison.

Both files have data-quality mask 127. The injection masks are 31 for H1 and 23 for L1: the latter records a continuous-wave hardware injection, while the bits indicating absence of transient hardware injections are set in both files. A bit value of one means no injection of that class GWOSC, technical details, “Injections for the event data releasers”. Requiring every injection mask to equal 31 would incorrectly reject the supplied Livingston record.

What the executable calculation computes

Write $\Delta t=1/4096\,\mathrm{s}$ and $N=131072$. The frequency spacing of the full record is $\Delta f=1/(N\Delta t)=1/32\,\mathrm{Hz}$. For real sampled data $d_j$, the site’s Fourier convention gives

$$D_k=\Delta t\sum_{j=0}^{N-1} d_j e^{+2\pi i jk/N}.$$

NumPy’s forward FFT has the opposite exponential sign, so the code’s dt*rfft(data) is $D_k^*$, not $D_k$ NumPy 2.0 manual, “Implementation details”. This difference is harmless when translated consistently; it is consequential when interpreting the sign of a recovered delay.

Let $H_k$ be the transform of the padded template and $S_k$ the one-sided strain-noise PSD. Over the selected positive-frequency bins $\mathcal B$, define

$$\begin{aligned} \sigma_h^2&=4\Delta f\sum_{k\in\mathcal B} \frac{|H_k|^2}{S_k},\\ z_\ell&=\frac{4\Delta f}{\sigma_h} \sum_{k\in\mathcal B} \frac{D_k^*H_k}{S_k} e^{+2\pi i k\ell/N}, \qquad \rho_\ell=|z_\ell|. \end{aligned}$$

The real and imaginary parts correlate against two orthogonal templates with equal norm: the supplied real waveform and a quadrature waveform whose positive-frequency transform is $-iH_k$. The latter shifts each retained frequency by a quarter cycle; conjugate negative-frequency components make the time-domain waveform real. Taking the modulus maximizes over their relative phase. This is a precisely defined two-dimensional template family, not a general model of the physical polarizations or phase dependence of arbitrary binary waveforms. For a noise-free copy $d_j=A h_{j-m}$, the expected peak is at lag $m$, with $\rho_m=|A|\sigma_h$. Changing the arbitrary overall template amplitude therefore leaves the normalized peak statistic unchanged.

The inverse FFT supplies the sum over lags efficiently. The code fills only the positive-frequency bins before calling ifft; mirroring the bins would discard the imaginary quadrature. Its prefactor compensates for the inverse FFT’s $1/N$ normalization. DC and Nyquist are omitted because the chosen band is far from both endpoints.

Three practical choices are equally important:

1. Estimate noise away from the event. Relative to the record start, use $[0,14)$ and $[18,32)$ seconds. The baseline divides these intervals into four-second segments with 50% overlap, subtracts each segment’s mean, applies a periodic Hann window, and averages 12 periodograms. No segment crosses the excluded interval. For a segment of length $M$ and window $w_j$, its interior-frequency periodogram is

   $$\widehat S_k= \frac{2\Delta t}{\sum_{j=0}^{M-1}w_j^2} \left|\sum_{j=0}^{M-1}w_j(d_j-\bar d) e^{-2\pi i jk/M}\right|^2.$$

   The denominator corrects for the window’s squared amplitude; the factor two combines positive and negative frequencies for real data. The PSD is linearly interpolated from its $0.25\,\mathrm{Hz}$ grid onto the full-record grid. Interpolation adds no independent spectral information.

2. Control endpoint discontinuities. Subtract the full strain record’s mean and apply a Tukey window with parameter $0.125$, tapering approximately two seconds at each edge. It equals one over the event waveform. Smooth the supplied template’s first $0.1\,\mathrm{s}$ and last $0.005\,\mathrm{s}$, then pad it to the record length. These operations reduce spectral leakage from the large low-frequency strain and the finite waveform; they do not make a finite record exactly stationary. The numerical result uses this specified finite-record approximation to the stationary-noise inner product.

3. Translate lag into the waveform’s reference time. The first template sample has time $t_0=-0.620404109\,\mathrm{s}$ relative to the release waveform’s peak reference. The padded array places that sample at index zero. A correlation lag $\ell$ therefore corresponds to reference time $\ell\Delta t-t_0$ after the strain record starts, not simply $\ell\Delta t$. Search only the declared interior interval $[16.2,16.6]$ seconds. Circular FFT correlation is not a trustworthy search near the record boundaries.

For background on the statistic and its noise assumptions, see Likelihood Ratios and Matched Filtering and Power Spectral Density and Nonstationarity.

Recovered peaks and sensitivity to analysis choices

The executed baseline, with band 35–350 Hz, gives:

Detector	Peak statistic $\rho$	Reference-time offset from GPS 1126259462
Hanford H1	19.890864	$0.422894343\,\mathrm{s}$
Livingston L1	14.175285	$0.415814265\,\mathrm{s}$

The Hanford peak is 29 samples later, or $7.080078125\,\mathrm{ms}$. Thus this calculation recovers the detector ordering and approximate separation illustrated in the GW150914 release, Fig. 1 caption. These digits identify a reproducible grid point; they are not an uncertainty estimate. The grid spacing is $0.244140625\,\mathrm{ms}$, and waveform, noise, and calibration uncertainties are not included here.

The same script changes one control at a time:

PSD segment	Band (Hz)	Averaged segments	H1 peak $\rho$	L1 peak $\rho$	H1 minus L1 (ms)
4 s, baseline	35–350	12	19.891	14.175	7.080
2 s	35–350	26	19.189	14.059	7.080
8 s	35–350	4	21.222	15.480	7.324
4 s	30–350	12	20.162	14.205	7.080
4 s	40–350	12	19.467	13.856	7.080
4 s	35–500	12	19.892	14.178	7.080

Longer PSD segments resolve frequency more finely but leave fewer averages in this short record. Their noisy inverse PSD can move the statistic substantially. The eight-second case is therefore not automatically a better answer, and this table is a sensitivity study, not a convergent sequence or confidence interval. The timing changes by at most one sample across these choices; the peak statistic is less stable. Longer off-source data and a justified treatment of noise-spectrum uncertainty would be necessary for a more precise analysis. Overlapping periodograms are also not independent observations.

Independent checks of the implementation

Before reading the event files, the executable tests the calculation on small arrays whose answers can be obtained another way:

• A sum of three sinusoids has an analytic noise-weighted norm. Both its amplitude factors and the one-sided PSD normalization must agree with that norm.
• An explicit positive-exponent discrete Fourier sum checks every lag of the FFT correlation, while a separate time-domain correlation checks its real part for a flat PSD. Parseval’s equality compares time- and frequency-domain energy.
• Array indexing inserts copies delayed by $+37$ and $-29$ samples, with amplitudes $0.7$ and $-1.3$. Recovering both signed delays and amplitudes checks the Fourier sign and distinguishes a negative waveform from a negative delay.
• The hand-written periodogram average is compared with the explicitly configured SciPy 1.13.1 welch implementation.

All six relative discrepancies in the reference run are below $3\times10^{-15}$, against a required tolerance $10^{-10}$. The rerun comparison allows an absolute difference $10^{-7}$ in the statistic and $10^{-10}\,\mathrm{s}$ in reported grid time. Those are numerical reproduction tolerances, not physical error bars. Passing them checks this finite calculation; it does not validate a source model, Gaussian-noise assumption, detection significance, or posterior distribution.

A later GW170817 inference project

Start by choosing a claim that the available products can test. The official GW170817 release provides versioned strain and links to released parameter-estimation samples; its persistent dataset identifier is DOI 10.7935/K5B8566F. A defensible open-data target is therefore:

Using one declared strain branch and detector set, recover selected released signal and posterior features—such as coalescence time, detector-frame chirp mass, and luminosity-distance structure—within predeclared numerical and Monte Carlo tolerances.

This target is narrower than reproducing the published false-alarm rate. Search significance requires the search configuration, background construction, trials, vetoes, and analyzed live time, not just an event-centered strain segment. The original discovery paper identifies the observed network, signal, and source interpretation Abbott et al. 2017, §§I–IV, but its rounded tables are not a substitute for released machine-readable targets.

Freeze the following before computation:

• event and product revision, detector set, GPS interval, sample rate, and before- or after-noise-subtraction branch;
• byte hashes, data-quality and hardware-injection flags, calibration notes, valid band, and event-time convention;
• conditioning operations, analysis segment, window, high-pass or band-pass filters, resampling method, and PSD data;
• waveform family and version, reference frequency, parameter frame, calibration treatment, and likelihood approximations;
• every prior and bound, sampling transformation, marginalization, sampler settings, stopping rule, random seeds, runtime, package lock, and numerical precision;
• reference posterior release, sample weights, comparison parameters, and acceptance tolerances.

The product selection is explained in Open Gravitational-Wave Data Products. In particular, the GW170817 page exposes distinct strain branches with different glitch and noise treatment. Mixing the PSD from one branch with strain from another changes the likelihood and is not a harmless file substitution.

From strain to a likelihood

For detector $a$, define the noise-weighted inner product over the declared band,

$$(x\mid y)_a =4\,\operatorname{Re} \int_{f_{\min}}^{f_{\max}} \frac{\widetilde x_a^*(f)\widetilde y_a(f)} {S_{n,a}(f)}\,df.$$

Under the stationary Gaussian-noise approximation, the network log likelihood is

$$\ln\mathcal L(d\mid\theta) =-\frac12\sum_a \bigl(d_a-h_a(\theta)\mid d_a-h_a(\theta)\bigr)_a +C,$$

where $h_a$ includes arrival-time delay, antenna response, calibration convention, and waveform parameters. This likelihood is a model for the conditioned data, not a statement that detector noise is globally stationary or Gaussian. The standard Bayesian parameter-estimation construction and its implementation checks are described by Veitch et al. 2015, §§II–IV.

A reproducible deterministic chain is:

1. decode the frozen files and verify detector, GPS grid, sample count, units, hash, and quality state;
2. apply the declared conditioning in a fixed order, accounting for filter transients and excised intervals;
3. estimate or load $S_{n,a}(f)$ from the declared data and record averaging, window normalization, overlap, frequency resolution, and interpolation;
4. construct $h_a(\theta)$ in the same detector/source-frame and time/phase conventions;
5. test the inner product and likelihood on analytic vectors and on a known software or waveform injection;
6. run the sampler only after the deterministic log likelihood agrees across two implementations or against an independent benchmark;
7. transform and compare samples using the same parameter definitions and weights as the reference release.

The open-source Bilby framework is one documented route for defining gravitational-wave likelihoods and sampler interfaces Ashton et al. 2019, §§2–4. Using it does not itself establish reproduction: the data branch, PSD, waveform, priors, marginalizations, and versioned dependencies still determine the result.

Checks before sampling

Sampling can conceal a wrong deterministic calculation. Complete these checks first:

Fourier and normalization check. Parseval consistency, one-sided PSD factors, FFT cadence, window normalization, and frequency-bin width must agree. Verify a sinusoid and Gaussian vector whose norm can be calculated independently.

Time and detector check. Inject a waveform with known geocentric time and sky direction, project it into each detector, and recover the expected relative arrival times and phases. A detector-name swap can leave a plausible network statistic while corrupting localization.

Signal-recovery check. For a zero-noise injection $d=h(\theta_0)$, the likelihood should be maximized at $\theta_0$ along identifiable directions. Degenerate directions may remain broad, and prior boundaries can move posterior summaries even when the likelihood is correct.

Branch and glitch check. Inspect the release notes and a time–frequency representation around all documented transients. A gate, subtraction, or altered interval must be represented in both the data model and the PSD choice; it cannot be hidden as generic “cleaning.”

Residual check. At a representative posterior point, compare the residual spectrum with the PSD over the analysis band and inspect detector-by-detector contributions to $\ln\mathcal L$. A network total can look stable while one detector is misnormalized.

These checks are independent only to the extent that they do not reuse the same faulty implementation. Generating and recovering an injection with the same projection or FFT routine is useful regression testing, but an analytic norm, alternate library, or independently decoded format is needed to challenge that routine.

Convergence and comparison

A single sampler run is not a convergence study. Vary the controls that can move the result:

Refinement axis	Stable quantities to inspect	Failure signature
Sample rate and frequency band	Network signal-to-noise ratio, log-likelihood differences, selected posterior quantiles	Drift from aliasing, omitted signal power, or an invalid release band
PSD duration, segment length, overlap, and window	Per-detector norms, residual spectrum, mass and distance summaries	Signal contamination, nonstationarity, or window-normalization error
Analysis duration and conditioning	Coalescence time, chirp mass, residuals	Edge transients, glitch leakage, or a mismatched strain branch
Waveform resolution or approximant	Match, likelihood, masses, inclination, distance	Modeling uncertainty larger than the claimed reproduction tolerance
Sampler resolution, stopping rule, and seed	Effective samples, evidence when relevant, quantiles and correlations	Multimodal loss, prior-boundary sticking, or stochastic undercoverage
Calibration draws or response variants	Per-detector phase/amplitude and distance	Instrument uncertainty omitted from an apparently precise posterior

Choose tolerances from independent replicate runs and numerical refinements, not from whether the published value falls inside after the fact. For posterior quantile $q$, a stochastic diagnostic is

$$z_q =\frac{\left|q_{\rm run}-q_{\rm ref}\right|} {\sqrt{\sigma_{{\rm MC,run}}^2+\sigma_{{\rm MC,ref}}^2}},$$

where both Monte Carlo errors must account for weights and correlation. This is not a complete scientific distance: also compare marginal and joint structure, tail probability, support, and physically important degeneracies. A normalized one-dimensional Wasserstein distance or an energy distance on standardized joint samples can supplement quantiles, but its threshold must be calibrated on same-model replicates.

For waveform or reconstructed-signal comparison, use the match

$$\mathcal M =\max_{t_c,\phi_c} \frac{(h_1\mid h_2)} {\sqrt{(h_1\mid h_1)(h_2\mid h_2)}}, \qquad \text{mismatch}=1-\mathcal M,$$

with PSD and band stated. A high match does not guarantee compatible distance or inclination posteriors; it measures agreement in one noise-weighted signal geometry.

Interpreting outcomes

The outcome should be classified, not compressed into “reproduced” or “failed”:

• Artifact failure: a hash, schema, time grid, flag, or release identity is wrong.
• Deterministic-method failure: norms, injections, projections, or likelihood values disagree before sampling.
• Stochastic failure: independent chains or seeds do not converge within the predeclared Monte Carlo envelope.
• Model difference: both computations converge, but a waveform, prior, PSD, calibration treatment, or data branch explains a stable posterior shift.
• Unavailable-evidence boundary: the public products cannot test a search background, private calibration state, or production detail required by the original claim.
• Compatible reproduction: frozen inputs pass integrity and deterministic checks, refinements are stable, and selected signal and posterior diagnostics meet the predeclared tolerance.

Compatibility is claim-specific. Recovering selected posterior marginals does not reproduce evidence values, sky localization, counterpart association, or population selection unless those were separately targeted and their products supplied.

The GW170817 protocol above is not an executed parameter-inference result. Its future realization must freeze the selected release, priors, likelihood, software, and reference samples and then report its own checks. The completed GW150914 calculation establishes only the narrower filtering result described earlier. Refresh either calculation after a relevant product correction, calibration revision, software change, or independent check that reveals a missing systematic.

Common pitfalls

Matching a paper’s rounded median. One matching number can occur despite the wrong posterior shape, prior, or detector contribution. Compare machine-readable distributions and deterministic likelihood checkpoints.

Changing several controls at once. If PSD, waveform, prior, and sampler all change, a posterior shift has no diagnosis. Establish a same-model baseline and introduce one consequential variation at a time.

Calling a same-code injection independent. It checks internal consistency, not a shared sign, unit, or normalization error. Add an analytic or alternate-implementation benchmark.

Exercises

1. Translate a correlation lag. The H1 maximum occurs at template-start lag 64727 samples. Using $4096\,\mathrm{Hz}$ and $t_0=-0.620404109\,\mathrm{s}$, recover its reference-time offset from GPS 1126259462. What mistake results from treating the lag as the waveform peak time?

Solution

The file begins 16 seconds before that GPS reference, so

$$\frac{64727}{4096}-(-0.620404109)-16 =0.422894343375\,\mathrm{s}.$$

Omitting $-t_0$ reports the first waveform sample’s alignment instead, shifting the claimed reference time earlier by about $0.6204\,\mathrm{s}$. The inter-detector difference can still look correct because the same error cancels in the subtraction. Correct relative timing alone does not establish the absolute reference convention.

2. Convert waveform mismatch to a likelihood scale. Two equal-norm, time- and phase-aligned waveforms each have optimal signal-to-noise ratio $\rho=20$ and match $\mathcal M=0.995$. In zero noise, compute their squared difference and the log-likelihood loss when one is treated as the model for the other.

Solution

For equal norms,

$$(h_1-h_2\mid h_1-h_2) =2\rho^2(1-\mathcal M) =2(20)^2(0.005)=4.$$

The Gaussian log-likelihood difference is

$$\Delta\ln\mathcal L=-\frac12(4)=-2.$$

Thus a match that looks very close to unity can still produce a meaningful likelihood penalty at appreciable signal-to-noise ratio. This calculation holds for the declared PSD, band, and maximizations; it does not quantify parameter bias by itself.

3. Compare a posterior quantile. A reference quantile is $q_{\rm ref}=1.20$ with Monte Carlo error $0.03$, while a rerun gives $q_{\rm run}=1.28$ with Monte Carlo error $0.04$. Evaluate $z_q$ and interpret it.

Solution

The combined stochastic scale is

$$\sqrt{0.03^2+0.04^2}=0.05,$$

so

$$z_q=\frac{0.08}{0.05}=1.6.$$

The shift does not exceed a hypothetical two-standard-error stochastic screen. That is not proof of reproduction: the screen must have been declared in advance, both errors must be reliable, multiple-comparison effects may matter, and joint posterior structure plus deterministic checks must also agree.

References

• Abbott, B. P., et al. “GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral.” Physical Review Letters 119 (2017): 161101. DOI. Open PDF.
• Ashton, Gregory, et al. “BILBY: A User-Friendly Bayesian Inference Library for Gravitational-Wave Astronomy.” Astrophysical Journal Supplement Series 241 (2019): 27. DOI. Open PDF.
• Gravitational Wave Open Science Center. “Data Release for Event GW150914.” Dataset DOI 10.7935/K5MW2F23. V2 strain files and supplied SXS:BBH:0305 waveform. Accessed September 23, 2026. Release. Dataset DOI. Release checksums.
• Gravitational Wave Open Science Center. “Data Release for Event GW170817.” Dataset DOI 10.7935/K5B8566F; revision history through 2019. Accessed August 18, 2026. Release. Dataset DOI.
• Gravitational Wave Open Science Center. “Technical Details for Advanced LIGO and Virgo Event Data Releases.” Accessed September 23, 2026. Technical notes.
• NumPy Developers. “Discrete Fourier Transform (numpy.fft).” NumPy 2.0 Manual, “Implementation details.” Documentation.
• SciPy Developers. “scipy.signal.welch.” SciPy 1.13.1 Manual. Documentation.
• Veitch, John, et al. “Parameter Estimation for Compact Binaries with Ground-Based Gravitational-Wave Observations Using the LALInference Software Library.” Physical Review D 91 (2015): 042003. DOI. Open PDF.