diff --git a/crates/rustitude-gluex/src/dalitz.rs b/crates/rustitude-gluex/src/dalitz.rs index 06f6b83..4289f4f 100644 --- a/crates/rustitude-gluex/src/dalitz.rs +++ b/crates/rustitude-gluex/src/dalitz.rs @@ -1,5 +1,5 @@ use rayon::prelude::*; -use rustitude_core::prelude::*; +use rustitude_core::{convert, prelude::*}; use crate::utils::Decay; @@ -38,22 +38,25 @@ impl Node for OmegaDalitz { let dalitz_t = (pip + pi0).m2(); let dalitz_u = (pim + pi0).m2(); - let m3pi = (F::TWO * pip.m()) + pi0.m(); - let dalitz_d = F::TWO * omega.m() * (omega.m() - m3pi); - let dalitz_sc = (F::ONE / F::THREE) * (omega.m2() + pip.m2() + pim.m2() + pi0.m2()); - let dalitz_x = F::fsqrt(F::THREE) * (dalitz_t - dalitz_u) / dalitz_d; - let dalitz_y = F::THREE * (dalitz_sc - dalitz_s) / dalitz_d; + let m3pi = (convert!(2.0, F) * pip.m()) + pi0.m(); + let dalitz_d = convert!(2.0, F) * omega.m() * (omega.m() - m3pi); + let dalitz_sc = + (F::one() / convert!(3.0, F)) * (omega.m2() + pip.m2() + pim.m2() + pi0.m2()); + let dalitz_x = F::sqrt(convert!(3.0, F)) * (dalitz_t - dalitz_u) / dalitz_d; + let dalitz_y = convert!(3.0, F) * (dalitz_sc - dalitz_s) / dalitz_d; let dalitz_z = dalitz_x * dalitz_x + dalitz_y * dalitz_y; - let dalitz_sin3theta = F::fsin(F::THREE * F::fasin(dalitz_y / F::fsqrt(dalitz_z))); + let dalitz_sin3theta = + F::sin(convert!(3.0, F) * F::asin(dalitz_y / F::sqrt(dalitz_z))); let pip_omega = pip.boost_along(&omega); let pim_omega = pim.boost_along(&omega); let pi_cross = pip_omega.momentum().cross(&pim_omega.momentum()); - let lambda = (F::FOUR / F::THREE) * F::fabs(pi_cross.dot(&pi_cross)) - / ((F::ONE / F::NINE) - * (omega.m2() - (F::TWO * pip.m() + pi0.m()).fpowi(2)).fpowi(2)); + let lambda = (convert!(4.0, F) / convert!(3.0, F)) + * F::abs(pi_cross.dot(&pi_cross)) + / ((F::one() / convert!(9.0, F)) + * (omega.m2() - (convert!(2.0, F) * pip.m() + pi0.m()).powi(2)).powi(2)); (dalitz_z, (dalitz_sin3theta, lambda)) }) @@ -69,13 +72,19 @@ impl Node for OmegaDalitz { let beta = parameters[1]; let gamma = parameters[2]; let delta = parameters[3]; - Ok(F::fsqrt(F::fabs( + Ok(F::sqrt(F::abs( lambda - * (F::ONE - + F::TWO * alpha * dalitz_z - + F::TWO * beta * dalitz_z.fpowf(F::THREE / F::TWO) * dalitz_sin3theta - + F::TWO * gamma * dalitz_z.fpowi(2) - + F::TWO * delta * dalitz_z.fpowf(F::FIVE / F::TWO) * dalitz_sin3theta), + * (F::one() + + convert!(2.0, F) * alpha * dalitz_z + + convert!(2.0, F) + * beta + * dalitz_z.powf(convert!(3.0, F) / convert!(2.0, F)) + * dalitz_sin3theta + + convert!(2.0, F) * gamma * dalitz_z.powi(2) + + convert!(2.0, F) + * delta + * dalitz_z.powf(convert!(5.0, F) / convert!(2.0, F)) + * dalitz_sin3theta), )) .into()) } diff --git a/crates/rustitude-gluex/src/harmonics.rs b/crates/rustitude-gluex/src/harmonics.rs index f1d582a..280197d 100644 --- a/crates/rustitude-gluex/src/harmonics.rs +++ b/crates/rustitude-gluex/src/harmonics.rs @@ -1,5 +1,5 @@ use rayon::prelude::*; -use rustitude_core::prelude::*; +use rustitude_core::{convert, prelude::*}; use sphrs::{ComplexSH, SHEval}; use crate::utils::{Decay, Frame, Sign, Wave}; @@ -68,25 +68,21 @@ impl Node for Zlm { .map(|event| { let (_, y, _, p) = self.decay.coordinates(self.frame, 0, event); let ylm = ComplexSH::Spherical.eval(self.wave.l(), self.wave.m(), &p); - let big_phi = F::fatan2( + let big_phi = F::atan2( y.dot(&event.eps), - event - .beam_p4 - .momentum() - .normalize() - .dot(&event.eps.cross(&y)), + event.beam_p4.direction().dot(&event.eps.cross(&y)), ); - let pgamma = event.eps.norm(); + let pgamma = event.eps_mag(); let phase = Complex::cis(-big_phi); let zlm = ylm * phase; match self.reflectivity { Sign::Positive => Complex::new( - F::fsqrt(F::ONE + pgamma) * zlm.re, - F::fsqrt(F::ONE - pgamma) * zlm.im, + F::sqrt(F::one() + pgamma) * zlm.re, + F::sqrt(F::one() - pgamma) * zlm.im, ), Sign::Negative => Complex::new( - F::fsqrt(F::ONE - pgamma) * zlm.re, - F::fsqrt(F::ONE + pgamma) * zlm.im, + F::sqrt(F::one() - pgamma) * zlm.re, + F::sqrt(F::one() + pgamma) * zlm.im, ), } }) @@ -122,24 +118,21 @@ impl Node for OnePS { .par_iter() .map(|event| { let (_, y, _, _) = self.decay.coordinates(self.frame, 0, event); - let pol_angle = event.eps[0].facos(); - let big_phi = y.dot(&event.eps).fatan2( - event - .beam_p4 - .momentum() - .normalize() - .dot(&event.eps.cross(&y)), + let pol_angle = F::acos(event.eps[0]); + let big_phi = F::atan2( + y.dot(&event.eps), + event.beam_p4.direction().dot(&event.eps.cross(&y)), ); - let pgamma = event.eps.norm(); + let pgamma = event.eps_mag(); let phase = Complex::cis(-(pol_angle + big_phi)); match self.reflectivity { Sign::Positive => Complex::new( - F::fsqrt(F::ONE + pgamma) * phase.re, - F::fsqrt(F::ONE - pgamma) * phase.im, + F::sqrt(F::one() + pgamma) * phase.re, + F::sqrt(F::one() - pgamma) * phase.im, ), Sign::Negative => Complex::new( - F::fsqrt(F::ONE - pgamma) * phase.re, - F::fsqrt(F::ONE + pgamma) * phase.im, + F::sqrt(F::one() - pgamma) * phase.re, + F::sqrt(F::one() + pgamma) * phase.im, ), } }) @@ -184,19 +177,22 @@ impl Node for TwoPS { let ylm_m = ComplexSH::Spherical .eval(self.wave.l(), -self.wave.m(), &p) .conj(); - let m_refl = F::convert_isize(if self.wave.m() % 2 == 0 { - self.reflectivity as isize - } else { - -(self.reflectivity as isize) - }); + let m_refl = convert!( + if self.wave.m() % 2 == 0 { + self.reflectivity as isize + } else { + -(self.reflectivity as isize) + }, + F + ); let big_theta = match self.wave.m().cmp(&0) { - std::cmp::Ordering::Less => F::ZERO, - std::cmp::Ordering::Equal => F::f(0.5), - std::cmp::Ordering::Greater => F::fsqrt(F::f(0.5)), + std::cmp::Ordering::Less => F::zero(), + std::cmp::Ordering::Equal => convert!(0.5, F), + std::cmp::Ordering::Greater => F::FRAC_1_SQRT_2(), }; - let wigner_d_lm0_m = ylm_m.scale(F::fsqrt( - F::FOUR * F::PI() - / (F::TWO * ::convert_usize(self.wave.l() as usize) + F::ONE), + let wigner_d_lm0_m = ylm_m.scale(F::sqrt( + convert!(4, F) * F::PI() + / (convert!(2, F) * convert!(self.wave.l(), F) + F::one()), )); ylm_p.scale(big_theta) - wigner_d_lm0_m.scale(m_refl) }) diff --git a/crates/rustitude-gluex/src/polarization.rs b/crates/rustitude-gluex/src/polarization.rs index 287bb35..beb93d7 100644 --- a/crates/rustitude-gluex/src/polarization.rs +++ b/crates/rustitude-gluex/src/polarization.rs @@ -1,6 +1,6 @@ use crate::utils::{self, Decay, Frame, Sign}; use rayon::prelude::*; -use rustitude_core::prelude::*; +use rustitude_core::{convert, prelude::*}; use sphrs::{ComplexSH, SHEval}; #[derive(Clone)] @@ -35,13 +35,13 @@ impl ThreePiPolFrac { match (decay_resonance, decay_isobar) { (Decay::ThreeBodyDecay(_), Decay::TwoBodyDecay(_)) => Self { beam_pol: match beam_pol { - Sign::Positive => F::ONE, - Sign::Negative => -F::ONE, + Sign::Positive => F::one(), + Sign::Negative => -F::one(), }, j_resonance, p_resonance: match p_resonance { - Sign::Positive => F::ONE, - Sign::Negative => -F::ONE, + Sign::Positive => F::one(), + Sign::Negative => -F::one(), }, i_resonance, l_resonance, @@ -59,7 +59,7 @@ impl ThreePiPolFrac { impl Node for ThreePiPolFrac { fn calculate(&self, parameters: &[F], event: &Event) -> Result, RustitudeError> { - Ok(self.data[event.index] * (F::ONE + self.beam_pol * parameters[0]) / F::FOUR) + Ok(self.data[event.index] * (F::one() + self.beam_pol * parameters[0]) / convert!(4, F)) } fn precalculate(&mut self, dataset: &Dataset) -> Result<(), RustitudeError> { @@ -85,37 +85,44 @@ impl Node for ThreePiPolFrac { self.decay_isobar.primary_p4(event).m(), self.decay_isobar.secondary_p4(event).m(), ); - let neg_res_hel_prod = - Complex::new(F::fcos(F::TWO * alpha), F::fsin(F::TWO * alpha)) - * self.beam_pol - * self.p_resonance - * (F::convert_u32((self.j_resonance % 2) * 2) / F::TWO); - let mut res = F::ZERO.c(); + let neg_res_hel_prod = Complex::new( + F::cos(convert!(2, F) * alpha), + F::sin(convert!(2, F) * alpha), + ) * self.beam_pol + * self.p_resonance + * convert!(self.j_resonance % 2, F); + let mut res = Complex::default(); for m_res in -(self.l_resonance as isize)..(self.l_resonance as isize) { - let mut term = F::ZERO.c(); + let mut term = Complex::default(); for m_iso in -(self.j_isobar as isize)..(self.j_isobar as isize) { let ylm = ComplexSH::Spherical.eval( self.j_isobar as i64, m_iso as i64, &p1_iso_coords, ); - let cg_neg = F::f(wigners::clebsch_gordan( - self.j_isobar, - self.l_resonance as i32, - m_iso as u32, - m_res as i32, - self.j_resonance, - -1, - )); - let cg_pos = F::f(wigners::clebsch_gordan( - self.j_isobar, - self.l_resonance as i32, - m_iso as u32, - m_res as i32, - self.j_resonance, - 1, - )); - term += ylm * (neg_res_hel_prod * cg_neg + cg_pos); + let cg_neg = convert!( + wigners::clebsch_gordan( + self.j_isobar, + self.l_resonance as i32, + m_iso as u32, + m_res as i32, + self.j_resonance, + -1, + ), + F + ); + let cg_pos = convert!( + wigners::clebsch_gordan( + self.j_isobar, + self.l_resonance as i32, + m_iso as u32, + m_res as i32, + self.j_resonance, + 1, + ), + F + ); + term += ylm * Complex::from(neg_res_hel_prod * cg_neg + cg_pos); } let ylm = ComplexSH::Spherical.eval( self.l_resonance as i64, @@ -125,7 +132,7 @@ impl Node for ThreePiPolFrac { term *= ylm; res += term; } - res *= F::f( + res *= convert!( wigners::clebsch_gordan( 1, 1, @@ -141,8 +148,9 @@ impl Node for ThreePiPolFrac { self.i_resonance as u32, (self.iz_daughters[0] + self.iz_daughters[1] + self.iz_daughters[2]) as i32, ), - ) * k.fpowi(self.l_resonance as i32) - * q.fpowi(self.j_isobar as i32); + F + ) * k.powi(self.l_resonance as i32) + * q.powi(self.j_isobar as i32); res }) .collect(); diff --git a/crates/rustitude-gluex/src/resonances.rs b/crates/rustitude-gluex/src/resonances.rs index 3982c87..8351ff4 100644 --- a/crates/rustitude-gluex/src/resonances.rs +++ b/crates/rustitude-gluex/src/resonances.rs @@ -1,9 +1,9 @@ use crate::utils; use crate::utils::Decay; -use nalgebra::{SMatrix, SVector}; +use nalgebra::{RealField, SMatrix, SVector}; use rayon::prelude::*; -use rustitude_core::prelude::*; +use rustitude_core::{convert, convert_array, convert_vec, prelude::*}; #[derive(Default, Clone)] pub struct BreitWigner { @@ -53,9 +53,9 @@ impl Node for BreitWigner { let g0 = parameters[1]; let f0 = utils::blatt_weisskopf(m0, m1, m2, self.l); let q0 = utils::breakup_momentum(m0, m1, m2); - let g = g0 * (m0 / m) * (q / q0) * (f.fpowi(2) / f0.fpowi(2)); - Ok(Complex::new(f * (m0 * g0 / F::PI()), F::ZERO) - / Complex::new(m0.fpowi(2) - m.fpowi(2), -F::ONE * m0 * g)) + let g = g0 * (m0 / m) * (q / q0) * (f.powi(2) / f0.powi(2)); + Ok(Complex::new(f * (m0 * g0 / F::PI()), F::zero()) + / Complex::new(m0.powi(2) - m.powi(2), -F::one() * m0 * g)) } fn parameters(&self) -> Vec { @@ -98,9 +98,9 @@ impl Node for Flatte { .map(|event| { let res_mass = self.decay.resonance_p4(event).m(); let br_mom_channel_1 = - utils::breakup_momentum_c(F::fsqrt(res_mass), self.m1s[0], self.m1s[1]); + utils::breakup_momentum_c(F::sqrt(res_mass), self.m1s[0], self.m1s[1]); let br_mom_channel_2 = - utils::breakup_momentum_c(F::fsqrt(res_mass), self.m2s[0], self.m2s[1]); + utils::breakup_momentum_c(F::sqrt(res_mass), self.m2s[0], self.m2s[1]); (res_mass, [br_mom_channel_1, br_mom_channel_2]) }) .collect(); @@ -114,16 +114,15 @@ impl Node for Flatte { fn calculate(&self, parameters: &[F], event: &Event) -> Result, RustitudeError> { let (res_mass, br_momenta) = self.data[event.index]; let gammas = [ - parameters[1].c() * br_momenta[0], - parameters[2].c() * br_momenta[1], + Complex::from(parameters[1]) * br_momenta[0], + Complex::from(parameters[2]) * br_momenta[1], ]; let gamma_low = gammas[self.low_channel]; let gamma_j = gammas[self.channel]; - let mass = parameters[0].c(); - Ok( - (mass * Complex::sqrt(gamma_low * gamma_j)) / (mass.powi(2) - res_mass.fpowi(2).c()) - - F::I * mass * (gammas[0] * gammas[1]), - ) + let mass = Complex::from(parameters[0]); + Ok((mass * Complex::sqrt(gamma_low * gamma_j)) + / (mass.powi(2) - Complex::from(res_mass.powi(2))) + - Complex::::i() * mass * (gammas[0] * gammas[1])) } } @@ -143,7 +142,7 @@ pub struct KMatrixConstants { l: usize, } -impl KMatrixConstants { +impl KMatrixConstants { fn c_matrix(&self, s: F) -> SMatrix, C, C> { SMatrix::from_diagonal(&SVector::from_fn(|i, _| { utils::rho(s, self.m1s[i], self.m2s[i]) / F::PI() @@ -154,11 +153,11 @@ impl KMatrixConstants { .ln() - utils::chi_plus(s, self.m1s[i], self.m2s[i]) / F::PI() * ((self.m2s[i] - self.m1s[i]) / (self.m1s[i] + self.m2s[i])) - * (self.m2s[i] / self.m1s[i]).fln() + * F::ln(self.m2s[i] / self.m1s[i]) })) } fn barrier_factor(s: F, m1: F, m2: F, mr: F, l: usize) -> F { - utils::blatt_weisskopf(s.fsqrt(), m1, m2, l) / utils::blatt_weisskopf(mr, m1, m2, l) + utils::blatt_weisskopf(F::sqrt(s), m1, m2, l) / utils::blatt_weisskopf(mr, m1, m2, l) } fn barrier_matrix(&self, s: F) -> SMatrix { SMatrix::from_fn(|i, a| { @@ -171,17 +170,17 @@ impl KMatrixConstants { SMatrix::from_fn(|i, j| { (0..R) .map(|a| { - (bf[(i, a)] - * bf[(j, a)] - * (self.g[(i, a)] * self.g[(j, a)] - + (self.c[(i, j)]) * (self.mrs[a].fpowi(2) - s))) - .c() - * self.pole_product_remainder(s, a) + Complex::from( + bf[(i, a)] + * bf[(j, a)] + * (self.g[(i, a)] * self.g[(j, a)] + + (self.c[(i, j)]) * (self.mrs[a].powi(2) - s)), + ) * self.pole_product_remainder(s, a) }) .sum::>() * self .adler_zero - .map_or(F::ONE, |az| (s - az.s_0) / az.s_norm) + .map_or(F::one(), |az| (s - az.s_0) / az.s_norm) }) } @@ -189,7 +188,7 @@ impl KMatrixConstants { (0..R) .filter_map(|a| { if a != a_i { - Some(self.mrs[a].fpowi(2) - s) + Some(self.mrs[a].powi(2) - s) } else { None } @@ -197,16 +196,7 @@ impl KMatrixConstants { .product() } fn pole_product(&self, s: F) -> F { - (0..R).map(|a| (self.mrs[a].fpowi(2) - s)).product() - } - - fn ikc_inv(&self, s: F, channel: usize) -> SVector, C> { - let i_mat = SMatrix::, C, C>::identity().scale(self.pole_product(s)); - let k_mat = self.k_matrix(s); - let c_mat = self.c_matrix(s); - let ikc_mat = i_mat + k_mat * c_mat; - let ikc_inv_mat = ikc_mat.try_inverse().unwrap(); - ikc_inv_mat.row(channel).transpose() + (0..R).map(|a| (self.mrs[a].powi(2) - s)).product() } fn p_vector( @@ -226,6 +216,17 @@ impl KMatrixConstants { ikc_inv_vec.dot(&Self::p_vector(betas, pvector_constants_mat)) } } +impl KMatrixConstants { + fn ikc_inv(&self, s: F, channel: usize) -> SVector, C> { + let i_mat = SMatrix::, C, C>::identity().scale(self.pole_product(s)); + let k_mat = self.k_matrix(s); + let c_mat = self.c_matrix(s); + let ikc_mat = i_mat + k_mat * c_mat; + let ikc_inv_mat = ikc_mat.try_inverse().unwrap(); + ikc_inv_mat.row(channel).transpose() + } +} + #[derive(Clone)] #[allow(clippy::type_complexity)] pub struct KMatrixF0 { @@ -235,32 +236,32 @@ pub struct KMatrixF0 { data: Vec<(SVector, 5>, SMatrix, 5, 5>)>, } #[rustfmt::skip] -impl KMatrixF0 { +impl KMatrixF0 { pub fn new(channel: usize, decay: Decay) -> Self { Self { channel, decay, constants: KMatrixConstants { - g: SMatrix::::from_vec(F::fv(vec![ + g: SMatrix::::from_vec(convert_vec!(vec![ 0.74987, -0.01257, 0.27536, -0.15102, 0.36103, 0.06401, 0.00204, 0.77413, 0.50999, 0.13112, -0.23417, -0.01032, 0.72283, 0.11934, 0.36792, 0.01270, 0.26700, 0.09214, 0.02742, -0.04025, -0.14242, 0.22780, 0.15981, 0.16272, -0.17397, - ])), - c: SMatrix::::from_vec(F::fv(vec![ + ], F)), + c: SMatrix::::from_vec(convert_vec!(vec![ 0.03728, 0.00000, -0.01398, -0.02203, 0.01397, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, -0.01398, 0.00000, 0.02349, 0.03101, -0.04003, -0.02203, 0.00000, 0.03101, -0.13769, -0.06722, 0.01397, 0.00000, -0.04003, -0.06722, -0.28401, - ])), - m1s: F::fa([0.1349768, 2.0 * 0.1349768, 0.493677, 0.547862, 0.547862]), - m2s: F::fa([0.1349768, 2.0 * 0.1349768, 0.497611, 0.547862, 0.95778]), - mrs: F::fa([0.51461, 0.90630, 1.23089, 1.46104, 1.69611]), + ], F)), + m1s: convert_array!([0.1349768, 2.0 * 0.1349768, 0.493677, 0.547862, 0.547862], F), + m2s: convert_array!([0.1349768, 2.0 * 0.1349768, 0.497611, 0.547862, 0.95778], F), + mrs: convert_array!([0.51461, 0.90630, 1.23089, 1.46104, 1.69611], F), adler_zero: Some(AdlerZero { - s_0: F::f(0.0091125), - s_norm: F::ONE, + s_0: convert!(0.0091125, F), + s_norm: F::one(), }), l: 0, }, @@ -269,7 +270,7 @@ impl KMatrixF0 { } } -impl Node for KMatrixF0 { +impl Node for KMatrixF0 { fn precalculate(&mut self, dataset: &Dataset) -> Result<(), RustitudeError> { self.data = dataset .events @@ -278,7 +279,7 @@ impl Node for KMatrixF0 { let s = self.decay.resonance_p4(event).m2(); let barrier_mat = self.constants.barrier_matrix(s); let pvector_constants = SMatrix::, 5, 5>::from_fn(|i, a| { - barrier_mat[(i, a)].c() + Complex::from(barrier_mat[(i, a)]) * self.constants.g[(i, a)] * self.constants.pole_product_remainder(s, a) }); @@ -326,27 +327,27 @@ pub struct KMatrixF2 { data: Vec<(SVector, 4>, SMatrix, 4, 4>)>, } #[rustfmt::skip] -impl KMatrixF2 { +impl KMatrixF2 { pub fn new(channel: usize, decay: Decay) -> Self { Self { channel, decay, constants: KMatrixConstants { - g: SMatrix::::from_vec(F::fv(vec![ + g: SMatrix::::from_vec(convert_vec!(vec![ 0.40033, 0.15479, -0.08900, -0.00113, 0.01820, 0.17300, 0.32393, 0.15256, -0.06709, 0.22941, -0.43133, 0.23721, -0.49924, 0.19295, 0.27975, -0.03987, - ])), - c: SMatrix::::from_vec(F::fv(vec![ + ], F)), + c: SMatrix::::from_vec(convert_vec!(vec![ -0.04319, 0.00000, 0.00984, 0.01028, 0.00000, 0.00000, 0.00000, 0.00000, 0.00984, 0.00000, -0.07344, 0.05533, 0.01028, 0.00000, 0.05533, -0.05183, - ])), - m1s: F::fa([0.1349768, 2.0 * 0.1349768, 0.493677, 0.547862]), - m2s: F::fa([0.1349768, 2.0 * 0.1349768, 0.497611, 0.547862]), - mrs: F::fa([1.15299, 1.48359, 1.72923, 1.96700]), + ], F)), + m1s: convert_array!([0.1349768, 2.0 * 0.1349768, 0.493677, 0.547862], F), + m2s: convert_array!([0.1349768, 2.0 * 0.1349768, 0.497611, 0.547862], F), + mrs: convert_array!([1.15299, 1.48359, 1.72923, 1.96700], F), adler_zero: None, l: 2, }, @@ -355,7 +356,7 @@ impl KMatrixF2 { } } -impl Node for KMatrixF2 { +impl Node for KMatrixF2 { fn precalculate(&mut self, dataset: &Dataset) -> Result<(), RustitudeError> { self.data = dataset .events @@ -364,7 +365,7 @@ impl Node for KMatrixF2 { let s = self.decay.resonance_p4(event).m2(); let barrier_mat = self.constants.barrier_matrix(s); let pvector_constants = SMatrix::, 4, 4>::from_fn(|i, a| { - barrier_mat[(i, a)].c() + Complex::from(barrier_mat[(i, a)]) * self.constants.g[(i, a)] * self.constants.pole_product_remainder(s, a) }); @@ -410,23 +411,23 @@ pub struct KMatrixA0 { data: Vec<(SVector, 2>, SMatrix, 2, 2>)>, } #[rustfmt::skip] -impl KMatrixA0 { +impl KMatrixA0 { pub fn new(channel: usize, decay: Decay) -> Self { Self { channel, decay, constants: KMatrixConstants { - g: SMatrix::::from_vec(F::fv(vec![ + g: SMatrix::::from_vec(convert_vec!(vec![ 0.43215, -0.28825, 0.19000, 0.43372 - ])), - c: SMatrix::::from_vec(F::fv(vec![ + ], F)), + c: SMatrix::::from_vec(convert_vec!(vec![ 0.00000, 0.00000, 0.00000, 0.00000 - ])), - m1s: F::fa([0.1349768, 0.493677]), - m2s: F::fa([0.547862, 0.497611]), - mrs: F::fa([0.95395, 1.26767]), + ], F)), + m1s: convert_array!([0.1349768, 0.493677], F), + m2s: convert_array!([0.547862, 0.497611], F), + mrs: convert_array!([0.95395, 1.26767], F), adler_zero: None, l: 0, }, @@ -435,7 +436,7 @@ impl KMatrixA0 { } } -impl Node for KMatrixA0 { +impl Node for KMatrixA0 { fn precalculate(&mut self, dataset: &Dataset) -> Result<(), RustitudeError> { self.data = dataset .events @@ -444,7 +445,7 @@ impl Node for KMatrixA0 { let s = self.decay.resonance_p4(event).m2(); let barrier_mat = self.constants.barrier_matrix(s); let pvector_constants = SMatrix::, 2, 2>::from_fn(|i, a| { - barrier_mat[(i, a)].c() + Complex::from(barrier_mat[(i, a)]) * self.constants.g[(i, a)] * self.constants.pole_product_remainder(s, a) }); @@ -484,24 +485,24 @@ pub struct KMatrixA2 { data: Vec<(SVector, 3>, SMatrix, 3, 2>)>, } #[rustfmt::skip] -impl KMatrixA2 { +impl KMatrixA2 { pub fn new(channel: usize, decay: Decay) -> Self { Self { channel, decay, constants: KMatrixConstants { - g: SMatrix::::from_vec(F::fv(vec![ + g: SMatrix::::from_vec(convert_vec!(vec![ 0.30073, 0.21426, -0.09162, 0.68567, 0.12543, 0.00184 - ])), - c: SMatrix::::from_vec(F::fv(vec![ + ], F)), + c: SMatrix::::from_vec(convert_vec!(vec![ -0.40184, 0.00033, -0.08707, 0.00033, -0.21416, -0.06193, -0.08707, -0.06193, -0.17435, - ])), - m1s: F::fa([0.1349768, 0.493677, 0.1349768]), - m2s: F::fa([0.547862, 0.497611, 0.95778]), - mrs: F::fa([1.30080, 1.75351]), + ], F)), + m1s: convert_array!([0.1349768, 0.493677, 0.1349768], F), + m2s: convert_array!([0.547862, 0.497611, 0.95778], F), + mrs: convert_array!([1.30080, 1.75351], F), adler_zero: None, l: 2, }, @@ -510,7 +511,7 @@ impl KMatrixA2 { } } -impl Node for KMatrixA2 { +impl Node for KMatrixA2 { fn precalculate(&mut self, dataset: &Dataset) -> Result<(), RustitudeError> { self.data = dataset .events @@ -519,7 +520,7 @@ impl Node for KMatrixA2 { let s = self.decay.resonance_p4(event).m2(); let barrier_mat = self.constants.barrier_matrix(s); let pvector_constants = SMatrix::, 3, 2>::from_fn(|i, a| { - barrier_mat[(i, a)].c() + Complex::from(barrier_mat[(i, a)]) * self.constants.g[(i, a)] * self.constants.pole_product_remainder(s, a) }); @@ -559,24 +560,24 @@ pub struct KMatrixRho { data: Vec<(SVector, 3>, SMatrix, 3, 2>)>, } #[rustfmt::skip] -impl KMatrixRho { +impl KMatrixRho { pub fn new(channel: usize, decay: Decay) -> Self { Self { channel, decay, constants: KMatrixConstants { - g: SMatrix::::from_vec(F::fv(vec![ + g: SMatrix::::from_vec(convert_vec!(vec![ 0.28023, 0.01806, 0.06501, 0.16318, 0.53879, 0.00495, - ])), - c: SMatrix::::from_vec(F::fv(vec![ + ], F)), + c: SMatrix::::from_vec(convert_vec!(vec![ -0.06948, 0.00000, 0.07958, 0.00000, 0.00000, 0.00000, 0.07958, 0.00000, -0.60000, - ])), - m1s: F::fa([0.1349768, 2.0 * 0.1349768, 0.493677]), - m2s: F::fa([0.1349768, 2.0 * 0.1349768, 0.497611]), - mrs: F::fa([0.71093, 1.58660]), + ], F)), + m1s: convert_array!([0.1349768, 2.0 * 0.1349768, 0.493677], F), + m2s: convert_array!([0.1349768, 2.0 * 0.1349768, 0.497611], F), + mrs: convert_array!([0.71093, 1.58660], F), adler_zero: None, l: 1, }, @@ -585,7 +586,7 @@ impl KMatrixRho { } } -impl Node for KMatrixRho { +impl Node for KMatrixRho { fn precalculate(&mut self, dataset: &Dataset) -> Result<(), RustitudeError> { self.data = dataset .events @@ -594,7 +595,7 @@ impl Node for KMatrixRho { let s = self.decay.resonance_p4(event).m2(); let barrier_mat = self.constants.barrier_matrix(s); let pvector_constants = SMatrix::, 3, 2>::from_fn(|i, a| { - barrier_mat[(i, a)].c() + Complex::from(barrier_mat[(i, a)]) * self.constants.g[(i, a)] * self.constants.pole_product_remainder(s, a) }); @@ -634,22 +635,22 @@ pub struct KMatrixPi1 { data: Vec<(SVector, 2>, SMatrix, 2, 1>)>, } #[rustfmt::skip] -impl KMatrixPi1 { +impl KMatrixPi1 { pub fn new(channel: usize, decay: Decay) -> Self { Self { channel, decay, constants: KMatrixConstants { - g: SMatrix::::from_vec(F::fv(vec![ + g: SMatrix::::from_vec(convert_vec!(vec![ 0.80564, 1.04595 - ])), - c: SMatrix::::from_vec(F::fv(vec![ + ], F)), + c: SMatrix::::from_vec(convert_vec!(vec![ 1.05000, 0.15163, 0.15163, -0.24611, - ])), - m1s: F::fa([0.1349768, 0.1349768]), - m2s: F::fa([0.547862, 0.95778]), - mrs: F::fa([1.38552]), + ], F)), + m1s: convert_array!([0.1349768, 0.1349768], F), + m2s: convert_array!([0.547862, 0.95778], F), + mrs: convert_array!([1.38552], F), adler_zero: None, l: 1, }, @@ -658,7 +659,7 @@ impl KMatrixPi1 { } } -impl Node for KMatrixPi1 { +impl Node for KMatrixPi1 { fn precalculate(&mut self, dataset: &Dataset) -> Result<(), RustitudeError> { self.data = dataset .events @@ -667,7 +668,7 @@ impl Node for KMatrixPi1 { let s = self.decay.resonance_p4(event).m2(); let barrier_mat = self.constants.barrier_matrix(s); let pvector_constants = SMatrix::, 2, 1>::from_fn(|i, a| { - barrier_mat[(i, a)].c() + Complex::from(barrier_mat[(i, a)]) * self.constants.g[(i, a)] * self.constants.pole_product_remainder(s, a) }); diff --git a/crates/rustitude-gluex/src/sdmes.rs b/crates/rustitude-gluex/src/sdmes.rs index 5b659ca..1cd9668 100644 --- a/crates/rustitude-gluex/src/sdmes.rs +++ b/crates/rustitude-gluex/src/sdmes.rs @@ -32,18 +32,15 @@ impl Node for TwoPiSDME { .par_iter() .map(|event| { let (_, y, _, p) = self.decay.coordinates(self.frame, 0, event); - let big_phi = y.dot(&event.eps).fatan2( - event - .beam_p4 - .momentum() - .normalize() - .dot(&event.eps.cross(&y)), + let big_phi = F::atan2( + y.dot(&event.eps), + event.beam_p4.direction().dot(&event.eps.cross(&y)), ); - let pgamma = event.eps.norm(); + let pgamma = event.eps_mag(); ( - p.theta_cos().fpowi(2), - F::fsin(p.theta()).fpowi(2), - F::fsin(F::TWO * p.theta()), + p.theta_cos().powi(2), + F::sin(p.theta()).powi(2), + F::sin(convert!(2, F) * p.theta()), p.phi(), big_phi, pgamma, @@ -65,23 +62,22 @@ impl Node for TwoPiSDME { let rho_102 = parameters[7]; let rho_1n12 = parameters[8]; - Ok(F::fsqrt(F::fabs( - (F::THREE / (F::FOUR * F::PI())) - * (F::f(0.5) * (F::ONE - rho_000) - + F::f(0.5) * (F::THREE * rho_000 - F::ONE) * cossqtheta - - F::SQRT_2() * rho_100 * sin2theta * F::fcos(phi) - - rho_1n10 * sinsqtheta * F::fcos(F::TWO * phi)) + Ok(Complex::from(F::sqrt(F::abs( + (convert!(3, F) / (convert!(4, F) * F::PI())) + * (convert!(0.5, F) * (F::one() - rho_000) + + convert!(0.5, F) * (convert!(3, F) * rho_000 - F::one()) * cossqtheta + - F::SQRT_2() * rho_100 * sin2theta * F::cos(phi) + - rho_1n10 * sinsqtheta * F::cos(convert!(2, F) * phi)) - pgamma - * F::fcos(F::TWO * big_phi) + * F::cos(convert!(2, F) * big_phi) * (rho_111 * sinsqtheta + rho_001 * cossqtheta - - F::SQRT_2() * rho_101 * sin2theta * F::fcos(phi) - - rho_1n11 * sinsqtheta * F::fcos(F::TWO * phi)) + - F::SQRT_2() * rho_101 * sin2theta * F::cos(phi) + - rho_1n11 * sinsqtheta * F::cos(convert!(2, F) * phi)) - pgamma - * F::fsin(F::TWO * big_phi) - * (F::SQRT_2() * rho_102 * sin2theta * F::fsin(phi) - + rho_1n12 * sinsqtheta * F::fsin(F::TWO * phi)), - )) - .c()) + * F::sin(convert!(2, F) * big_phi) + * (F::SQRT_2() * rho_102 * sin2theta * F::sin(phi) + + rho_1n12 * sinsqtheta * F::sin(convert!(2, F) * phi)), + )))) } fn parameters(&self) -> Vec { @@ -129,26 +125,19 @@ impl Node for ThreePiSDME { let res_p4 = self.decay.resonance_p4(event); let p1_res_p4 = self.decay.primary_p4(event).boost_along(&res_p4); let p2_res_p4 = self.decay.primary_p4(event).boost_along(&res_p4); - let norm = p1_res_p4 - .momentum() - .cross(&p2_res_p4.momentum()) - .normalize(); + let norm = p1_res_p4.momentum().cross(&p2_res_p4.momentum()).unit(); let (_, y, _, p) = self .frame .coordinates_from_boosted_vec(self.decay, &norm, event); - let big_phi = F::fatan2( + let big_phi = F::atan2( y.dot(&event.eps), - event - .beam_p4 - .momentum() - .normalize() - .dot(&event.eps.cross(&y)), + event.beam_p4.direction().dot(&event.eps.cross(&y)), ); - let pgamma = event.eps.norm(); + let pgamma = event.eps_mag(); ( - p.theta_cos().fpowi(2), - F::fsin(p.theta()).fpowi(2), - F::fsin(F::TWO * p.theta()), + p.theta_cos().powi(2), + F::sin(p.theta()).powi(2), + F::sin(convert!(2, F) * p.theta()), p.phi(), big_phi, pgamma, @@ -170,23 +159,22 @@ impl Node for ThreePiSDME { let rho_102 = parameters[7]; let rho_1n12 = parameters[8]; - Ok(F::fsqrt(F::fabs( - (F::THREE / (F::FOUR * F::PI())) - * (F::f(0.5) * (F::ONE - rho_000) - + F::f(0.5) * (F::THREE * rho_000 - F::ONE) * cossqtheta - - F::SQRT_2() * rho_100 * sin2theta * F::fcos(phi) - - rho_1n10 * sinsqtheta * F::fcos(F::TWO * phi)) + Ok(Complex::from(F::sqrt(F::abs( + (convert!(3, F) / (convert!(4, F) * F::PI())) + * (convert!(0.5, F) * (F::one() - rho_000) + + convert!(0.5, F) * (convert!(3, F) * rho_000 - F::one()) * cossqtheta + - F::SQRT_2() * rho_100 * sin2theta * F::cos(phi) + - rho_1n10 * sinsqtheta * F::cos(convert!(2, F) * phi)) - pgamma - * F::fcos(F::TWO * big_phi) + * F::cos(convert!(2, F) * big_phi) * (rho_111 * sinsqtheta + rho_001 * cossqtheta - - F::SQRT_2() * rho_101 * sin2theta * F::fcos(phi) - - rho_1n11 * sinsqtheta * F::fcos(F::TWO * phi)) + - F::SQRT_2() * rho_101 * sin2theta * F::cos(phi) + - rho_1n11 * sinsqtheta * F::cos(convert!(2, F) * phi)) - pgamma - * F::fsin(F::TWO * big_phi) - * (F::SQRT_2() * rho_102 * sin2theta * F::fsin(phi) - + rho_1n12 * sinsqtheta * F::fsin(F::TWO * phi)), - )) - .c()) + * F::sin(convert!(2, F) * big_phi) + * (F::SQRT_2() * rho_102 * sin2theta * F::sin(phi) + + rho_1n12 * sinsqtheta * F::sin(convert!(2, F) * phi)), + )))) } fn parameters(&self) -> Vec { @@ -232,18 +220,15 @@ impl Node for VecRadiativeSDME { .par_iter() .map(|event| { let (_, y, _, p) = self.decay.coordinates(self.frame, 0, event); - let big_phi = y.dot(&event.eps).fatan2( - event - .beam_p4 - .momentum() - .normalize() - .dot(&event.eps.cross(&y)), + let big_phi = F::atan2( + y.dot(&event.eps), + event.beam_p4.direction().dot(&event.eps.cross(&y)), ); - let pgamma = event.eps.norm(); + let pgamma = event.eps_mag(); ( - p.theta_cos().fpowi(2), - F::fsin(p.theta()).fpowi(2), - F::fsin(F::TWO * p.theta()), + p.theta_cos().powi(2), + F::sin(p.theta()).powi(2), + F::sin(convert!(2, F) * p.theta()), p.phi(), big_phi, pgamma, @@ -265,23 +250,24 @@ impl Node for VecRadiativeSDME { let rho_102 = parameters[7]; let rho_1n12 = parameters[8]; - Ok(F::fsqrt(F::fabs( - (F::THREE / (F::EIGHT * F::PI())) - * (F::ONE - sinsqtheta * F::f(0.5) * (F::ONE - rho_000) - rho_000 * cossqtheta - + rho_1n10 * sinsqtheta * F::fcos(F::TWO * phi) - + F::SQRT_2() * rho_100 * sin2theta * F::fcos(phi) + Ok(Complex::from(F::sqrt(F::abs( + (convert!(3, F) / (convert!(8, F) * F::PI())) + * (F::one() + - sinsqtheta * convert!(0.5, F) * (F::one() - rho_000) + - rho_000 * cossqtheta + + rho_1n10 * sinsqtheta * F::cos(convert!(2, F) * phi) + + F::SQRT_2() * rho_100 * sin2theta * F::cos(phi) - pgamma - * F::fcos(F::TWO * big_phi) - * (F::TWO * rho_111 + * F::cos(convert!(2, F) * big_phi) + * (convert!(2, F) * rho_111 + (rho_001 - rho_111) * sinsqtheta - + rho_1n11 * sinsqtheta * F::fcos(F::TWO * phi) - + F::SQRT_2() * rho_101 * sin2theta * F::fcos(phi)) + + rho_1n11 * sinsqtheta * F::cos(convert!(2, F) * phi) + + F::SQRT_2() * rho_101 * sin2theta * F::cos(phi)) + pgamma - * F::fsin(F::TWO * big_phi) - * (rho_1n12 * sinsqtheta * F::fsin(F::TWO * phi) - + F::SQRT_2() * rho_102 * sin2theta * F::fsin(phi))), - )) - .c()) + * F::sin(convert!(2, F) * big_phi) + * (rho_1n12 * sinsqtheta * F::sin(convert!(2, F) * phi) + + F::SQRT_2() * rho_102 * sin2theta * F::sin(phi))), + )))) } fn parameters(&self) -> Vec { diff --git a/crates/rustitude-gluex/src/utils.rs b/crates/rustitude-gluex/src/utils.rs index be7ea27..074f940 100644 --- a/crates/rustitude-gluex/src/utils.rs +++ b/crates/rustitude-gluex/src/utils.rs @@ -1,32 +1,30 @@ use std::{fmt::Display, num::ParseIntError, str::FromStr}; use factorial::Factorial; -use rustitude_core::prelude::*; +use rustitude_core::{convert, prelude::*}; use sphrs::Coordinates; use thiserror::Error; pub fn breakup_momentum(m0: F, m1: F, m2: F) -> F { - F::fsqrt(F::fabs( - m0.fpowi(4) + m1.fpowi(4) + m2.fpowi(4) - - F::TWO - * (m0.fpowi(2) * m1.fpowi(2) - + m0.fpowi(2) * m2.fpowi(2) - + m1.fpowi(2) * m2.fpowi(2)), - )) / (F::TWO * m0) + F::sqrt(F::abs( + m0.powi(4) + m1.powi(4) + m2.powi(4) + - convert!(2, F) + * (m0.powi(2) * m1.powi(2) + m0.powi(2) * m2.powi(2) + m1.powi(2) * m2.powi(2)), + )) / (convert!(2, F) * m0) } /// Computes the ([`Complex`]) breakup momentum of a particle with mass `m0` decaying into two particles /// with masses `m1` and `m2`. pub fn breakup_momentum_c(m0: F, m1: F, m2: F) -> Complex { - rho(m0.fpowi(2), m1, m2) * m0 / F::TWO + rho(m0.powi(2), m1, m2) * m0 / convert!(2, F) } pub fn chi_plus(s: F, m1: F, m2: F) -> Complex { - (F::ONE - ((m1 + m2) * (m1 + m2)) / s).c() + Complex::from(F::one() - ((m1 + m2) * (m1 + m2)) / s) } pub fn chi_minus(s: F, m1: F, m2: F) -> Complex { - (F::ONE - ((m1 - m2) * (m1 - m2)) / s).c() + Complex::from(F::one() - ((m1 - m2) * (m1 - m2)) / s) } pub fn rho(s: F, m1: F, m2: F) -> Complex { @@ -35,18 +33,22 @@ pub fn rho(s: F, m1: F, m2: F) -> Complex { pub fn blatt_weisskopf(m0: F, m1: F, m2: F, l: usize) -> F { let q = breakup_momentum(m0, m1, m2); - let z = q.fpowi(2) / F::f(0.1973).fpowi(2); + let z = q.powi(2) / convert!(0.1973, F).powi(2); match l { - 0 => F::ONE, - 1 => F::fsqrt((F::TWO * z) / (z + F::ONE)), - 2 => F::fsqrt((F::f(13.0) * z.fpowi(2)) / ((z - F::THREE).fpowi(2) + F::NINE * z)), - 3 => F::fsqrt( - (F::f(277.0) * z.fpowi(3)) - / (z * (z - F::f(15.0)).fpowi(2) + F::NINE * (F::TWO * z - F::FIVE).fpowi(2)), + 0 => F::one(), + 1 => F::sqrt((convert!(2, F) * z) / (z + F::one())), + 2 => F::sqrt( + (convert!(13.0, F) * z.powi(2)) / ((z - convert!(3, F)).powi(2) + convert!(9, F) * z), ), - 4 => F::fsqrt( - (F::f(12746.0) * z.fpowi(4)) / (z.fpowi(2) - F::f(45.0) * z + F::f(105.0)).fpowi(2) - + F::f(25.0) * z * (F::TWO * z - F::f(21.0)).fpowi(2), + 3 => F::sqrt( + (convert!(277.0, F) * z.powi(3)) + / (z * (z - convert!(15.0, F)).powi(2) + + convert!(9, F) * (convert!(2, F) * z - convert!(5, F)).powi(2)), + ), + 4 => F::sqrt( + (convert!(12746.0, F) * z.powi(4)) + / (z.powi(2) - convert!(45.0, F) * z + convert!(105.0, F)).powi(2) + + convert!(25.0, F) * z * (convert!(2, F) * z - convert!(21.0, F)).powi(2), ), l => panic!("L = {l} is not yet implemented"), } @@ -59,19 +61,23 @@ pub fn blatt_weisskopf(m0: F, m1: F, m2: F, l: usize) -> F { /// `m2`, the absolute value of this function can be safely assumed to be equal to its value. pub fn blatt_weisskopf_c(m0: F, m1: F, m2: F, l: usize) -> Complex { let q = breakup_momentum_c(m0, m1, m2); - let z = q.powi(2) / F::f(0.1973).fpowi(2); + let z = q.powi(2) / convert!(0.1973, F).powi(2); match l { - 0 => F::ONE.c(), - 1 => Complex::sqrt((F::TWO.c() * z) / (z + F::ONE)), - 2 => Complex::sqrt((z.powi(2) * F::f(13.0)) / ((z - F::THREE).powi(2) + z * F::NINE)), + 0 => Complex::from(F::one()), + 1 => Complex::sqrt((Complex::from(convert!(2, F)) * z) / (z + F::one())), + 2 => Complex::sqrt( + (z.powi(2) * convert!(13.0, F)) / ((z - convert!(3, F)).powi(2) + z * convert!(9, F)), + ), 3 => Complex::sqrt( - (z.powi(3) * F::f(277.0)) - / (z * (z - F::f(15.0)).powi(2) + (z * F::TWO - F::FIVE).powi(2)) - * F::NINE, + (z.powi(3) * convert!(277.0, F)) + / (z * (z - convert!(15.0, F)).powi(2) + + (z * convert!(2, F) - convert!(5, F)).powi(2)) + * convert!(9, F), ), 4 => Complex::sqrt( - (z.powi(4) * F::f(12746.0)) / (z.powi(2) - z * F::f(45.0) + F::f(105.0)).powi(2) - + z * F::f(25.0) * (z * F::TWO - F::f(21.0)).powi(2), + (z.powi(4) * convert!(12746.0, F)) + / (z.powi(2) - z * convert!(45.0, F) + convert!(105.0, F)).powi(2) + + z * convert!(25.0, F) * (z * convert!(2, F) - convert!(21.0, F)).powi(2), ), l => panic!("L = {l} is not yet implemented"), } @@ -82,22 +88,24 @@ pub fn small_wigner_d_matrix(beta: F, j: usize, m: isize, n: isize) -> let jmm = (j as i32 - m as i32) as u32; let jpn = (j as i32 + n as i32) as u32; let jmn = (j as i32 - n as i32) as u32; - let prefactor = F::fsqrt(F::convert_u32( + let prefactor = F::sqrt(convert!( jpm.factorial() * jmm.factorial() * jpn.factorial() * jmn.factorial(), + F )); let s_min = isize::max(0, n - m) as usize; let s_max = isize::min(jpn as isize, jmm as isize) as usize; let sum: F = (s_min..=s_max) .map(|s| { - (F::fpowi(-F::ONE, m as i32 - n as i32 + s as i32) - * (F::fcos(beta / F::TWO) - .fpowi(2 * (j as i32) + n as i32 - m as i32 - 2 * (s as i32))) - * (F::fsin(beta / F::TWO).fpowi(m as i32 - n as i32 + 2 * s as i32))) - / F::convert_u32( + (F::powi(-F::one(), m as i32 - n as i32 + s as i32) + * (F::cos(beta / convert!(2, F)) + .powi(2 * (j as i32) + n as i32 - m as i32 - 2 * (s as i32))) + * (F::sin(beta / convert!(2, F)).powi(m as i32 - n as i32 + 2 * s as i32))) + / convert!( (jpm - s as u32).factorial() * (s as u32).factorial() * ((m - n + s as isize) as u32).factorial() * (jmm - s as u32).factorial(), + F ) }) .sum(); @@ -112,9 +120,9 @@ pub fn wigner_d_matrix( m: isize, n: isize, ) -> Complex { - Complex::cis(-(F::convert_isize(m)) * alpha) + Complex::cis(convert!(-m, F) * alpha) * small_wigner_d_matrix(beta, j, m, n) - * Complex::cis(-(F::convert_isize(n)) * gamma) + * Complex::cis(convert!(-n, F) * gamma) } #[derive(Clone, Copy, Default, PartialEq)] @@ -223,7 +231,7 @@ impl FromStr for Frame { } } -pub fn coordinates( +pub fn coordinates( x: &Vector3, y: &Vector3, z: &Vector3, @@ -245,18 +253,14 @@ impl Frame { let other_res_vec = other_p4.boost_along(&resonance_p4).momentum(); let (x, y, z) = match self { Frame::Helicity => { - let z = -recoil_res_vec.normalize(); - let y = beam_res_vec.cross(&z).normalize(); + let z = -recoil_res_vec.unit(); + let y = beam_res_vec.cross(&z).unit(); let x = y.cross(&z); (x, y, z) } Frame::GottfriedJackson => { - let z = beam_res_vec.normalize(); - let y = event - .beam_p4 - .momentum() - .cross(&(-recoil_res_vec)) - .normalize(); + let z = beam_res_vec.unit(); + let y = event.beam_p4.momentum().cross(&(-recoil_res_vec)).unit(); let x = y.cross(&z); (x, y, z) } @@ -274,18 +278,14 @@ impl Frame { let recoil_res_vec = event.recoil_p4.boost_along(&resonance_p4).momentum(); let (x, y, z) = match self { Frame::Helicity => { - let z = -recoil_res_vec.normalize(); - let y = beam_res_vec.cross(&z).normalize(); + let z = -recoil_res_vec.unit(); + let y = beam_res_vec.cross(&z).unit(); let x = y.cross(&z); (x, y, z) } Frame::GottfriedJackson => { - let z = beam_res_vec.normalize(); - let y = event - .beam_p4 - .momentum() - .cross(&(-recoil_res_vec)) - .normalize(); + let z = beam_res_vec.unit(); + let y = event.beam_p4.momentum().cross(&(-recoil_res_vec)).unit(); let x = y.cross(&z); (x, y, z) }