Internal function improvements (#195)

* Remove extra particle copy

This was an unnecessary extra copy that was being made (legacy from
when we allowed the internal method to be called directly).

* More consistent internal reconstruction interfaces

Rename all of the internal reconstruction interfaces to ! functions as
they mutate the input particle vector.

Make _ee_genkt_algorithm! take particles as a positional parameter
for consistency reasons with the other functions.

Improve docstrings.

* Format fixes

* Improved documentation

Co-authored-by: Mateusz Jakub Fila <[email protected]>

* Additional docstring fixes

---------

Co-authored-by: Mateusz Jakub Fila <[email protected]>

Author: graeme-a-stewart

src/EEAlgorithm.jl

@@ -201,11 +201,12 @@ Run an e+e- reconstruction algorithm on a set of initial particles.
 - `particles::AbstractVector{T}`: A vector of particles to be clustered.
 - `algorithm::JetAlgorithm.Algorithm`: The jet algorithm to use.
 - `p::Union{Real, Nothing} = nothing`: The power parameter for the algorithm.
-  Must be specified for EEKt algorithm. Other algorithms will ignore this value.
-- `R = 4.0`: The jet radius parameter. Not required / ignored for the Durham
+  This is not required for the `Durham` algorithm, but must be specified for the
+  `EEKt`` algorithm.
+- `R = 4.0`: The jet radius parameter. Not required and ignored for the `Durham`
   algorithm.
-- `recombine`: The recombination scheme to use.
-- `preprocess`: Preprocessing function for input particles.
+- `recombine = addjets`: The recombination scheme to use.
+- `preprocess = nothing`: Preprocessing function for input particles.
 
 # Returns
 - The result of the jet clustering as a `ClusterSequence` object.
@@ -214,10 +215,10 @@ Run an e+e- reconstruction algorithm on a set of initial particles.
 This is the public interface to the e+e- jet clustering algorithm. The function
 will check for consistency between the algorithm and the power parameter as
 needed. It will then prepare the internal EDM particles for the clustering
-itself, and call the actual reconstruction method `_ee_genkt_algorithm`.
+itself, and call the actual reconstruction method `_ee_genkt_algorithm!`.
 
-If the algorithm is Durham, `R` is nominally set to 4.
-If the algorithm is EEkt, power `p` must be specified.
+If the algorithm is Durham, `R` is nominally set to 4. If the algorithm is EEkt,
+power `p` must also be specified.
 """
 function ee_genkt_algorithm(particles::AbstractVector{T}; algorithm::JetAlgorithm.Algorithm,
                             p::Union{Real, Nothing} = nothing, R = 4.0, recombine = addjets,
@@ -261,21 +262,41 @@ function ee_genkt_algorithm(particles::AbstractVector{T}; algorithm::JetAlgorith
     end
 
     # Now call the actual reconstruction method, tuned for our internal EDM
-    _ee_genkt_algorithm(particles = recombination_particles, p = p, R = R,
-                        algorithm = algorithm,
-                        recombine = recombine)
+    _ee_genkt_algorithm!(recombination_particles; p = p, R = R,
+                         algorithm = algorithm,
+                         recombine = recombine)
 end
 
 """
-    _ee_genkt_algorithm(; particles::AbstractVector{EEJet},
+    _ee_genkt_algorithm!(particles::AbstractVector{EEJet};
                         algorithm::JetAlgorithm.Algorithm, p::Real, R = 4.0,
                         recombine = addjets)
 
-This function is the actual implementation of the e+e- jet clustering algorithm.
+This function is the internal implementation of the e+e- jet clustering
+algorithm. It takes a vector of `EEJet` `particles` representing the input
+particles and reconstructs jets based on the specified parameters.
+
+Users of the package should use the `ee_genkt_algorithm` function as their
+entry point to this jet reconstruction.
+
+# Arguments
+- `particles::AbstractVector{EEJet}`: A vector of `EEJet` particles used
+  as input for jet reconstruction. This vector must supply the correct
+  `cluster_hist_index` values and will be *mutated* as part of the returned
+  `ClusterSequence`.
+- `algorithm::JetAlgorithm.Algorithm`: The jet reconstruction algorithm to use.
+- `p::Real`: The power to which the transverse momentum (`pt`) of each particle
+  is raised.
+- `R = 4.0`: The jet radius parameter.
+- `recombine = addjets`: The recombination function used to merge two jets.
+
+# Returns
+- `clusterseq`: The resulting `ClusterSequence` object representing the
+  reconstructed jets.
 """
-function _ee_genkt_algorithm(; particles::AbstractVector{EEJet},
-                             algorithm::JetAlgorithm.Algorithm, p::Real, R = 4.0,
-                             recombine = addjets)
+function _ee_genkt_algorithm!(particles::AbstractVector{EEJet};
+                              algorithm::JetAlgorithm.Algorithm, p::Real, R::Real = 4.0,
+                              recombine = addjets)
     # Bounds
     N::Int = length(particles)
 

src/PlainAlgo.jl

@@ -194,18 +194,24 @@ end
 
 Perform pp jet reconstruction using the plain algorithm.
 
+The power value maps to specific pp jet reconstruction algorithms, but can be
+omitted when the algorithm implies the power value to use. It must be specified
+for the `GenKt` algorithm.
+
 # Arguments
 - `particles::AbstractVector{T}`: A vector of particles used for jet
    reconstruction, any array of particles, which supports suitable 4-vector
    methods, viz. pt2(), phi(), rapidity(), px(), py(), pz(), energy(), can be
    used for each element.
 - `algorithm::JetAlgorithm.Algorithm`: The jet algorithm to use.
-- `p::Union{Real, Nothing} = nothing`: The power value used for jet reconstruction.
-   Must be specified for GenKt algorithm. Other algorithms will ignore this value.
+- `p::Union{Real, Nothing} = nothing`: The power value used for jet
+   reconstruction. Must be specified for GenKt algorithm. Other algorithms will
+   ignore this value.
 - `R = 1.0`: The radius parameter used for jet reconstruction.
 - `recombine::Function = addjets`: The recombination function used to combine
   particles into a new jet.
-- `preprocess::Function = nothing`: A function to preprocess the input particles.
+- `preprocess::Function = nothing`: A function to preprocess the input
+  particles.
 
 **Note** for the `particles` argument, the 4-vector methods need to exist in the
 JetReconstruction package namespace.
@@ -214,12 +220,13 @@ This code will use the `k_t` algorithm types, operating in `(rapidity, φ)`
 space.
 
 # Returns
-- `Vector{PseudoJet}`: A vector of reconstructed jets.
+- `clusterseq`: The resulting `ClusterSequence` object representing the
+  reconstructed jets.
 
 # Example
 ```julia
-jets = plain_jet_reconstruct(particles; algorithm = JetAlgorithm.GenKt, p = -1, R = 0.4)
 jets = plain_jet_reconstruct(particles; algorithm = JetAlgorithm.Kt, R = 1.0)
+jets = plain_jet_reconstruct(particles; algorithm = JetAlgorithm.GenKt, p = -0.5, R = 0.4)
 ```
 """
 function plain_jet_reconstruct(particles::AbstractVector{T};
@@ -260,45 +267,40 @@ function plain_jet_reconstruct(particles::AbstractVector{T};
     end
 
     # Now call the actual reconstruction method, tuned for our internal EDM
-    _plain_jet_reconstruct(recombination_particles; algorithm = algorithm, p = p, R = R,
-                           recombine = recombine)
+    _plain_jet_reconstruct!(recombination_particles; algorithm = algorithm, p = p, R = R,
+                            recombine = recombine)
 end
 
 """
-    _plain_jet_reconstruct(particles::AbstractVector{PseudoJet};
+    _plain_jet_reconstruct!(particles::AbstractVector{PseudoJet};
                            algorithm::JetAlgorithm.Algorithm, p::Real, R = 1.0,
                            recombine = addjets)
 
 This is the internal implementation of jet reconstruction using the plain
-algorithm. It takes a vector of `particles` representing the input particles and
-reconstructs jets based on the specified parameters. Here the particles must be
-of type `PseudoJet`.
+algorithm. It takes a vector of `PseudoJet` `particles` representing the input
+particles and reconstructs jets based on the specified parameters.
 
 Users of the package should use the `plain_jet_reconstruct` function as their
 entry point to this jet reconstruction.
 
-The power value maps to specific pp jet reconstruction algorithms: -1 = AntiKt,
-0 = Cambridge/Aachen, 1 = Inclusive Kt. Floating point values are allowed for
-generalised k_t algorithm. The algorithm parameter must be consistent with the
-power parameter.
-
 # Arguments
-- `particles::AbstractVector{PseudoJet}`: A vector of `PseudoJet` objects
-  representing the input particles.
+- `particles::AbstractVector{PseudoJet}`: A vector of `PseudoJet` particles used
+  as input for jet reconstruction. This vector must supply the correct
+  `cluster_hist_index` values and will be *mutated* as part of the returned
+  `ClusterSequence`.
 - `algorithm::JetAlgorithm.Algorithm`: The jet reconstruction algorithm to use.
-- `p::Real`: The power to which the transverse momentum (`pt`) of each particle is
-  raised.
+- `p::Real`: The power to which the transverse momentum (`pt`) of each particle
+  is raised.
 - `R = 1.0`: The jet radius parameter.
-- `recombine`: The recombination function used to merge two jets. Default is `+`
-  (additive recombination).
+- `recombine = addjets`: The recombination scheme to use.
 
 # Returns
 - `clusterseq`: The resulting `ClusterSequence` object representing the
   reconstructed jets.
 """
-function _plain_jet_reconstruct(particles::AbstractVector{PseudoJet};
-                                algorithm::JetAlgorithm.Algorithm, p::Real, R = 1.0,
-                                recombine = addjets)
+function _plain_jet_reconstruct!(particles::AbstractVector{PseudoJet};
+                                 algorithm::JetAlgorithm.Algorithm, p::Real, R = 1.0,
+                                 recombine = addjets)
     # Bounds
     N::Int = length(particles)
     # Parameters
@@ -319,9 +321,6 @@ function _plain_jet_reconstruct(particles::AbstractVector{PseudoJet};
 
     # Setup the initial history and get the total energy
     history, Qtot = initial_history(particles)
-    # 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)
 

src/TiledAlgoLL.jl

@@ -312,7 +312,8 @@ If both are given they must be consistent or an exception is thrown.
 - `preprocess::Function = nothing`: A function to preprocess the input particles.
 
 ## Returns
-- `Vector{PseudoJet}`: A vector of reconstructed jets.
+- `ClusterSequence`: The resulting `ClusterSequence` object representing the
+  reconstructed jets.
 
 ## Example
 ```julia
@@ -354,22 +355,24 @@ function tiled_jet_reconstruct(particles::AbstractVector{T};
         end
     end
 
-    _tiled_jet_reconstruct(recombination_particles; algorithm = algorithm, p = p, R = R,
-                           recombine = recombine)
+    _tiled_jet_reconstruct!(recombination_particles; algorithm = algorithm, p = p, R = R,
+                            recombine = recombine)
 end
 
 """
-    _tiled_jet_reconstruct(particles::AbstractVector{PseudoJet};
+    _tiled_jet_reconstruct!(particles::AbstractVector{PseudoJet};
                            algorithm::JetAlgorithm.Algorithm,
                            p::Real, R = 1.0, recombine = addjets)
 
-Main jet reconstruction algorithm entry point for reconstructing jets once preprocessing
-of data types are done. The algorithm parameter must be consistent with the
-power parameter.
+Main jet internal reconstruction algorithm entry point for reconstructing jets
+once preprocessing of data types are done. The algorithm parameter must be
+consistent with the power parameter.
 
 ## Arguments
-- `particles::AbstractVector{PseudoJet}`: A vector of `PseudoJet` particles used as input for jet
-  reconstruction.
+- `particles::AbstractVector{PseudoJet}`: A vector of `PseudoJet` particles used
+  as input for jet reconstruction. This vector must supply the correct
+  `cluster_hist_index` values and will be *mutated* as part of the returned
+  `ClusterSequence`.
 - `algorithm::JetAlgorithm.Algorithm`: The jet reconstruction algorithm to use.
 - `p::Real`: The power parameter for the jet reconstruction algorithm, thus
   switching between different algorithms.
@@ -378,16 +381,17 @@ power parameter.
   pseudojets.
 
 ## Returns
-- `Vector{PseudoJet}`: A vector of reconstructed jets.
+- `clusterseq`: The resulting `ClusterSequence` object representing the
+  reconstructed jets.
 
 ## Example
 ```julia
-_tiled_jet_reconstruct(particles::Vector{PseudoJet}; algorithm = JetAlgorithm.Kt, p = 1, R = 0.4)
+_tiled_jet_reconstruct!(particles::Vector{PseudoJet}; algorithm = JetAlgorithm.Kt, p = 1, R = 0.4)
 ```
 """
-function _tiled_jet_reconstruct(particles::AbstractVector{PseudoJet};
-                                algorithm::JetAlgorithm.Algorithm,
-                                p::Real, R = 1.0, recombine = addjets)
+function _tiled_jet_reconstruct!(particles::AbstractVector{PseudoJet};
+                                 algorithm::JetAlgorithm.Algorithm,
+                                 p::Real, R = 1.0, recombine = addjets)
     # Bounds
     N::Int = length(particles)
 
@@ -407,17 +411,8 @@ function _tiled_jet_reconstruct(particles::AbstractVector{PseudoJet};
     # 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[]
-    sizehint!(jets, N * 2)
-    resize!(jets, N)
-
-    # Copy input data into the jets container
-    copyto!(jets, particles)
-
     # Setup the initial history and get the total energy
-    history, Qtot = initial_history(jets)
+    history, Qtot = initial_history(particles)
 
     # Now get the tiling setup
     _eta = Vector{Float64}(undef, length(particles))
@@ -428,7 +423,8 @@ function _tiled_jet_reconstruct(particles::AbstractVector{PseudoJet};
     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(algorithm, p, R, RecoStrategy.N2Tiled, particles, history,
+                                 Qtot)
 
     # Tiled jets is a structure that has additional variables for tracking which tile a jet is in
     tiledjets = similar(clusterseq.jets, TiledJet)