Adding support for SEOBNRv5EHM

dingo/asimov/dingo.ini

@@ -203,9 +203,9 @@ prior-dict-updates = {
 {%- assign p = priors['eccentricity'] %}
    eccentricity = {{ p['type'] }}(name='eccentricity', minimum={{ p['minimum'] }}, maximum={{ p['maximum'] }}),
 {%- endif %}
-{%- if priors.keys() contains "mean_anomaly" %}
-{%- assign p = priors['mean_anomaly'] %}
-   mean_anomaly = {{ p['type'] }}(name='mean_anomaly', minimum={{ p['minimum'] }}, maximum={{ p['maximum'] }}, boundary={{ p['boundary'] | default: 'periodic' }} ),
+{%- if priors.keys() contains "mean_per_ano" %}
+{%- assign p = priors['mean_per_ano'] %}
+   mean_per_ano = {{ p['type'] }}(name='mean_per_ano', minimum={{ p['minimum'] }}, maximum={{ p['maximum'] }}, boundary={{ p['boundary'] | default: 'periodic' }} ),
 {%- endif %}
 {%- if priors.keys() contains "luminosity distance" %}
 {%- assign p = priors['luminosity distance'] %}

dingo/gw/likelihood.py

@@ -14,7 +14,7 @@
     DecimateWaveformsAndASDS,
     create_mask_based_on_frequency_update,
 )
