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
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120

use crate::pseudojet::PseudoJet;

use std::cmp::min;

use noisy_float::prelude::*;

/// Distance measure between pseudojets used for clustering
pub trait Distance {
    /// Distance between pseudojets
    fn distance(&self, p1: &PseudoJet, p2: &PseudoJet) -> N64;
    /// Distance to the beam axis
    fn beam_distance(&self, p1: &PseudoJet) -> N64;
}

pub struct AntiKt {
    r2: N64,
}

/// anti-kt distance measure with radius parameter `r`
pub fn anti_kt(r: N64) -> AntiKt {
    AntiKt { r2: r * r }
}

/// anti-kt distance measure with radius parameter `r`
pub fn anti_kt_f(r: f64) -> AntiKt {
    anti_kt(n64(r))
}

impl Distance for AntiKt {
    fn distance(&self, p1: &PseudoJet, p2: &PseudoJet) -> N64 {
        min(p1.inv_pt2(), p2.inv_pt2()) * p1.delta_r2(p2) / self.r2
    }

    fn beam_distance(&self, p1: &PseudoJet) -> N64 {
        p1.inv_pt2()
    }
}

pub struct Kt {
    r2: N64,
}

/// kt distance measure with radius parameter `r`
pub fn kt(r: N64) -> Kt {
    Kt { r2: r * r }
}

/// kt distance measure with radius parameter `r`
pub fn kt_f(r: f64) -> Kt {
    kt(n64(r))
}

impl Distance for Kt {
    fn distance(&self, p1: &PseudoJet, p2: &PseudoJet) -> N64 {
        min(p1.pt2(), p2.pt2()) * p1.delta_r2(p2) / self.r2
    }

    fn beam_distance(&self, p1: &PseudoJet) -> N64 {
        p1.pt2()
    }
}

pub struct CambridgeAachen {
    r2: N64,
}

/// Cambridge/Aachen distance measure with radius parameter `r`
pub fn cambridge_aachen(r: N64) -> CambridgeAachen {
    CambridgeAachen { r2: r * r }
}

/// Cambridge/Aachen distance measure with radius parameter `r`
pub fn cambridge_aachen_f(r: f64) -> CambridgeAachen {
    cambridge_aachen(n64(r))
}

impl Distance for CambridgeAachen {
    fn distance(&self, p1: &PseudoJet, p2: &PseudoJet) -> N64 {
        p1.delta_r2(p2) / self.r2
    }

    fn beam_distance(&self, _p1: &PseudoJet) -> N64 {
        n64(1.)
    }
}

pub struct GenKt {
    r2: N64,
    p: N64,
}

/// Generalised kt distance measure with radius parameter `r` and exponent `p`
pub fn gen_kt(r: N64, p: N64) -> GenKt {
    GenKt { r2: r * r, p }
}

/// Generalised kt distance measure with radius parameter `r` and exponent `p`
pub fn gen_kt_f(r: f64, p: f64) -> GenKt {
    gen_kt(n64(r), n64(p))
}

impl Distance for GenKt {
    fn distance(&self, p1: &PseudoJet, p2: &PseudoJet) -> N64 {
        min(p1.pt2().powf(self.p), p2.pt2().powf(self.p)) * p1.delta_r2(p2) / self.r2
    }

    fn beam_distance(&self, p1: &PseudoJet) -> N64 {
        p1.pt2().powf(self.p)
    }
}

impl<T: Distance> Distance for &T {
    fn distance(&self, p1: &PseudoJet, p2: &PseudoJet) -> N64 {
        (*self).distance(p1, p2)
    }

    fn beam_distance(&self, p1: &PseudoJet) -> N64 {
        (*self).beam_distance(p1)
    }
}