Skip to content

Commit

Permalink
feat: update all amplitudes and methods to be compatible with the lat…
Browse files Browse the repository at this point in the history 
…est changes in the `Field` trait
  • Loading branch information
@denehoffman
denehoffman committed Aug 17, 2024
1 parent a843c3e commit 1814640
Show file tree
Hide file tree
Showing 6 changed files with 312 additions and 312 deletions.
41 changes: 25 additions & 16 deletions crates/rustitude-gluex/src/dalitz.rs
View file
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
use rayon::prelude::*;
use rustitude_core::prelude::*;
use rustitude_core::{convert, prelude::*};

use crate::utils::Decay;

Expand Down Expand Up @@ -38,22 +38,25 @@ impl<F: Field> Node<F> for OmegaDalitz<F> {
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))
})
Expand All @@ -69,13 +72,19 @@ impl<F: Field> Node<F> for OmegaDalitz<F> {
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())
}
Expand Down
66 changes: 31 additions & 35 deletions crates/rustitude-gluex/src/harmonics.rs
View file
Original file line number Diff line number Diff line change
@@ -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};
Expand Down Expand Up @@ -68,25 +68,21 @@ impl<F: Field + num::Float> Node<F> for Zlm<F> {
.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,
),
}
})
Expand Down Expand Up @@ -122,24 +118,21 @@ impl<F: Field> Node<F> for OnePS<F> {
.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,
),
}
})
Expand Down Expand Up @@ -184,19 +177,22 @@ impl<F: Field> Node<F> for TwoPS<F> {
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 * <F as Field>::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)
})
Expand Down
74 changes: 41 additions & 33 deletions crates/rustitude-gluex/src/polarization.rs
View file
Original file line number Diff line number Diff line change
@@ -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)]
Expand Down Expand Up @@ -35,13 +35,13 @@ impl<F: Field> ThreePiPolFrac<F> {
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,
Expand All @@ -59,7 +59,7 @@ impl<F: Field> ThreePiPolFrac<F> {

impl<F: Field> Node<F> for ThreePiPolFrac<F> {
fn calculate(&self, parameters: &[F], event: &Event<F>) -> Result<Complex<F>, 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<F>) -> Result<(), RustitudeError> {
Expand All @@ -85,37 +85,44 @@ impl<F: Field> Node<F> for ThreePiPolFrac<F> {
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,
Expand All @@ -125,7 +132,7 @@ impl<F: Field> Node<F> for ThreePiPolFrac<F> {
term *= ylm;
res += term;
}
res *= F::f(
res *= convert!(
wigners::clebsch_gordan(
1,
1,
Expand All @@ -141,8 +148,9 @@ impl<F: Field> Node<F> for ThreePiPolFrac<F> {
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();
Expand Down
Loading

0 comments on commit 1814640

Please sign in to comment.