-from dingo.gw.waveform_generator import WaveformGenerator
+from dingo.gw.waveform_generator import WaveformGenerator, wfg_utils
 from dingo.gw.domains import (
     Domain,
     UniformFrequencyDomain,
@@ -382,6 +382,18 @@ def _log_likelihood_phase_grid_mode_decomposed(self, theta, phases=None):
         pol_m = self.signal_m({**theta, "phase": 0})
         pol_m = {k: pol["waveform"] for k, pol in pol_m.items()}
 
+        # For TEOB, generate_hplus_hcross_m applies a phase-dependent time shift
+        # bookkept in wfg._deferred_timeshift_data. The analytical precomputation
+        # below assumes mu(phi) = sum_m mu_m(0) * exp(-i m phi), which is only
+        # valid when that timeshift is absent. When it's present we fall back to
+        # a per-phase resum so the phase-grid likelihood stays consistent with
+        # the full direct likelihood.
+        dts = getattr(self.waveform_generator, "_deferred_timeshift_data", None)
+        if dts is not None:
+            return self._log_likelihood_phase_grid_mode_decomposed_with_dts(
+                pol_m, d, phases, dts
+            )
+
         # Step 2: Precompute complex inner products (mu, mu) and (d, mu) for the
         # individual modes m.
         min_idx = self.data_domain.min_idx
@@ -475,6 +487,56 @@ def _log_likelihood_phase_grid_mode_decomposed(self, theta, phases=None):
 
         return log_likelihoods
 
+    def _log_likelihood_phase_grid_mode_decomposed_with_dts(
+        self, pol_m, d, phases, dts
+    ):
+        """Variant of _log_likelihood_phase_grid_mode_decomposed for approximants
+        (currently TEOB) whose mode decomposition carries a phase-dependent time
+        shift recorded in `dts` (= wfg._deferred_timeshift_data).
+
+        We cannot use the analytical precomputation trick because the timeshift
+        depends on the argmax of the TD-resummed h+, which varies with phase.
+        Instead we resum modes per phase and apply the exp(-2 pi i f dt(phi))
+        correction before computing inner products. The geocenter-level time
+        shift commutes with detector projection, so we can apply it directly to
+        the per-detector signals returned by signal_m.
+        """
+        m_vals = sorted(pol_m.keys())
+        ifos = list(d.keys())
+        freqs = self.data_domain()
+
+        log_likelihoods = np.empty(len(phases))
+        for idx, phi in enumerate(phases):
+            # Analytical resum of per-detector mode contributions.
+            mu = {
+                ifo: sum(
+                    pol_m[m][ifo] * np.exp(-1j * m * phi) for m in m_vals
+                )
+                for ifo in ifos
+            }
+            # Deferred-timeshift correction (no-op when phi == 0.0).
+            if phi != 0.0:
+                target_phase = dts["phase_ref"] + phi
+                epoch_target = wfg_utils.compute_epoch_from_resized_td_modes(
+                    dts["resized_td_modes"], dts["iota"], target_phase,
+                    dts["delta_t"],
+                )
+                dt_target = 1.0 / dts["delta_f"] + epoch_target
+                dt_correction = dt_target - dts["dt_ref"]
+                shift = np.exp(-2j * np.pi * dt_correction * freqs)
+                for ifo in ifos:
+                    mu[ifo] = mu[ifo] * shift
+
+            rho2opt = sum(
+                inner_product(mu_ifo, mu_ifo) for mu_ifo in mu.values()
+            )
+            kappa2 = sum(
+                inner_product(d_ifo, mu_ifo)
+                for d_ifo, mu_ifo in zip(d.values(), mu.values())
+            )
+            log_likelihoods[idx] = self.log_Zn + kappa2 - 0.5 * rho2opt
+        return log_likelihoods
+
     def _log_likelihood_phase_marginalized(self, theta):
         """
 

dingo/gw/result.py

@@ -701,6 +701,10 @@ def get_pesummary_samples(self, num_processes=1):
             * The spin angles phi_jl and theta_jn are transformed to account for a
             difference in phase definition.
             * Some columns are dropped: delta_log_prob_target, log_prob
+            * Aligned-spin columns ``chi_1``/``chi_2`` are renamed to
+            ``spin_1z``/``spin_2z`` so that PESummary's standard-name
+            dictionary recognises them and derives ``a_1``, ``a_2``,
+            ``chi_eff``, etc.
         """
         if hasattr(self, "_pesummary_samples"):
             return self._pesummary_samples
@@ -727,18 +731,27 @@ def get_pesummary_samples(self, num_processes=1):
             if "time" in col:
                 samples.loc[:, col] += self.t_ref
 
+        # Rename aligned-spin columns to match PESummary's standard names.
+        # PESummary's standard_names.py maps spin_1z/spin_2z but not chi_1/chi_2,
+        # so without this rename pesummary cannot derive a_1, chi_eff, chi_p, etc.
+        samples.rename(
+            columns={"chi_1": "spin_1z", "chi_2": "spin_2z"}, inplace=True
+        )
+
         spin_conversion_phase_old = self.base_metadata["dataset_settings"][
             "waveform_generator"
         ].get("spin_conversion_phase")
 
-        # Redefine phase parameter to be consistent with Bilby. COMMENTED BECAUSE SLOW
-        samples = change_spin_conversion_phase(
-            samples,
-            self.f_ref,
-            spin_conversion_phase_old,
-            None,
-            num_processes=num_processes,
-        )
+        # Redefine phase parameter to be consistent with Bilby.
+        # Only apply spin conversion when precessing spin parameters are present.
+        if "phi_jl" in samples.columns:
+            samples = change_spin_conversion_phase(
+                samples,
+                self.f_ref,
+                spin_conversion_phase_old,
+                None,
+                num_processes=num_processes,
+            )
 
         self._pesummary_samples = samples
 
@@ -750,7 +763,9 @@ def pesummary_prior(self):
 
         By convention, Dingo stores all times *relative* to a reference time, typically
         the trigger time for an event. The prior returned here corrects for that offset to
-        be consistent with other codes.
+        be consistent with other codes. Aligned-spin keys ``chi_1``/``chi_2`` are
+        renamed to ``spin_1z``/``spin_2z`` to match the renaming applied in
+        :meth:`get_pesummary_samples`.
         """
         prior = copy.deepcopy(self.prior)
         for p in prior:
@@ -760,4 +775,7 @@ def pesummary_prior(self):
                     prior[p].minimum += self.t_ref
                 except AttributeError:
                     continue
+        for old, new in (("chi_1", "spin_1z"), ("chi_2", "spin_2z")):
+            if old in prior and new not in prior:
+                prior[new] = prior.pop(old)
         return prior