laddu_amplitudes/
breit_wigner.rs

1use serde::{Deserialize, Serialize};
2
3use laddu_core::{
4    amplitudes::{Amplitude, AmplitudeID, ParameterLike},
5    data::Event,
6    resources::{Cache, ParameterID, Parameters, Resources},
7    utils::{
8        functions::{blatt_weisskopf, breakup_momentum},
9        variables::{Mass, Variable},
10    },
11    Complex, DVector, Float, LadduError, PI,
12};
13
14#[cfg(feature = "python")]
15use laddu_python::{
16    amplitudes::{PyAmplitude, PyParameterLike},
17    utils::variables::PyMass,
18};
19#[cfg(feature = "python")]
20use pyo3::prelude::*;
21
22/// A relativistic Breit-Wigner [`Amplitude`], parameterized as follows:
23/// ```math
24/// I_{\ell}(m; m_0, \Gamma_0, m_1, m_2) =  \frac{1}{\pi}\frac{m_0 \Gamma_0 B_{\ell}(m, m_1, m_2)}{(m_0^2 - m^2) - \imath m_0 \Gamma}
25/// ```
26/// where
27/// ```math
28/// \Gamma = \Gamma_0 \frac{m_0}{m} \frac{q(m, m_1, m_2)}{q(m_0, m_1, m_2)} \left(\frac{B_{\ell}(m, m_1, m_2)}{B_{\ell}(m_0, m_1, m_2)}\right)^2
29/// ```
30/// is the relativistic width correction, $`q(m_a, m_b, m_c)`$ is the breakup momentum of a particle with mass $`m_a`$ decaying into two particles with masses $`m_b`$ and $`m_c`$, $`B_{\ell}(m_a, m_b, m_c)`$ is the Blatt-Weisskopf barrier factor for the same decay assuming particle $`a`$ has angular momentum $`\ell`$, $`m_0`$ is the mass of the resonance, $`\Gamma_0`$ is the nominal width of the resonance, $`m_1`$ and $`m_2`$ are the masses of the decay products, and $`m`$ is the "input" mass.
31#[derive(Clone, Serialize, Deserialize)]
32pub struct BreitWigner {
33    name: String,
34    mass: ParameterLike,
35    width: ParameterLike,
36    pid_mass: ParameterID,
37    pid_width: ParameterID,
38    l: usize,
39    daughter_1_mass: Mass,
40    daughter_2_mass: Mass,
41    resonance_mass: Mass,
42}
43impl BreitWigner {
44    /// Construct a [`BreitWigner`] with the given name, mass, width, and angular momentum (`l`).
45    /// This uses the given `resonance_mass` as the "input" mass and two daughter masses of the
46    /// decay products to determine phase-space and Blatt-Weisskopf factors.
47    pub fn new(
48        name: &str,
49        mass: ParameterLike,
50        width: ParameterLike,
51        l: usize,
52        daughter_1_mass: &Mass,
53        daughter_2_mass: &Mass,
54        resonance_mass: &Mass,
55    ) -> Box<Self> {
56        Self {
57            name: name.to_string(),
58            mass,
59            width,
60            pid_mass: ParameterID::default(),
61            pid_width: ParameterID::default(),
62            l,
63            daughter_1_mass: daughter_1_mass.clone(),
64            daughter_2_mass: daughter_2_mass.clone(),
65            resonance_mass: resonance_mass.clone(),
66        }
67        .into()
68    }
69}
70
71#[typetag::serde]
72impl Amplitude for BreitWigner {
73    fn register(&mut self, resources: &mut Resources) -> Result<AmplitudeID, LadduError> {
74        self.pid_mass = resources.register_parameter(&self.mass);
75        self.pid_width = resources.register_parameter(&self.width);
76        resources.register_amplitude(&self.name)
77    }
78
79    fn compute(&self, parameters: &Parameters, event: &Event, _cache: &Cache) -> Complex<Float> {
80        let mass = self.resonance_mass.value(event);
81        let mass0 = parameters.get(self.pid_mass).abs();
82        let width0 = parameters.get(self.pid_width).abs();
83        let mass1 = self.daughter_1_mass.value(event);
84        let mass2 = self.daughter_2_mass.value(event);
85        let q0 = breakup_momentum(mass0, mass1, mass2);
86        let q = breakup_momentum(mass, mass1, mass2);
87        let f0 = blatt_weisskopf(mass0, mass1, mass2, self.l);
88        let f = blatt_weisskopf(mass, mass1, mass2, self.l);
89        let width = width0 * (mass0 / mass) * (q / q0) * (f / f0).powi(2);
90        let n = Float::sqrt(mass0 * width0 / PI);
91        let d = Complex::new(mass0.powi(2) - mass.powi(2), -(mass0 * width));
92        Complex::from(f * n) / d
93    }
94
95    fn compute_gradient(
96        &self,
97        parameters: &Parameters,
98        event: &Event,
99        cache: &Cache,
100        gradient: &mut DVector<Complex<Float>>,
101    ) {
102        let mut indices = Vec::with_capacity(2);
103        if let ParameterID::Parameter(index) = self.pid_mass {
104            indices.push(index)
105        }
106        if let ParameterID::Parameter(index) = self.pid_width {
107            indices.push(index)
108        }
109        self.central_difference_with_indices(&indices, parameters, event, cache, gradient)
110    }
111}
112
113/// An relativistic Breit-Wigner Amplitude
114///
115/// This Amplitude represents a relativistic Breit-Wigner with known angular momentum
116///
117/// Parameters
118/// ----------
119/// name : str
120///     The Amplitude name
121/// mass : laddu.ParameterLike
122///     The mass of the resonance
123/// width : laddu.ParameterLike
124///     The (nonrelativistic) width of the resonance
125/// l : int
126///     The total orbital momentum (:math:`l > 0`)
127/// daughter_1_mass : laddu.Mass
128///     The mass of the first decay product
129/// daughter_2_mass : laddu.Mass
130///     The mass of the second decay product
131/// resonance_mass: laddu.Mass
132///     The total mass of the resonance
133///
134/// Returns
135/// -------
136/// laddu.Amplitude
137///     An Amplitude which can be registered by a laddu.Manager
138///
139/// See Also
140/// --------
141/// laddu.Manager
142///
143#[cfg(feature = "python")]
144#[pyfunction(name = "BreitWigner")]
145pub fn py_breit_wigner(
146    name: &str,
147    mass: PyParameterLike,
148    width: PyParameterLike,
149    l: usize,
150    daughter_1_mass: &PyMass,
151    daughter_2_mass: &PyMass,
152    resonance_mass: &PyMass,
153) -> PyAmplitude {
154    PyAmplitude(BreitWigner::new(
155        name,
156        mass.0,
157        width.0,
158        l,
159        &daughter_1_mass.0,
160        &daughter_2_mass.0,
161        &resonance_mass.0,
162    ))
163}
164
165#[cfg(test)]
166mod tests {
167    use std::sync::Arc;
168
169    use super::*;
170    use approx::assert_relative_eq;
171    use laddu_core::{data::test_dataset, parameter, Manager, Mass};
172
173    #[test]
174    fn test_bw_evaluation() {
175        let mut manager = Manager::default();
176        let amp = BreitWigner::new(
177            "bw",
178            parameter("mass"),
179            parameter("width"),
180            2,
181            &Mass::new([2]),
182            &Mass::new([3]),
183            &Mass::new([2, 3]),
184        );
185        let aid = manager.register(amp).unwrap();
186
187        let dataset = Arc::new(test_dataset());
188        let expr = aid.into();
189        let model = manager.model(&expr);
190        let evaluator = model.load(&dataset);
191
192        let result = evaluator.evaluate(&[1.5, 0.3]);
193
194        assert_relative_eq!(result[0].re, 1.45856917, epsilon = Float::EPSILON.sqrt());
195        assert_relative_eq!(result[0].im, 1.4107341, epsilon = Float::EPSILON.sqrt());
196    }
197
198    #[test]
199    fn test_bw_gradient() {
200        let mut manager = Manager::default();
201        let amp = BreitWigner::new(
202            "bw",
203            parameter("mass"),
204            parameter("width"),
205            2,
206            &Mass::new([2]),
207            &Mass::new([3]),
208            &Mass::new([2, 3]),
209        );
210        let aid = manager.register(amp).unwrap();
211
212        let dataset = Arc::new(test_dataset());
213        let expr = aid.into();
214        let model = manager.model(&expr);
215        let evaluator = model.load(&dataset);
216        dbg!(Mass::new([2, 3]).value_on(&dataset));
217
218        let result = evaluator.evaluate_gradient(&[1.7, 0.3]);
219        assert_relative_eq!(result[0][0].re, -2.410585, epsilon = Float::EPSILON.cbrt());
220        assert_relative_eq!(result[0][0].im, -1.8880913, epsilon = Float::EPSILON.cbrt());
221        assert_relative_eq!(result[0][1].re, 1.0467031, epsilon = Float::EPSILON.cbrt());
222        assert_relative_eq!(result[0][1].im, 1.3683612, epsilon = Float::EPSILON.cbrt());
223    }
224}
laddu_amplitudes/breit_wigner.rs

laddu_amplitudes/
breit_wigner.rs
