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/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
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/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],
diff --git a/src/ClusterSequence.jl b/src/ClusterSequence.jl
index 9bdf653..a6d0515 100644
--- a/src/ClusterSequence.jl
+++ b/src/ClusterSequence.jl
@@ -148,6 +148,13 @@ ClusterSequence(algorithm::JetAlgorithm.Algorithm, p::Real, R::Float64, strategy
                        Qtot)
 end
 
+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 d9c502b..e7709f6 100644
--- a/src/EEAlgorithm.jl
+++ b/src/EEAlgorithm.jl
@@ -242,12 +242,15 @@ 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[]
+        # 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)
         for i in eachindex(particles)
             push!(recombination_particles,
-                  EEjet(px(particles[i]), py(particles[i]), pz(particles[i]),
-                        energy(particles[i])))
+                  EEjet{ParticleType}(px(particles[i]), py(particles[i]), pz(particles[i]),
+                                      energy(particles[i])))
         end
     end
 
@@ -264,15 +267,18 @@ 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)
 
     # R squared
     R2 = R^2
 
+    # Numerical type?
+    ParticleType = T
+
     # Constant factor for the dij metric and the beam distance function
     if algorithm == JetAlgorithm.Durham
         dij_factor = 2.0
@@ -289,14 +295,15 @@ 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
     history, Qtot = initial_history(particles)
 
-    clusterseq = ClusterSequence(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 c61f44f..4e85d6d 100644
--- a/src/EEjet.jl
+++ b/src/EEjet.jl
@@ -1,35 +1,123 @@
 """
-    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 reconstruction 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)
+"""
+    Base.eltype(::Type{EEjet{T}}) where T
+
+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(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}
+
+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)
 
+"""
+    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
@@ -87,15 +175,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.
+
+# 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 raised to the power of 2p for this jet.
 
-mutable struct EERecoJet
+# 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
diff --git a/src/PlainAlgo.jl b/src/PlainAlgo.jl
index fabcb91..b648478 100644
--- a/src/PlainAlgo.jl
+++ b/src/PlainAlgo.jl
@@ -230,13 +230,16 @@ 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[]
+        # 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)
         for i in eachindex(particles)
             push!(recombination_particles,
@@ -278,23 +281,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 for this reconstruction
+    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
     # 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,8 +310,10 @@ 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,
-                                 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 3ecfd43..a06524c 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,100 @@ 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.
-
-"""
-mutable struct PseudoJet <: FourMomentum
-    px::Float64
-    py::Float64
-    pz::Float64
-    E::Float64
+- `_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{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,
+    Base.eltype(::Type{PseudoJet{T}}) where T
+
+Return the element type of the `PseudoJet` struct.
+"""
+Base.eltype(::Type{PseudoJet{T}}) where {T} = T
+
+"""
+    PseudoJet{T}(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.
+Constructs a parametrised 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{T}(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.
+A PseudoJet object, with type parameter `T`.
+"""
+PseudoJet(px::T, py::T,
+pz::T, E::T) where {T <: Real} = PseudoJet{T}(px, py, pz, E, 0, px^2 + py^2)
+
 """
-PseudoJet(px::Real, py::Real,
-pz::Real, E::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
+
+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),
+                                                                                    U(pz),
+                                                                                    U(E), 0,
+                                                                                    U(px^2 +
+                                                                                      py^2))
 
 import Base.show
 """
diff --git a/src/TiledAlgoLL.jl b/src/TiledAlgoLL.jl
index ff3b470..efe09f9 100644
--- a/src/TiledAlgoLL.jl
+++ b/src/TiledAlgoLL.jl
@@ -374,18 +374,22 @@ 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[]
+        # 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)
         for i in eachindex(particles)
             push!(recombination_particles,
-                  PseudoJet(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
 
@@ -421,9 +425,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 +443,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 for this reconstruction
+    ParticleType = T
+
     # 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 +463,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 +471,9 @@ 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)
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)