rustitude_gluex/
polarization.rs

1use crate::utils::{self, Decay, Frame, Sign};
2use rayon::prelude::*;
3use rustitude_core::{convert, prelude::*};
4use sphrs::{ComplexSH, SHEval};
5
6#[derive(Clone)]
7pub struct ThreePiPolFrac<F> {
8    beam_pol: F,
9    j_resonance: u32,
10    p_resonance: F,
11    i_resonance: usize,
12    l_resonance: u32,
13    j_isobar: u32,
14    i_isobar: usize,
15    iz_daughters: [usize; 3],
16    decay_resonance: Decay,
17    decay_isobar: Decay,
18    data: Vec<Complex<F>>,
19}
20
21impl<F: Field> ThreePiPolFrac<F> {
22    #[allow(clippy::too_many_arguments)]
23    pub fn new(
24        beam_pol: Sign,
25        j_resonance: u32,
26        p_resonance: Sign,
27        i_resonance: usize,
28        l_resonance: u32,
29        j_isobar: u32,
30        i_isobar: usize,
31        iz_daughters: [usize; 3],
32        decay_resonance: Decay,
33        decay_isobar: Decay,
34    ) -> Self {
35        match (decay_resonance, decay_isobar) {
36            (Decay::ThreeBodyDecay(_), Decay::TwoBodyDecay(_)) => Self {
37                beam_pol: match beam_pol {
38                    Sign::Positive => F::one(),
39                    Sign::Negative => -F::one(),
40                },
41                j_resonance,
42                p_resonance: match p_resonance {
43                    Sign::Positive => F::one(),
44                    Sign::Negative => -F::one(),
45                },
46                i_resonance,
47                l_resonance,
48                j_isobar,
49                i_isobar,
50                iz_daughters,
51                decay_resonance,
52                decay_isobar,
53                data: Vec::default(),
54            },
55            _ => unimplemented!(),
56        }
57    }
58}
59
60impl<F: Field> Node<F> for ThreePiPolFrac<F> {
61    fn calculate(&self, parameters: &[F], event: &Event<F>) -> Result<Complex<F>, RustitudeError> {
62        Ok(self.data[event.index] * (F::one() + self.beam_pol * parameters[0]) / convert!(4, F))
63    }
64
65    fn precalculate(&mut self, dataset: &Dataset<F>) -> Result<(), RustitudeError> {
66        self.data = dataset
67            .events
68            .par_iter()
69            .map(|event| {
70                let isobar_p4 = self.decay_isobar.resonance_p4(event);
71                let resonance_p4 = self.decay_resonance.resonance_p4(event);
72                let alpha = event.recoil_p4.phi();
73                let (_, _, _, p3_res_coords) =
74                    self.decay_resonance.coordinates(Frame::Helicity, 2, event);
75                let p1_iso_p4 = self.decay_isobar.primary_p4(event).boost_along(&isobar_p4);
76                let (_, _, _, p1_iso_coords) =
77                    Frame::Helicity.coordinates(self.decay_resonance, &p1_iso_p4, event);
78                let k = utils::breakup_momentum(
79                    resonance_p4.m(),
80                    isobar_p4.m(),
81                    self.decay_resonance.tertiary_p4(event).m(),
82                );
83                let q = utils::breakup_momentum(
84                    isobar_p4.m(),
85                    self.decay_isobar.primary_p4(event).m(),
86                    self.decay_isobar.secondary_p4(event).m(),
87                );
88                let neg_res_hel_prod = Complex::new(
89                    F::cos(convert!(2, F) * alpha),
90                    F::sin(convert!(2, F) * alpha),
91                ) * self.beam_pol
92                    * self.p_resonance
93                    * convert!(self.j_resonance % 2, F);
94                let mut res = Complex::default();
95                for m_res in -(self.l_resonance as isize)..(self.l_resonance as isize) {
96                    let mut term = Complex::default();
97                    for m_iso in -(self.j_isobar as isize)..(self.j_isobar as isize) {
98                        let ylm = ComplexSH::Spherical.eval(
99                            self.j_isobar as i64,
100                            m_iso as i64,
101                            &p1_iso_coords,
102                        );
103                        let cg_neg = convert!(
104                            wigners::clebsch_gordan(
105                                self.j_isobar,
106                                self.l_resonance as i32,
107                                m_iso as u32,
108                                m_res as i32,
109                                self.j_resonance,
110                                -1,
111                            ),
112                            F
113                        );
114                        let cg_pos = convert!(
115                            wigners::clebsch_gordan(
116                                self.j_isobar,
117                                self.l_resonance as i32,
118                                m_iso as u32,
119                                m_res as i32,
120                                self.j_resonance,
121                                1,
122                            ),
123                            F
124                        );
125                        term += ylm * neg_res_hel_prod * cg_neg + cg_pos;
126                    }
127                    let ylm = ComplexSH::Spherical.eval(
128                        self.l_resonance as i64,
129                        m_res as i64,
130                        &p3_res_coords,
131                    );
132                    term *= ylm;
133                    res += term;
134                }
135                res *= convert!(
136                    wigners::clebsch_gordan(
137                        1,
138                        1,
139                        self.iz_daughters[0] as u32,
140                        self.iz_daughters[1] as i32,
141                        self.i_isobar as u32,
142                        (self.iz_daughters[0] + self.iz_daughters[1]) as i32,
143                    ) * wigners::clebsch_gordan(
144                        self.i_isobar as u32,
145                        1,
146                        (self.iz_daughters[0] + self.iz_daughters[1]) as u32,
147                        self.iz_daughters[2] as i32,
148                        self.i_resonance as u32,
149                        (self.iz_daughters[0] + self.iz_daughters[1] + self.iz_daughters[2]) as i32,
150                    ),
151                    F
152                ) * k.powi(self.l_resonance as i32)
153                    * q.powi(self.j_isobar as i32);
154                res
155            })
156            .collect();
157        Ok(())
158    }
159
160    fn parameters(&self) -> Vec<String> {
161        std::vec!["polarization fraction".to_string()]
162    }
163}
rustitude_gluex/polarization.rs

rustitude_gluex/
polarization.rs
