Make algorithm argument mandatory in jet_reconstruct

README.md

@@ -27,7 +27,7 @@ algorithm and generalised $`k_\text{T}`$ for $`e^+e^-`$.
 The simplest interface is to call:
 
 ```julia
-cs = jet_reconstruct(particles::AbstractVector{T}; algorithm = JetAlgorithm.AntiKt, R = 1.0, [p = -1,] [strategy = RecoStrategy.Best])
+cs = jet_reconstruct(particles::AbstractVector{T}; algorithm = JetAlgorithm.AntiKt, R = 1.0, [p = nothing,] [strategy = RecoStrategy.Best])
 ```
 
 - `particles` - a one dimensional array (vector) of input particles for the clustering
@@ -42,21 +42,11 @@ cs = jet_reconstruct(particles::AbstractVector{T}; algorithm = JetAlgorithm.Anti
   - `JetAlgorithm.Durham` the $e^+e-$ $k_\text{T}$ algorithm, also known as the Durham algorithm
   - `JetAlgorithm.EEKt` the $e^+e-$ generalised $k_\text{T}$ algorithm
 - `R` - the cone size parameter; no particles more geometrically distance than `R` will be merged (default 1.0; note this parameter is ignored for the Durham algorithm)
+- `p` - algorithm power which must be provided for generalised algorithms, for other algorithms this parameter is ignored
 - `strategy` - the algorithm strategy to adopt, as described below (default `RecoStrategy.Best`)
 
 The object returned is a `ClusterSequence`, which internally tracks all merge steps.
 
-Alternatively, *for pp reconstruction*, one can swap the `algorithm=...`
-parameter for the value of `p`, the transverse momentum power used in the
-$d_{ij}$ metric for deciding on closest jets, as $k^{2p}_\text{T}$. Different
-values of $p$ then correspond to different reconstruction algorithms:
-
-- `-1` gives anti-$`{k}_\text{T}`$ clustering (default)
-- `0` gives Cambridge/Aachen
-- `1` gives inclusive $k_\text{T}$
-
-Note, for the `GenKt` and `EEKt` algorithms the `p` value *must* also be given to specify the algorithm fully.
-
 To obtain the final inclusive jets, use the `inclusive_jets` method:
 
 ```julia

docs/src/index.md

@@ -31,8 +31,7 @@ jet_reconstruct(particles; algorithm = JetAlgorithm.GenKt, R = 0.4,
 ```
 
 Where `particles` is a collection of 4-vector objects (see [Input Particle
-Types](@ref)) to reconstruct and the algorithm is either given explicitly or
-implied by the power value.
+Types](@ref)) to reconstruct and the algorithm is given explicitly.
 
 For the case of generalised ``k_T`` (for ``pp`` and ``e^+e^-``) both the
 algorithm (`GenKt`, `EEKt`) and `p` are needed.
@@ -58,12 +57,6 @@ Each known algorithm is referenced using a `JetAlgorithm` scoped enum value.
 | ``e^+e-`` ``k_\text{T}`` / Durham | `JetAlgorithm.Durham` | `R` value ignored and can be omitted |
 | generalised ``e^+e-`` ``k_\text{T}`` | `JetAlgorithm.EEKt` | For ``e^+e^-``, value of `p` must also be specified |
 
