From c1d88c3ca8d0f75ee92a37b3993e0c8ad36ff13a Mon Sep 17 00:00:00 2001
From: Graeme A Stewart <graeme.andrew.stewart@cern.ch>
Date: Fri, 22 Nov 2024 15:56:58 +0100
Subject: [PATCH 01/16] Update PseudoJet to be type T <: Real

Allowing for different Float values, as well as more exotic types, if needed
---
 src/Pseudojet.jl | 72 ++++++++++++++++++++++++------------------------
 1 file changed, 36 insertions(+), 36 deletions(-)

diff --git a/src/Pseudojet.jl b/src/Pseudojet.jl
index 3ecfd43..0d52302 100644
--- a/src/Pseudojet.jl
+++ b/src/Pseudojet.jl
@@ -18,7 +18,7 @@ const _invalid_phi = -100.0
 const _invalid_rap = -1.e200
 
 """
-    mutable struct PseudoJet <: FourMomentum
+    mutable struct PseudoJet{T <: Real} <: FourMomentum
 
 The `PseudoJet` struct represents a pseudojet, a four-momentum object used in
 jet reconstruction algorithms. Additional information for the link back into the
@@ -26,70 +26,70 @@ history of the clustering is stored in the `_cluster_hist_index` field. There is
 caching of the more expensive calculations for rapidity and azimuthal angle.
 
 # Fields
-- `px::Float64`: The x-component of the momentum.
-- `py::Float64`: The y-component of the momentum.
-- `pz::Float64`: The z-component of the momentum.
-- `E::Float64`: The energy component of the momentum.
+- `px::T`: The x-component of the momentum.
+- `py::T`: The y-component of the momentum.
+- `pz::T`: The z-component of the momentum.
+- `E::T`: The energy component of the momentum.
 - `_cluster_hist_index::Int`: The index of the cluster history.
-- `_pt2::Float64`: The squared transverse momentum.
-- `_inv_pt2::Float64`: The inverse squared transverse momentum.
-- `_rap::Float64`: The rapidity.
-- `_phi::Float64`: The azimuthal angle.
+- `_pt2::T`: The squared transverse momentum.
+- `_inv_pt2::T`: The inverse squared transverse momentum.
+- `_rap::T`: The rapidity.
+- `_phi::T`: The azimuthal angle.
 
 """
-mutable struct PseudoJet <: FourMomentum
-    px::Float64
-    py::Float64
-    pz::Float64
-    E::Float64
+mutable struct PseudoJet{T <: Real} <: FourMomentum
+    px::T
+    py::T
+    pz::T
+    E::T
     _cluster_hist_index::Int
-    _pt2::Float64
-    _inv_pt2::Float64
-    _rap::Float64
-    _phi::Float64
+    _pt2::T
+    _inv_pt2::T
+    _rap::T
+    _phi::T
 end
 
 """
-    PseudoJet(px::Real, py::Real, pz::Real, E::Real,
+    PseudoJet(px::T, py::T, pz::T, E::T,
         _cluster_hist_index::Int,
-        pt2::Real)
+        pt2::T) where {T <: Real}
 
 Constructs a PseudoJet object with the given momentum components and energy and
 history index.
 
 # Arguments
-- `px::Real`: The x-component of the momentum.
-- `py::Real`: The y-component of the momentum.
-- `pz::Real`: The z-component of the momentum.
-- `E::Real`: The energy.
+- `px::T`: The x-component of the momentum.
+- `py::T`: The y-component of the momentum.
+- `pz::T`: The z-component of the momentum.
+- `E::T`: The energy.
 - `_cluster_hist_index::Int`: The cluster history index.
-- `pt2::Real`: The transverse momentum squared.
+- `pt2::T`: The transverse momentum squared.
 
 # Returns
 A `PseudoJet` object.
 """
-PseudoJet(px::Real, py::Real, pz::Real, E::Real,
+PseudoJet(px::T, py::T, pz::T, E::T,
 _cluster_hist_index::Int,
-pt2::Real) = PseudoJet(px,
-                       py, pz, E, _cluster_hist_index,
-                       pt2, 1.0 / pt2, _invalid_rap, _invalid_phi)
+pt2::T) where {T <: Real} = PseudoJet{T}(px,
+                                         py, pz, E, _cluster_hist_index,
+                                         pt2, 1.0 / pt2, _invalid_rap, _invalid_phi)
 
 """
-    PseudoJet(px::Real, py::Real, pz::Real, E::Real)
+    PseudoJet(px::T, py::T, pz::T, E::T) where T <: Real
 
 Constructs a PseudoJet object with the given momentum components and energy.
 
 # Arguments
-- `px::Real`: The x-component of the momentum.
-- `py::Real`: The y-component of the momentum.
-- `pz::Real`: The z-component of the momentum.
-- `E::Real`: The energy.
+- `px::T`: The x-component of the momentum.
+- `py::T`: The y-component of the momentum.
+- `pz::T`: The z-component of the momentum.
+- `E::T`: The energy.
 
 # Returns
 A PseudoJet object.
 """
-PseudoJet(px::Real, py::Real,
-pz::Real, E::Real) = PseudoJet(px, py, pz, E, 0, px^2 + py^2)
+PseudoJet(px::T, py::T,
+pz::T, E::T) where {T <: Real} = PseudoJet(px, py, pz, E, 0, px^2 + py^2)
 
 import Base.show
 """

From 84f8796cf1734bd7177aff3f519f6a6272462665 Mon Sep 17 00:00:00 2001
From: Graeme A Stewart <graeme.andrew.stewart@cern.ch>
Date: Fri, 22 Nov 2024 16:04:58 +0100
Subject: [PATCH 02/16] Update EEJet to be type T <: Real

Allowing for different Float values, as well as more exotic types,
if needed.
---
 src/EEjet.jl | 73 +++++++++++++++++++++++++++++++++-------------------
 1 file changed, 47 insertions(+), 26 deletions(-)

diff --git a/src/EEjet.jl b/src/EEjet.jl
index c61f44f..c9b2c41 100644
--- a/src/EEjet.jl
+++ b/src/EEjet.jl
@@ -1,34 +1,39 @@
 """
-    struct EEjet
+    struct EEjet{T <: Real} <: FourMomentum
 
-The `EEjet` struct is a 4-momentum object used for the e+e jet reconstruction routines.
+The `EEjet` struct is a 4-momentum object used for the e+e jet reconstruction
+routines. Internal fields are used to track the reconstuction and to cache
+values needed during the execution of the algorithm.
 
 # Fields
-- `px::Float64`: The x-component of the jet momentum.
-- `py::Float64`: The y-component of the jet momentum.
-- `pz::Float64`: The z-component of the jet momentum.
-- `E::Float64`: The energy of the jet.
+- `px::T`: The x-component of the jet momentum.
+- `py::T`: The y-component of the jet momentum.
+- `pz::T`: The z-component of the jet momentum.
+- `E::T`: The energy of the jet.
 - `_cluster_hist_index::Int`: The index of the cluster histogram.
-- `_p2::Float64`: The squared momentum of the jet.
-- `_inv_p::Float64`: The inverse momentum of the jet.
+- `_p2::T`: The squared momentum of the jet.
+- `_inv_p::T`: The inverse momentum of the jet.
+
+# Type Parameters
+- `T <: Real`: The type of the numerical values.
 """
-mutable struct EEjet <: FourMomentum
-    px::Float64
-    py::Float64
-    pz::Float64
-    E::Float64
-    _p2::Float64
-    _inv_p::Float64
+mutable struct EEjet{T <: Real} <: FourMomentum
+    px::T
+    py::T
+    pz::T
+    E::T
+    _p2::T
+    _inv_p::T
     _cluster_hist_index::Int
 end
 
-function EEjet(px::Real, py::Real, pz::Real, E::Real, _cluster_hist_index::Int)
+function EEjet(px::T, py::T, pz::T, E::T, _cluster_hist_index::Int) where {T <: Real}
     @muladd p2 = px * px + py * py + pz * pz
     inv_p = @fastmath 1.0 / sqrt(p2)
-    EEjet(px, py, pz, E, p2, inv_p, _cluster_hist_index)
+    EEjet{T}(px, py, pz, E, p2, inv_p, _cluster_hist_index)
 end
 
-EEjet(px::Real, py::Real, pz::Real, E::Real) = EEjet(px, py, pz, E, 0)
+EEjet(px::T, py::T, pz::T, E::T) where {T <: Real} = EEjet{T}(px, py, pz, E, 0)
 
 EEjet(pj::PseudoJet) = EEjet(px(pj), py(pj), pz(pj), energy(pj), cluster_hist_index(pj))
 
@@ -87,15 +92,31 @@ function show(io::IO, eej::EEjet)
           " cluster_hist_index: ", eej._cluster_hist_index, ")")
 end
 
-# Optimised reconstruction struct for e+e jets
+"""
+    mutable struct EERecoJet{T <: Real}
+
+Optimised struct for e+e jets reconstruction, to be used with StructArrays.
 
