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
121
122
123
124
125
126
127
128
129
130

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 {
        let dy = p1.rap() - p2.rap();
        let dphi = p1.phi() - p2.phi();

        let delta_sq = dy*dy + dphi*dphi;

        min(p1.inv_pt2(), p2.inv_pt2()) * delta_sq/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 {
        let dy = p1.rap() - p2.rap();
        let dphi = p1.phi() - p2.phi();

        let delta_sq = dy*dy + dphi*dphi;

        min(p1.pt2(), p2.pt2()) * delta_sq/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 {
        let dy = p1.rap() - p2.rap();
        let dphi = p1.phi() - p2.phi();

        let delta_sq = dy*dy + dphi*dphi;

        delta_sq/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 {
        let dy = p1.rap() - p2.rap();
        let dphi = p1.phi() - p2.phi();

        let delta_sq = dy*dy + dphi*dphi;

        min(p1.pt2().powf(self.p), p2.pt2().powf(self.p)) * delta_sq/self.r2
    }

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