laddu_core/kinematics/
frame_axes.rs

1use serde::{Deserialize, Serialize};
2
3use super::{
4    decay_angles::DecayAngles,
5    rest_frame::RestFrame,
6    support::{checked_boost_vector, unit_vector},
7};
8use crate::{
9    quantum::Frame,
10    vectors::{Vec3, Vec4},
11    LadduError, LadduResult,
12};
13
14/// Orthonormal axes used to project decay momenta into an angular-analysis frame.
15#[derive(Clone, Copy, Debug, PartialEq, Serialize, Deserialize)]
16pub struct FrameAxes {
17    x: Vec3,
18    y: Vec3,
19    z: Vec3,
20}
21
22impl Default for FrameAxes {
23    fn default() -> Self {
24        Self {
25            x: Vec3::x(),
26            y: Vec3::y(),
27            z: Vec3::z(),
28        }
29    }
30}
31
32impl FrameAxes {
33    /// Construct normalized right-handed axes after validating orthogonality.
34    pub fn new(x: Vec3, y: Vec3, z: Vec3) -> LadduResult<Self> {
35        const ORTHOGONALITY_TOL: f64 = 1.0e-12;
36        const HANDEDNESS_TOL: f64 = 1.0e-12;
37
38        let x = unit_vector(x, "x axis")?;
39        let y = unit_vector(y, "y axis")?;
40        let z = unit_vector(z, "z axis")?;
41        if x.dot(&y).abs() > ORTHOGONALITY_TOL
42            || x.dot(&z).abs() > ORTHOGONALITY_TOL
43            || y.dot(&z).abs() > ORTHOGONALITY_TOL
44        {
45            return Err(LadduError::Custom(
46                "frame axes must be mutually orthogonal".to_string(),
47            ));
48        }
49        if x.cross(&y).dot(&z) < -HANDEDNESS_TOL {
50            return Err(LadduError::Custom(
51                "frame axes must form a right-handed basis".to_string(),
52            ));
53        }
54        Ok(Self { x, y, z })
55    }
56
57    /// Construct right-handed axes from a `z` axis and a plane normal.
58    pub fn from_z_and_plane_normal(z: Vec3, plane_normal: Vec3) -> LadduResult<Self> {
59        let z = unit_vector(z, "z axis")?;
60        let plane_normal = unit_vector(plane_normal, "frame-plane normal")?;
61        let x = unit_vector(plane_normal.cross(&z), "x axis")?;
62        let y = unit_vector(z.cross(&x), "y axis")?;
63        Self::new(x, y, z)
64    }
65
66    /// Construct decay-frame axes for a reaction root from caller-selected production momenta.
67    ///
68    /// `reference`, `parent`, and `spectator` are lab-frame four-momenta. `system_boost` is the
69    /// boost into the frame where the production plane is defined. This keeps named event topology
70    /// handling outside the frame helper while still sharing the convention-sensitive geometry.
71    pub fn from_decay_frame(
72        frame: Frame,
73        reference: Vec4,
74        parent: Vec4,
75        spectator: Vec4,
76        system_boost: Vec3,
77    ) -> LadduResult<Self> {
78        checked_boost_vector(system_boost, "production system")?;
79        let reference = reference.boost(&system_boost);
80        let parent = parent.boost(&system_boost);
81        let spectator = spectator.boost(&system_boost);
82
83        let parent_rest = RestFrame::new(parent)?;
84        let reference_in_parent = parent_rest.transform(reference).vec3();
85        let spectator_in_parent = parent_rest.transform(spectator).vec3();
86
87        let plane_normal = unit_vector(
88            reference_in_parent.cross(&(-spectator_in_parent)),
89            "production-plane normal",
90        )?;
91
92        let z = match frame {
93            Frame::Helicity => unit_vector(-spectator_in_parent, "decay-frame z axis")?,
94            Frame::GottfriedJackson => {
95                unit_vector(reference_in_parent, "Gottfried-Jackson z axis")?
96            }
97            Frame::Adair => {
98                return Err(LadduError::Custom(
99                    "Adair frame construction is not implemented yet".to_string(),
100                ));
101            }
102        };
103
104        Self::from_z_and_plane_normal(z, plane_normal)
105    }
106
107    /// Construct daughter axes in the current rest frame.
108    ///
109    /// The daughter direction defines the new `z` axis. The current `z` axis and daughter
110    /// direction define the rotation plane, which makes this helper reusable with any parent-axis
111    /// convention.
112    pub fn for_daughter(self, daughter_momentum: Vec3) -> LadduResult<Self> {
113        let z = unit_vector(daughter_momentum, "daughter z axis")?;
114        let mut plane_normal = self.z.cross(&z);
115        if plane_normal.mag2() <= f64::EPSILON * f64::EPSILON {
116            plane_normal = self.y;
117        }
118        Self::from_z_and_plane_normal(z, plane_normal)
119    }
120
121    /// Return the unit `x` axis.
122    pub const fn x(self) -> Vec3 {
123        self.x
124    }
125
126    /// Return the unit `y` axis.
127    pub const fn y(self) -> Vec3 {
128        self.y
129    }
130
131    /// Return the unit `z` axis.
132    pub const fn z(self) -> Vec3 {
133        self.z
134    }
135
136    /// Project a vector onto these frame axes.
137    pub fn components(self, vector: Vec3) -> Vec3 {
138        Vec3::new(
139            vector.dot(&self.x),
140            vector.dot(&self.y),
141            vector.dot(&self.z),
142        )
143    }
144
145    /// Compute spherical decay angles for a vector expressed in the same rest frame as the axes.
146    pub fn angles(self, vector: Vec3) -> LadduResult<DecayAngles> {
147        let components = self.components(vector);
148        DecayAngles::from_components(components)
149    }
150}
151
152#[cfg(test)]
153mod tests {
154    use approx::assert_relative_eq;
155
156    use super::*;
157
158    fn transverse_momenta() -> (Vec4, Vec4, Vec4) {
159        (
160            Vec4::new(0.0, 0.0, 5.0, 5.0),
161            Vec4::new(1.0, 0.0, 0.0, 2.0),
162            Vec4::new(-1.0, 0.0, 0.0, 2.0),
163        )
164    }
165
166    fn transverse_axes(frame: Frame) -> FrameAxes {
167        let (reference, parent, spectator) = transverse_momenta();
168        FrameAxes::from_decay_frame(frame, reference, parent, spectator, Vec3::zero()).unwrap()
169    }
170
171    fn assert_vec3_close(actual: Vec3, expected: Vec3) {
172        assert_relative_eq!(actual.x, expected.x);
173        assert_relative_eq!(actual.y, expected.y);
174        assert_relative_eq!(actual.z, expected.z);
175    }
176
177    #[test]
178    fn helicity_frame_axes_are_fixed_for_transverse_production() {
179        let axes = transverse_axes(Frame::Helicity);
180
181        assert_vec3_close(axes.x(), Vec3::new(0.0, 0.0, -1.0));
182        assert_vec3_close(axes.y(), Vec3::new(0.0, 1.0, 0.0));
183        assert_vec3_close(axes.z(), Vec3::new(1.0, 0.0, 0.0));
184    }
185
186    #[test]
187    fn gottfried_jackson_frame_uses_boosted_reference_axis() {
188        let axes = transverse_axes(Frame::GottfriedJackson);
189        let sqrt_three_over_two = 3.0_f64.sqrt() / 2.0;
190
191        assert_vec3_close(axes.x(), Vec3::new(sqrt_three_over_two, 0.0, 0.5));
192        assert_vec3_close(axes.y(), Vec3::new(0.0, 1.0, 0.0));
193        assert_vec3_close(axes.z(), Vec3::new(-0.5, 0.0, sqrt_three_over_two));
194    }
195
196    #[test]
197    fn rest_frame_boosts_parent_to_zero_momentum() {
198        let parent = Vec4::new(1.0, 2.0, 3.0, 5.0);
199        let rest_parent = RestFrame::new(parent).unwrap().transform(parent);
200
201        assert_relative_eq!(rest_parent.px(), 0.0);
202        assert_relative_eq!(rest_parent.py(), 0.0);
203        assert_relative_eq!(rest_parent.pz(), 0.0);
204    }
205
206    #[test]
207    fn daughter_axes_support_cascade_frames() {
208        let child_axes = FrameAxes::default()
209            .for_daughter(Vec3::new(1.0, 0.0, 0.0))
210            .unwrap();
211
212        assert_vec3_close(child_axes.x(), Vec3::new(0.0, 0.0, -1.0));
213        assert_vec3_close(child_axes.y(), Vec3::new(0.0, 1.0, 0.0));
214        assert_vec3_close(child_axes.z(), Vec3::new(1.0, 0.0, 0.0));
215    }
216
217    #[test]
218    fn projected_decay_angles_pin_azimuth_sign() {
219        let axes = transverse_axes(Frame::Helicity);
220        let angles = axes.angles(Vec3::new(0.0, 0.0, 1.0)).unwrap();
221
222        assert_relative_eq!(angles.costheta(), 0.0);
223        assert_relative_eq!(angles.phi(), std::f64::consts::PI);
224    }
225
226    #[test]
227    fn frame_axes_reject_degenerate_production_plane() {
228        let err = FrameAxes::from_decay_frame(
229            Frame::Helicity,
230            Vec4::new(0.0, 0.0, 5.0, 5.0),
231            Vec4::new(0.0, 0.0, 1.0, 2.0),
232            Vec4::new(0.0, 0.0, -1.0, 2.0),
233            Vec3::zero(),
234        );
235
236        assert!(err.is_err());
237    }
238}
laddu_core/kinematics/frame_axes.rs

laddu_core/kinematics/
frame_axes.rs
