Skip to content

Commit 24a35e5

Browse files
Jet utilities (#141) (#144)
Add functions to calculate the momentum fraction and the kt scale for jets. - pt_fraction() for transverse momentum fraction - kt_scale() These are implemented generically for both jet types, though may be primary interesting for pp events (PseudoJet). deltar is added as a helper that returns the jet separation in pseudorapidity-phi. This is distinct from deltaR which operates in rapidity-phi space. These are probably not the greatest names, but are not in the public API. deltaR is moved to the new source file, JetUtils.jl. There is a small tidy up of the logic for calculating ϕ where the check for ϕ>2π can never be true. Added functions to convert jets to LorentzVector and LorentzVectorCyl. This required some initial refactoring of jet methods to be general, rather than PseudoJet specific, which will be continued. (cherry picked from commit c7cfa9e)
1 parent a71556e commit 24a35e5

9 files changed

+214
-87
lines changed

docs/make.jl

+3-2
Original file line numberDiff line numberDiff line change
@@ -13,11 +13,12 @@ makedocs(sitename = "JetReconstruction.jl",
1313
pages = [
1414
"Home" => "index.md",
1515
"Examples" => "examples.md",
16-
"EDM4hep" => "EDM4hep.md",
17-
"Visualisation" => "visualisation.md",
1816
"Particle Inputs" => "particles.md",
1917
"Reconstruction Strategies" => "strategy.md",
2018
"Substructure" => "substructure.md",
19+
"Jet Helpers" => "helpers.md",
20+
"EDM4hep" => "EDM4hep.md",
21+
"Visualisation" => "visualisation.md",
2122
"Reference Docs" => Any["Public API" => "lib/public.md",
2223
"Internal API" => "lib/internal.md"],
2324
"Extras" => Any["Serialisation" => "extras/serialisation.md"]

docs/src/helpers.md

+25
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,25 @@
1+
# Jet Helper Functions
2+
3+
These functions are provided as a convenient way to work with the results of jet
4+
reconstruction.
5+
6+
## Jet Pair Helpers
7+
8+
- [`pt_fraction`](@ref)
9+
10+
Returns the transverse momentum fraction in the softer of the two jets.
11+
12+
- [`kt_scale`](@ref)
13+
14+
Returns the transverse momentum scale between two jets, the product of the
15+
smaller ``p_T`` and the angular separation in the the η-ϕ plane.
16+
17+
## Conversion Functions
18+
19+
- [`lorentzvector`](@ref)
20+
21+
Return a cartesian `LorentzVector` from a jet.
22+
23+
- [`lorentzvector_cyl`](@ref)
24+
25+
Return a cylindrical `LorentzVectorCyl` from a jet.

src/EEjet.jl

-2
Original file line numberDiff line numberDiff line change
@@ -49,8 +49,6 @@ phi(eej::EEjet) = begin
4949
phi = pt2(eej) == 0.0 ? 0.0 : atan(eej.py, eej.px)
5050
if phi < 0.0
5151
phi += 2π
52-
elseif phi >= 2π
53-
phi -= 2π # can happen if phi=-|eps<1e-15|?
5452
end
5553
phi
5654
end

src/JetReconstruction.jl

+2-1
Original file line numberDiff line numberDiff line change
@@ -39,7 +39,8 @@ energy(p::LorentzVectorCyl) = LorentzVectorHEP.energy(p)
3939
# Pseudojet and EEjet types
4040
include("Pseudojet.jl")
4141
include("EEjet.jl")
42-
export PseudoJet, EEjet
42+
include("JetUtils.jl")
43+
export PseudoJet, EEjet, pt_fraction, kt_scale, lorentzvector, lorentzvector_cyl
4344

4445
# Jet reconstruction strategies and algorithms (enums!)
4546
include("AlgorithmStrategyEnums.jl")

src/JetUtils.jl

+137
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,137 @@
1+
# Generic utility functions for jet structs
2+
3+
"""
4+
mag(jet::T) where {T <: FourMomentum}
5+
6+
Return the magnitude of the momentum of a jet, `|p|`.
7+
8+
# Returns
9+
The magnitude of the jet.
10+
"""
11+
mag(jet::T) where {T <: FourMomentum} = sqrt(muladd(jet.px, jet.px, jet.py^2) + jet.pz^2)
12+
13+
"""
14+
CosTheta(jet::T) where {T <: FourMomentum}
15+
16+
Compute the cosine of the angle between the momentum vector of `jet` and the z-axis.
17+
18+
# Returns
19+
- The cosine of the angle between the jet and the z-axis.
20+
"""
21+
@inline function CosTheta(jet::T) where {T <: FourMomentum}
22+
fZ = jet.pz
23+
ptot = mag(jet)
24+
return ifelse(ptot == 0.0, 1.0, fZ / ptot)
25+
end
26+
27+
"""
28+
eta(jet::T) where {T <: FourMomentum}
29+
30+
Compute the pseudorapidity (η) of a jet.
31+
32+
# Returns
33+
- The pseudorapidity (η) of the jet.
34+
"""
35+
function eta(jet::T) where {T <: FourMomentum}
36+
cosTheta = CosTheta(jet)
37+
(cosTheta^2 < 1.0) && return -0.5 * log((1.0 - cosTheta) / (1.0 + cosTheta))
38+
fZ = jet.pz
39+
iszero(fZ) && return 0.0
40+
# Warning("PseudoRapidity","transverse momentum = 0! return +/- 10e10");
41+
fZ > 0.0 && return 10e10
42+
return -10e10
43+
end
44+
45+
"""
46+
const η = eta
47+
48+
Alias for the pseudorapidity function, `eta`.
49+
"""
50+
const η = eta
51+
52+
"""
53+
deltaR(jet1::T, jet2::T) where T <: FourMomentum
54+
55+
Function to calculate the distance in the y-ϕ plane between two jets `jet1` and
56+
`jet2` (that is using *rapidity* and *azimuthal angle*).
57+
58+
# Returns
59+
- The Euclidean distance in the y-ϕ plane for the two jets.
60+
"""
61+
function deltaR(jet1::T, jet2::T) where {T <: FourMomentum}
62+
δy = rapidity(jet1) - rapidity(jet2)
63+
δϕ = phi(jet1) - phi(jet2)
64+
δϕ = abs(δϕ) > π ? 2π - abs(δϕ) : δϕ
65+
66+
return sqrt(δy^2 + δϕ^2)
67+
end
68+
69+
"""
70+
deltar(jet1::T, jet2::T) where T <: FourMomentum
71+
72+
Function to calculate the distance in the η-ϕ plane between two jets `jet1` and
73+
`jet2` (that is, using the *pseudorapidity* and *azimuthal angle*).
74+
75+
# Returns
76+
- The Euclidean distance in the η-ϕ plane for the two jets.
77+
"""
78+
function deltar(jet1::T, jet2::T) where {T <: FourMomentum}
79+
δη = eta(jet1) - eta(jet2)
80+
δϕ = phi(jet1) - phi(jet2)
81+
δϕ = abs(δϕ) > π ? 2π - abs(δϕ) : δϕ
82+
83+
return sqrt(δη^2 + δϕ^2)
84+
end
85+
86+
"""
87+
pt_fraction(jet1::T, jet2::T) where T <: FourMomentum
88+
89+
Computes the transverse momentum fraction of the softer of two jets.
90+
91+
# Returns
92+
- The transverse momentum fraction of the softer of the two jets.
93+
"""
94+
function pt_fraction(jet1::T, jet2::T) where {T <: FourMomentum}
95+
pt1 = JetReconstruction.pt(jet1)
96+
pt2 = JetReconstruction.pt(jet2)
97+
return min(pt1, pt2) / (pt1 + pt2)
98+
end
99+
100+
"""
101+
kt_scale(jet1::T, jet2::T) where {T <: FourMomentum}
102+
103+
Computes the transverse momentum scale as the product of the minimum pt and
104+
the angular separation in the η-ϕ plane (using *pseudorapidity*).
105+
106+
# Returns
107+
- The transverse momentum scale of the two jets.
108+
"""
109+
function kt_scale(jet1::T, jet2::T) where {T <: FourMomentum}
110+
pt1 = JetReconstruction.pt(jet1)
111+
pt2 = JetReconstruction.pt(jet2)
112+
return min(pt1, pt2) * deltar(jet1, jet2)
113+
end
114+
115+
"""
116+
lorentzvector_cyl(jet::T) where T <: FourMomentum -> LorentzVectorHEP.LorentzVectorCyl
117+
118+
Return a cylindrical `LorentzVectorCyl` from a jet.
119+
"""
120+
function lorentzvector_cyl(jet::T) where {T <: FourMomentum}
121+
return LorentzVectorHEP.fromPtEtaPhiE(JetReconstruction.pt(jet),
122+
JetReconstruction.eta(jet),
123+
JetReconstruction.phi(jet),
124+
JetReconstruction.energy(jet))
125+
end
126+
127+
"""
128+
lorentzvector(jet::T) where {T <: FourMomentum} -> -> LorentzVector
129+
130+
Return a cartesian `LorentzVector` from a jet.
131+
"""
132+
function lorentzvector(jet::T) where {T <: FourMomentum}
133+
return LorentzVector(JetReconstruction.energy(jet),
134+
JetReconstruction.px(jet),
135+
JetReconstruction.py(jet),
136+
JetReconstruction.pz(jet))
137+
end

src/Pseudojet.jl

-60
Original file line numberDiff line numberDiff line change
@@ -167,8 +167,6 @@ _set_rap_phi!(p::PseudoJet) = begin
167167
p._phi = p._pt2 == 0.0 ? 0.0 : atan(p.py, p.px)
168168
if p._phi < 0.0
169169
p._phi += 2π
170-
elseif p._phi >= 2π
171-
p._phi -= 2π # can happen if phi=-|eps<1e-15|?
172170
end
173171

174172
if p.E == abs(p.pz) && iszero(p._pt2)
@@ -278,64 +276,6 @@ Calculate the invariant mass squared (m^2) of a PseudoJet.
278276
"""
279277
m2(p::PseudoJet) = (p.E + p.pz) * (p.E - p.pz) - p._pt2
280278

281-
"""
282-
mag(p::PseudoJet)
283-
284-
Return the magnitude of the momentum of a `PseudoJet`, `|p|`.
285-
286-
# Arguments
287-
- `p::PseudoJet`: The `PseudoJet` object for which to compute the magnitude.
288-
289-
# Returns
290-
The magnitude of the `PseudoJet` object.
291-
"""
292-
mag(p::PseudoJet) = sqrt(muladd(p.px, p.px, p.py^2) + p.pz^2)
293-
294-
"""
295-
CosTheta(p::PseudoJet)
296-
297-
Compute the cosine of the angle between the momentum vector `p` and the z-axis.
298-
299-
# Arguments
300-
- `p::PseudoJet`: The PseudoJet object representing the momentum vector.
301-
302-
# Returns
303-
- The cosine of the angle between `p` and the z-axis.
304-
"""
305-
@inline function CosTheta(p::PseudoJet)
306-
fZ = p.pz
307-
ptot = mag(p)
308-
return ifelse(ptot == 0.0, 1.0, fZ / ptot)
309-
end
310-
311-
"""
312-
eta(p::PseudoJet)
313-
314-
Compute the pseudorapidity (η) of a PseudoJet.
315-
316-
# Arguments
317-
- `p::PseudoJet`: The PseudoJet object for which to compute the pseudorapidity.
318-
319-
# Returns
320-
- The pseudorapidity (η) of the PseudoJet.
321-
"""
322-
function eta(p::PseudoJet)
323-
cosTheta = CosTheta(p)
324-
(cosTheta^2 < 1.0) && return -0.5 * log((1.0 - cosTheta) / (1.0 + cosTheta))
325-
fZ = p.pz
326-
iszero(fZ) && return 0.0
327-
# Warning("PseudoRapidity","transverse momentum = 0! return +/- 10e10");
328-
fZ > 0.0 && return 10e10
329-
return -10e10
330-
end
331-
332-
"""
333-
const η = eta
334-
335-
Alias for the pseudorapidity function, `eta`.
336-
"""
337-
const η = eta
338-
339279
"""
340280
m(p::PseudoJet)
341281

src/Substructure.jl

+1-22
Original file line numberDiff line numberDiff line change
@@ -1,25 +1,4 @@
1-
"""
2-
deltaR(jet1, jet2) -> Float64
3-
4-
Function to calculate the distance in the y-ϕ plane between two jets `jet1` and `jet2`
5-
6-
# Arguments
7-
- `jet1::PseudoJet`: The first jet.
8-
- `jet2::PseudoJet`: The second jet.
9-
10-
# Returns
11-
- `Float64`: The Euclidean distance in the y-ϕ plane for the two jets.
12-
"""
13-
function deltaR(jet1::PseudoJet, jet2::PseudoJet)
14-
y1, phi1 = rapidity(jet1), phi(jet1)
15-
y2, phi2 = rapidity(jet2), phi(jet2)
16-
17-
d_y = y1 - y2
18-
d_phi = phi1 - phi2
19-
d_phi = abs(d_phi) > π ? 2π - abs(d_phi) : d_phi
20-
21-
return sqrt(d_y^2 + d_phi^2)
22-
end
1+
# Substructure specific functions for jet grooming and filtering
232

243
"""
254
recluster(jet, clusterseq; R = 1.0, algorithm = JetAlgorithm.CA) -> ClusterSequence

test/runtests.jl

+3
Original file line numberDiff line numberDiff line change
@@ -6,6 +6,9 @@ function main()
66
# Algorithm/power consistency checks
77
include("test-algpower-consistency.jl")
88

9+
# Jet utilities tests
10+
include("test-jet-utils.jl")
11+
912
# jet_reconstruct() interface check
1013
include("test-jet_reconstruct-interface.jl")
1114

test/test-jet-utils.jl

+43
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,43 @@
1+
# Some basic tests of common utilities for jet types
2+
3+
include("common.jl")
4+
5+
pj1 = PseudoJet(1.0, 1.3, -2.3, 25.0)
6+
pj2 = PseudoJet(1.0, 0.3, -1.3, 17.0)
7+
8+
eej1 = EEjet(1.0, 1.3, -2.3, 25.0)
9+
eej2 = EEjet(-1.0, 3.2, -1.2, 39.0)
10+
11+
@testset "Common jet utilities" begin
12+
@test JetReconstruction.pt_fraction(pj1, pj2) 0.3889609897118418
13+
@test JetReconstruction.pt_fraction(eej1, eej2) 0.32850184248668773
14+
15+
jr_kt_scale = JetReconstruction.kt_scale(pj1, pj2)
16+
lvhep_kt_scale = min(JetReconstruction.pt(pj1), JetReconstruction.pt(pj2)) *
17+
deltar(JetReconstruction.lorentzvector_cyl(pj1),
18+
JetReconstruction.lorentzvector_cyl(pj2))
19+
@test jr_kt_scale lvhep_kt_scale
20+
21+
# Test conversions
22+
lv_pj1 = JetReconstruction.lorentzvector(pj1)
23+
lv_eej1 = JetReconstruction.lorentzvector(eej1)
24+
@test LorentzVectorHEP.energy(lv_pj1) JetReconstruction.energy(pj1)
25+
JetReconstruction.energy(eej1) LorentzVectorHEP.energy(lv_eej1)
26+
@test LorentzVectorHEP.px(lv_pj1) JetReconstruction.px(pj1)
27+
JetReconstruction.px(eej1) LorentzVectorHEP.px(lv_eej1)
28+
@test LorentzVectorHEP.py(lv_pj1) JetReconstruction.py(pj1)
29+
JetReconstruction.py(eej1) LorentzVectorHEP.py(lv_eej1)
30+
@test LorentzVectorHEP.pz(lv_pj1) JetReconstruction.pz(pj1)
31+
JetReconstruction.pz(eej1) LorentzVectorHEP.pz(lv_eej1)
32+
33+
lvc_pj1 = JetReconstruction.lorentzvector_cyl(pj1)
34+
lvc_eej1 = JetReconstruction.lorentzvector_cyl(eej1)
35+
@test LorentzVectorHEP.pt(lvc_pj1) JetReconstruction.pt(eej1)
36+
JetReconstruction.pt(eej1) LorentzVectorHEP.pt(lvc_eej1)
37+
@test LorentzVectorHEP.eta(lvc_pj1) JetReconstruction.eta(eej1)
38+
JetReconstruction.eta(eej1) LorentzVectorHEP.eta(lvc_eej1)
39+
@test LorentzVectorHEP.phi(lvc_pj1) JetReconstruction.phi(eej1)
40+
JetReconstruction.phi(eej1) LorentzVectorHEP.phi(lvc_eej1)
41+
@test LorentzVectorHEP.mass(lvc_pj1) JetReconstruction.mass(eej1)
42+
JetReconstruction.mass(eej1) LorentzVectorHEP.mass(lvc_eej1)
43+
end

0 commit comments

Comments
 (0)