Make algorithm argument mandatory in jet_reconstruct (#161)

* make `algorithm` argument mandatory in `jet_reconstruct`

* Small improvements to README

---------

Co-authored-by: Graeme A Stewart <[email protected]>
Co-authored-by: Jerry Ling <[email protected]>

Author: 3 people

README.md

@@ -27,35 +27,25 @@ 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)
 ```
 
 - `particles` - a one dimensional array (vector) of input particles for the clustering
   - Any type that supplies the methods `pt2()`, `phi()`, `rapidity()`, `px()`, `py()`, `pz()`, `energy()` can be used
   - These methods have to be defined in the namespace of this package, i.e., `JetReconstruction.pt2(::T)`
   - The `PseudoJet` or `EEJet` types from this package, a 4-vector from `LorentzVectorHEP`, or a `ReconstructedParticle` from [EDM4hep](https://github.com/peremato/EDM4hep.jl) are suitable (and have the appropriate definitions)
 - `algorithm` is the name of the jet algorithm to be used (from the `JetAlgorithm` enum)
-  - `JetAlgorithm.AntiKt` anti-$`{k}_\text{T}`$ clustering (default)
+  - `JetAlgorithm.AntiKt` anti-$`{k}_\text{T}`$ clustering
   - `JetAlgorithm.CA` Cambridge/Aachen clustering
   - `JetAlgorithm.Kt` inclusive $k_\text{T}$
   - `JetAlgorithm.GenKt` generalised $k_\text{T}$ (which also requires specification of `p`)
   - `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
+  - `JetAlgorithm.EEKt` the $e^+e-$ generalised $k_\text{T}$ algorithm (which also requires specification of `p`)
 - `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)
-- `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.
+For a more complete description of all possible parameters please [refer to the documentation](https://juliahep.github.io/JetReconstruction.jl/stable/#Reconstruction-Interface).
 
 To obtain the final inclusive jets, use the `inclusive_jets` method:
 

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)