-#### ``pp`` Algorithms
-
-For the three ``pp`` algorithms with fixed `p` values, the `p` value can be
-given instead of the algorithm name. However, this should be considered
-*deprecated* and will be removed in a future release.
-
 ### Strategy
 
 Generally one does not need to manually specify a strategy, but [Algorithm

examples/constituents/jetreco-constituents-nb.jl

@@ -26,7 +26,7 @@ md"Event to pick"
 event_no = 1
 
 # ╔═╡ 2a899d67-71f3-4fe0-8104-7633a44a06a8
-cluster_seq = jet_reconstruct(events[event_no], p = 1, R = 1.0)
+cluster_seq = jet_reconstruct(events[event_no]; algorithm = JetAlgorithm.Kt, R = 1.0)
 
 # ╔═╡ 043cc484-e537-409c-8aa2-e4904b5dc283
 md"Retrieve the exclusive pj_jets, but as `PseudoJet` types"

examples/constituents/jetreco-constituents.jl

@@ -15,7 +15,7 @@ events = read_final_state_particles(input_file)
 # Event to pick
 event_no = 1
 
-cluster_seq = jet_reconstruct(events[event_no], p = 1, R = 1.0)
+cluster_seq = jet_reconstruct(events[event_no]; algorithm = JetAlgorithm.Kt, R = 1.0)
 
 # Retrieve the exclusive pj_jets, but as `PseudoJet` types
 pj_jets = inclusive_jets(cluster_seq; ptmin = 5.0, T = PseudoJet)

examples/instrumented-jetreco.jl

@@ -76,18 +76,16 @@ function profile_code(events::Vector{Vector{T}}, profile, nsamples; R = 0.4, p =
 end
 
 """
-    allocation_stats(events::Vector{Vector{T}}; distance::Real = 0.4,
-                          p::Union{Real, Nothing} = nothing,
-                          algorithm::Union{JetAlgorithm.Algorithm, Nothing} = nothing,
+    allocation_stats(events::Vector{Vector{T}}; algorithm::JetAlgorithm.Algorithm,
+                          distance::Real = 0.4, p::Union{Real, Nothing} = nothing,
                           strategy::RecoStrategy.Strategy,
+                          recombine = RecombinationMethods[RecombinationScheme.EScheme],
                           ptmin::Real = 5.0) where {T <: JetReconstruction.FourMomentum}
-
 Take a memory allocation profile of the jet reconstruction code, printing the
 output.
 """
-function allocation_stats(events::Vector{Vector{T}}; distance::Real = 0.4,
-                          p::Union{Real, Nothing} = nothing,
-                          algorithm::Union{JetAlgorithm.Algorithm, Nothing} = nothing,
+function allocation_stats(events::Vector{Vector{T}}; algorithm::JetAlgorithm.Algorithm,
+                          distance::Real = 0.4, p::Union{Real, Nothing} = nothing,
                           strategy::RecoStrategy.Strategy,
                           recombine = RecombinationMethods[RecombinationScheme.EScheme],
                           ptmin::Real = 5.0) where {T <: JetReconstruction.FourMomentum}
@@ -103,8 +101,8 @@ end
 
 """
     benchmark_jet_reco(events::Vector{Vector{T}};
+                            algorithm::JetAlgorithm.Algorithm,
                             distance::Real = 0.4,
-                            algorithm::Union{JetAlgorithm.Algorithm, Nothing} = nothing,
                             p::Union{Real, Nothing} = nothing,
                             ptmin::Real = 5.0,
                             dcut = nothing,
@@ -125,8 +123,8 @@ print summary statistics on the runtime.
 - `dump_cs`: If `true`, dump the cluster sequence for each event.
 """
 function benchmark_jet_reco(events::Vector{Vector{T}};
+                            algorithm::JetAlgorithm.Algorithm,
                             distance::Real = 0.4,
-                            algorithm::Union{JetAlgorithm.Algorithm, Nothing} = nothing,
                             p::Union{Real, Nothing} = nothing,
                             strategy::RecoStrategy.Strategy,
                             recombine = RecombinationMethods[RecombinationScheme.EScheme],
@@ -143,9 +141,8 @@ function benchmark_jet_reco(events::Vector{Vector{T}};
         jet_collection = FinalJets[]
     end
 
-    # Set consistent algorithm and power
-    (p, algorithm) = JetReconstruction.get_algorithm_power_consistency(p = p,
-                                                                       algorithm = algorithm)
+    # Set consistent algorithm power
+    p = JetReconstruction.get_algorithm_power(p = p, algorithm = algorithm)
     @info "Jet reconstruction will use $(algorithm) with power $(p)"
 
     # Now setup timers and run the loop

examples/jetreco.jl

@@ -24,8 +24,8 @@ happens inside the JetReconstruction package itself.
 Final jets can be serialised if the "dump" option is given
 """
 function jet_process(events::Vector{Vector{T}};
+                     algorithm::JetAlgorithm.Algorithm,
                      distance::Real = 0.4,
-                     algorithm::Union{JetAlgorithm.Algorithm, Nothing} = nothing,
                      p::Union{Real, Nothing} = nothing,
                      ptmin::Real = 5.0,
                      dcut = nothing,
@@ -34,9 +34,8 @@ function jet_process(events::Vector{Vector{T}};
                      dump::Union{String, Nothing} = nothing) where {T <:
                                                                     JetReconstruction.FourMomentum}
 
-    # Set consistent algorithm and power
-    (p, algorithm) = JetReconstruction.get_algorithm_power_consistency(p = p,
-                                                                       algorithm = algorithm)
+    # Set consistent algorithm power
+    p = JetReconstruction.get_algorithm_power(p = p, algorithm = algorithm)
     @info "Jet reconstruction will use $(algorithm) with power $(p)"
 
     # A friendly label for the algorithm and final jet selection

examples/substructure/jet-grooming.jl

@@ -8,7 +8,7 @@ events = read_final_state_particles(input_file)
 # Event to pick
 event_no = 1
 
-cluster_seq = jet_reconstruct(events[event_no], p = 0, R = 1.0)
+cluster_seq = jet_reconstruct(events[event_no]; algorithm = JetAlgorithm.CA, R = 1.0)
 jets = inclusive_jets(cluster_seq; ptmin = 5.0, T = PseudoJet)
 
 filter = (radius = 0.3, hardest_jets = 3)

examples/substructure/jet-tagging.jl

@@ -8,7 +8,7 @@ events = read_final_state_particles(input_file)
 # Event to pick
 event_no = 1
 
-cluster_seq = jet_reconstruct(events[event_no], p = 0, R = 1.0)
+cluster_seq = jet_reconstruct(events[event_no]; algorithm = JetAlgorithm.CA, R = 1.0)
 jets = inclusive_jets(cluster_seq; ptmin = 5.0, T = PseudoJet)
 
 MDtagger = (mu = 0.67, y = 0.09)

examples/visualisation/animate-reconstruction.jl

@@ -56,10 +56,8 @@ function main()
         @warn "Neither algorithm nor power specified, defaulting to AntiKt"
         args[:algorithm] = JetAlgorithm.AntiKt
     end
-    # Set consistent algorithm and power
-    (p, algorithm) = JetReconstruction.get_algorithm_power_consistency(p = args[:power],
-                                                                       algorithm = args[:algorithm])
-    cs = jet_reconstruct(events[1], R = args[:distance], p = p, algorithm = algorithm,
+    cs = jet_reconstruct(events[1], R = args[:distance], p = args[:power],
+                         algorithm = args[:algorithm],
                          strategy = args[:strategy])
 
     animatereco(cs, args[:output]; azimuth = (1.8, 3.0), elevation = 0.5,

examples/visualisation/visualise-jets-nb.jl

@@ -30,7 +30,7 @@ md"""# Jet Reconstruction Visualisation
 
 This Pluto script visualises the result of a jet reconstruction process. 
 
-Use the sliders below to change the reconstructed event, the algorithm and the jet radius parameter."""
+Use the sliders below to change the reconstructed event, the GenKt algorithm power and the jet radius parameter."""
 
 # ╔═╡ 4e569f74-570b-4b30-9ea7-9cbc420f50f8
 md"Event number:"
@@ -60,7 +60,7 @@ events::Vector{Vector{PseudoJet}} = read_final_state_particles(input_file,
                                                                skipevents = event_no);
 
 # ╔═╡ 2a899d67-71f3-4fe0-8104-7633a44a06a8
-cs = jet_reconstruct(events[1], p = power, R = radius)
+cs = jet_reconstruct(events[1]; algorithm = JetAlgorithm.GenKt, p = power, R = radius)
 
 # ╔═╡ b5fd4e96-d073-4e5f-8de1-41addaa0dc3d
 jetreco_vis = jetsplot(events[1], cs; Module = GLMakie)

examples/visualisation/visualise-jets.jl

@@ -65,11 +65,8 @@ function main()
     events::Vector{Vector{PseudoJet}} = read_final_state_particles(args[:file],
                                                                    maxevents = args[:event],
                                                                    skipevents = args[:event])
-
-    (p,
-    algorithm) = JetReconstruction.get_algorithm_power_consistency(p = args[:power],
-                                                                   algorithm = args[:algorithm])
-    cs = jet_reconstruct(events[1], R = args[:distance], p = p, algorithm = algorithm,
+    cs = jet_reconstruct(events[1], R = args[:distance], p = args[:power],
+                         algorithm = args[:algorithm],
                          strategy = args[:strategy])
 
     plt = jetsplot(events[1], cs; Module = CairoMakie)

src/AlgorithmStrategyEnums.jl

@@ -38,88 +38,49 @@ A constant array that contains the jet algorithms for which power is variable.
 """
 const varpower_algorithms = [JetAlgorithm.GenKt, JetAlgorithm.EEKt]
 
-"""
-    power2algorithm
-
-A dictionary that maps power values to corresponding jet algorithm used for pp
-jet reconstruction, where these are fixed.
-"""
-const power2algorithm = Dict(-1 => JetAlgorithm.AntiKt,
-                             0 => JetAlgorithm.CA,
-                             1 => JetAlgorithm.Kt)
-
 """
     algorithm2power
 
-A dictionary that maps pp algorithm names to their corresponding power values.
-
-The dictionary is created by iterating over the `power2algorithm` dictionary and swapping the keys and values.
+A dictionary that maps algorithm names to their corresponding power values.
 """
-const algorithm2power = Dict((i.second, i.first) for i in power2algorithm)
+const algorithm2power = Dict(JetAlgorithm.AntiKt => -1,
+                             JetAlgorithm.CA => 0,
+                             JetAlgorithm.Kt => 1,
+                             JetAlgorithm.Durham => 1)
 
 """
-    get_algorithm_power_consistency(; p::Union{Real, Nothing}, algorithm::Union{JetAlgorithm.Algorithm, Nothing})
+    get_algorithm_power(; algorithm::JetAlgorithm.Algorithm, p::Union{Real, Nothing}) -> Real
 
-Get the algorithm and power consistency correct
+Get the algorithm power
 
-This function checks the consistency between the algorithm and power parameters.
-If the algorithm is specified, it checks if the power parameter is consistent
-with the algorithm's known power. If the power parameter is not specified, it
-sets the power parameter based on the algorithm. If neither the algorithm nor
-the power parameter is specified, it throws an `ArgumentError`.
+This function returns appropriate power value for the specified jet algorithm if the algorithm
+is a fixed power algorithm (like `AntiKt`, `CA`, `Kt`, or `Durham`).
+If the algorithm is generalized (like `GenKt` or `EEKt`), it requires a power value
+to be specified and the function returns the same value.
 
 # Arguments
 - `p::Union{Real, Nothing}`: The power value.
-- `algorithm::Union{JetAlgorithm.Algorithm, Nothing}`: The algorithm.
+- `algorithm::JetAlgorithm.Algorithm`: The algorithm.
 
 # Returns
-A named tuple of the consistent power and algorithm values.
+Resolved algorithm power value.
 
 # Throws
-- `ArgumentError`: If the algorithm and power are inconsistent or if neither the
-  algorithm nor the power is specified.
+- `ArgumentError`: If no power is specified for a generalized algorithm
 
 """
-function get_algorithm_power_consistency(; p::Union{Real, Nothing},
-                                         algorithm::Union{JetAlgorithm.Algorithm, Nothing})
+function get_algorithm_power(; algorithm::JetAlgorithm.Algorithm, p::Union{Real, Nothing})
     # For the case where an algorithm has a variable power value
     # we need to check that the power value and algorithm are both specified
     if algorithm in varpower_algorithms
         if isnothing(p)
             throw(ArgumentError("Power must be specified for algorithm $algorithm"))
         end
         # And then we just return what these are...
-        return (p = p, algorithm = algorithm)
-    end
-
-    # Some algorithms have fixed power values
-    if algorithm == JetAlgorithm.Durham
-        return (p = 1, algorithm = algorithm)
+        return p
     end
-
-    # Otherwise we check the consistency between the algorithm and power
-    if !isnothing(algorithm)
-        power_from_alg = algorithm2power[algorithm]
-        if !isnothing(p) && p != power_from_alg
-            throw(ArgumentError("Algorithm and power are inconsistent"))
-        end
-        return (p = power_from_alg, algorithm = algorithm)
-    else
-        if isnothing(p)
-            throw(ArgumentError("Either algorithm or power must be specified"))
-        end
-        # Set the algorithm from the power
-        algorithm_from_power = power2algorithm[p]
-        return (p = p, algorithm = algorithm_from_power)
-    end
-end
-
-"""Allow a check for algorithm and power consistency"""
-function check_algorithm_power_consistency(; p::Union{Real, Nothing},
-                                           algorithm::Union{JetAlgorithm.Algorithm,
-                                                            Nothing})
-    get_algorithm_power_consistency(p = p, algorithm = algorithm)
-    return true
+    # Otherwise the algorithm has a fixed power value
+    return algorithm2power[algorithm]
 end
 
 """

src/EEAlgorithm.jl

@@ -191,20 +191,19 @@ Copy the contents of slot `i` in the `eereco` array to slot `j`.
 end
 
 """
-    ee_genkt_algorithm(particles::AbstractVector{T}; p = -1, R = 4.0,
-                       algorithm::JetAlgorithm.Algorithm = JetAlgorithm.Durham,
-                       recombine = addjets, preprocess = nothing) where {T}
+    ee_genkt_algorithm(particles::AbstractVector{T}; algorithm::JetAlgorithm.Algorithm,
+                       p::Union{Real, Nothing} = nothing, R = 4.0, recombine = addjets,
+                       preprocess = nothing) where {T}
 
 Run an e+e- reconstruction algorithm on a set of initial particles.
 
 # Arguments
 - `particles::AbstractVector{T}`: A vector of particles to be clustered.
-- `p = 1`: The power parameter for the algorithm. Not required / ignored for
-  the Durham algorithm when it is set to 1.
+- `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
   algorithm.
-- `algorithm::JetAlgorithm.Algorithm = JetAlgorithm.Durham`: The specific jet
-  algorithm to use.
 - `recombine`: The recombination scheme to use.
 - `preprocess`: Preprocessing function for input particles.
 
@@ -217,17 +216,15 @@ 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`.
 
-If the algorithm is Durham, `p` is set to 1 and `R` is nominally set to 4.
-
-Note that unlike `pp` reconstruction the algorithm has to be specified
-explicitly.
+If the algorithm is Durham, `R` is nominally set to 4.
+If the algorithm is EEkt, power `p` must be specified.
 """
-function ee_genkt_algorithm(particles::AbstractVector{T}; p = 1,
-                            algorithm::JetAlgorithm.Algorithm = JetAlgorithm.Durham,
-                            R = 4.0, recombine = addjets, preprocess = nothing) where {T}
+function ee_genkt_algorithm(particles::AbstractVector{T}; algorithm::JetAlgorithm.Algorithm,
+                            p::Union{Real, Nothing} = nothing, R = 4.0, recombine = addjets,
+                            preprocess = nothing) where {T}
 
-    # Check for consistency between algorithm and power
-    (p, algorithm) = get_algorithm_power_consistency(p = p, algorithm = algorithm)
+    # Check for consistency algorithm power
+    p = get_algorithm_power(p = p, algorithm = algorithm)
 
     # Integer p if possible
     p = (round(p) == p) ? Int(p) : p
@@ -270,14 +267,14 @@ function ee_genkt_algorithm(particles::AbstractVector{T}; p = 1,
 end
 
 """
-    _ee_genkt_algorithm(; particles::AbstractVector{EEJet}, p = 1, R = 4.0,
-                       algorithm::JetAlgorithm.Algorithm = JetAlgorithm.Durham,
-                       recombine = addjets)
+    _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.
 """
-function _ee_genkt_algorithm(; particles::AbstractVector{EEJet}, p = 1, R = 4.0,
-                             algorithm::JetAlgorithm.Algorithm = JetAlgorithm.Durham,
+function _ee_genkt_algorithm(; particles::AbstractVector{EEJet},
+                             algorithm::JetAlgorithm.Algorithm, p::Real, R = 4.0,
                              recombine = addjets)
     # Bounds
     N::Int = length(particles)