-mutable struct EERecoJet
+# Fields
+- `index::Int`: The index of the jet.
+- `nni::Int`: The nearest neighbour index.
+- `nndist::T`: The distance to the nearest neighbour.
+- `dijdist::T`: The distance between jets.
+- `nx::T`: The x-component of the jet's momentum.
+- `ny::T`: The y-component of the jet's momentum.
+- `nz::T`: The z-component of the jet's momentum.
+- `E2p::T`: The energy-to-momentum ratio of the jet.
+
+# Type Parameters
+- `T <: Real`: The type of the numerical values.
+"""
+mutable struct EERecoJet{T <: Real}
     index::Int
     nni::Int
-    nndist::Float64
-    dijdist::Float64
-    nx::Float64
-    ny::Float64
-    nz::Float64
-    E2p::Float64
+    nndist::T
+    dijdist::T
+    nx::T
+    ny::T
+    nz::T
+    E2p::T
 end

From b00d370332884d2da7b2754452d9416f04bb475d Mon Sep 17 00:00:00 2001
From: Graeme A Stewart <graeme.andrew.stewart@cern.ch>
Date: Mon, 25 Nov 2024 14:10:34 +0100
Subject: [PATCH 03/16] Minor corrections to docstrings

---
 src/EEjet.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/EEjet.jl b/src/EEjet.jl
index c9b2c41..b673c01 100644
--- a/src/EEjet.jl
+++ b/src/EEjet.jl
@@ -2,7 +2,7 @@
     struct EEjet{T <: Real} <: FourMomentum
 
 The `EEjet` struct is a 4-momentum object used for the e+e jet reconstruction
-routines. Internal fields are used to track the reconstuction and to cache
+routines. Internal fields are used to track the reconstruction and to cache
 values needed during the execution of the algorithm.
 
 # Fields
@@ -105,7 +105,7 @@ Optimised struct for e+e jets reconstruction, to be used with StructArrays.
 - `nx::T`: The x-component of the jet's momentum.
 - `ny::T`: The y-component of the jet's momentum.
 - `nz::T`: The z-component of the jet's momentum.
-- `E2p::T`: The energy-to-momentum ratio of the jet.
+- `E2p::T`: The energy raised to the power of 2p for this jet.
 
 # Type Parameters
 - `T <: Real`: The type of the numerical values.

From 4022fb9a852445d7bfdf7e1d8317777a2b152cc8 Mon Sep 17 00:00:00 2001
From: Graeme A Stewart <graeme.andrew.stewart@cern.ch>
Date: Mon, 25 Nov 2024 14:48:09 +0100
Subject: [PATCH 04/16] Fix types for EEJet

---
 src/EEAlgorithm.jl | 5 ++++-
 src/EEjet.jl       | 4 ++--
 2 files changed, 6 insertions(+), 3 deletions(-)

diff --git a/src/EEAlgorithm.jl b/src/EEAlgorithm.jl
index d9c502b..7c317e6 100644
--- a/src/EEAlgorithm.jl
+++ b/src/EEAlgorithm.jl
@@ -273,6 +273,9 @@ function _ee_genkt_algorithm(; particles::Vector{EEjet}, p = 1, R = 4.0,
     # R squared
     R2 = R^2
 
+    # Numerical type?
+    ParticleType = typeof(particles[1].E)
+
     # Constant factor for the dij metric and the beam distance function
     if algorithm == JetAlgorithm.Durham
         dij_factor = 2.0
@@ -289,7 +292,7 @@ function _ee_genkt_algorithm(; particles::Vector{EEjet}, p = 1, R = 4.0,
     # For optimised reconstruction generate an SoA containing the necessary
     # jet information and populate it accordingly
     # We need N slots for this array
-    eereco = StructArray{EERecoJet}(undef, N)
+    eereco = StructArray{EERecoJet{ParticleType}}(undef, N)
     fill_reco_array!(eereco, particles, R2, p)
 
     # Setup the initial history and get the total energy
diff --git a/src/EEjet.jl b/src/EEjet.jl
index b673c01..3d3bb68 100644
--- a/src/EEjet.jl
+++ b/src/EEjet.jl
@@ -27,13 +27,13 @@ mutable struct EEjet{T <: Real} <: FourMomentum
     _cluster_hist_index::Int
 end
 
-function EEjet(px::T, py::T, pz::T, E::T, _cluster_hist_index::Int) where {T <: Real}
+function EEjet(px::T, py::T, pz::T, E::T, _cluster_hist_index::Integer) where {T <: Real}
     @muladd p2 = px * px + py * py + pz * pz
     inv_p = @fastmath 1.0 / sqrt(p2)
     EEjet{T}(px, py, pz, E, p2, inv_p, _cluster_hist_index)
 end
 
-EEjet(px::T, py::T, pz::T, E::T) where {T <: Real} = EEjet{T}(px, py, pz, E, 0)
+EEjet(px::T, py::T, pz::T, E::T) where {T <: Real} = EEjet(px, py, pz, E, 0)
 
 EEjet(pj::PseudoJet) = EEjet(px(pj), py(pj), pz(pj), energy(pj), cluster_hist_index(pj))
 

From cc365b9f77222ee87c67d72faa84af3cb5458f91 Mon Sep 17 00:00:00 2001
From: Graeme A Stewart <graeme.andrew.stewart@cern.ch>
Date: Mon, 25 Nov 2024 22:07:12 +0100
Subject: [PATCH 05/16] Fix types for EE reco

Ensure that types are stable for EEReco, meaning that they are
EEjet{T} for a specific T, not just generic EEjet
---
 src/ClusterSequence.jl |  6 ++++++
 src/EEAlgorithm.jl     | 15 ++++++++++-----
 src/EEjet.jl           |  1 +
 3 files changed, 17 insertions(+), 5 deletions(-)

diff --git a/src/ClusterSequence.jl b/src/ClusterSequence.jl
index 9bdf653..6b46bb5 100644
--- a/src/ClusterSequence.jl
+++ b/src/ClusterSequence.jl
@@ -148,6 +148,12 @@ ClusterSequence(algorithm::JetAlgorithm.Algorithm, p::Real, R::Float64, strategy
                        Qtot)
 end
 
+ClusterSequence{T}(algorithm::JetAlgorithm.Algorithm, p::Real, R::Float64, strategy::RecoStrategy.Strategy, jets::Vector{T}, history, Qtot) where {T <: FourMomentum} = begin
+    ClusterSequence{T}(algorithm, Float64(p), R, strategy, jets, length(jets), history,
+                       Qtot)
+end
+
+
 """
     add_step_to_history!(clusterseq::ClusterSequence, parent1, parent2, jetp_index, dij)
 
diff --git a/src/EEAlgorithm.jl b/src/EEAlgorithm.jl
index 7c317e6..61a0b3b 100644
--- a/src/EEAlgorithm.jl
+++ b/src/EEAlgorithm.jl
@@ -3,6 +3,8 @@
 const large_distance = 16.0 # = 4^2
 const large_dij = 1.0e6
 
+using Infiltrator
+
 """
     angular_distance(eereco, i, j) -> Float64
 
@@ -242,11 +244,12 @@ function ee_genkt_algorithm(particles::AbstractArray{T, 1}; p = 1, R = 4.0,
         recombination_particles = copy(particles)
         sizehint!(recombination_particles, length(particles) * 2)
     else
-        recombination_particles = EEjet[]
+        ParticleType = typeof(px(particles[1]))
+        recombination_particles = EEjet{ParticleType}[]
         sizehint!(recombination_particles, length(particles) * 2)
         for i in eachindex(particles)
             push!(recombination_particles,
-                  EEjet(px(particles[i]), py(particles[i]), pz(particles[i]),
+                  EEjet{ParticleType}(px(particles[i]), py(particles[i]), pz(particles[i]),
                         energy(particles[i])))
         end
     end
@@ -264,9 +267,9 @@ end
 
 This function is the actual implementation of the e+e- jet clustering algorithm.
 """
-function _ee_genkt_algorithm(; particles::Vector{EEjet}, p = 1, R = 4.0,
+function _ee_genkt_algorithm(; particles::Vector{EEjet{T}}, p = 1, R = 4.0,
                              algorithm::JetAlgorithm.Algorithm = JetAlgorithm.Durham,
-                             recombine = +)
+                             recombine = +) where {T <: Real}
     # Bounds
     N::Int = length(particles)
 
