1use crate::utils;
2use crate::utils::Decay;
3
4use nalgebra::{RealField, SMatrix, SVector};
5use rayon::prelude::*;
6use rustitude_core::{convert, convert_array, convert_vec, prelude::*};
7
8#[derive(Default, Clone)]
9pub struct BreitWigner<F: Field> {
10 l: usize,
11 decay: Decay,
12 m: Vec<F>,
13 m1: Vec<F>,
14 m2: Vec<F>,
15 q: Vec<F>,
16 f: Vec<F>,
17}
18impl<F: Field> BreitWigner<F> {
19 pub fn new(l: usize, decay: Decay) -> Self {
20 Self {
21 l,
22 decay,
23 ..Default::default()
24 }
25 }
26}
27impl<F: Field> Node<F> for BreitWigner<F> {
28 fn precalculate(&mut self, dataset: &Dataset<F>) -> Result<(), RustitudeError> {
29 (self.m, (self.m1, (self.m2, (self.q, self.f)))) = dataset
30 .events
31 .par_iter()
32 .map(|event| {
33 let p1 = self.decay.primary_p4(event);
34 let p2 = self.decay.secondary_p4(event);
35 let m = (p1 + p2).m();
36 let m1 = p1.m();
37 let m2 = p2.m();
38 let q = utils::breakup_momentum(m, m1, m2);
39 let f = utils::blatt_weisskopf(m, m1, m2, self.l);
40 (m, (m1, (m2, (q, f))))
41 })
42 .unzip();
43 Ok(())
44 }
45
46 fn calculate(&self, parameters: &[F], event: &Event<F>) -> Result<Complex<F>, RustitudeError> {
47 let m = self.m[event.index];
48 let m1 = self.m1[event.index];
49 let m2 = self.m2[event.index];
50 let q = self.q[event.index];
51 let f = self.f[event.index];
52 let m0 = parameters[0];
53 let g0 = parameters[1];
54 let f0 = utils::blatt_weisskopf(m0, m1, m2, self.l);
55 let q0 = utils::breakup_momentum(m0, m1, m2);
56 let g = g0 * (m0 / m) * (q / q0) * (f.powi(2) / f0.powi(2));
57 Ok(Complex::new(f * (m0 * g0 / F::PI()), F::zero())
58 / Complex::new(m0.powi(2) - m.powi(2), -F::one() * m0 * g))
59 }
60
61 fn parameters(&self) -> Vec<String> {
62 vec!["mass".to_string(), "width".to_string()]
63 }
64}
65
66#[derive(Clone)]
67pub struct Flatte<F: Field> {
68 channel: usize,
69 m1s: [F; 2],
70 m2s: [F; 2],
71 decay: Decay,
72 low_channel: usize,
73 data: Vec<(F, [Complex<F>; 2])>,
74}
75impl<F: Field> Flatte<F> {
76 pub fn new(channel: usize, m1s: [F; 2], m2s: [F; 2], decay: Decay) -> Self {
77 let low_channel = if (m1s[0] + m1s[1]) < (m2s[0] + m2s[1]) {
78 0
79 } else {
80 1
81 };
82 Self {
83 channel,
84 m1s,
85 m2s,
86 decay,
87 low_channel,
88 data: Vec::default(),
89 }
90 }
91}
92
93impl<F: Field> Node<F> for Flatte<F> {
94 fn precalculate(&mut self, dataset: &Dataset<F>) -> Result<(), RustitudeError> {
95 self.data = dataset
96 .events
97 .par_iter()
98 .map(|event| {
99 let res_mass = self.decay.resonance_p4(event).m();
100 let br_mom_channel_1 =
101 utils::breakup_momentum_c(F::sqrt(res_mass), self.m1s[0], self.m1s[1]);
102 let br_mom_channel_2 =
103 utils::breakup_momentum_c(F::sqrt(res_mass), self.m2s[0], self.m2s[1]);
104 (res_mass, [br_mom_channel_1, br_mom_channel_2])
105 })
106 .collect();
107 Ok(())
108 }
109
110 fn parameters(&self) -> Vec<String> {
111 vec!["mass".to_string(), "g1".to_string(), "g2".to_string()]
112 }
113
114 fn calculate(&self, parameters: &[F], event: &Event<F>) -> Result<Complex<F>, RustitudeError> {
115 let (res_mass, br_momenta) = self.data[event.index];
116 let gammas = [
117 Complex::from(parameters[1]) * br_momenta[0],
118 Complex::from(parameters[2]) * br_momenta[1],
119 ];
120 let gamma_low = gammas[self.low_channel];
121 let gamma_j = gammas[self.channel];
122 let mass = Complex::from(parameters[0]);
123 Ok((mass * Complex::sqrt(gamma_low * gamma_j))
124 / (mass.powi(2) - Complex::from(res_mass.powi(2)))
125 - Complex::<F>::i() * mass * (gammas[0] * gammas[1]))
126 }
127}
128
129#[derive(Clone, Copy)]
130pub struct AdlerZero<F: Field> {
131 pub s_0: F,
132 pub s_norm: F,
133}
134#[derive(Clone)]
135pub struct KMatrixConstants<F: Field, const C: usize, const R: usize> {
136 g: SMatrix<F, C, R>,
137 c: SMatrix<F, C, C>,
138 m1s: [F; C],
139 m2s: [F; C],
140 mrs: [F; R],
141 adler_zero: Option<AdlerZero<F>>,
142 l: usize,
143}
144
145impl<F: Field + 'static, const C: usize, const R: usize> KMatrixConstants<F, C, R> {
146 fn c_matrix(&self, s: F) -> SMatrix<Complex<F>, C, C> {
147 SMatrix::from_diagonal(&SVector::from_fn(|i, _| {
148 utils::rho(s, self.m1s[i], self.m2s[i]) / F::PI()
149 * ((utils::chi_plus(s, self.m1s[i], self.m2s[i])
150 + utils::rho(s, self.m1s[i], self.m2s[i]))
151 / (utils::chi_plus(s, self.m1s[i], self.m2s[i])
152 - utils::rho(s, self.m1s[i], self.m2s[i])))
153 .ln()
154 - utils::chi_plus(s, self.m1s[i], self.m2s[i]) / F::PI()
155 * ((self.m2s[i] - self.m1s[i]) / (self.m1s[i] + self.m2s[i]))
156 * F::ln(self.m2s[i] / self.m1s[i])
157 }))
158 }
159 fn barrier_factor(s: F, m1: F, m2: F, mr: F, l: usize) -> F {
160 utils::blatt_weisskopf(F::sqrt(s), m1, m2, l) / utils::blatt_weisskopf(mr, m1, m2, l)
161 }
162 fn barrier_matrix(&self, s: F) -> SMatrix<F, C, R> {
163 SMatrix::from_fn(|i, a| {
164 Self::barrier_factor(s, self.m1s[i], self.m2s[i], self.mrs[a], self.l)
165 })
166 }
167
168 fn k_matrix(&self, s: F) -> SMatrix<Complex<F>, C, C> {
169 let bf = self.barrier_matrix(s);
170 SMatrix::from_fn(|i, j| {
171 (0..R)
172 .map(|a| {
173 Complex::from(
174 bf[(i, a)]
175 * bf[(j, a)]
176 * (self.g[(i, a)] * self.g[(j, a)]
177 + (self.c[(i, j)]) * (self.mrs[a].powi(2) - s)),
178 ) * self.pole_product_remainder(s, a)
179 })
180 .sum::<Complex<F>>()
181 * self
182 .adler_zero
183 .map_or(F::one(), |az| (s - az.s_0) / az.s_norm)
184 })
185 }
186
187 fn pole_product_remainder(&self, s: F, a_i: usize) -> F {
188 (0..R)
189 .filter_map(|a| {
190 if a != a_i {
191 Some(self.mrs[a].powi(2) - s)
192 } else {
193 None
194 }
195 })
196 .product()
197 }
198 fn pole_product(&self, s: F) -> F {
199 (0..R).map(|a| (self.mrs[a].powi(2) - s)).product()
200 }
201
202 fn p_vector(
203 betas: &SVector<Complex<F>, R>,
204 pvector_constants: &SMatrix<Complex<F>, C, R>,
205 ) -> SVector<Complex<F>, C> {
206 SVector::<Complex<F>, C>::from_fn(|j, _| {
207 (0..R).map(|a| betas[a] * pvector_constants[(j, a)]).sum()
208 })
209 }
210
211 pub fn calculate_k_matrix(
212 betas: &SVector<Complex<F>, R>,
213 ikc_inv_vec: &SVector<Complex<F>, C>,
214 pvector_constants_mat: &SMatrix<Complex<F>, C, R>,
215 ) -> Complex<F> {
216 ikc_inv_vec.dot(&Self::p_vector(betas, pvector_constants_mat))
217 }
218}
219impl<F: Field + RealField + 'static, const C: usize, const R: usize> KMatrixConstants<F, C, R> {
220 fn ikc_inv(&self, s: F, channel: usize) -> SVector<Complex<F>, C> {
221 let i_mat = SMatrix::<Complex<F>, C, C>::identity().scale(self.pole_product(s));
222 let k_mat = self.k_matrix(s);
223 let c_mat = self.c_matrix(s);
224 let ikc_mat = i_mat + k_mat * c_mat;
225 let ikc_inv_mat = ikc_mat.try_inverse().unwrap();
226 ikc_inv_mat.row(channel).transpose()
227 }
228}
229
230#[derive(Clone)]
231#[allow(clippy::type_complexity)]
232pub struct KMatrixF0<F: Field> {
233 channel: usize,
234 decay: Decay,
235 constants: KMatrixConstants<F, 5, 5>,
236 data: Vec<(SVector<Complex<F>, 5>, SMatrix<Complex<F>, 5, 5>)>,
237}
238#[rustfmt::skip]
239impl<F: Field + 'static> KMatrixF0<F> {
240 pub fn new(channel: usize, decay: Decay) -> Self {
241 Self {
242 channel,
243 decay,
244 constants: KMatrixConstants {
245 g: SMatrix::<F, 5, 5>::from_vec(convert_vec!(vec![
246 0.74987, -0.01257, 0.27536, -0.15102, 0.36103,
247 0.06401, 0.00204, 0.77413, 0.50999, 0.13112,
248 -0.23417, -0.01032, 0.72283, 0.11934, 0.36792,
249 0.01270, 0.26700, 0.09214, 0.02742, -0.04025,
250 -0.14242, 0.22780, 0.15981, 0.16272, -0.17397,
251 ], F)),
252 c: SMatrix::<F, 5, 5>::from_vec(convert_vec!(vec![
253 0.03728, 0.00000, -0.01398, -0.02203, 0.01397,
254 0.00000, 0.00000, 0.00000, 0.00000, 0.00000,
255 -0.01398, 0.00000, 0.02349, 0.03101, -0.04003,
256 -0.02203, 0.00000, 0.03101, -0.13769, -0.06722,
257 0.01397, 0.00000, -0.04003, -0.06722, -0.28401,
258 ], F)),
259 m1s: convert_array!([0.1349768, 2.0 * 0.1349768, 0.493677, 0.547862, 0.547862], F),
260 m2s: convert_array!([0.1349768, 2.0 * 0.1349768, 0.497611, 0.547862, 0.95778], F),
261 mrs: convert_array!([0.51461, 0.90630, 1.23089, 1.46104, 1.69611], F),
262 adler_zero: Some(AdlerZero {
263 s_0: convert!(0.0091125, F),
264 s_norm: F::one(),
265 }),
266 l: 0,
267 },
268 data: Vec::default(),
269 }
270 }
271}
272
273impl<F: Field + RealField> Node<F> for KMatrixF0<F> {
274 fn precalculate(&mut self, dataset: &Dataset<F>) -> Result<(), RustitudeError> {
275 self.data = dataset
276 .events
277 .par_iter()
278 .map(|event| {
279 let s = self.decay.resonance_p4(event).m2();
280 let barrier_mat = self.constants.barrier_matrix(s);
281 let pvector_constants = SMatrix::<Complex<F>, 5, 5>::from_fn(|i, a| {
282 Complex::from(barrier_mat[(i, a)])
283 * self.constants.g[(i, a)]
284 * self.constants.pole_product_remainder(s, a)
285 });
286 (self.constants.ikc_inv(s, self.channel), pvector_constants)
287 })
288 .collect();
289 Ok(())
290 }
291 fn calculate(&self, parameters: &[F], event: &Event<F>) -> Result<Complex<F>, RustitudeError> {
292 let betas = SVector::<Complex<F>, 5>::new(
293 Complex::new(parameters[0], parameters[1]),
294 Complex::new(parameters[2], parameters[3]),
295 Complex::new(parameters[4], parameters[5]),
296 Complex::new(parameters[6], parameters[7]),
297 Complex::new(parameters[8], parameters[9]),
298 );
299 let (ikc_inv_vec, pvector_constants_mat) = self.data[event.index];
300 Ok(KMatrixConstants::calculate_k_matrix(
301 &betas,
302 &ikc_inv_vec,
303 &pvector_constants_mat,
304 ))
305 }
306 fn parameters(&self) -> Vec<String> {
307 vec![
308 "f0_500 re".to_string(),
309 "f0_500 im".to_string(),
310 "f0_980 re".to_string(),
311 "f0_980 im".to_string(),
312 "f0_1370 re".to_string(),
313 "f0_1370 im".to_string(),
314 "f0_1500 re".to_string(),
315 "f0_1500 im".to_string(),
316 "f0_1710 re".to_string(),
317 "f0_1710 im".to_string(),
318 ]
319 }
320}
321#[derive(Clone)]
322#[allow(clippy::type_complexity)]
323pub struct KMatrixF2<F: Field> {
324 channel: usize,
325 decay: Decay,
326 constants: KMatrixConstants<F, 4, 4>,
327 data: Vec<(SVector<Complex<F>, 4>, SMatrix<Complex<F>, 4, 4>)>,
328}
329#[rustfmt::skip]
330impl<F: Field + 'static> KMatrixF2<F> {
331 pub fn new(channel: usize, decay: Decay) -> Self {
332 Self {
333 channel,
334 decay,
335 constants: KMatrixConstants {
336 g: SMatrix::<F, 4, 4>::from_vec(convert_vec!(vec![
337 0.40033, 0.15479, -0.08900, -0.00113,
338 0.01820, 0.17300, 0.32393, 0.15256,
339 -0.06709, 0.22941, -0.43133, 0.23721,
340 -0.49924, 0.19295, 0.27975, -0.03987,
341 ], F)),
342 c: SMatrix::<F, 4, 4>::from_vec(convert_vec!(vec![
343 -0.04319, 0.00000, 0.00984, 0.01028,
344 0.00000, 0.00000, 0.00000, 0.00000,
345 0.00984, 0.00000, -0.07344, 0.05533,
346 0.01028, 0.00000, 0.05533, -0.05183,
347 ], F)),
348 m1s: convert_array!([0.1349768, 2.0 * 0.1349768, 0.493677, 0.547862], F),
349 m2s: convert_array!([0.1349768, 2.0 * 0.1349768, 0.497611, 0.547862], F),
350 mrs: convert_array!([1.15299, 1.48359, 1.72923, 1.96700], F),
351 adler_zero: None,
352 l: 2,
353 },
354 data: Vec::default()
355 }
356 }
357}
358
359impl<F: Field + RealField> Node<F> for KMatrixF2<F> {
360 fn precalculate(&mut self, dataset: &Dataset<F>) -> Result<(), RustitudeError> {
361 self.data = dataset
362 .events
363 .par_iter()
364 .map(|event| {
365 let s = self.decay.resonance_p4(event).m2();
366 let barrier_mat = self.constants.barrier_matrix(s);
367 let pvector_constants = SMatrix::<Complex<F>, 4, 4>::from_fn(|i, a| {
368 Complex::from(barrier_mat[(i, a)])
369 * self.constants.g[(i, a)]
370 * self.constants.pole_product_remainder(s, a)
371 });
372 (self.constants.ikc_inv(s, self.channel), pvector_constants)
373 })
374 .collect();
375 Ok(())
376 }
377 fn calculate(&self, parameters: &[F], event: &Event<F>) -> Result<Complex<F>, RustitudeError> {
378 let betas = SVector::<Complex<F>, 4>::new(
379 Complex::new(parameters[0], parameters[1]),
380 Complex::new(parameters[2], parameters[3]),
381 Complex::new(parameters[4], parameters[5]),
382 Complex::new(parameters[6], parameters[7]),
383 );
384 let (ikc_inv_vec, pvector_constants_mat) = self.data[event.index];
385 Ok(KMatrixConstants::calculate_k_matrix(
386 &betas,
387 &ikc_inv_vec,
388 &pvector_constants_mat,
389 ))
390 }
391 fn parameters(&self) -> Vec<String> {
392 vec![
393 "f2_1270 re".to_string(),
394 "f2_1270 im".to_string(),
395 "f2_1525 re".to_string(),
396 "f2_1525 im".to_string(),
397 "f2_1810 re".to_string(),
398 "f2_1810 im".to_string(),
399 "f2_1950 re".to_string(),
400 "f2_1950 im".to_string(),
401 ]
402 }
403}
404
405#[derive(Clone)]
406#[allow(clippy::type_complexity)]
407pub struct KMatrixA0<F: Field> {
408 channel: usize,
409 decay: Decay,
410 constants: KMatrixConstants<F, 2, 2>,
411 data: Vec<(SVector<Complex<F>, 2>, SMatrix<Complex<F>, 2, 2>)>,
412}
413#[rustfmt::skip]
414impl<F: Field + 'static> KMatrixA0<F> {
415 pub fn new(channel: usize, decay: Decay) -> Self {
416 Self {
417 channel,
418 decay,
419 constants: KMatrixConstants {
420 g: SMatrix::<F, 2, 2>::from_vec(convert_vec!(vec![
421 0.43215, -0.28825,
422 0.19000, 0.43372
423 ], F)),
424 c: SMatrix::<F, 2, 2>::from_vec(convert_vec!(vec![
425 0.00000, 0.00000,
426 0.00000, 0.00000
427 ], F)),
428 m1s: convert_array!([0.1349768, 0.493677], F),
429 m2s: convert_array!([0.547862, 0.497611], F),
430 mrs: convert_array!([0.95395, 1.26767], F),
431 adler_zero: None,
432 l: 0,
433 },
434 data: Vec::default()
435 }
436 }
437}
438
439impl<F: Field + RealField> Node<F> for KMatrixA0<F> {
440 fn precalculate(&mut self, dataset: &Dataset<F>) -> Result<(), RustitudeError> {
441 self.data = dataset
442 .events
443 .par_iter()
444 .map(|event| {
445 let s = self.decay.resonance_p4(event).m2();
446 let barrier_mat = self.constants.barrier_matrix(s);
447 let pvector_constants = SMatrix::<Complex<F>, 2, 2>::from_fn(|i, a| {
448 Complex::from(barrier_mat[(i, a)])
449 * self.constants.g[(i, a)]
450 * self.constants.pole_product_remainder(s, a)
451 });
452 (self.constants.ikc_inv(s, self.channel), pvector_constants)
453 })
454 .collect();
455 Ok(())
456 }
457 fn calculate(&self, parameters: &[F], event: &Event<F>) -> Result<Complex<F>, RustitudeError> {
458 let betas = SVector::<Complex<F>, 2>::new(
459 Complex::new(parameters[0], parameters[1]),
460 Complex::new(parameters[2], parameters[3]),
461 );
462 let (ikc_inv_vec, pvector_constants_mat) = self.data[event.index];
463 Ok(KMatrixConstants::calculate_k_matrix(
464 &betas,
465 &ikc_inv_vec,
466 &pvector_constants_mat,
467 ))
468 }
469 fn parameters(&self) -> Vec<String> {
470 vec![
471 "a0_980 re".to_string(),
472 "a0_980 im".to_string(),
473 "a0_1450 re".to_string(),
474 "a0_1450 im".to_string(),
475 ]
476 }
477}
478
479#[derive(Clone)]
480#[allow(clippy::type_complexity)]
481pub struct KMatrixA2<F: Field> {
482 channel: usize,
483 decay: Decay,
484 constants: KMatrixConstants<F, 3, 2>,
485 data: Vec<(SVector<Complex<F>, 3>, SMatrix<Complex<F>, 3, 2>)>,
486}
487#[rustfmt::skip]
488impl<F: Field + 'static> KMatrixA2<F> {
489 pub fn new(channel: usize, decay: Decay) -> Self {
490 Self {
491 channel,
492 decay,
493 constants: KMatrixConstants {
494 g: SMatrix::<F, 3, 2>::from_vec(convert_vec!(vec![
495 0.30073, 0.21426, -0.09162,
496 0.68567, 0.12543, 0.00184
497 ], F)),
498 c: SMatrix::<F, 3, 3>::from_vec(convert_vec!(vec![
499 -0.40184, 0.00033, -0.08707,
500 0.00033, -0.21416, -0.06193,
501 -0.08707, -0.06193, -0.17435,
502 ], F)),
503 m1s: convert_array!([0.1349768, 0.493677, 0.1349768], F),
504 m2s: convert_array!([0.547862, 0.497611, 0.95778], F),
505 mrs: convert_array!([1.30080, 1.75351], F),
506 adler_zero: None,
507 l: 2,
508 },
509 data: Vec::default()
510 }
511 }
512}
513
514impl<F: Field + RealField> Node<F> for KMatrixA2<F> {
515 fn precalculate(&mut self, dataset: &Dataset<F>) -> Result<(), RustitudeError> {
516 self.data = dataset
517 .events
518 .par_iter()
519 .map(|event| {
520 let s = self.decay.resonance_p4(event).m2();
521 let barrier_mat = self.constants.barrier_matrix(s);
522 let pvector_constants = SMatrix::<Complex<F>, 3, 2>::from_fn(|i, a| {
523 Complex::from(barrier_mat[(i, a)])
524 * self.constants.g[(i, a)]
525 * self.constants.pole_product_remainder(s, a)
526 });
527 (self.constants.ikc_inv(s, self.channel), pvector_constants)
528 })
529 .collect();
530 Ok(())
531 }
532 fn calculate(&self, parameters: &[F], event: &Event<F>) -> Result<Complex<F>, RustitudeError> {
533 let betas = SVector::<Complex<F>, 2>::new(
534 Complex::new(parameters[0], parameters[1]),
535 Complex::new(parameters[2], parameters[3]),
536 );
537 let (ikc_inv_vec, pvector_constants_mat) = self.data[event.index];
538 Ok(KMatrixConstants::calculate_k_matrix(
539 &betas,
540 &ikc_inv_vec,
541 &pvector_constants_mat,
542 ))
543 }
544 fn parameters(&self) -> Vec<String> {
545 vec![
546 "a2_1320 re".to_string(),
547 "a2_1320 im".to_string(),
548 "a2_1700 re".to_string(),
549 "a2_1700 im".to_string(),
550 ]
551 }
552}
553
554#[derive(Clone)]
555#[allow(clippy::type_complexity)]
556pub struct KMatrixRho<F: Field> {
557 channel: usize,
558 decay: Decay,
559 constants: KMatrixConstants<F, 3, 2>,
560 data: Vec<(SVector<Complex<F>, 3>, SMatrix<Complex<F>, 3, 2>)>,
561}
562#[rustfmt::skip]
563impl<F: Field + 'static> KMatrixRho<F> {
564 pub fn new(channel: usize, decay: Decay) -> Self {
565 Self {
566 channel,
567 decay,
568 constants: KMatrixConstants {
569 g: SMatrix::<F, 3, 2>::from_vec(convert_vec!(vec![
570 0.28023, 0.01806, 0.06501,
571 0.16318, 0.53879, 0.00495,
572 ], F)),
573 c: SMatrix::<F, 3, 3>::from_vec(convert_vec!(vec![
574 -0.06948, 0.00000, 0.07958,
575 0.00000, 0.00000, 0.00000,
576 0.07958, 0.00000, -0.60000,
577 ], F)),
578 m1s: convert_array!([0.1349768, 2.0 * 0.1349768, 0.493677], F),
579 m2s: convert_array!([0.1349768, 2.0 * 0.1349768, 0.497611], F),
580 mrs: convert_array!([0.71093, 1.58660], F),
581 adler_zero: None,
582 l: 1,
583 },
584 data: Vec::default(),
585 }
586 }
587}
588
589impl<F: Field + RealField> Node<F> for KMatrixRho<F> {
590 fn precalculate(&mut self, dataset: &Dataset<F>) -> Result<(), RustitudeError> {
591 self.data = dataset
592 .events
593 .par_iter()
594 .map(|event| {
595 let s = self.decay.resonance_p4(event).m2();
596 let barrier_mat = self.constants.barrier_matrix(s);
597 let pvector_constants = SMatrix::<Complex<F>, 3, 2>::from_fn(|i, a| {
598 Complex::from(barrier_mat[(i, a)])
599 * self.constants.g[(i, a)]
600 * self.constants.pole_product_remainder(s, a)
601 });
602 (self.constants.ikc_inv(s, self.channel), pvector_constants)
603 })
604 .collect();
605 Ok(())
606 }
607 fn calculate(&self, parameters: &[F], event: &Event<F>) -> Result<Complex<F>, RustitudeError> {
608 let betas = SVector::<Complex<F>, 2>::new(
609 Complex::new(parameters[0], parameters[1]),
610 Complex::new(parameters[2], parameters[3]),
611 );
612 let (ikc_inv_vec, pvector_constants_mat) = self.data[event.index];
613 Ok(KMatrixConstants::calculate_k_matrix(
614 &betas,
615 &ikc_inv_vec,
616 &pvector_constants_mat,
617 ))
618 }
619 fn parameters(&self) -> Vec<String> {
620 vec![
621 "rho_770 re".to_string(),
622 "rho_770 im".to_string(),
623 "rho_1700 re".to_string(),
624 "rho_1700 im".to_string(),
625 ]
626 }
627}
628
629#[derive(Clone)]
630#[allow(clippy::type_complexity)]
631pub struct KMatrixPi1<F: Field> {
632 channel: usize,
633 decay: Decay,
634 constants: KMatrixConstants<F, 2, 1>,
635 data: Vec<(SVector<Complex<F>, 2>, SMatrix<Complex<F>, 2, 1>)>,
636}
637#[rustfmt::skip]
638impl<F: Field + 'static> KMatrixPi1<F> {
639 pub fn new(channel: usize, decay: Decay) -> Self {
640 Self {
641 channel,
642 decay,
643 constants: KMatrixConstants {
644 g: SMatrix::<F, 2, 1>::from_vec(convert_vec!(vec![
645 0.80564, 1.04595
646 ], F)),
647 c: SMatrix::<F, 2, 2>::from_vec(convert_vec!(vec![
648 1.05000, 0.15163,
649 0.15163, -0.24611,
650 ], F)),
651 m1s: convert_array!([0.1349768, 0.1349768], F),
652 m2s: convert_array!([0.547862, 0.95778], F),
653 mrs: convert_array!([1.38552], F),
654 adler_zero: None,
655 l: 1,
656 },
657 data: Vec::default()
658 }
659 }
660}
661
662impl<F: Field + RealField> Node<F> for KMatrixPi1<F> {
663 fn precalculate(&mut self, dataset: &Dataset<F>) -> Result<(), RustitudeError> {
664 self.data = dataset
665 .events
666 .par_iter()
667 .map(|event| {
668 let s = self.decay.resonance_p4(event).m2();
669 let barrier_mat = self.constants.barrier_matrix(s);
670 let pvector_constants = SMatrix::<Complex<F>, 2, 1>::from_fn(|i, a| {
671 Complex::from(barrier_mat[(i, a)])
672 * self.constants.g[(i, a)]
673 * self.constants.pole_product_remainder(s, a)
674 });
675 (self.constants.ikc_inv(s, self.channel), pvector_constants)
676 })
677 .collect();
678 Ok(())
679 }
680 fn calculate(&self, parameters: &[F], event: &Event<F>) -> Result<Complex<F>, RustitudeError> {
681 let betas = SVector::<Complex<F>, 1>::new(Complex::new(parameters[0], parameters[1]));
682 let (ikc_inv_vec, pvector_constants_mat) = self.data[event.index];
683 Ok(KMatrixConstants::calculate_k_matrix(
684 &betas,
685 &ikc_inv_vec,
686 &pvector_constants_mat,
687 ))
688 }
689 fn parameters(&self) -> Vec<String> {
690 vec!["pi1_1600 re".to_string(), "pi1_1600 im".to_string()]
691 }
692}