rustitude_gluex/
utils.rs

1use std::{fmt::Display, num::ParseIntError, str::FromStr};
2
3use factorial::Factorial;
4use rustitude_core::{convert, prelude::*};
5use sphrs::Coordinates;
6use thiserror::Error;
7
8pub fn breakup_momentum<F: Field>(m0: F, m1: F, m2: F) -> F {
9    F::sqrt(F::abs(
10        m0.powi(4) + m1.powi(4) + m2.powi(4)
11            - convert!(2, F)
12                * (m0.powi(2) * m1.powi(2) + m0.powi(2) * m2.powi(2) + m1.powi(2) * m2.powi(2)),
13    )) / (convert!(2, F) * m0)
14}
15
16/// Computes the ([`Complex<F>`]) breakup momentum of a particle with mass `m0` decaying into two particles
17/// with masses `m1` and `m2`.
18pub fn breakup_momentum_c<F: Field>(m0: F, m1: F, m2: F) -> Complex<F> {
19    rho(m0.powi(2), m1, m2) * m0 / convert!(2, F)
20}
21
22pub fn chi_plus<F: Field>(s: F, m1: F, m2: F) -> Complex<F> {
23    Complex::from(F::one() - ((m1 + m2) * (m1 + m2)) / s)
24}
25
26pub fn chi_minus<F: Field>(s: F, m1: F, m2: F) -> Complex<F> {
27    Complex::from(F::one() - ((m1 - m2) * (m1 - m2)) / s)
28}
29
30pub fn rho<F: Field>(s: F, m1: F, m2: F) -> Complex<F> {
31    Complex::sqrt(chi_plus(s, m1, m2) * chi_minus(s, m1, m2))
32}
33
34pub fn blatt_weisskopf<F: Field>(m0: F, m1: F, m2: F, l: usize) -> F {
35    let q = breakup_momentum(m0, m1, m2);
36    let z = q.powi(2) / convert!(0.1973, F).powi(2);
37    match l {
38        0 => F::one(),
39        1 => F::sqrt((convert!(2, F) * z) / (z + F::one())),
40        2 => F::sqrt(
41            (convert!(13.0, F) * z.powi(2)) / ((z - convert!(3, F)).powi(2) + convert!(9, F) * z),
42        ),
43        3 => F::sqrt(
44            (convert!(277.0, F) * z.powi(3))
45                / (z * (z - convert!(15.0, F)).powi(2)
46                    + convert!(9, F) * (convert!(2, F) * z - convert!(5, F)).powi(2)),
47        ),
48        4 => F::sqrt(
49            (convert!(12746.0, F) * z.powi(4))
50                / (z.powi(2) - convert!(45.0, F) * z + convert!(105.0, F)).powi(2)
51                + convert!(25.0, F) * z * (convert!(2, F) * z - convert!(21.0, F)).powi(2),
52        ),
53        l => panic!("L = {l} is not yet implemented"),
54    }
55}
56
57/// Computes the ([`Complex<F>`]) Blatt-Weisskopf barrier factor representing the energy required for a particle
58/// with mass `m0` to decay into two particles with masses `m1` and `m2` and angular momentum `l`.
59///
60/// In applications where `m0` is expected to be above the mass threshold to produce `m1` and
61/// `m2`, the absolute value of this function can be safely assumed to be equal to its value.
62pub fn blatt_weisskopf_c<F: Field>(m0: F, m1: F, m2: F, l: usize) -> Complex<F> {
63    let q = breakup_momentum_c(m0, m1, m2);
64    let z = q.powi(2) / convert!(0.1973, F).powi(2);
65    match l {
66        0 => Complex::from(F::one()),
67        1 => Complex::sqrt((Complex::from(convert!(2, F)) * z) / (z + F::one())),
68        2 => Complex::sqrt(
69            (z.powi(2) * convert!(13.0, F)) / ((z - convert!(3, F)).powi(2) + z * convert!(9, F)),
70        ),
71        3 => Complex::sqrt(
72            (z.powi(3) * convert!(277.0, F))
73                / (z * (z - convert!(15.0, F)).powi(2)
74                    + (z * convert!(2, F) - convert!(5, F)).powi(2))
75                * convert!(9, F),
76        ),
77        4 => Complex::sqrt(
78            (z.powi(4) * convert!(12746.0, F))
79                / (z.powi(2) - z * convert!(45.0, F) + convert!(105.0, F)).powi(2)
80                + z * convert!(25.0, F) * (z * convert!(2, F) - convert!(21.0, F)).powi(2),
81        ),
82        l => panic!("L = {l} is not yet implemented"),
83    }
84}
85
86pub fn small_wigner_d_matrix<F: Field>(beta: F, j: usize, m: isize, n: isize) -> F {
87    let jpm = (j as i32 + m as i32) as u32;
88    let jmm = (j as i32 - m as i32) as u32;
89    let jpn = (j as i32 + n as i32) as u32;
90    let jmn = (j as i32 - n as i32) as u32;
91    let prefactor = F::sqrt(convert!(
92        jpm.factorial() * jmm.factorial() * jpn.factorial() * jmn.factorial(),
93        F
94    ));
95    let s_min = isize::max(0, n - m) as usize;
96    let s_max = isize::min(jpn as isize, jmm as isize) as usize;
97    let sum: F = (s_min..=s_max)
98        .map(|s| {
99            (F::powi(-F::one(), m as i32 - n as i32 + s as i32)
100                * (F::cos(beta / convert!(2, F))
101                    .powi(2 * (j as i32) + n as i32 - m as i32 - 2 * (s as i32)))
102                * (F::sin(beta / convert!(2, F)).powi(m as i32 - n as i32 + 2 * s as i32)))
103                / convert!(
104                    (jpm - s as u32).factorial()
105                        * (s as u32).factorial()
106                        * ((m - n + s as isize) as u32).factorial()
107                        * (jmm - s as u32).factorial(),
108                    F
109                )
110        })
111        .sum();
112    prefactor * sum
113}
114
115pub fn wigner_d_matrix<F: Field>(
116    alpha: F,
117    beta: F,
118    gamma: F,
119    j: usize,
120    m: isize,
121    n: isize,
122) -> Complex<F> {
123    Complex::cis(convert!(-m, F) * alpha)
124        * small_wigner_d_matrix(beta, j, m, n)
125        * Complex::cis(convert!(-n, F) * gamma)
126}
127
128#[derive(Clone, Copy, Default, PartialEq)]
129#[rustfmt::skip]
130pub enum Wave {
131    #[default]
132    S,
133    S0,
134    Pn1, P0, P1, P,
135    Dn2, Dn1, D0, D1, D2, D,
136    Fn3, Fn2, Fn1, F0, F1, F2, F3, F,
137}
138
139#[rustfmt::skip]
140impl Wave {
141    pub fn new(l: usize, m: isize) -> Self {
142        match l {
143            0 => Self::S0,
144            1 => match m {
145                -1 => Self::Pn1,
146                0 => Self::P0,
147                1 => Self::P1,
148                _ => panic!("|m = {m}| > (l = {l})")
149            }
150            2 => match m {
151                -2 => Self::Dn2,
152                -1 => Self::Dn1,
153                0 => Self::D0,
154                1 => Self::D1,
155                2 => Self::D2,
156                _ => panic!("|m = {m}| > (l = {l})")
157            }
158            3 => match m {
159                -3 => Self::Fn3,
160                -2 => Self::Fn2,
161                -1 => Self::Fn1,
162                0 => Self::F0,
163                1 => Self::F1,
164                2 => Self::F2,
165                3 => Self::F3,
166                _ => panic!("|m = {m}| > (l = {l})")
167            }
168            _ => panic!("(l = {l}) > 3 is not yet implemented!")
169        }
170    }
171    pub fn l(&self) -> i64 {
172        match self {
173            Self::S0 | Self::S => 0,
174            Self::Pn1 | Self::P0 | Self::P1 | Self::P => 1,
175            Self::Dn2 | Self::Dn1 | Self::D0 | Self::D1 | Self::D2 | Self::D => 2,
176            Self::Fn3 | Self::Fn2 | Self::Fn1 | Self::F0 | Self::F1 | Self::F2 | Self::F3 | Self::F => 3,
177        }
178    }
179    pub fn m(&self) -> i64 {
180        match self {
181            Self::S | Self::P | Self::D | Self::F => 0,
182            Self::S0 | Self::P0 | Self::D0 | Self::F0 => 0,
183            Self::Pn1 | Self::Dn1 | Self::Fn1 => -1,
184            Self::P1 | Self::D1 | Self::F1 => 1,
185            Self::Dn2 | Self::Fn2 => -2,
186            Self::D2 | Self::F2 => 2,
187            Self::Fn3 => -3,
188            Self::F3 => 3,
189        }
190    }
191}
192
193impl Display for Wave {
194    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
195        let l_string = match self.l() {
196            0 => "S",
197            1 => "P",
198            2 => "D",
199            3 => "F",
200            _ => unimplemented!(),
201        };
202        write!(f, "{} {:+}", l_string, self.m())
203    }
204}
205
206#[derive(Copy, Clone, PartialEq)]
207pub enum Frame {
208    Helicity,
209    GottfriedJackson,
210}
211
212#[derive(Debug, PartialEq, Eq, Error)]
213#[error("Unknown frame: {0}")]
214pub struct ParseFrameError(String);
215
216impl From<ParseFrameError> for RustitudeError {
217    fn from(value: ParseFrameError) -> Self {
218        RustitudeError::ParseError(value.to_string())
219    }
220}
221
222impl FromStr for Frame {
223    type Err = ParseFrameError;
224
225    fn from_str(s: &str) -> Result<Self, Self::Err> {
226        match s.to_lowercase().as_ref() {
227            "helicity" | "hx" => Ok(Frame::Helicity),
228            "gottfried-jackson" | "gj" => Ok(Frame::GottfriedJackson),
229            _ => Err(ParseFrameError(s.to_string())),
230        }
231    }
232}
233
234pub fn coordinates<F: Field + 'static>(
235    x: &Vector3<F>,
236    y: &Vector3<F>,
237    z: &Vector3<F>,
238    p: &Vector3<F>,
239) -> Coordinates<F> {
240    Coordinates::cartesian(p.dot(x), p.dot(y), p.dot(z))
241}
242
243impl Frame {
244    pub fn coordinates<F: Field>(
245        &self,
246        decay: Decay,
247        other_p4: &FourMomentum<F>,
248        event: &Event<F>,
249    ) -> (Vector3<F>, Vector3<F>, Vector3<F>, Coordinates<F>) {
250        let resonance_p4 = decay.resonance_p4(event);
251        let beam_res_vec = event.beam_p4.boost_along(&resonance_p4).momentum();
252        let recoil_res_vec = event.recoil_p4.boost_along(&resonance_p4).momentum();
253        let other_res_vec = other_p4.boost_along(&resonance_p4).momentum();
254        let (x, y, z) = match self {
255            Frame::Helicity => {
256                let z = -recoil_res_vec.unit();
257                let y = beam_res_vec.cross(&z).unit();
258                let x = y.cross(&z);
259                (x, y, z)
260            }
261            Frame::GottfriedJackson => {
262                let z = beam_res_vec.unit();
263                let y = event.beam_p4.momentum().cross(&(-recoil_res_vec)).unit();
264                let x = y.cross(&z);
265                (x, y, z)
266            }
267        };
268        (x, y, z, coordinates(&x, &y, &z, &other_res_vec))
269    }
270    pub fn coordinates_from_boosted_vec<F: Field>(
271        &self,
272        decay: Decay,
273        other_res_vec: &Vector3<F>,
274        event: &Event<F>,
275    ) -> (Vector3<F>, Vector3<F>, Vector3<F>, Coordinates<F>) {
276        let resonance_p4 = decay.resonance_p4(event);
277        let beam_res_vec = event.beam_p4.boost_along(&resonance_p4).momentum();
278        let recoil_res_vec = event.recoil_p4.boost_along(&resonance_p4).momentum();
279        let (x, y, z) = match self {
280            Frame::Helicity => {
281                let z = -recoil_res_vec.unit();
282                let y = beam_res_vec.cross(&z).unit();
283                let x = y.cross(&z);
284                (x, y, z)
285            }
286            Frame::GottfriedJackson => {
287                let z = beam_res_vec.unit();
288                let y = event.beam_p4.momentum().cross(&(-recoil_res_vec)).unit();
289                let x = y.cross(&z);
290                (x, y, z)
291            }
292        };
293        (x, y, z, coordinates(&x, &y, &z, other_res_vec))
294    }
295}
296
297#[derive(Copy, Clone, PartialEq)]
298pub enum Sign {
299    Positive = 1,
300    Negative = -1,
301}
302
303#[derive(Debug, PartialEq, Eq, Error)]
304#[error("Unknown sign: {0}")]
305pub struct ParseSignError(String);
306
307impl From<ParseSignError> for RustitudeError {
308    fn from(value: ParseSignError) -> Self {
309        RustitudeError::ParseError(value.to_string())
310    }
311}
312
313impl FromStr for Sign {
314    type Err = ParseSignError;
315
316    fn from_str(s: &str) -> Result<Self, Self::Err> {
317        match s.to_lowercase().as_ref() {
318            "positive" | "pos" | "p" | "+" | "plus" | "+1" => Ok(Sign::Positive),
319            "negative" | "neg" | "n" | "-" | "minus" | "m" | "-1" => Ok(Sign::Negative),
320            _ => Err(ParseSignError(s.to_string())),
321        }
322    }
323}
324
325impl Display for Sign {
326    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
327        match self {
328            Sign::Positive => write!(f, "+"),
329            Sign::Negative => write!(f, "-"),
330        }
331    }
332}
333
334#[derive(Clone, Copy)]
335pub enum Decay {
336    TwoBodyDecay([usize; 2]),
337    ThreeBodyDecay([usize; 3]),
338}
339#[derive(Debug, PartialEq, Eq, Error)]
340pub enum ParseDecayError {
341    #[error("Invalid format: {0}")]
342    InvalidFormat(String),
343    #[error("Invalid number: {0}")]
344    InvalidNumber(#[from] ParseIntError),
345    #[error("Invalid number of final states: {0}")]
346    InvalidLength(usize),
347}
348
349impl From<ParseDecayError> for RustitudeError {
350    fn from(value: ParseDecayError) -> Self {
351        RustitudeError::ParseError(value.to_string())
352    }
353}
354
355impl FromStr for Decay {
356    type Err = ParseDecayError;
357
358    fn from_str(s: &str) -> Result<Self, Self::Err> {
359        let cleaned = s.replace(|c: char| c.is_whitespace() || c == '[' || c == ']', "");
360        let parts: Vec<&str> = cleaned.split(',').collect();
361        match parts.len() {
362            2 => {
363                let values = [
364                    parts[0].parse().map_err(ParseDecayError::InvalidNumber)?,
365                    parts[1].parse().map_err(ParseDecayError::InvalidNumber)?,
366                ];
367                Ok(Decay::TwoBodyDecay(values))
368            }
369            3 => {
370                let values = [
371                    parts[0].parse().map_err(ParseDecayError::InvalidNumber)?,
372                    parts[1].parse().map_err(ParseDecayError::InvalidNumber)?,
373                    parts[2].parse().map_err(ParseDecayError::InvalidNumber)?,
374                ];
375                Ok(Decay::ThreeBodyDecay(values))
376            }
377            _ => Err(ParseDecayError::InvalidLength(parts.len())),
378        }
379    }
380}
381
382impl Default for Decay {
383    fn default() -> Self {
384        Self::TwoBodyDecay([0, 1])
385    }
386}
387
388impl Decay {
389    pub fn resonance_p4<F: Field>(&self, event: &Event<F>) -> FourMomentum<F> {
390        match self {
391            Decay::TwoBodyDecay(inds) => inds.iter().map(|&i| event.daughter_p4s[i]).sum(),
392            Decay::ThreeBodyDecay(inds) => inds.iter().map(|&i| event.daughter_p4s[i]).sum(),
393        }
394    }
395    pub fn daughter_p4<'a, F: Field>(
396        &self,
397        index: usize,
398        event: &'a Event<F>,
399    ) -> &'a FourMomentum<F> {
400        match self {
401            Decay::TwoBodyDecay(inds) => &event.daughter_p4s[inds[index]],
402            Decay::ThreeBodyDecay(inds) => &event.daughter_p4s[inds[index]],
403        }
404    }
405    pub fn primary_p4<'a, F: Field>(&self, event: &'a Event<F>) -> &'a FourMomentum<F> {
406        match self {
407            Decay::TwoBodyDecay(inds) => &event.daughter_p4s[inds[0]],
408            Decay::ThreeBodyDecay(inds) => &event.daughter_p4s[inds[0]],
409        }
410    }
411    pub fn secondary_p4<'a, F: Field>(&self, event: &'a Event<F>) -> &'a FourMomentum<F> {
412        match self {
413            Decay::TwoBodyDecay(inds) => &event.daughter_p4s[inds[1]],
414            Decay::ThreeBodyDecay(inds) => &event.daughter_p4s[inds[1]],
415        }
416    }
417    pub fn tertiary_p4<'a, F: Field>(&self, event: &'a Event<F>) -> &'a FourMomentum<F> {
418        match self {
419            Decay::TwoBodyDecay(_) => panic!(),
420            Decay::ThreeBodyDecay(inds) => &event.daughter_p4s[inds[2]],
421        }
422    }
423    pub fn coordinates<F: Field>(
424        &self,
425        frame: Frame,
426        index: usize,
427        event: &Event<F>,
428    ) -> (Vector3<F>, Vector3<F>, Vector3<F>, Coordinates<F>) {
429        frame.coordinates(*self, self.daughter_p4(index, event), event)
430    }
431}
rustitude_gluex/utils.rs

rustitude_gluex/
utils.rs
