deep_causality_physics/nuclear/
quantities.rs

1/*
2 * SPDX-License-Identifier: MIT
3 * Copyright (c) 2023 - 2026. The DeepCausality Authors and Contributors. All Rights Reserved.
4 */
5
6use crate::error::PhysicsError;
7
8/// Amount of Substance (Moles).
9#[derive(Debug, Clone, Copy, PartialEq, PartialOrd, Default)]
10pub struct AmountOfSubstance(f64);
11
12impl AmountOfSubstance {
13    pub fn new(val: f64) -> Result<Self, PhysicsError> {
14        if val < 0.0 {
15            return Err(PhysicsError::PhysicalInvariantBroken(
16                "Negative AmountOfSubstance".into(),
17            ));
18        }
19        Ok(Self(val))
20    }
21    pub fn new_unchecked(val: f64) -> Self {
22        Self(val)
23    }
24    pub fn value(&self) -> f64 {
25        self.0
26    }
27}
28impl From<AmountOfSubstance> for f64 {
29    fn from(val: AmountOfSubstance) -> Self {
30        val.0
31    }
32}
33
34/// Half-Life (Seconds).
35#[derive(Debug, Clone, Copy, PartialEq, PartialOrd, Default)]
36pub struct HalfLife(f64);
37
38impl HalfLife {
39    /// Creates a new `HalfLife` instance.
40    ///
41    /// # Errors
42    /// Returns `PhysicsError` if `val <= 0.0`.
43    pub fn new(val: f64) -> Result<Self, PhysicsError> {
44        if !val.is_finite() {
45            return Err(PhysicsError::PhysicalInvariantBroken(format!(
46                "HalfLife must be finite: {}",
47                val
48            )));
49        }
50        if val <= 0.0 {
51            return Err(PhysicsError::PhysicalInvariantBroken(
52                "HalfLife must be positive (zero implies infinite decay rate)".into(),
53            ));
54        }
55        Ok(Self(val))
56    }
57    pub fn new_unchecked(val: f64) -> Self {
58        Self(val)
59    }
60    pub fn value(&self) -> f64 {
61        self.0
62    }
63}
64impl From<HalfLife> for f64 {
65    fn from(val: HalfLife) -> Self {
66        val.0
67    }
68}
69
70/// Radioactivity (Becquerels).
71#[derive(Debug, Clone, Copy, PartialEq, PartialOrd, Default)]
72pub struct Activity(f64);
73
74impl Activity {
75    pub fn new(val: f64) -> Result<Self, PhysicsError> {
76        if val < 0.0 {
77            return Err(PhysicsError::PhysicalInvariantBroken(
78                "Negative Activity".into(),
79            ));
80        }
81        Ok(Self(val))
82    }
83    pub fn new_unchecked(val: f64) -> Self {
84        Self(val)
85    }
86    pub fn value(&self) -> f64 {
87        self.0
88    }
89}
90impl From<Activity> for f64 {
91    fn from(val: Activity) -> Self {
92        val.0
93    }
94}
95
96/// Energy Density (Joules per cubic meter).
97#[derive(Debug, Clone, Copy, PartialEq, PartialOrd, Default)]
98pub struct EnergyDensity(f64);
99
100impl EnergyDensity {
101    pub fn new(val: f64) -> Result<Self, PhysicsError> {
102        if val < 0.0 {
103            return Err(PhysicsError::PhysicalInvariantBroken(
104                "Negative EnergyDensity".into(),
105            ));
106        }
107        Ok(Self(val))
108    }
109    pub fn new_unchecked(val: f64) -> Self {
110        Self(val)
111    }
112    pub fn value(&self) -> f64 {
113        self.0
114    }
115}
116impl From<EnergyDensity> for f64 {
117    fn from(val: EnergyDensity) -> Self {
118        val.0
119    }
120}
121
122// =============================================================================
123// Lund String Fragmentation Types
124// =============================================================================
125
126use core::ops::{Add, Sub};
127
128/// Lorentz 4-momentum (E, px, py, pz) in natural units (c = 1).
129///
130/// Components are stored in (+---) signature convention (particle physics).
131/// Energy E is component 0, spatial momentum (px, py, pz) are components 1-3.
132#[derive(Debug, Clone, Copy, Default, PartialEq)]
133pub struct FourMomentum {
134    /// Energy component E (GeV)
135    e: f64,
136    /// x-component of 3-momentum px (GeV/c)
137    px: f64,
138    /// y-component of 3-momentum py (GeV/c)
139    py: f64,
140    /// z-component of 3-momentum pz (GeV/c)
141    pz: f64,
142}
143
144impl FourMomentum {
145    /// Creates a new 4-momentum.
146    pub const fn new(e: f64, px: f64, py: f64, pz: f64) -> Self {
147        Self { e, px, py, pz }
148    }
149
150    /// Creates a 4-momentum from mass and 3-momentum.
151    ///
152    /// Computes E = √(m² + |p|²)
153    pub fn from_mass_and_momentum(mass: f64, px: f64, py: f64, pz: f64) -> Self {
154        let p_sq = px * px + py * py + pz * pz;
155        let e = (mass * mass + p_sq).sqrt();
156        Self { e, px, py, pz }
157    }
158
159    /// Creates a 4-momentum for a particle at rest.
160    pub const fn at_rest(mass: f64) -> Self {
161        Self {
162            e: mass,
163            px: 0.0,
164            py: 0.0,
165            pz: 0.0,
166        }
167    }
168
169    /// Energy component E.
170    pub const fn e(&self) -> f64 {
171        self.e
172    }
173
174    /// x-component of 3-momentum.
175    pub const fn px(&self) -> f64 {
176        self.px
177    }
178
179    /// y-component of 3-momentum.
180    pub const fn py(&self) -> f64 {
181        self.py
182    }
183
184    /// z-component of 3-momentum.
185    pub const fn pz(&self) -> f64 {
186        self.pz
187    }
188
189    /// Lorentz-invariant mass squared: m² = E² - |p|²
190    pub fn invariant_mass_squared(&self) -> f64 {
191        self.e * self.e - self.px * self.px - self.py * self.py - self.pz * self.pz
192    }
193
194    /// Lorentz-invariant mass: m = √(E² - |p|²)
195    ///
196    /// Returns 0.0 for massless or spacelike 4-vectors.
197    pub fn invariant_mass(&self) -> f64 {
198        let m_sq = self.invariant_mass_squared();
199        if m_sq > 0.0 { m_sq.sqrt() } else { 0.0 }
200    }
201
202    /// Magnitude of 3-momentum: |p| = √(px² + py² + pz²)
203    pub fn momentum_magnitude(&self) -> f64 {
204        (self.px * self.px + self.py * self.py + self.pz * self.pz).sqrt()
205    }
206
207    /// Transverse momentum: pT = √(px² + py²)
208    pub fn transverse_momentum(&self) -> f64 {
209        (self.px * self.px + self.py * self.py).sqrt()
210    }
211
212    /// Transverse mass: mT = √(m² + pT²)
213    pub fn transverse_mass(&self) -> f64 {
214        let m_sq = self.invariant_mass_squared();
215        let pt_sq = self.px * self.px + self.py * self.py;
216        (m_sq + pt_sq).sqrt()
217    }
218
219    /// Rapidity: y = 0.5 * ln((E + pz) / (E - pz))
220    ///
221    /// Returns 0.0 if denominator is near zero.
222    pub fn rapidity(&self) -> f64 {
223        let denom = self.e - self.pz;
224        if denom.abs() < 1e-10 {
225            return 0.0;
226        }
227        0.5 * ((self.e + self.pz) / denom).ln()
228    }
229
230    /// Pseudorapidity: η = -ln(tan(θ/2))
231    ///
232    /// where θ is the polar angle from the beam axis (z).
233    pub fn pseudorapidity(&self) -> f64 {
234        let p = self.momentum_magnitude();
235        if p.abs() < 1e-10 {
236            return 0.0;
237        }
238        0.5 * ((p + self.pz) / (p - self.pz)).ln()
239    }
240
241    /// Azimuthal angle φ in the transverse plane.
242    pub fn phi(&self) -> f64 {
243        self.py.atan2(self.px)
244    }
245
246    /// Lightcone plus component: p+ = E + pz
247    pub fn lightcone_plus(&self) -> f64 {
248        self.e + self.pz
249    }
250
251    /// Lightcone minus component: p- = E - pz
252    pub fn lightcone_minus(&self) -> f64 {
253        self.e - self.pz
254    }
255
256    /// Lorentz boost along z-axis.
257    pub fn boost_z(&self, beta: f64) -> Self {
258        let gamma = 1.0 / (1.0 - beta * beta).sqrt();
259        let e_new = gamma * (self.e - beta * self.pz);
260        let pz_new = gamma * (self.pz - beta * self.e);
261        Self {
262            e: e_new,
263            px: self.px,
264            py: self.py,
265            pz: pz_new,
266        }
267    }
268}
269
270impl Add for FourMomentum {
271    type Output = Self;
272
273    fn add(self, rhs: Self) -> Self::Output {
274        Self {
275            e: self.e + rhs.e,
276            px: self.px + rhs.px,
277            py: self.py + rhs.py,
278            pz: self.pz + rhs.pz,
279        }
280    }
281}
282
283impl Sub for FourMomentum {
284    type Output = Self;
285
286    fn sub(self, rhs: Self) -> Self::Output {
287        Self {
288            e: self.e - rhs.e,
289            px: self.px - rhs.px,
290            py: self.py - rhs.py,
291            pz: self.pz - rhs.pz,
292        }
293    }
294}
295
296/// A produced hadron from string fragmentation.
297#[derive(Debug, Clone, PartialEq)]
298pub struct Hadron {
299    /// PDG Monte Carlo particle ID
300    pdg_id: i32,
301    /// 4-momentum in the lab frame
302    momentum: FourMomentum,
303}
304
305impl Hadron {
306    /// Creates a new hadron.
307    pub const fn new(pdg_id: i32, momentum: FourMomentum) -> Self {
308        Self { pdg_id, momentum }
309    }
310
311    /// PDG Monte Carlo particle ID.
312    pub const fn pdg_id(&self) -> i32 {
313        self.pdg_id
314    }
315
316    /// 4-momentum in the lab frame.
317    pub const fn momentum(&self) -> FourMomentum {
318        self.momentum
319    }
320
321    /// Invariant mass of the hadron.
322    pub fn mass(&self) -> f64 {
323        self.momentum.invariant_mass()
324    }
325
326    /// Energy of the hadron.
327    pub fn energy(&self) -> f64 {
328        self.momentum.e()
329    }
330
331    /// Transverse momentum.
332    pub fn pt(&self) -> f64 {
333        self.momentum.transverse_momentum()
334    }
335
336    /// Rapidity.
337    pub fn rapidity(&self) -> f64 {
338        self.momentum.rapidity()
339    }
340}
341
342/// Configuration parameters for Lund String Fragmentation.
343///
344/// The Lund symmetric fragmentation function is:
345/// $$ f(z) = \frac{1}{z} (1-z)^a \exp\left(-\frac{b \cdot m_T^2}{z}\right) $$
346///
347/// where z is the lightcone momentum fraction and mT is the transverse mass.
348#[derive(Debug, Clone, PartialEq)]
349pub struct LundParameters {
350    /// String tension κ (GeV/fm) - energy per unit length of the color flux tube
351    kappa: f64,
352    /// Lund 'a' parameter - controls the shape of the fragmentation function
353    lund_a: f64,
354    /// Lund 'b' parameter (GeV⁻²) - controls mass dependence
355    lund_b: f64,
356    /// Transverse momentum width σ_pT (GeV) - Gaussian width for pT generation
357    sigma_pt: f64,
358    /// Strange quark suppression factor s/u (typically 0.2-0.4)
359    strange_suppression: f64,
360    /// Diquark suppression factor qq/q (typically 0.05-0.15)
361    diquark_suppression: f64,
362    /// Spin-1 meson suppression relative to spin-0 (vector/pseudoscalar ratio)
363    vector_meson_fraction: f64,
364    /// Minimum invariant mass to continue fragmentation (GeV)
365    min_invariant_mass: f64,
366}
367
368impl Default for LundParameters {
369    /// Default parameters tuned to LEP e+e- data.
370    ///
371    /// Based on PYTHIA 8 Monash tune.
372    fn default() -> Self {
373        Self {
374            kappa: 1.0,                  // GeV/fm
375            lund_a: 0.68,                // dimensionless
376            lund_b: 0.98,                // GeV⁻²
377            sigma_pt: 0.36,              // GeV
378            strange_suppression: 0.30,   // s/u ratio
379            diquark_suppression: 0.10,   // qq/q ratio
380            vector_meson_fraction: 0.50, // V/(V+PS)
381            min_invariant_mass: 0.5,     // GeV
382        }
383    }
384}
385
386impl LundParameters {
387    /// Creates new Lund parameters.
388    #[allow(clippy::too_many_arguments)]
389    pub const fn new(
390        kappa: f64,
391        lund_a: f64,
392        lund_b: f64,
393        sigma_pt: f64,
394        strange_suppression: f64,
395        diquark_suppression: f64,
396        vector_meson_fraction: f64,
397        min_invariant_mass: f64,
398    ) -> Self {
399        Self {
400            kappa,
401            lund_a,
402            lund_b,
403            sigma_pt,
404            strange_suppression,
405            diquark_suppression,
406            vector_meson_fraction,
407            min_invariant_mass,
408        }
409    }
410
411    /// String tension κ (GeV/fm).
412    pub const fn kappa(&self) -> f64 {
413        self.kappa
414    }
415
416    /// Lund 'a' parameter.
417    pub const fn lund_a(&self) -> f64 {
418        self.lund_a
419    }
420
421    /// Lund 'b' parameter (GeV⁻²).
422    pub const fn lund_b(&self) -> f64 {
423        self.lund_b
424    }
425
426    /// Transverse momentum width σ_pT (GeV).
427    pub const fn sigma_pt(&self) -> f64 {
428        self.sigma_pt
429    }
430
431    /// Strange quark suppression factor s/u.
432    pub const fn strange_suppression(&self) -> f64 {
433        self.strange_suppression
434    }
435
436    /// Diquark suppression factor qq/q.
437    pub const fn diquark_suppression(&self) -> f64 {
438        self.diquark_suppression
439    }
440
441    /// Vector meson fraction V/(V+PS).
442    pub const fn vector_meson_fraction(&self) -> f64 {
443        self.vector_meson_fraction
444    }
445
446    /// Minimum invariant mass to continue fragmentation (GeV).
447    pub const fn min_invariant_mass(&self) -> f64 {
448        self.min_invariant_mass
449    }
450}
deep_causality_physics/nuclear/quantities.rs

deep_causality_physics/nuclear/
quantities.rs