@@ -298,7 +301,7 @@ function _ee_genkt_algorithm(; particles::Vector{EEjet}, p = 1, R = 4.0,
     # Setup the initial history and get the total energy
     history, Qtot = initial_history(particles)
 
-    clusterseq = ClusterSequence(algorithm, p, R, RecoStrategy.N2Plain, particles, history,
+    clusterseq = ClusterSequence{EEjet{ParticleType}}(algorithm, p, R, RecoStrategy.N2Plain, particles, history,
                                  Qtot)
 
     # Run over initial pairs of jets to find nearest neighbours
@@ -307,6 +310,8 @@ function _ee_genkt_algorithm(; particles::Vector{EEjet}, p = 1, R = 4.0,
     # Only for debugging purposes...
     # ee_check_consistency(clusterseq, clusterseq_index, N, nndist, nndij, nni, "Start")
 
+    #@infiltrate
+
     # Now we can start the main loop
     iter = 0
     while N != 0
diff --git a/src/EEjet.jl b/src/EEjet.jl
index 3d3bb68..9b945ff 100644
--- a/src/EEjet.jl
+++ b/src/EEjet.jl
@@ -34,6 +34,7 @@ function EEjet(px::T, py::T, pz::T, E::T, _cluster_hist_index::Integer) where {T
 end
 
 EEjet(px::T, py::T, pz::T, E::T) where {T <: Real} = EEjet(px, py, pz, E, 0)
+EEjet{T}(px::T, py::T, pz::T, E::T) where {T <: Real} = EEjet(px, py, pz, E, 0)
 
 EEjet(pj::PseudoJet) = EEjet(px(pj), py(pj), pz(pj), energy(pj), cluster_hist_index(pj))
 

From b64a483316235533d5b6bf777df05a5a6de1c626 Mon Sep 17 00:00:00 2001
From: Graeme A Stewart <graeme.andrew.stewart@cern.ch>
Date: Tue, 26 Nov 2024 08:29:41 +0100
Subject: [PATCH 06/16] Fixes for type variation from CLI

Allow instrumented-jetreco to have a type argument that makes the particle
collection have that specific type, allowing for switches to different Floats

N.B, Currently Float32 is slower than Float64
---
 examples/instrumented-jetreco.jl | 9 +++++++--
 src/EEjet.jl                     | 5 ++++-
 2 files changed, 11 insertions(+), 3 deletions(-)

diff --git a/examples/instrumented-jetreco.jl b/examples/instrumented-jetreco.jl
index df471e7..dc46672 100644
--- a/examples/instrumented-jetreco.jl
+++ b/examples/instrumented-jetreco.jl
@@ -244,6 +244,11 @@ function parse_command_line(args)
         arg_type = RecoStrategy.Strategy
         default = RecoStrategy.Best
 
+        "--type", "-T"
+        help = """Numerical type to use for the reconstruction (Float32, Float64)"""
+        arg_type = Symbol
+        default = :Float64
+
         "--nsamples", "-m"
         help = "Number of measurement points to acquire."
         arg_type = Int
@@ -296,9 +301,9 @@ function main()
     global_logger(logger)
     # Try to read events into the correct type!
     if JetReconstruction.is_ee(args[:algorithm])
-        jet_type = EEjet
+        jet_type = EEjet{eval(args[:type])}
     else
-        jet_type = PseudoJet
+        jet_type = PseudoJet{eval(args[:type])}
     end
     events::Vector{Vector{jet_type}} = read_final_state_particles(args[:file],
                                                                   maxevents = args[:maxevents],
diff --git a/src/EEjet.jl b/src/EEjet.jl
index 9b945ff..ccadba7 100644
--- a/src/EEjet.jl
+++ b/src/EEjet.jl
@@ -33,8 +33,11 @@ function EEjet(px::T, py::T, pz::T, E::T, _cluster_hist_index::Integer) where {T
     EEjet{T}(px, py, pz, E, p2, inv_p, _cluster_hist_index)
 end
 
+# Constructor with type T passes through without history index
 EEjet(px::T, py::T, pz::T, E::T) where {T <: Real} = EEjet(px, py, pz, E, 0)
-EEjet{T}(px::T, py::T, pz::T, E::T) where {T <: Real} = EEjet(px, py, pz, E, 0)
+
+# Constructor with type U does type conversion before initialising
+EEjet{U}(px::T, py::T, pz::T, E::T) where {T <: Real, U <: Real} = EEjet(U(px), U(py), U(pz), U(E), 0)
 
 EEjet(pj::PseudoJet) = EEjet(px(pj), py(pj), pz(pj), energy(pj), cluster_hist_index(pj))
 

From dc58c70bfb4d178597e45a7d6cd3358f8d91151b Mon Sep 17 00:00:00 2001
From: Graeme A Stewart <graeme.andrew.stewart@cern.ch>
Date: Tue, 26 Nov 2024 08:48:46 +0100
Subject: [PATCH 07/16] Remove Infiltrator temporary

---
 Project.toml       | 4 ++--
 src/EEAlgorithm.jl | 4 ----
 2 files changed, 2 insertions(+), 6 deletions(-)

diff --git a/Project.toml b/Project.toml
index 3b2bab7..52cc055 100644
--- a/Project.toml
+++ b/Project.toml
@@ -16,12 +16,12 @@ MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 StructArrays = "09ab397b-f2b6-538f-b94a-2f83cf4a842a"
 
 [weakdeps]
-Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 EDM4hep = "eb32b910-dde9-4347-8fce-cd6be3498f0c"
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 
 [extensions]
-JetVisualisation = "Makie"
 EDM4hepJets = "EDM4hep"
+JetVisualisation = "Makie"
 
 [compat]
 Accessors = "0.1.36"
diff --git a/src/EEAlgorithm.jl b/src/EEAlgorithm.jl
index 61a0b3b..7142940 100644
--- a/src/EEAlgorithm.jl
+++ b/src/EEAlgorithm.jl
@@ -3,8 +3,6 @@
 const large_distance = 16.0 # = 4^2
 const large_dij = 1.0e6
 
-using Infiltrator
-
 """
     angular_distance(eereco, i, j) -> Float64
 
@@ -310,8 +308,6 @@ function _ee_genkt_algorithm(; particles::Vector{EEjet{T}}, p = 1, R = 4.0,
     # Only for debugging purposes...
     # ee_check_consistency(clusterseq, clusterseq_index, N, nndist, nndij, nni, "Start")
 
-    #@infiltrate
-
     # Now we can start the main loop
     iter = 0
     while N != 0

From d5f47e37d61ffbe8468243ef45a4301455bcac97 Mon Sep 17 00:00:00 2001
From: Graeme A Stewart <graeme.andrew.stewart@cern.ch>
Date: Tue, 26 Nov 2024 08:49:50 +0100
Subject: [PATCH 08/16] Fixes for pp reconstruction to use typed PseudoJets

It is very important to pass the fully typed PseudoJet{T} around, not a generic PseudoJet
---
 src/PlainAlgo.jl   | 24 ++++++++++++++----------
 src/Pseudojet.jl   | 17 +++++++++++++----
 src/TiledAlgoLL.jl | 20 ++++++++++++--------
 3 files changed, 39 insertions(+), 22 deletions(-)

diff --git a/src/PlainAlgo.jl b/src/PlainAlgo.jl
index fabcb91..6bb9681 100644
--- a/src/PlainAlgo.jl
+++ b/src/PlainAlgo.jl
@@ -230,13 +230,14 @@ function plain_jet_reconstruct(particles::AbstractArray{T, 1}; p::Union{Real, No
     # Integer p if possible
     p = (round(p) == p) ? Int(p) : p
 
-    if T == PseudoJet
+    if T isa PseudoJet
         # recombination_particles will become part of the cluster sequence, so size it for
         # the starting particles and all N recombinations
         recombination_particles = copy(particles)
         sizehint!(recombination_particles, length(particles) * 2)
     else
-        recombination_particles = PseudoJet[]
+        ParticleType = typeof(px(particles[1]))
+        recombination_particles = PseudoJet{ParticleType}[]
         sizehint!(recombination_particles, length(particles) * 2)
         for i in eachindex(particles)
             push!(recombination_particles,
@@ -278,23 +279,26 @@ generalised k_t algorithm.
 - `clusterseq`: The resulting `ClusterSequence` object representing the
   reconstructed jets.
 """
-function _plain_jet_reconstruct(; particles::Vector{PseudoJet}, p = -1, R = 1.0,
+function _plain_jet_reconstruct(; particles::Vector{PseudoJet{T}}, p = -1, R = 1.0,
                                 algorithm::JetAlgorithm.Algorithm = JetAlgorithm.AntiKt,
-                                recombine = +)
+                                recombine = +) where T <: Real
     # Bounds
     N::Int = length(particles)
     # Parameters
     R2 = R^2
 
+    # Numerical type?
+    ParticleType = typeof(particles[1].E)
+
     # Optimised compact arrays for determining the next merge step
     # We make sure these arrays are type stable - have seen issues where, depending on the values
     # returned by the methods, they can become unstable and performance degrades
-    kt2_array::Vector{Float64} = pt2.(particles) .^ p
-    phi_array::Vector{Float64} = phi.(particles)
-    rapidity_array::Vector{Float64} = rapidity.(particles)
+    kt2_array::Vector{ParticleType} = pt2.(particles) .^ p
+    phi_array::Vector{ParticleType} = phi.(particles)
+    rapidity_array::Vector{ParticleType} = rapidity.(particles)
     nn::Vector{Int} = Vector(1:N) # nearest neighbours
-    nndist::Vector{Float64} = fill(float(R2), N) # geometric distances to the nearest neighbour
-    nndij::Vector{Float64} = zeros(N) # dij metric distance
+    nndist::Vector{ParticleType} = fill(float(R2), N) # geometric distances to the nearest neighbour
+    nndij::Vector{ParticleType} = zeros(N) # dij metric distance
 
     # Maps index from the compact array to the clusterseq jet vector
     clusterseq_index::Vector{Int} = collect(1:N)
@@ -304,7 +308,7 @@ function _plain_jet_reconstruct(; particles::Vector{PseudoJet}, p = -1, R = 1.0,
     # Current implementation mutates the particles vector, so need to copy it
     # for the cluster sequence (there is too much copying happening, so this
     # needs to be rethought and reoptimised)
-    clusterseq = ClusterSequence(algorithm, p, R, RecoStrategy.N2Plain, particles, history,
+    clusterseq = ClusterSequence{PseudoJet{ParticleType}}(algorithm, p, R, RecoStrategy.N2Plain, particles, history,
                                  Qtot)
 
     # Initialize nearest neighbours
diff --git a/src/Pseudojet.jl b/src/Pseudojet.jl
index 0d52302..d47f2df 100644
--- a/src/Pseudojet.jl
+++ b/src/Pseudojet.jl
@@ -69,10 +69,16 @@ history index.
 A `PseudoJet` object.
 """
 PseudoJet(px::T, py::T, pz::T, E::T,
-_cluster_hist_index::Int,
-pt2::T) where {T <: Real} = PseudoJet{T}(px,
-                                         py, pz, E, _cluster_hist_index,
-                                         pt2, 1.0 / pt2, _invalid_rap, _invalid_phi)
+    _cluster_hist_index::Int,
+    pt2::T) where {T <: Real} = PseudoJet{T}(px,
+                                            py, pz, E, _cluster_hist_index,
+                                             pt2, 1.0 / pt2, _invalid_rap, _invalid_phi)
+
+PseudoJet{T}(px::T, py::T, pz::T, E::T,
+            _cluster_hist_index::Int,
+            pt2::T) where {T <: Real} = PseudoJet{T}(px,
+                                                    py, pz, E, _cluster_hist_index,
+                                                    pt2, 1.0 / pt2, _invalid_rap, _invalid_phi)
 
 """
     PseudoJet(px::T, py::T, pz::T, E::T) where T <: Real
@@ -90,6 +96,9 @@ A PseudoJet object.
 """
 PseudoJet(px::T, py::T,
 pz::T, E::T) where {T <: Real} = PseudoJet(px, py, pz, E, 0, px^2 + py^2)
+PseudoJet{T}(px::T, py::T,
+pz::T, E::T) where {T <: Real} = PseudoJet{T}(px, py, pz, E, 0, px^2 + py^2)
+
 
 import Base.show
 """
diff --git a/src/TiledAlgoLL.jl b/src/TiledAlgoLL.jl
index ff3b470..df6390e 100644
--- a/src/TiledAlgoLL.jl
+++ b/src/TiledAlgoLL.jl
@@ -374,17 +374,18 @@ function tiled_jet_reconstruct(particles::AbstractArray{T, 1}; p::Union{Real, No
     (p, algorithm) = get_algorithm_power_consistency(p = p, algorithm = algorithm)
 
     # If we have PseudoJets, we can just call the main algorithm...
-    if T == PseudoJet
+    if T isa PseudoJet
         # recombination_particles will become part of the cluster sequence, so size it for
         # the starting particles and all N recombinations
         recombination_particles = copy(particles)
         sizehint!(recombination_particles, length(particles) * 2)
     else
-        recombination_particles = PseudoJet[]
+        ParticleType = typeof(px(particles[1]))
+        recombination_particles = PseudoJet{ParticleType}[]
         sizehint!(recombination_particles, length(particles) * 2)
         for i in eachindex(particles)
             push!(recombination_particles,
-                  PseudoJet(px(particles[i]), py(particles[i]), pz(particles[i]),
+                  PseudoJet{ParticleType}(px(particles[i]), py(particles[i]), pz(particles[i]),
                             energy(particles[i])))
         end
     end
@@ -421,9 +422,9 @@ of data types are done.
 tiled_jet_reconstruct(particles::Vector{PseudoJet}; p = 1, R = 1.0, recombine = +)
 ```
 """
-function _tiled_jet_reconstruct(particles::Vector{PseudoJet}; p::Real = -1, R = 1.0,
+function _tiled_jet_reconstruct(particles::Vector{PseudoJet{T}}; p::Real = -1, R = 1.0,
                                 algorithm::JetAlgorithm.Algorithm = JetAlgorithm.AntiKt,
-                                recombine = +)
+                                recombine = +) where T <: Real
     # Bounds
     N::Int = length(particles)
 
@@ -439,13 +440,16 @@ function _tiled_jet_reconstruct(particles::Vector{PseudoJet}; p::Real = -1, R =
     R2::Float64 = R * R
     p = (round(p) == p) ? Int(p) : p # integer p if possible
 
+    # Numerical type?
+    ParticleType = typeof(particles[1].E)
+
     # This will be used quite deep inside loops, so declare it here so that
     # memory (de)allocation gets done only once
     tile_union = Vector{Int}(undef, 3 * _n_tile_neighbours)
 
     # Container for pseudojets, sized for all initial particles, plus all of the
     # merged jets that can be created during reconstruction
-    jets = PseudoJet[]
+    jets = PseudoJet{ParticleType}[]
     sizehint!(jets, N * 2)
     resize!(jets, N)
 
@@ -456,7 +460,7 @@ function _tiled_jet_reconstruct(particles::Vector{PseudoJet}; p::Real = -1, R =
     history, Qtot = initial_history(jets)
 
     # Now get the tiling setup
-    _eta = Vector{Float64}(undef, length(particles))
+    _eta = Vector{ParticleType}(undef, length(particles))
     for ijet in 1:length(particles)
         _eta[ijet] = rapidity(particles[ijet])
     end
@@ -464,7 +468,7 @@ function _tiled_jet_reconstruct(particles::Vector{PseudoJet}; p::Real = -1, R =
     tiling = Tiling(setup_tiling(_eta, R))
 
     # ClusterSequence is the struct that holds the state of the reconstruction
-    clusterseq = ClusterSequence(algorithm, p, R, RecoStrategy.N2Tiled, jets, history, Qtot)
+    clusterseq = ClusterSequence{PseudoJet{ParticleType}}(algorithm, p, R, RecoStrategy.N2Tiled, jets, history, Qtot)
 
     # Tiled jets is a structure that has additional variables for tracking which tile a jet is in
     tiledjets = similar(clusterseq.jets, TiledJet)

From c51faa1faec10096b4eb3cfa0a4418df132eb018 Mon Sep 17 00:00:00 2001
From: Graeme A Stewart <graeme.andrew.stewart@cern.ch>
Date: Tue, 26 Nov 2024 10:58:26 +0100
Subject: [PATCH 09/16] Add quick little benchmark script

---
 examples/benchmark.sh | 11 +++++++++++
 1 file changed, 11 insertions(+)
 create mode 100755 examples/benchmark.sh

diff --git a/examples/benchmark.sh b/examples/benchmark.sh
new file mode 100755
index 0000000..1bf6c50
--- /dev/null
+++ b/examples/benchmark.sh
@@ -0,0 +1,11 @@
+#! /bin/sh
+#
+# Quick and dirty set of benchmarks for the most important cases
+echo "pp 14TeV Tiled"
+julia --project instrumented-jetreco.jl --algorithm=AntiKt -R 0.4 ../test/data/events.pp13TeV.hepmc3.gz -S N2Tiled -m 16
+
+echo "pp 14 TeV Plain"
+julia --project instrumented-jetreco.jl --algorithm=AntiKt -R 0.4 ../test/data/events.pp13TeV.hepmc3.gz -S N2Plain -m 16
+
+echo "ee H Durham"
+julia --project instrumented-jetreco.jl --algorithm=Durham ../test/data/events.eeH.hepmc3.gz -m 16

From a42e19ef1e7e3d366d2f08cf7531f235e5fa6811 Mon Sep 17 00:00:00 2001
From: Graeme A Stewart <graeme.andrew.stewart@cern.ch>
Date: Tue, 26 Nov 2024 11:08:47 +0100
Subject: [PATCH 10/16] Define eltype() for PseudoJet and EEjet

Returns the parametrised element type
---
 src/EEAlgorithm.jl | 2 +-
 src/EEjet.jl       | 8 ++++++++
 src/PlainAlgo.jl   | 4 ++--
 src/Pseudojet.jl   | 7 +++++++
 src/TiledAlgoLL.jl | 4 ++--
 5 files changed, 20 insertions(+), 5 deletions(-)

diff --git a/src/EEAlgorithm.jl b/src/EEAlgorithm.jl
index 7142940..589578e 100644
--- a/src/EEAlgorithm.jl
+++ b/src/EEAlgorithm.jl
@@ -275,7 +275,7 @@ function _ee_genkt_algorithm(; particles::Vector{EEjet{T}}, p = 1, R = 4.0,
     R2 = R^2
 
     # Numerical type?
-    ParticleType = typeof(particles[1].E)
+    ParticleType = eltype(particles[1])
 
     # Constant factor for the dij metric and the beam distance function
     if algorithm == JetAlgorithm.Durham
diff --git a/src/EEjet.jl b/src/EEjet.jl
index ccadba7..112edff 100644
--- a/src/EEjet.jl
+++ b/src/EEjet.jl
@@ -27,6 +27,14 @@ mutable struct EEjet{T <: Real} <: FourMomentum
     _cluster_hist_index::Int
 end
 
+
+"""
+    Base.eltype(::Type{EEjet{T}}) where T
+
+Return the element type of the `EEjet` struct.
+"""
+Base.eltype(::Type{EEjet{T}}) where T = T
+
 function EEjet(px::T, py::T, pz::T, E::T, _cluster_hist_index::Integer) where {T <: Real}
     @muladd p2 = px * px + py * py + pz * pz
     inv_p = @fastmath 1.0 / sqrt(p2)
diff --git a/src/PlainAlgo.jl b/src/PlainAlgo.jl
index 6bb9681..cbc396d 100644
--- a/src/PlainAlgo.jl
+++ b/src/PlainAlgo.jl
@@ -287,8 +287,8 @@ function _plain_jet_reconstruct(; particles::Vector{PseudoJet{T}}, p = -1, R = 1
     # Parameters
     R2 = R^2
 
-    # Numerical type?
-    ParticleType = typeof(particles[1].E)
+    # Numerical type for this reconstruction
+    ParticleType = eltype(particles[1])
 
     # Optimised compact arrays for determining the next merge step
     # We make sure these arrays are type stable - have seen issues where, depending on the values
diff --git a/src/Pseudojet.jl b/src/Pseudojet.jl
index d47f2df..8f87511 100644
--- a/src/Pseudojet.jl
+++ b/src/Pseudojet.jl
@@ -49,6 +49,13 @@ mutable struct PseudoJet{T <: Real} <: FourMomentum
     _phi::T
 end
 
+"""
+    Base.eltype(::Type{PseudoJet{T}}) where T
+
+Return the element type of the `PseudoJet` struct.
+"""
+Base.eltype(::Type{PseudoJet{T}}) where T = T
+
 """
     PseudoJet(px::T, py::T, pz::T, E::T,
         _cluster_hist_index::Int,
diff --git a/src/TiledAlgoLL.jl b/src/TiledAlgoLL.jl
index df6390e..34c53de 100644
--- a/src/TiledAlgoLL.jl
+++ b/src/TiledAlgoLL.jl
@@ -440,8 +440,8 @@ function _tiled_jet_reconstruct(particles::Vector{PseudoJet{T}}; p::Real = -1, R
     R2::Float64 = R * R
     p = (round(p) == p) ? Int(p) : p # integer p if possible
 
-    # Numerical type?
-    ParticleType = typeof(particles[1].E)
+    # Numerical type for this reconstruction
+    ParticleType = eltype(particles[1])
 
     # This will be used quite deep inside loops, so declare it here so that
     # memory (de)allocation gets done only once

From 53ca16b84ff0a3b292fab6902f847d6def0c75f2 Mon Sep 17 00:00:00 2001
From: Graeme A Stewart <graeme.andrew.stewart@cern.ch>
Date: Tue, 26 Nov 2024 11:30:19 +0100
Subject: [PATCH 11/16] Add more documentation for constructors

Also support type switching constructor for PseudoJet
---
 src/EEjet.jl     | 78 ++++++++++++++++++++++++++++++++++++++++++++++--
 src/Pseudojet.jl | 29 +++++++++++++++---
 2 files changed, 101 insertions(+), 6 deletions(-)

diff --git a/src/EEjet.jl b/src/EEjet.jl
index 112edff..9f6ec68 100644
--- a/src/EEjet.jl
+++ b/src/EEjet.jl
@@ -35,20 +35,94 @@ Return the element type of the `EEjet` struct.
 """
 Base.eltype(::Type{EEjet{T}}) where T = T
 
+
+"""
+    EEjet(px::T, py::T, pz::T, E::T, _cluster_hist_index::Integer) where {T <: Real}
+
+Constructs an `EEjet` object with the given momentum components `px`, `py`,
+`pz`, energy `E`, and cluster histogram index `_cluster_hist_index`.
+
+The constructed EEjet object will be parametrised by the type `T`.
+
+# Arguments
+- `px::T`: The x-component of the momentum.
+- `py::T`: The y-component of the momentum.
+- `pz::T`: The z-component of the momentum.
+- `E::T`: The energy of the jet.
+- `_cluster_hist_index::Integer`: The index of the cluster histogram.
+
+# Returns
+- The initialised `EEjet` object.
+
+# Note
+- `T` must be a subtype of `Real`.
+- The `@muladd` macro is used to perform fused multiply-add operations for
+  computing `p2`.
+- The `@fastmath` macro is used to allow the compiler to perform optimizations
+  for computing `inv_p`.
+"""
 function EEjet(px::T, py::T, pz::T, E::T, _cluster_hist_index::Integer) where {T <: Real}
     @muladd p2 = px * px + py * py + pz * pz
     inv_p = @fastmath 1.0 / sqrt(p2)
     EEjet{T}(px, py, pz, E, p2, inv_p, _cluster_hist_index)
 end
 
-# Constructor with type T passes through without history index
+"""
+    EEjet(px::T, py::T, pz::T, E::T) where {T <: Real}
+
+Constructs an `EEjet` object with the given momentum components `px`, `py`,
+`pz`, energy `E`, and the cluster histogram index set to zero.
+
+The constructed EEjet object will be parametrised by the type `T`.
+
+# Arguments
+- `px::T`: The x-component of the momentum.
+- `py::T`: The y-component of the momentum.
+- `pz::T`: The z-component of the momentum.
+- `E::T`: The energy of the jet.
+
+# Returns
+- The initialised `EEjet` object.
+"""
 EEjet(px::T, py::T, pz::T, E::T) where {T <: Real} = EEjet(px, py, pz, E, 0)
 
-# Constructor with type U does type conversion before initialising
+"""
+    EEjet{U}(px::T, py::T, pz::T, E::T) where {T <: Real, U <: Real}
+
+Constructs an `EEjet` object with conversion of the given momentum components
+(`px`, `py`, `pz`) and energy (`E`) from type `T` to type `U`.
+
+# Arguments
+- `px::T`: The x-component of the momentum.
+- `py::T`: The y-component of the momentum.
+- `pz::T`: The z-component of the momentum.
+- `E::T`: The energy.
+
+# Type Parameters
+- `T <: Real`: The type of the input momentum components and energy.
+- `U <: Real`: The type to which the input values will be converted
+
+# Returns
+An `EEjet` object with the momentum components and energy parametrised to type
+`U`.
+"""
 EEjet{U}(px::T, py::T, pz::T, E::T) where {T <: Real, U <: Real} = EEjet(U(px), U(py), U(pz), U(E), 0)
 
+
+"""
+    EEjet(pj::PseudoJet) -> EEjet
+
+Constructs an `EEjet` object from a given `PseudoJet` object `pj`.
+
+# Arguments
+- `pj::PseudoJet`: A `PseudoJet` object used to create the `EEjet`.
+
+# Returns
+- An `EEjet` object initialized with the same properties of the given `PseudoJet`.
+"""
 EEjet(pj::PseudoJet) = EEjet(px(pj), py(pj), pz(pj), energy(pj), cluster_hist_index(pj))
 
+
 p2(eej::EEjet) = eej._p2
 pt2(eej::EEjet) = eej.px^2 + eej.py^2
 const kt2 = pt2
diff --git a/src/Pseudojet.jl b/src/Pseudojet.jl
index 8f87511..b3a5a27 100644
--- a/src/Pseudojet.jl
+++ b/src/Pseudojet.jl
@@ -35,7 +35,6 @@ caching of the more expensive calculations for rapidity and azimuthal angle.
 - `_inv_pt2::T`: The inverse squared transverse momentum.
 - `_rap::T`: The rapidity.
 - `_phi::T`: The azimuthal angle.
-
 """
 mutable struct PseudoJet{T <: Real} <: FourMomentum
     px::T
@@ -81,6 +80,25 @@ PseudoJet(px::T, py::T, pz::T, E::T,
                                             py, pz, E, _cluster_hist_index,
                                              pt2, 1.0 / pt2, _invalid_rap, _invalid_phi)
 
+"""
+    PseudoJet{T}(px::T, py::T, pz::T, E::T,
+        _cluster_hist_index::Int,
+        pt2::T) where {T <: Real}
+
+Constructs a parametrised PseudoJet object with the given momentum components
+and energy and history index.
+
+# Arguments
+- `px::T`: The x-component of the momentum.
+- `py::T`: The y-component of the momentum.
+- `pz::T`: The z-component of the momentum.
+- `E::T`: The energy.
+- `_cluster_hist_index::Int`: The cluster history index.
+- `pt2::T`: The transverse momentum squared.
+
+# Returns
+A `PseudoJet` object.
+"""
 PseudoJet{T}(px::T, py::T, pz::T, E::T,
             _cluster_hist_index::Int,
             pt2::T) where {T <: Real} = PseudoJet{T}(px,
@@ -99,12 +117,15 @@ Constructs a PseudoJet object with the given momentum components and energy.
 - `E::T`: The energy.
 
 # Returns
-A PseudoJet object.
+A PseudoJet object, with type parameter `T`.
 """
 PseudoJet(px::T, py::T,
 pz::T, E::T) where {T <: Real} = PseudoJet(px, py, pz, E, 0, px^2 + py^2)
-PseudoJet{T}(px::T, py::T,
-pz::T, E::T) where {T <: Real} = PseudoJet{T}(px, py, pz, E, 0, px^2 + py^2)
+
+"""
+    PseudoJet{U}(px::T, py::T, pz::T, E::T) where T <: Real, U <: Real
+"""
+PseudoJet{U}(px::T, py::T, pz::T, E::T) where {T <: Real, U <: Real} = PseudoJet{U}(U(px), U(py), U(pz), U(E), 0, U(px^2 + py^2))
 
 
 import Base.show

From 3bbe3fbc3b3d68bc787a365b9679159b55f5e93e Mon Sep 17 00:00:00 2001
From: Graeme A Stewart <graeme.andrew.stewart@cern.ch>
Date: Tue, 26 Nov 2024 11:53:45 +0100
Subject: [PATCH 12/16] Formatting!

---
 src/ClusterSequence.jl |  5 +++--
 src/EEAlgorithm.jl     |  7 ++++---
 src/EEjet.jl           |  9 +++------
 src/PlainAlgo.jl       |  8 +++++---
 src/Pseudojet.jl       | 26 +++++++++++++++-----------
 src/TiledAlgoLL.jl     | 11 +++++++----
 6 files changed, 37 insertions(+), 29 deletions(-)

diff --git a/src/ClusterSequence.jl b/src/ClusterSequence.jl
index 6b46bb5..a6d0515 100644
--- a/src/ClusterSequence.jl
+++ b/src/ClusterSequence.jl
@@ -148,12 +148,13 @@ ClusterSequence(algorithm::JetAlgorithm.Algorithm, p::Real, R::Float64, strategy
                        Qtot)
 end
 
-ClusterSequence{T}(algorithm::JetAlgorithm.Algorithm, p::Real, R::Float64, strategy::RecoStrategy.Strategy, jets::Vector{T}, history, Qtot) where {T <: FourMomentum} = begin
+function ClusterSequence{T}(algorithm::JetAlgorithm.Algorithm, p::Real, R::Float64,
+                            strategy::RecoStrategy.Strategy, jets::Vector{T}, history,
+                            Qtot) where {T <: FourMomentum}
     ClusterSequence{T}(algorithm, Float64(p), R, strategy, jets, length(jets), history,
                        Qtot)
 end
 
-
 """
     add_step_to_history!(clusterseq::ClusterSequence, parent1, parent2, jetp_index, dij)
 
diff --git a/src/EEAlgorithm.jl b/src/EEAlgorithm.jl
index 589578e..12192ae 100644
--- a/src/EEAlgorithm.jl
+++ b/src/EEAlgorithm.jl
@@ -248,7 +248,7 @@ function ee_genkt_algorithm(particles::AbstractArray{T, 1}; p = 1, R = 4.0,
         for i in eachindex(particles)
             push!(recombination_particles,
                   EEjet{ParticleType}(px(particles[i]), py(particles[i]), pz(particles[i]),
-                        energy(particles[i])))
+                                      energy(particles[i])))
         end
     end
 
@@ -299,8 +299,9 @@ function _ee_genkt_algorithm(; particles::Vector{EEjet{T}}, p = 1, R = 4.0,
     # Setup the initial history and get the total energy
     history, Qtot = initial_history(particles)
 
-    clusterseq = ClusterSequence{EEjet{ParticleType}}(algorithm, p, R, RecoStrategy.N2Plain, particles, history,
-                                 Qtot)
+    clusterseq = ClusterSequence{EEjet{ParticleType}}(algorithm, p, R, RecoStrategy.N2Plain,
+                                                      particles, history,
+                                                      Qtot)
 
     # Run over initial pairs of jets to find nearest neighbours
     get_angular_nearest_neighbours!(eereco, algorithm, dij_factor)
diff --git a/src/EEjet.jl b/src/EEjet.jl
index 9f6ec68..4e85d6d 100644
--- a/src/EEjet.jl
+++ b/src/EEjet.jl
@@ -27,14 +27,12 @@ mutable struct EEjet{T <: Real} <: FourMomentum
     _cluster_hist_index::Int
 end
 
-
 """
     Base.eltype(::Type{EEjet{T}}) where T
 
 Return the element type of the `EEjet` struct.
 """
-Base.eltype(::Type{EEjet{T}}) where T = T
-
+Base.eltype(::Type{EEjet{T}}) where {T} = T
 
 """
     EEjet(px::T, py::T, pz::T, E::T, _cluster_hist_index::Integer) where {T <: Real}
@@ -106,8 +104,8 @@ Constructs an `EEjet` object with conversion of the given momentum components
 An `EEjet` object with the momentum components and energy parametrised to type
 `U`.
 """
-EEjet{U}(px::T, py::T, pz::T, E::T) where {T <: Real, U <: Real} = EEjet(U(px), U(py), U(pz), U(E), 0)
-
+EEjet{U}(px::T, py::T, pz::T, E::T) where {T <: Real, U <: Real} = EEjet(U(px), U(py),
+                                                                         U(pz), U(E), 0)
 
 """
     EEjet(pj::PseudoJet) -> EEjet
@@ -122,7 +120,6 @@ Constructs an `EEjet` object from a given `PseudoJet` object `pj`.
 """
 EEjet(pj::PseudoJet) = EEjet(px(pj), py(pj), pz(pj), energy(pj), cluster_hist_index(pj))
 
-
 p2(eej::EEjet) = eej._p2
 pt2(eej::EEjet) = eej.px^2 + eej.py^2
 const kt2 = pt2
diff --git a/src/PlainAlgo.jl b/src/PlainAlgo.jl
index cbc396d..8b81c49 100644
--- a/src/PlainAlgo.jl
+++ b/src/PlainAlgo.jl
@@ -281,7 +281,7 @@ generalised k_t algorithm.
 """
 function _plain_jet_reconstruct(; particles::Vector{PseudoJet{T}}, p = -1, R = 1.0,
                                 algorithm::JetAlgorithm.Algorithm = JetAlgorithm.AntiKt,
-                                recombine = +) where T <: Real
+                                recombine = +) where {T <: Real}
     # Bounds
     N::Int = length(particles)
     # Parameters
@@ -308,8 +308,10 @@ function _plain_jet_reconstruct(; particles::Vector{PseudoJet{T}}, p = -1, R = 1
     # Current implementation mutates the particles vector, so need to copy it
     # for the cluster sequence (there is too much copying happening, so this
     # needs to be rethought and reoptimised)
-    clusterseq = ClusterSequence{PseudoJet{ParticleType}}(algorithm, p, R, RecoStrategy.N2Plain, particles, history,
-                                 Qtot)
+    clusterseq = ClusterSequence{PseudoJet{ParticleType}}(algorithm, p, R,
+                                                          RecoStrategy.N2Plain, particles,
+                                                          history,
+                                                          Qtot)
 
     # Initialize nearest neighbours
     @simd for i in 1:N
diff --git a/src/Pseudojet.jl b/src/Pseudojet.jl
index b3a5a27..3ce0ee6 100644
--- a/src/Pseudojet.jl
+++ b/src/Pseudojet.jl
@@ -53,7 +53,7 @@ end
 
 Return the element type of the `PseudoJet` struct.
 """
-Base.eltype(::Type{PseudoJet{T}}) where T = T
+Base.eltype(::Type{PseudoJet{T}}) where {T} = T
 
 """
     PseudoJet(px::T, py::T, pz::T, E::T,
@@ -75,10 +75,10 @@ history index.
 A `PseudoJet` object.
 """
 PseudoJet(px::T, py::T, pz::T, E::T,
-    _cluster_hist_index::Int,
-    pt2::T) where {T <: Real} = PseudoJet{T}(px,
-                                            py, pz, E, _cluster_hist_index,
-                                             pt2, 1.0 / pt2, _invalid_rap, _invalid_phi)
+_cluster_hist_index::Int,
+pt2::T) where {T <: Real} = PseudoJet{T}(px,
+                                         py, pz, E, _cluster_hist_index,
+                                         pt2, 1.0 / pt2, _invalid_rap, _invalid_phi)
 
 """
     PseudoJet{T}(px::T, py::T, pz::T, E::T,
@@ -100,10 +100,10 @@ and energy and history index.
 A `PseudoJet` object.
 """
 PseudoJet{T}(px::T, py::T, pz::T, E::T,
-            _cluster_hist_index::Int,
-            pt2::T) where {T <: Real} = PseudoJet{T}(px,
-                                                    py, pz, E, _cluster_hist_index,
-                                                    pt2, 1.0 / pt2, _invalid_rap, _invalid_phi)
+_cluster_hist_index::Int,
+pt2::T) where {T <: Real} = PseudoJet{T}(px,
+                                         py, pz, E, _cluster_hist_index,
+                                         pt2, 1.0 / pt2, _invalid_rap, _invalid_phi)
 
 """
     PseudoJet(px::T, py::T, pz::T, E::T) where T <: Real
@@ -125,8 +125,12 @@ pz::T, E::T) where {T <: Real} = PseudoJet(px, py, pz, E, 0, px^2 + py^2)
 """
     PseudoJet{U}(px::T, py::T, pz::T, E::T) where T <: Real, U <: Real
 """
-PseudoJet{U}(px::T, py::T, pz::T, E::T) where {T <: Real, U <: Real} = PseudoJet{U}(U(px), U(py), U(pz), U(E), 0, U(px^2 + py^2))
-
+PseudoJet{U}(px::T, py::T, pz::T, E::T) where {T <: Real, U <: Real} = PseudoJet{U}(U(px),
+                                                                                    U(py),
+                                                                                    U(pz),
+                                                                                    U(E), 0,
+                                                                                    U(px^2 +
+                                                                                      py^2))
 
 import Base.show
 """
diff --git a/src/TiledAlgoLL.jl b/src/TiledAlgoLL.jl
index 34c53de..dd578c4 100644
--- a/src/TiledAlgoLL.jl
+++ b/src/TiledAlgoLL.jl
@@ -385,8 +385,9 @@ function tiled_jet_reconstruct(particles::AbstractArray{T, 1}; p::Union{Real, No
         sizehint!(recombination_particles, length(particles) * 2)
         for i in eachindex(particles)
             push!(recombination_particles,
-                  PseudoJet{ParticleType}(px(particles[i]), py(particles[i]), pz(particles[i]),
-                            energy(particles[i])))
+                  PseudoJet{ParticleType}(px(particles[i]), py(particles[i]),
+                                          pz(particles[i]),
+                                          energy(particles[i])))
         end
     end
 
@@ -424,7 +425,7 @@ tiled_jet_reconstruct(particles::Vector{PseudoJet}; p = 1, R = 1.0, recombine =
 """
 function _tiled_jet_reconstruct(particles::Vector{PseudoJet{T}}; p::Real = -1, R = 1.0,
                                 algorithm::JetAlgorithm.Algorithm = JetAlgorithm.AntiKt,
-                                recombine = +) where T <: Real
+                                recombine = +) where {T <: Real}
     # Bounds
     N::Int = length(particles)
 
@@ -468,7 +469,9 @@ function _tiled_jet_reconstruct(particles::Vector{PseudoJet{T}}; p::Real = -1, R
     tiling = Tiling(setup_tiling(_eta, R))
 
     # ClusterSequence is the struct that holds the state of the reconstruction
-    clusterseq = ClusterSequence{PseudoJet{ParticleType}}(algorithm, p, R, RecoStrategy.N2Tiled, jets, history, Qtot)
+    clusterseq = ClusterSequence{PseudoJet{ParticleType}}(algorithm, p, R,
+                                                          RecoStrategy.N2Tiled, jets,
+                                                          history, Qtot)
 
     # Tiled jets is a structure that has additional variables for tracking which tile a jet is in
     tiledjets = similar(clusterseq.jets, TiledJet)

From f52aa36456f47b0365e3227e5f20653fca6db241 Mon Sep 17 00:00:00 2001
From: Graeme A Stewart <graeme.andrew.stewart@cern.ch>
Date: Fri, 6 Dec 2024 16:27:20 +0100
Subject: [PATCH 13/16] Remove non-parametrised outer constructor

At the full constructor level the type parameter really should be known,
the use case for the non-explicit parametrised version is weak, so
it is removed. (This was also causing a docstring warning.)
---
 src/Pseudojet.jl | 41 +++++++++++++++--------------------------
 1 file changed, 15 insertions(+), 26 deletions(-)

diff --git a/src/Pseudojet.jl b/src/Pseudojet.jl
index 3ce0ee6..a06524c 100644
--- a/src/Pseudojet.jl
+++ b/src/Pseudojet.jl
@@ -55,31 +55,6 @@ Return the element type of the `PseudoJet` struct.
 """
 Base.eltype(::Type{PseudoJet{T}}) where {T} = T
 
-"""
-    PseudoJet(px::T, py::T, pz::T, E::T,
-        _cluster_hist_index::Int,
-        pt2::T) where {T <: Real}
-
-Constructs a PseudoJet object with the given momentum components and energy and
-history index.
-
-# Arguments
-- `px::T`: The x-component of the momentum.
-- `py::T`: The y-component of the momentum.
-- `pz::T`: The z-component of the momentum.
-- `E::T`: The energy.
-- `_cluster_hist_index::Int`: The cluster history index.
-- `pt2::T`: The transverse momentum squared.
-
-# Returns
-A `PseudoJet` object.
-"""
-PseudoJet(px::T, py::T, pz::T, E::T,
-_cluster_hist_index::Int,
-pt2::T) where {T <: Real} = PseudoJet{T}(px,
-                                         py, pz, E, _cluster_hist_index,
-                                         pt2, 1.0 / pt2, _invalid_rap, _invalid_phi)
-
 """
     PseudoJet{T}(px::T, py::T, pz::T, E::T,
         _cluster_hist_index::Int,
@@ -120,10 +95,24 @@ Constructs a PseudoJet object with the given momentum components and energy.
 A PseudoJet object, with type parameter `T`.
 """
 PseudoJet(px::T, py::T,
-pz::T, E::T) where {T <: Real} = PseudoJet(px, py, pz, E, 0, px^2 + py^2)
+pz::T, E::T) where {T <: Real} = PseudoJet{T}(px, py, pz, E, 0, px^2 + py^2)
 
 """
     PseudoJet{U}(px::T, py::T, pz::T, E::T) where T <: Real, U <: Real
+
+Constructs a PseudoJet object with the given momentum components and energy,
+parameterised by `U`.
+
+# Arguments
+- `px::T`: The x-component of the momentum.
+- `py::T`: The y-component of the momentum.
+- `pz::T`: The z-component of the momentum.
+- `E::T`: The energy.
+
+- U <: Real: The type parameter for the PseudoJet object.
+
+# Returns
+A PseudoJet object, with type parameter `U`.
 """
 PseudoJet{U}(px::T, py::T, pz::T, E::T) where {T <: Real, U <: Real} = PseudoJet{U}(U(px),
                                                                                     U(py),

From 5487316c383d58a6417b62488ac8a818f2363a35 Mon Sep 17 00:00:00 2001
From: Graeme A Stewart <graeme.andrew.stewart@cern.ch>
Date: Fri, 6 Dec 2024 16:38:15 +0100
Subject: [PATCH 14/16] Make use of the known jet type parameter

The particle type is know in the function signature of the underlying reconstruction
---
 src/EEAlgorithm.jl | 4 +++-
 src/PlainAlgo.jl   | 4 +++-
 src/TiledAlgoLL.jl | 4 +++-
 3 files changed, 9 insertions(+), 3 deletions(-)

diff --git a/src/EEAlgorithm.jl b/src/EEAlgorithm.jl
index 12192ae..e7709f6 100644
--- a/src/EEAlgorithm.jl
+++ b/src/EEAlgorithm.jl
@@ -242,6 +242,8 @@ function ee_genkt_algorithm(particles::AbstractArray{T, 1}; p = 1, R = 4.0,
         recombination_particles = copy(particles)
         sizehint!(recombination_particles, length(particles) * 2)
     else
+        # We don't really know what element type we have here, so we need to
+        # drill down to a component to get that underlying type
         ParticleType = typeof(px(particles[1]))
         recombination_particles = EEjet{ParticleType}[]
         sizehint!(recombination_particles, length(particles) * 2)
@@ -275,7 +277,7 @@ function _ee_genkt_algorithm(; particles::Vector{EEjet{T}}, p = 1, R = 4.0,
     R2 = R^2
 
     # Numerical type?
-    ParticleType = eltype(particles[1])
+    ParticleType = T
 
     # Constant factor for the dij metric and the beam distance function
     if algorithm == JetAlgorithm.Durham
diff --git a/src/PlainAlgo.jl b/src/PlainAlgo.jl
index 8b81c49..b648478 100644
--- a/src/PlainAlgo.jl
+++ b/src/PlainAlgo.jl
@@ -236,6 +236,8 @@ function plain_jet_reconstruct(particles::AbstractArray{T, 1}; p::Union{Real, No
         recombination_particles = copy(particles)
         sizehint!(recombination_particles, length(particles) * 2)
     else
+        # We don't really know what element type we have here, so we need to
+        # drill down to a component to get that underlying type
         ParticleType = typeof(px(particles[1]))
         recombination_particles = PseudoJet{ParticleType}[]
         sizehint!(recombination_particles, length(particles) * 2)
@@ -288,7 +290,7 @@ function _plain_jet_reconstruct(; particles::Vector{PseudoJet{T}}, p = -1, R = 1
     R2 = R^2
 
     # Numerical type for this reconstruction
-    ParticleType = eltype(particles[1])
+    ParticleType = T
 
     # Optimised compact arrays for determining the next merge step
     # We make sure these arrays are type stable - have seen issues where, depending on the values
diff --git a/src/TiledAlgoLL.jl b/src/TiledAlgoLL.jl
index dd578c4..efe09f9 100644
--- a/src/TiledAlgoLL.jl
+++ b/src/TiledAlgoLL.jl
@@ -380,6 +380,8 @@ function tiled_jet_reconstruct(particles::AbstractArray{T, 1}; p::Union{Real, No
         recombination_particles = copy(particles)
         sizehint!(recombination_particles, length(particles) * 2)
     else
+        # We don't really know what element type we have here, so we need to
+        # drill down to a component to get that underlying type
         ParticleType = typeof(px(particles[1]))
         recombination_particles = PseudoJet{ParticleType}[]
         sizehint!(recombination_particles, length(particles) * 2)
@@ -442,7 +444,7 @@ function _tiled_jet_reconstruct(particles::Vector{PseudoJet{T}}; p::Real = -1, R
     p = (round(p) == p) ? Int(p) : p # integer p if possible
 
     # Numerical type for this reconstruction
-    ParticleType = eltype(particles[1])
+    ParticleType = T
 
     # This will be used quite deep inside loops, so declare it here so that
     # memory (de)allocation gets done only once

From faf67cf432f9a40c04a61a7e31bf3551e8a5085d Mon Sep 17 00:00:00 2001
From: Graeme A Stewart <graeme.andrew.stewart@cern.ch>
Date: Fri, 6 Dec 2024 16:50:32 +0100
Subject: [PATCH 15/16] Fix for jetreco.jl jet types

We explicitly use Float64 here, because this is deliberately the simple case.
---
 examples/jetreco.jl | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/examples/jetreco.jl b/examples/jetreco.jl
index 82a6a11..b80c062 100644
--- a/examples/jetreco.jl
+++ b/examples/jetreco.jl
@@ -140,17 +140,19 @@ function main()
     logger = ConsoleLogger(stdout, Logging.Info)
     global_logger(logger)
     # Try to read events into the correct type!
-    if JetReconstruction.is_ee(args[:algorithm])
-        jet_type = EEjet
+    # If we don't have an algorithm we default to PseudoJet
+    if !isnothing(args[:algorithm])
+        JetReconstruction.is_ee(args[:algorithm])
+        jet_type = EEjet{Float64}
     else
-        jet_type = PseudoJet
+        jet_type = PseudoJet{Float64}
     end
     events::Vector{Vector{jet_type}} = read_final_state_particles(args[:file],
                                                                   maxevents = args[:maxevents],
                                                                   skipevents = args[:skip],
                                                                   T = jet_type)
     if isnothing(args[:algorithm]) && isnothing(args[:power])
-        @warn "Neither algorithm nor power specified, defaulting to AntiKt"
+        @warn "Neither algorithm nor power specified, defaulting to pp event AntiKt"
         args[:algorithm] = JetAlgorithm.AntiKt
     end
     jet_process(events, distance = args[:distance], algorithm = args[:algorithm],

From 9f96393af1f67da0c3901eaa4867cdc576722966 Mon Sep 17 00:00:00 2001
From: Graeme A Stewart <graeme.andrew.stewart@cern.ch>
Date: Sat, 14 Dec 2024 17:06:45 +0100
Subject: [PATCH 16/16] Adapt tests to allow variable float types

Numerical type can be passed as a parameter to the test struct
This is passed down to read the events in the required precision.
Add some tests for e+e- and pp in Float32
---
 src/Utils.jl                    |  3 ++-
 test/_common.jl                 | 40 +++++++++++++++++++++++++++------
 test/runtests.jl                |  3 +++
 test/test-f32-reconstruction.jl | 20 +++++++++++++++++
 4 files changed, 58 insertions(+), 8 deletions(-)
 create mode 100644 test/test-f32-reconstruction.jl

diff --git a/src/Utils.jl b/src/Utils.jl
index 24448bf..ae82a57 100644
--- a/src/Utils.jl
+++ b/src/Utils.jl
@@ -21,7 +21,8 @@ the particles of a particular event. In particular T can be `PseudoJet` or
 a `LorentzVector` type. Note, if T is not `PseudoJet`, the order of the
 arguments in the constructor must be `(t, x, y, z)`.
 """
-function read_final_state_particles(fname; maxevents = -1, skipevents = 0, T = PseudoJet)
+function read_final_state_particles(fname; maxevents = -1, skipevents = 0,
+                                    T = PseudoJet{Float64})
     f = open(fname)
     if endswith(fname, ".gz")
         @debug "Reading gzipped file $fname"
diff --git a/test/_common.jl b/test/_common.jl
index 7628a5b..04ffb60 100644
--- a/test/_common.jl
+++ b/test/_common.jl
@@ -45,7 +45,7 @@ a vector of `FinalJet` objects, e.g.,
 selector = (cs) -> exclusive_jets(cs; njets = 2)
 ```
 """
-struct ComparisonTest
+struct ComparisonTest{T <: Real}
     events_file::AbstractString
     fastjet_outputs::AbstractString
     algorithm::JetAlgorithm.Algorithm
@@ -54,14 +54,24 @@ struct ComparisonTest
     R::Real
     selector::Any
     selector_name::AbstractString
+    numtype::Type{T}
+end
+
+"Define the element type based off the type parameter"
+Base.eltype(::Type{ComparisonTest{T}}) where {T} = T
+
+"""Constructor where there is no type given"""
+function ComparisonTest(events_file, fastjet_outputs, algorithm, strategy, power, R,
+                        selector, selector_name)
+    ComparisonTest(events_file, fastjet_outputs, algorithm, strategy, power, R, selector,
+                   selector_name, Float64)
 end
 
 """Constructor where there is no selector_name given"""
 function ComparisonTest(events_file, fastjet_outputs, algorithm, strategy, power, R,
                         selector)
-    selector_name = ""
     ComparisonTest(events_file, fastjet_outputs, algorithm, strategy, power, R, selector,
-                   selector_name)
+                   "", Float64)
 end
 
 """Read JSON file with fastjet jets in it"""
@@ -94,7 +104,13 @@ end
 
 function run_reco_test(test::ComparisonTest; testname = nothing)
     # Read the input events
-    events = JetReconstruction.read_final_state_particles(test.events_file)
+    if JetReconstruction.is_pp(test.algorithm)
+        events = JetReconstruction.read_final_state_particles(test.events_file;
+                                                              T = PseudoJet{test.numtype})
+    else
+        events = JetReconstruction.read_final_state_particles(test.events_file;
+                                                              T = EEjet{test.numtype})
+    end
     # Read the fastjet results
     fastjet_jets = read_fastjet_outputs(test.fastjet_outputs)
     sort_jets!(fastjet_jets)
@@ -117,6 +133,16 @@ function run_reco_test(test::ComparisonTest; testname = nothing)
         end
     end
 
+    # For Float64 we use a tolerance of 1e-6
+    # For Float32 we use a tolerance of 1e-2, which is very low
+    #  but the problem is that we are comparing with FastJet in double
+    #  precision and in an iterative algorithm, this can lead to accumulating
+    #  differences that end up quite large
+    tol = 1e-6
+    if test.numtype == Float32
+        tol = 1e-2
+    end
+
     @testset "$testname" begin
         # Test each event in turn...
         for (ievt, event) in enumerate(jet_collection)
@@ -131,9 +157,9 @@ function run_reco_test(test::ComparisonTest; testname = nothing)
                         # the momentum
                         # Sometimes phi could be in the range [-π, π], but FastJet always is [0, 2π]
                         normalised_phi = jet.phi < 0.0 ? jet.phi + 2π : jet.phi
-                        @test jet.rap≈fastjet_jets[ievt]["jets"][ijet]["rap"] atol=1e-7
-                        @test normalised_phi≈fastjet_jets[ievt]["jets"][ijet]["phi"] atol=1e-7
-                        @test jet.pt≈fastjet_jets[ievt]["jets"][ijet]["pt"] rtol=1e-6
+                        @test jet.rap≈fastjet_jets[ievt]["jets"][ijet]["rap"] rtol=tol
+                        @test normalised_phi≈fastjet_jets[ievt]["jets"][ijet]["phi"] rtol=tol
+                        @test jet.pt≈fastjet_jets[ievt]["jets"][ijet]["pt"] rtol=tol
                     end
                 end
             end
diff --git a/test/runtests.jl b/test/runtests.jl
index 3dab8d2..6c290d0 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -41,6 +41,9 @@ function main()
     include("test-pp-reconstruction.jl")
     include("test-ee-reconstruction.jl")
 
+    # Also test reconstruction with Float32
+    include("test-f32-reconstruction.jl")
+
     # Compare inputting data in PseudoJet with using a LorentzVector
     do_test_compare_types(RecoStrategy.N2Plain, algname = pp_algorithms[-1], power = -1)
     do_test_compare_types(RecoStrategy.N2Tiled, algname = pp_algorithms[-1], power = -1)
diff --git a/test/test-f32-reconstruction.jl b/test/test-f32-reconstruction.jl
new file mode 100644
index 0000000..e03ada5
--- /dev/null
+++ b/test/test-f32-reconstruction.jl
@@ -0,0 +1,20 @@
+# Tests of reconstruction algorithms in Float32
+
+include("common.jl")
+
+durham_njets3_f32 = ComparisonTest(events_file_ee,
+                                   joinpath(@__DIR__, "data",
+                                            "jet-collections-fastjet-njets3-Durham-eeH.json.gz"),
+                                   JetAlgorithm.Durham, RecoStrategy.N2Plain, 1, 4.0,
+                                   (cs) -> exclusive_jets(cs; njets = 3),
+                                   "exclusive njets Float32", Float32)
+
+run_reco_test(durham_njets3_f32)
+
+antikt_pcut_f32 = ComparisonTest(events_file_pp,
+                                 joinpath(@__DIR__, "data",
+                                          "jet-collections-fastjet-inclusive-AntiKt.json.gz"),
+                                 JetAlgorithm.AntiKt, RecoStrategy.N2Tiled, -1, 0.4,
+                                 (cs) -> inclusive_jets(cs; ptmin = 5.0),
+                                 "inclusive-float32", Float32)
+run_reco_test(antikt_pcut_f32)