reaction_kinematics

57class KinematicsResult(Mapping[str, _V]):
58    """
59    A read-only, dict-like kinematics result: ``result["theta3_lab"]`` works
60    exactly like a plain dict of arrays/lists, plus a ``.units`` dict (str ->
61    ``pint.Unit``) giving the unit of each key.
62
63    Values are plain numbers/arrays, not pint.Quantity — they work unmodified
64    with numpy, matplotlib, pandas, etc. Check ``.units[key]`` when you need to
65    know or convert the unit; wrap a value yourself (``ureg.Quantity(result[key],
66    result.units[key])``) if you want pint-native arithmetic for a specific step.
67    """
68
69    def __init__(self, data: dict[str, _V], units: dict[str, pint.Unit]) -> None:
70        self._data = data
71        self.units = units
72
73    def __getitem__(self, key: str) -> _V:
74        return self._data[key]
75
76    def __iter__(self) -> Iterator[str]:
77        return iter(self._data)
78
79    def __len__(self) -> int:
80        return len(self._data)
81
82    def __repr__(self) -> str:
83        return f"{type(self).__name__}({self._data!r}, units={self.units!r})"

122class Reaction:
123    """
124    Defines a two-body nuclear reaction: projectile + target → ejectile + recoil.
125
126    All energy-dependent methods accept a ``beam_energy`` parameter (beam kinetic
127    energy, MeV by default) and an ``energy_unit`` keyword that governs both that
128    input and every energy- and momentum-valued output; likewise ``angle_value``/
129    ``angle_unit`` govern every angle-valued input and output.
130
131    The three ``kinematics_*`` methods below return a ``KinematicsResult``: a
132    read-only, dict-like object (``result["theta3_lab"]`` works like a plain
133    dict of arrays/lists) with a ``.units`` dict (str -> ``pint.Unit``) giving
134    the unit of each key. ``Reaction.output_units(angle_unit=..., energy_unit=...)``
135    gives that same ``{key: pint.Unit}`` mapping without needing to run a
136    computation first. See the ``Returns`` section of each method below for
137    per-key detail. Internal computation is cached per energy so repeated
138    calls at the same energy are efficient.
139
140    Parameters
141    ----------
142    mass1 : str, MassInput, or float
143        Projectile mass, or a full reaction string in ``"target(beam,ejectile)recoil"``
144        notation (e.g. ``"3H(p,n)3He"``). If notation is given, ``mass2``–``mass4``
145        must be omitted.
146    mass2, mass3, mass4 : str, MassInput, or float, optional
147        Target, ejectile, and recoil masses. Required when ``mass1`` is not a
148        reaction notation string. Strings like ``"p"``, ``"12C"``, ``"alpha"``
149        are looked up in the mass table (no ``mass_unit`` needed — the table is
150        already in the right units). Floats require ``mass_unit``.
151    mass_unit : str or EnergyUnit, optional
152        Unit for numeric masses (e.g. ``"MeV"``, ``"keV"``). Only used for
153        ``float``/``int`` mass arguments; ignored for isotope-string arguments.
154
155    Attributes
156    ----------
157    n_cm_grid_points : int
158        Total number of CM angle grid points (default 1001). Changing this
159        invalidates the cache.
160
161    Examples
162    --------
163    >>> rxn = Reaction("p", "3H", "n", "3He")
164    >>> rxn = Reaction("3H(p,n)3He")  # equivalent
165    >>> rxn.q_value
166    -0.763...
167    >>> data = rxn.kinematics_table_at_beam_energy(1.2)
168    >>> result = rxn.kinematics_at_beam_energy_and_angle(1.2, "theta3_lab", 30)
169    >>> branches = rxn.kinematics_curve_at_angle(np.linspace(1.0, 5.0, 200), 30)
170    """
171
172    def __init__(
173        self,
174        mass1: MassArg,
175        mass2: MassArg | None = None,
176        mass3: MassArg | None = None,
177        mass4: MassArg | None = None,
178        *,
179        mass_unit: str | EnergyUnit | None = None,
180    ) -> None:
181        if isinstance(mass1, str) and "(" in mass1:
182            if any(m is not None for m in (mass2, mass3, mass4)):
183                raise ValueError(
184                    "Cannot mix reaction notation string with separate mass arguments."
185                )
186            _m1, _m2, _m3, _m4 = _parse_reaction_notation(mass1)
187        else:
188            if mass2 is None or mass3 is None or mass4 is None:
189                raise ValueError("Provide either a reaction notation string or all four masses.")
190            _m1, _m2, _m3, _m4 = mass1, mass2, mass3, mass4
191        self._m1 = _parse_mass(_m1, mass_unit)
192        self._m2 = _parse_mass(_m2, mass_unit)
193        self._m3 = _parse_mass(_m3, mass_unit)
194        self._m4 = _parse_mass(_m4, mass_unit)
195        self.__n_cm_grid_points: int = 1001
196        self._cached_ek: float | None = None
197        self._table: dict[str, list[float]] | None = None
198        # per-energy kinematic state, populated by _bind
199        self._nogo: bool = False
200        self._pcmp: float | None = None
201        self._thesinh: float | None = None
202        self._thecosh: float | None = None
203        self._e03: float | None = None
204        self._e04: float | None = None
205        # lab-frame energy extrema
206        self._emax3: float | None = None
207        self._emin3: float | None = None
208        self._emax4: float | None = None
209        self._emin4: float | None = None
210        # max lab angles and associated quantities (None = no forward maximum)
211        self._theta3max: float | None = None
212        self._theta4max: float | None = None
213        self._e3atmaxang: float | None = None
214        self._e4atmaxang: float | None = None
215        self._cmcos3max: float | None = None
216        self._cmcos4max: float | None = None
217
218    @property
219    def n_cm_grid_points(self) -> int:
220        return self.__n_cm_grid_points
221
222    @n_cm_grid_points.setter
223    def n_cm_grid_points(self, val: int) -> None:
224        if not isinstance(val, int):
225            raise TypeError(f"n_cm_grid_points must be an int, got {type(val).__name__}")
226        if val < 2:
227            raise ValueError(f"n_cm_grid_points={val} must be at least 2")
228        self.__n_cm_grid_points = val
229        self._cached_ek = None
230        self._table = None
231
232    @property
233    def q_value(self) -> float:
234        """Q-value of the reaction in MeV: Q = (beam + target) - (ejectile + recoil)."""
235        return self._m1 + self._m2 - self._m3 - self._m4
236
237    def _bind(self, ek_mev: float) -> None:
238        """Compute and cache kinematic state for the given beam energy."""
239        if ek_mev == self._cached_ek:
240            return
241        self._cached_ek = ek_mev
242        self._table = None
243        self._compute(ek_mev)
244
245    def _compute(self, ek_mev: float) -> None:
246        self._nogo = False
247        self._pcmp = self._thesinh = self._thecosh = self._e03 = self._e04 = None
248
249        # Mandelstam s
250        s = (self._m1 + self._m2) ** 2 + 2.0 * self._m2 * ek_mev
251        if s <= 0.0:
252            self._nogo = True
253            return
254
255        # initial CM momentum
256        pcm2 = (s - self._m1**2 - self._m2**2) ** 2 - 4.0 * self._m1**2 * self._m2**2
257        if pcm2 < 0:
258            self._nogo = True
259            return
260        pcm = math.sqrt(pcm2 / (4.0 * s))
261
262        # CM rapidity → boost parameters
263        acmratio = (math.sqrt(self._m2**2 + pcm**2) + pcm) / self._m2
264        cmrap = math.log(acmratio)
265        self._thesinh = math.sinh(cmrap)
266        self._thecosh = math.cosh(cmrap)
267
268        # final-state CM momentum
269        pcmp2 = (s - self._m3**2 - self._m4**2) ** 2 - 4.0 * self._m3**2 * self._m4**2
270        if pcmp2 < 0:
271            self._nogo = True
272            return
273        self._pcmp = math.sqrt(pcmp2 / (4.0 * s))
274
275        # CM total energies of outgoing particles
276        self._e03 = math.sqrt(self._pcmp**2 + self._m3**2)
277        self._e04 = math.sqrt(self._pcmp**2 + self._m4**2)
278
279        # lab-frame energy extrema
280        self._emax3 = self._e03 * self._thecosh + self._pcmp * self._thesinh - self._m3
281        self._emin3 = self._e03 * self._thecosh - self._pcmp * self._thesinh - self._m3
282        self._emax4 = self._e04 * self._thecosh + self._pcmp * self._thesinh - self._m4
283        self._emin4 = self._e04 * self._thecosh - self._pcmp * self._thesinh - self._m4
284
285        # reset max-angle quantities to sentinel values
286        self._theta3max = None
287        self._theta4max = None
288        self._e3atmaxang = None
289        self._e4atmaxang = None
290        self._cmcos3max = None
291        self._cmcos4max = None
292
293        # max ejectile lab angle (only exists when pcmp < m3 * sinh(rapidity))
294        thetatest3: float | None = None
295        if self._m3 > 0.0:
296            thetatest3 = self._pcmp / (self._m3 * self._thesinh)
297            if thetatest3 < 1.0:
298                self._theta3max = math.asin(thetatest3)
299                patmax = (self._e03 * math.cos(self._theta3max) * self._thesinh) / (
300                    1.0 + thetatest3**2 * self._thesinh**2
301                )
302                eatmax = math.sqrt(patmax**2 + self._m3**2)
303                self._e3atmaxang = eatmax - self._m3
304                self._cmcos3max = (eatmax - self._e03 * self._thecosh) / (
305                    self._pcmp * self._thesinh
306                )
307
308        # elastic case: forward-angle symmetry forces theta3max = 90°
309        if (self._m1 + self._m2) == (self._m3 + self._m4) and thetatest3 is not None:
310            if abs(thetatest3 - 1.0) < 1e-3:
311                self._theta3max = math.pi / 2.0
312                self._cmcos3max = -1.0
313                self._e3atmaxang = (
314                    self._e03 * self._thecosh
315                    + self._cmcos3max * self._pcmp * self._thesinh
316                    - self._m3
317                )
318
319        # max recoil lab angle
320        thetatest4: float | None = None
321        if self._m4 > 0.0:
322            thetatest4 = self._pcmp / (self._m4 * self._thesinh)
323            if thetatest4 < 1.0:
324                self._theta4max = math.asin(thetatest4)
325                patmax = (self._e04 * math.cos(self._theta4max) * self._thesinh) / (
326                    1.0 + thetatest4**2 * self._thesinh**2
327                )
328                eatmax = math.sqrt(patmax**2 + self._m4**2)
329                self._e4atmaxang = eatmax - self._m4
330                self._cmcos4max = (eatmax - self._e04 * self._thecosh) / (
331                    self._pcmp * self._thesinh
332                )
333
334        # elastic case: theta4max = 90°
335        if (self._m1 + self._m2) == (self._m3 + self._m4) and thetatest4 is not None:
336            if abs(thetatest4 - 1.0) < 1e-3:
337                self._theta4max = math.pi / 2.0
338                self._cmcos4max = 1.0
339                self._e4atmaxang = (
340                    self._e04 * self._thecosh
341                    - self._cmcos4max * self._pcmp * self._thesinh
342                    - self._m4
343                )
344
345    def _kinematics_at_coscm(self, coscm: float) -> dict[str, float]:
346        if (
347            self._pcmp is None
348            or self._thecosh is None
349            or self._e03 is None
350            or self._thesinh is None
351            or self._e04 is None
352        ):
353            raise ValueError("Kinematic state not computed — call _bind first")
354
355        sincm = math.sqrt(max(0.0, 1.0 - coscm**2))
356
357        ppar3 = self._pcmp * self._thecosh * coscm + self._e03 * self._thesinh
358        pperp3 = self._pcmp * sincm
359        ptot3 = math.hypot(ppar3, pperp3)
360
361        ppar4 = -self._pcmp * self._thecosh * coscm + self._e04 * self._thesinh
362        pperp4 = self._pcmp * sincm
363        ptot4 = math.hypot(ppar4, pperp4)
364
365        e3 = self._e03 * self._thecosh + coscm * self._pcmp * self._thesinh - self._m3
366        e4 = self._e04 * self._thecosh - coscm * self._pcmp * self._thesinh - self._m4
367
368        # dOmega_lab/dOmega_cm for each ejectile, i.e. the factor that converts a
369        # lab-frame differential cross section to the cm-frame one:
370        # dsigma/dOmega_cm = dsigma/dOmega_lab * jacobianN_lab
371        #
372        # Both are d(cos_theta_lab)/d(coscm), coscm being cos_theta_cm of particle 3.
373        # ppar4's sign convention (see above) flips the sign of the thecosh term
374        # relative to particle 3's formula.
375        if ptot3 > 0:
376            jacobian3_lab = (self._pcmp**2 / ptot3**2) * (
377                coscm * ppar3 / ptot3 + sincm * pperp3 * self._thecosh / ptot3
378            )
379        else:
380            jacobian3_lab = math.nan
381        if ptot4 > 0:
382            jacobian4_lab = (self._pcmp**2 / ptot4**2) * (
383                coscm * ppar4 / ptot4 - sincm * pperp4 * self._thecosh / ptot4
384            )
385        else:
386            jacobian4_lab = math.nan
387
388        return {
389            "cos_theta_cm": coscm,
390            "theta_cm": math.acos(coscm),
391            "theta3_lab": math.acos(ppar3 / ptot3) if ptot3 > 0 else 0.0,
392            "theta4_lab": math.acos(ppar4 / ptot4) if ptot4 > 0 else 0.0,
393            "energy3_lab": e3,
394            "energy4_lab": e4,
395            "velocity3_lab": ptot3 / (e3 + self._m3),
396            "velocity4_lab": ptot4 / (e4 + self._m4),
397            "momentum3_lab": ptot3,
398            "momentum4_lab": ptot4,
399            "jacobian3_lab": jacobian3_lab,
400            "jacobian4_lab": jacobian4_lab,
401        }
402
403    def _build_table(self) -> None:
404        """Build the interpolation table over the full CM angle grid."""
405        keys = [
406            "cos_theta_cm",
407            "theta_cm",
408            "theta3_lab",
409            "theta4_lab",
410            "energy3_lab",
411            "energy4_lab",
412            "velocity3_lab",
413            "velocity4_lab",
414            "momentum3_lab",
415            "momentum4_lab",
416            "jacobian3_lab",
417            "jacobian4_lab",
418        ]
419        table: dict[str, list[float]] = {k: [] for k in keys}
420        for coscm in np.linspace(-1.0, 1.0, self.n_cm_grid_points):
421            row = self._kinematics_at_coscm(coscm)
422            for k in keys:
423                table[k].append(row[k])
424        self._table = table
425
426    def kinematics_table_at_beam_energy(
427        self,
428        beam_energy: float,
429        *,
430        angle_unit: AngleUnit = AngleUnit.deg,
431        energy_unit: EnergyUnit = EnergyUnit.MeV,
432    ) -> KinematicsResult[npt.NDArray[np.float64]]:
433        """
434        Compute full kinematics over a CM angle grid.
435
436        Values are plain arrays (work unmodified with numpy/matplotlib/etc.).
437        The returned ``KinematicsResult`` also carries a ``.units`` dict (str ->
438        ``pint.Unit``) giving the unit of each key.
439
440        Parameters
441        ----------
442        beam_energy : float
443            Beam kinetic energy, in ``energy_unit``.
444        angle_unit : AngleUnit, optional
445            Unit for angle outputs ``theta_cm``, ``theta3_lab``, ``theta4_lab``
446            (default degrees). Does not affect ``jacobianN_lab``, which is a
447            dimensionless ratio of solid angles regardless of angle unit.
448        energy_unit : EnergyUnit, optional
449            Unit of ``beam_energy`` (default MeV). Also governs the unit of every
450            energy- and momentum-valued output (``energy3_lab``, ``energy4_lab``,
451            ``momentum3_lab``, ``momentum4_lab``), exactly as ``angle_unit`` governs
452            every angle-valued output.
453
454        Returns
455        -------
456        KinematicsResult[np.ndarray]
457            ``"cos_theta_cm"`` : dimensionless, in [-1, 1].
458            ``"theta_cm"``, ``"theta3_lab"``, ``"theta4_lab"`` : ``angle_unit``
459            (default degrees).
460            ``"energy3_lab"``, ``"energy4_lab"`` : ``energy_unit`` (default MeV).
461            ``"velocity3_lab"``, ``"velocity4_lab"`` : dimensionless, as a fraction
462            of the speed of light ``c``.
463            ``"momentum3_lab"``, ``"momentum4_lab"`` : ``energy_unit``/``c``
464            (natural units, default MeV/``c``).
465            ``"jacobian3_lab"``, ``"jacobian4_lab"`` : dimensionless.
466            ``jacobianN_lab`` is ``dOmegaN_lab/dOmega_cm``, i.e. the factor that
467            converts a lab-frame differential cross section to the cm-frame one:
468            ``dsigma/dOmega_cm = dsigma/dOmegaN_lab * jacobianN_lab``.
469
470            ``result.units`` carries this same mapping as ``{key: pint.Unit}``.
471            ``Reaction.output_units(angle_unit=..., energy_unit=...)`` gives the
472            same thing without needing to run a computation first.
473
474        Raises
475        ------
476        ValueError
477            If the reaction is kinematically forbidden at this energy.
478        """
479        angle_unit = AngleUnit.from_any(angle_unit)
480        energy_unit = EnergyUnit.from_any(energy_unit)
481        ek_mev = _parse_energy(beam_energy, energy_unit)
482        self._bind(ek_mev)
483        if self._nogo:
484            raise ValueError(f"Reaction kinematically forbidden at beam_energy={ek_mev} MeV")
485        keys = [
486            "cos_theta_cm",
487            "theta_cm",
488            "theta3_lab",
489            "theta4_lab",
490            "energy3_lab",
491            "energy4_lab",
492            "velocity3_lab",
493            "velocity4_lab",
494            "momentum3_lab",
495            "momentum4_lab",
496            "jacobian3_lab",
497            "jacobian4_lab",
498        ]
499        rows = [
500            self._kinematics_at_coscm(coscm)
501            for coscm in np.linspace(-1.0, 1.0, self.n_cm_grid_points)
502        ]
503        result = {k: np.array([row[k] for row in rows]) for k in keys}
504        for k in keys:
505            if k in _ANGLE_KEYS:
506                result[k] = result[k] / angle_unit.value
507            elif k in _ENERGY_KEYS:
508                result[k] = result[k] / energy_unit.value
509        return KinematicsResult(result, _result_units(keys, angle_unit, energy_unit))
510
511    def kinematics_at_beam_energy_and_angle(
512        self,
513        beam_energy: float,
514        angle_name: str,
515        angle_value: float,
516        *,
517        angle_unit: AngleUnit = AngleUnit.deg,
518        energy_unit: EnergyUnit = EnergyUnit.MeV,
519        duplicate_tol: float = 1e-6,
520    ) -> KinematicsResult[list[float]]:
521        """
522        Interpolate kinematic quantities at a fixed beam energy and kinematic variable value.
523
524        Always returns lists to handle multi-valued cases (e.g. two ejectile
525        energies at the same lab angle). Values are plain lists (work unmodified
526        with numpy/pandas/etc.). The returned ``KinematicsResult`` also carries a
527        ``.units`` dict (str -> ``pint.Unit``) giving the unit of each key.
528
529        Parameters
530        ----------
531        beam_energy : float
532            Beam kinetic energy, in ``energy_unit``.
533        angle_name : str
534            Independent variable name, can be one of ``"theta3_lab"``, ``"theta4_lab"``,
535            ``"theta_cm"``, ``"cos_theta_cm"``.
536        angle_value : float
537            Value to evaluate at. For angle keys (``theta*``), interpreted in
538            ``angle_unit`` (default degrees). For ``"cos_theta_cm"``, treated
539            as a dimensionless cosine — ``angle_unit`` is ignored.
540        angle_unit : AngleUnit, optional
541            Unit of ``angle_value`` for angle keys (default degrees).
542        energy_unit : EnergyUnit, optional
543            Unit of ``beam_energy`` (default MeV). Also governs the unit of every
544            energy- and momentum-valued output and of ``duplicate_tol``, exactly
545            as ``angle_unit`` governs every angle-valued output.
546        duplicate_tol : float, optional
547            Tolerance for merging near-duplicate solutions (default 1e-6), in
548            ``energy_unit``. Solutions within this ``energy3_lab`` difference are
549            treated as the same physical solution.
550
551        Returns
552        -------
553        KinematicsResult[list[float]]
554            Full dict of all kinematic variables, each a list of solutions sorted
555            descending by ``energy3_lab``. Angle keys are in ``angle_unit``, energy
556            and momentum keys in ``energy_unit`` (both default degrees/MeV) — see
557            ``Reaction.output_units()`` for the full mapping, or just read
558            ``.units[key]`` off the returned result directly.
559
560        Raises
561        ------
562        ValueError
563            If ``beam_energy`` or ``angle_value`` is not finite, if the reaction is
564            kinematically forbidden at this energy, or if ``angle_value`` is outside
565            the physical range.
566
567        Examples
568        --------
569        >>> rxn = Reaction("p", "3H", "n", "3He")
570        >>> rxn.kinematics_at_beam_energy_and_angle(1.2, "theta3_lab", 30)
571        {'theta3_lab': [...], 'energy3_lab': [...], ...}
572        """
573        angle_unit = AngleUnit.from_any(angle_unit)
574        energy_unit = EnergyUnit.from_any(energy_unit)
575        output = self._kinematics_at_beam_energy_and_angle_raw(
576            beam_energy,
577            angle_name,
578            angle_value,
579            angle_unit=angle_unit,
580            energy_unit=energy_unit,
581            duplicate_tol=duplicate_tol,
582        )
583        return KinematicsResult(output, _result_units(output.keys(), angle_unit, energy_unit))
584
585    def _kinematics_at_beam_energy_and_angle_raw(
586        self,
587        beam_energy: float,
588        angle_name: str,
589        angle_value: float,
590        *,
591        angle_unit: AngleUnit,
592        energy_unit: EnergyUnit,
593        duplicate_tol: float,
594    ) -> dict[str, list[float]]:
595        """Plain-float core of kinematics_at_beam_energy_and_angle, reused internally
596        by kinematics_curve_at_angle."""
597        # Keep the value/unit as the caller wrote it for error messages — the
598        # internal, canonical-radians angle_value below is not what they typed.
599        is_angle = angle_name.startswith("theta")
600        angle_value_display = f"{angle_value} {angle_unit.name}" if is_angle else f"{angle_value}"
601        if is_angle:
602            angle_value = angle_value * angle_unit.value
603
604        if not math.isfinite(angle_value):
605            raise ValueError(f"{angle_name}={angle_value_display} is not a finite number")
606
607        duplicate_tol_mev = duplicate_tol * energy_unit.value
608
609        ek_mev = _parse_energy(beam_energy, energy_unit)
610        self._bind(ek_mev)
611        if self._nogo:
612            raise ValueError(f"Reaction kinematically forbidden at beam_energy={ek_mev} MeV")
613
614        if self._table is None:
615            self._build_table()
616        assert self._table is not None
617
618        keys = list(self._table.keys())
619        xs = self._table[angle_name]
620
621        solutions = []
622
623        exact_idx = np.where(np.isclose(xs, angle_value, atol=1e-12))[0]
624        if len(exact_idx) > 0:
625            for i in exact_idx:
626                solutions.append({k: self._table[k][i] for k in keys})
627        else:
628            found = False
629            for i in range(len(xs) - 1):
630                x0, x1 = xs[i], xs[i + 1]
631                if (x0 - angle_value) * (x1 - angle_value) <= 0 and x0 != x1:
632                    found = True
633                    t = (angle_value - x0) / (x1 - x0)
634                    solutions.append(
635                        {
636                            k: self._table[k][i] + t * (self._table[k][i + 1] - self._table[k][i])
637                            for k in keys
638                        }
639                    )
640            if not found:
641                raise ValueError(f"{angle_name}={angle_value_display} outside physical range")
642
643        unique: list = []
644        for sol in solutions:
645            if not any(
646                abs(sol["energy3_lab"] - u["energy3_lab"]) < duplicate_tol_mev for u in unique
647            ):
648                unique.append(sol)
649        unique.sort(key=lambda s: s["energy3_lab"], reverse=True)
650
651        output = {k: [s[k] for s in unique] for k in keys}
652        for k in keys:
653            if k in _ANGLE_KEYS:
654                output[k] = [v / angle_unit.value for v in output[k]]
655            elif k in _ENERGY_KEYS:
656                output[k] = [v / energy_unit.value for v in output[k]]
657        return output
658
659    def kinematics_curve_at_angle(
660        self,
661        beam_energy_array: Iterable[float],
662        theta3_lab: float,
663        *,
664        angle_unit: AngleUnit = AngleUnit.deg,
665        energy_unit: EnergyUnit = EnergyUnit.MeV,
666    ) -> list[KinematicsResult[npt.NDArray[np.float64]]]:
667        """
668        Compute ejectile kinematics at a fixed ejectile lab angle (``theta3_lab``)
669        over a range of beam energies.
670
671        Returns two branches (high- and low-energy) as a list of two
672        ``KinematicsResult``s. Each contains plain arrays indexed by beam energy,
673        with ``NaN`` where that branch does not exist, plus a ``.units`` dict
674        (str -> ``pint.Unit``) computed from the same ``angle_unit``/``energy_unit``
675        as the data. Branch 0 is always the higher-energy solution.
676
677        Parameters
678        ----------
679        beam_energy_array : array-like
680            Beam energies to sweep, in ``energy_unit``.
681        theta3_lab : float
682            Fixed ejectile lab angle (``theta3_lab`` in the output dict), in
683            ``angle_unit``.
684        angle_unit : AngleUnit, optional
685            Unit of ``theta3_lab`` input and ``theta4_lab`` output (default degrees).
686        energy_unit : EnergyUnit, optional
687            Unit of ``beam_energy_array`` values (default MeV). Also governs the
688            unit of every energy- and momentum-valued output (``beam_energy_lab``,
689            ``energy3_lab``, ``energy4_lab``, ``momentum3_lab``, ``momentum4_lab``),
690            exactly as ``angle_unit`` governs ``theta4_lab``.
691
692        Returns
693        -------
694        list[KinematicsResult[np.ndarray]]
695            List of two results, each with keys ``"beam_energy_lab"``, ``"energy3_lab"``,
696            ``"energy4_lab"``, ``"theta4_lab"``, ``"velocity3_lab"``, ``"velocity4_lab"``,
697            ``"momentum3_lab"``, ``"momentum4_lab"``, ``"jacobian3_lab"``, ``"jacobian4_lab"``.
698            Angle keys are in ``angle_unit``, energy and momentum keys in
699            ``energy_unit`` (both default degrees/MeV) — see ``Reaction.output_units()``
700            for the full mapping, or just read ``.units[key]`` off either result.
701
702        Examples
703        --------
704        >>> rxn = Reaction("p", "3H", "n", "3He")
705        >>> branches = rxn.kinematics_curve_at_angle(np.linspace(1.0, 5.0, 200), 30)
706        >>> for b in branches:
707        ...     plt.plot(b["beam_energy_lab"], b["energy3_lab"])
708        """
709        angle_unit = AngleUnit.from_any(angle_unit)
710        energy_unit = EnergyUnit.from_any(energy_unit)
711        theta_rad = theta3_lab * angle_unit.value
712
713        keys = [
714            "energy3_lab",
715            "energy4_lab",
716            "theta4_lab",
717            "velocity3_lab",
718            "velocity4_lab",
719            "momentum3_lab",
720            "momentum4_lab",
721            "jacobian3_lab",
722            "jacobian4_lab",
723        ]
724        branches = [
725            {"beam_energy_lab": [], **{k: [] for k in keys}},
726            {"beam_energy_lab": [], **{k: [] for k in keys}},
727        ]
728
729        for ek in beam_energy_array:
730            ek_mev = _parse_energy(ek, energy_unit)
731            try:
732                row = self._kinematics_at_beam_energy_and_angle_raw(
733                    ek_mev,
734                    "theta3_lab",
735                    theta_rad,
736                    angle_unit=AngleUnit.rad,
737                    energy_unit=EnergyUnit.MeV,
738                    duplicate_tol=1e-6,
739                )
740            except ValueError:
741                solutions = []
742            else:
743                n = len(row["energy3_lab"])
744                solutions = [{k: row[k][i] for k in keys} for i in range(n)]
745
746            for i, branch in enumerate(branches):
747                branch["beam_energy_lab"].append(ek_mev)
748                sol = solutions[i] if i < len(solutions) else None
749                for k in keys:
750                    branch[k].append(sol[k] if sol is not None else float("nan"))
751
752        result = []
753        for branch in branches:
754            converted = {}
755            for k, v in branch.items():
756                arr = np.array(v)
757                if k in _ANGLE_KEYS:
758                    arr = arr / angle_unit.value
759                elif k in _ENERGY_KEYS:
760                    arr = arr / energy_unit.value
761                converted[k] = arr
762            result.append(
763                KinematicsResult(
764                    converted, _result_units(converted.keys(), angle_unit, energy_unit)
765                )
766            )
767        return result
768
769    @staticmethod
770    def output_units(
771        *,
772        angle_unit: AngleUnit | str = AngleUnit.deg,
773        energy_unit: EnergyUnit | str = EnergyUnit.MeV,
774    ) -> dict[str, pint.Unit]:
775        """
776        The ``pint.Unit`` for every key that can appear in the dicts returned by
777        ``kinematics_table_at_beam_energy``, ``kinematics_at_beam_energy_and_angle``,
778        and ``kinematics_curve_at_angle`` — the same mapping each of those methods
779        attaches as ``result.units``, available here without running a computation
780        first.
781
782        Pass the same ``angle_unit``/``energy_unit`` you'd pass to one of those
783        methods to get a matching ``{key: pint.Unit}`` map — not every key appears
784        in every method's output, so look up only the keys actually present in the
785        result you're inspecting. ``velocity3_lab``/``velocity4_lab`` (fraction of
786        the speed of light) and ``jacobian3_lab``/``jacobian4_lab`` (ratio of solid
787        angles, dOmega_lab/dOmega_cm) are both genuinely dimensionless, not just
788        unlabeled.
789
790        Examples
791        --------
792        >>> str(Reaction.output_units()["theta3_lab"])
793        'degree'
794        >>> f"{Reaction.output_units(energy_unit='keV')['momentum3_lab']:~}"
795        'keV / c'
796        """
797        angle_unit = AngleUnit.from_any(angle_unit)
798        energy_unit = EnergyUnit.from_any(energy_unit)
799        # Built from the same _result_unit used to actually tag KinematicsResult.units,
800        # so this preview can't drift out of sync with real results.
801        return _result_units(_ALL_RESULT_KEYS, angle_unit, energy_unit)