1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105

use noisy_float::prelude::*;

/// A basic four-vector
///
/// The zero component is the energy/time component. The remainder are
/// the spatial components
#[derive(PartialEq, Eq, PartialOrd, Ord, Debug, Clone, Copy, Default)]
pub struct FourVector {
    pt: N64,
    p: [N64; 4],
}

impl FourVector {
    /// Construct a new four-vector
    pub fn new() -> Self {
        Self::default()
    }

    /// The euclidean norm \sqrt{\sum v_\mu^2} with \mu = 0,1,2,3
    pub fn euclid_norm(&self) -> N64 {
        self.euclid_norm_sq().sqrt()
    }
    /// The square \sum v_\mu^2 with \mu = 0,1,2,3 of the euclidean norm
    pub fn euclid_norm_sq(&self) -> N64 {
        self.p.iter().map(|e| *e * *e).sum()
    }

    /// The spatial norm \sqrt{\sum v_i^2} with i = 1,2,3
    pub fn spatial_norm(&self) -> N64 {
        self.spatial_norm_sq().sqrt()
    }

    /// The square \sum v_i^2 with i = 1,2,3 of the spatial norm
    pub fn spatial_norm_sq(&self) -> N64 {
        self.p.iter().skip(1).map(|e| *e * *e).sum()
    }

    /// The scalar transverse momentum
    pub fn pt(&self) -> N64 {
        self.pt
    }

    const fn len() -> usize {
        4
    }

    fn update_pt(&mut self) {
        self.pt = (self.p[1] * self.p[1] + self.p[2] * self.p[2]).sqrt();
    }
}

impl std::convert::From<[N64; 4]> for FourVector {
    fn from(p: [N64; 4]) -> FourVector {
        let mut res = FourVector {
            p,
            pt: std::default::Default::default(),
        };
        res.update_pt();
        res
    }
}

impl std::ops::Index<usize> for FourVector {
    type Output = N64;

    fn index(&self, i: usize) -> &Self::Output {
        &self.p[i]
    }
}

impl std::ops::AddAssign for FourVector {
    fn add_assign(&mut self, rhs: FourVector) {
        for i in 0..Self::len() {
            self.p[i] += rhs[i]
        }
        self.update_pt();
    }
}

impl std::ops::SubAssign for FourVector {
    fn sub_assign(&mut self, rhs: FourVector) {
        for i in 0..Self::len() {
            self.p[i] -= rhs[i]
        }
        self.update_pt();
    }
}

impl std::ops::Add for FourVector {
    type Output = Self;

    fn add(mut self, rhs: FourVector) -> Self::Output {
        self += rhs;
        self
    }
}

impl std::ops::Sub for FourVector {
    type Output = Self;

    fn sub(mut self, rhs: FourVector) -> Self::Output {
        self -= rhs;
        self
    }
}