diff --git a/README.md b/README.md index cf6b1249..4ca0996c 100644 --- a/README.md +++ b/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) ``` - `particles` - a one dimensional array (vector) of input particles for the clustering @@ -35,27 +35,17 @@ cs = jet_reconstruct(particles::AbstractVector{T}; algorithm = JetAlgorithm.Anti - 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: diff --git a/docs/src/index.md b/docs/src/index.md index 662a4540..f8eb4854 100644 --- a/docs/src/index.md +++ b/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 diff --git a/examples/constituents/jetreco-constituents-nb.jl b/examples/constituents/jetreco-constituents-nb.jl index e159c3fd..cc5b60a6 100644 --- a/examples/constituents/jetreco-constituents-nb.jl +++ b/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" diff --git a/examples/constituents/jetreco-constituents.jl b/examples/constituents/jetreco-constituents.jl index 50778156..50988d26 100644 --- a/examples/constituents/jetreco-constituents.jl +++ b/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) diff --git a/examples/instrumented-jetreco.jl b/examples/instrumented-jetreco.jl index cb60e7e4..1bc46c34 100644 --- a/examples/instrumented-jetreco.jl +++ b/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 diff --git a/examples/jetreco.jl b/examples/jetreco.jl index 37724a90..d52b231c 100644 --- a/examples/jetreco.jl +++ b/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 diff --git a/examples/substructure/jet-grooming.jl b/examples/substructure/jet-grooming.jl index 9dfb524f..9ff46dea 100644 --- a/examples/substructure/jet-grooming.jl +++ b/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) diff --git a/examples/substructure/jet-tagging.jl b/examples/substructure/jet-tagging.jl index aac9ecaf..23fa0be1 100644 --- a/examples/substructure/jet-tagging.jl +++ b/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) diff --git a/examples/visualisation/animate-reconstruction.jl b/examples/visualisation/animate-reconstruction.jl index 11534435..4bfd486a 100644 --- a/examples/visualisation/animate-reconstruction.jl +++ b/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, diff --git a/examples/visualisation/visualise-jets-nb.jl b/examples/visualisation/visualise-jets-nb.jl index 237c1012..09af6c8f 100644 --- a/examples/visualisation/visualise-jets-nb.jl +++ b/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) diff --git a/examples/visualisation/visualise-jets.jl b/examples/visualisation/visualise-jets.jl index 45c90aa6..fcb6338b 100644 --- a/examples/visualisation/visualise-jets.jl +++ b/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) diff --git a/src/AlgorithmStrategyEnums.jl b/src/AlgorithmStrategyEnums.jl index a0e6cec1..64fe8fa7 100644 --- a/src/AlgorithmStrategyEnums.jl +++ b/src/AlgorithmStrategyEnums.jl @@ -38,50 +38,38 @@ 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 @@ -89,37 +77,10 @@ function get_algorithm_power_consistency(; p::Union{Real, Nothing}, 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 """ diff --git a/src/EEAlgorithm.jl b/src/EEAlgorithm.jl index f9bb96d8..27dac346 100644 --- a/src/EEAlgorithm.jl +++ b/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) diff --git a/src/GenericAlgo.jl b/src/GenericAlgo.jl index b6cd3043..3e5177c8 100644 --- a/src/GenericAlgo.jl +++ b/src/GenericAlgo.jl @@ -1,7 +1,7 @@ """ - jet_reconstruct(particles::AbstractVector; p::Union{Real, Nothing} = nothing, - algorithm::Union{JetAlgorithm.Algorithm, Nothing} = nothing, - R = 1.0, recombine = addjets, preprocess = nothing, + jet_reconstruct(particles::AbstractVector; algorithm::JetAlgorithm.Algorithm, + p::Union{Real, Nothing} = nothing, R = 1.0, + recombine = addjets, preprocess = nothing, strategy::RecoStrategy.Strategy = RecoStrategy.Best) Reconstructs jets from a collection of particles using a specified algorithm and @@ -10,11 +10,10 @@ strategy. # Arguments - `particles::AbstractVector`: A collection of particles used for jet reconstruction. +- `algorithm::JetAlgorithm.Algorithm`: The algorithm to use for jet reconstruction. - `p::Union{Real, Nothing} = nothing`: The power value used for the distance - measure for generalised k_T, which maps to a particular reconstruction - algorithm (-1 = AntiKt, 0 = Cambridge/Aachen, 1 = Kt). -- `algorithm::Union{JetAlgorithm.Algorithm, Nothing} = nothing`: The algorithm - to use for jet reconstruction. + measure for generalised k_T algorithms (GenKt, EEKt). Other algorithms will + ignore this value. - `R = 1.0`: The jet radius parameter. - `recombine = addjets`: The recombination scheme used for combining particles. - `preprocess = nothing`: The function to preprocess the particles before @@ -24,8 +23,9 @@ strategy. strategy to use. `RecoStrategy.Best` makes a dynamic decision based on the number of starting particles. -Note that one of `p` or `algorithm` must be specified, with `algorithm` -preferred. +Note that `p` must be specified for `GenKt` and `EEKt` algorithms, +other algorithms will ignore its value. +When an algorithm has no `R` dependence the `R` parameter is ignored. # Returns A cluster sequence object containing the reconstructed jets and the merging @@ -49,20 +49,9 @@ input particles and which returns a new particle, usually matching a non-standard recombination scheme, e.g., massless particles for ``p_T`` or ``p_T^2`` recombination. `nothing` means no preprocessing is done. -## Consistency of `p`, `algorithm` and `R` arguments -If an algorithm is explicitly specified the `p` value should be consistent with -it or `nothing`. If the algorithm is one where `p` can vary, then it has to be -given, along with the algorithm. - -If the `p` parameter is passed and `algorithm=nothing`, then pp-type -reconstruction is implied (i.e., AntiKt, CA, Kt or GenKt will be used, depending -on the value of `p`). - -When an algorithm has no `R` dependence the `R` parameter is ignored. - # Example ```julia -jet_reconstruct(particles; p = -1, R = 0.4) +jet_reconstruct(particles; algorithm = JetAlgorithm.AntiKt, R = 0.4) jet_reconstruct(particles; algorithm = JetAlgorithm.Kt, R = 1.0) jet_reconstruct(particles; algorithm = JetAlgorithm.Durham) jet_reconstruct(particles; algorithm = JetAlgorithm.GenKt, p = 0.5, R = 1.0) @@ -70,18 +59,11 @@ jet_reconstruct(particles; algorithm = JetAlgorithm.AntiKt, R = 1.0, preprocess recombine = addjets_ptscheme) ``` """ -function jet_reconstruct(particles::AbstractVector; p::Union{Real, Nothing} = nothing, - algorithm::Union{JetAlgorithm.Algorithm, Nothing} = nothing, - R = 1.0, recombine = addjets, preprocess = nothing, +function jet_reconstruct(particles::AbstractVector; algorithm::JetAlgorithm.Algorithm, + p::Union{Real, Nothing} = nothing, R = 1.0, + recombine = addjets, preprocess = nothing, strategy::RecoStrategy.Strategy = RecoStrategy.Best) - - # Either map to the fixed algorithm corresponding to the strategy - # or to an optimal choice based on the density of initial particles - if (algorithm === nothing) && (p === nothing) - throw(ArgumentError("Please specify either an algorithm or a power value (power only for pp-type reconstruction)")) - end - - if (algorithm === nothing) || is_pp(algorithm) + if is_pp(algorithm) # We assume a pp reconstruction if strategy == RecoStrategy.Best # The breakpoint of ~80 is determined empirically on e+e- -> H and 0.5 TeV pp -> 5 GeV jets @@ -95,6 +77,8 @@ function jet_reconstruct(particles::AbstractVector; p::Union{Real, Nothing} = no end elseif is_ee(algorithm) alg = ee_genkt_algorithm + else + throw(ErrorException("Invalid algorithm neither pp nor ee: $(algorithm)")) end # Now call the chosen algorithm, passing through the other parameters diff --git a/src/PlainAlgo.jl b/src/PlainAlgo.jl index 2e2db535..465ff123 100644 --- a/src/PlainAlgo.jl +++ b/src/PlainAlgo.jl @@ -187,9 +187,10 @@ Base.@propagate_inbounds function upd_nn_step!(i, j, k, N, Nn, kt2_array, rapidi end """ - plain_jet_reconstruct(particles::AbstractVector{T}; p::Union{Real, Nothing} = -1, - algorithm::Union{JetAlgorithm.Algorithm, Nothing} = nothing, - R = 1.0, recombine = addjets, preprocess = nothing) where {T} + plain_jet_reconstruct(particles::AbstractVector{T}; + algorithm::JetAlgorithm.Algorithm, + p::Union{Real, Nothing} = nothing, R = 1.0, + recombine = addjets, preprocess = nothing) where {T} Perform pp jet reconstruction using the plain algorithm. @@ -198,10 +199,10 @@ Perform pp jet reconstruction using the plain algorithm. 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. -- `p::Union{Real, Nothing} = -1`: The power value used for jet reconstruction. -- `algorithm::Union{JetAlgorithm, Nothing} = nothing`: The explicit jet - algorithm to use. -- `R::Float64 = 1.0`: The radius parameter used for jet reconstruction. +- `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. +- `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. @@ -212,9 +213,6 @@ JetReconstruction package namespace. This code will use the `k_t` algorithm types, operating in `(rapidity, φ)` space. -It is not necessary to specify both the `algorithm` and the `p` (power) value. -If both are given they must be consistent or an exception is thrown. - # Returns - `Vector{PseudoJet}`: A vector of reconstructed jets. @@ -224,12 +222,13 @@ jets = plain_jet_reconstruct(particles; p = -1, R = 0.4) jets = plain_jet_reconstruct(particles; algorithm = JetAlgorithm.Kt, R = 1.0) ``` """ -function plain_jet_reconstruct(particles::AbstractVector{T}; p::Union{Real, Nothing} = -1, - algorithm::Union{JetAlgorithm.Algorithm, Nothing} = nothing, - R = 1.0, recombine = addjets, preprocess = nothing) where {T} +function plain_jet_reconstruct(particles::AbstractVector{T}; + algorithm::JetAlgorithm.Algorithm, + p::Union{Real, Nothing} = nothing, R = 1.0, + recombine = addjets, preprocess = nothing) where {T} - # Check for consistency between algorithm and power - (p, algorithm) = get_algorithm_power_consistency(p = p, algorithm = algorithm) + # Get consistent algorithm power + p = get_algorithm_power(p = p, algorithm = algorithm) # Integer p if possible p = (round(p) == p) ? Int(p) : p @@ -261,15 +260,14 @@ function plain_jet_reconstruct(particles::AbstractVector{T}; p::Union{Real, Noth end # Now call the actual reconstruction method, tuned for our internal EDM - _plain_jet_reconstruct(particles = recombination_particles, p = p, R = R, - algorithm = algorithm, + _plain_jet_reconstruct(recombination_particles; algorithm = algorithm, p = p, R = R, recombine = recombine) end """ - _plain_jet_reconstruct(; particles::AbstractVector{PseudoJet}, p = -1, - algorithm::JetAlgorithm.Algorithm = JetAlgorithm.AntiKt, - R = 1.0, recombine = addjets) + _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 @@ -285,12 +283,12 @@ generalised k_t algorithm. The algorithm parameter must be consistent with the power parameter. # Arguments -- `particles`: A vector of `PseudoJet` objects representing the input particles. -- `p = -1`: The power to which the transverse momentum (`pt`) of each particle is +- `particles::AbstractVector{PseudoJet}`: A vector of `PseudoJet` objects + representing the input particles. +- `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 = 1.0`: The jet radius parameter. -- `algorithm::JetAlgorithm.Algorithm = JetAlgorithm.AntiKt`: The jet reconstruction - algorithm to use. - `recombine`: The recombination function used to merge two jets. Default is `+` (additive recombination). @@ -298,9 +296,9 @@ power parameter. - `clusterseq`: The resulting `ClusterSequence` object representing the reconstructed jets. """ -function _plain_jet_reconstruct(; particles::AbstractVector{PseudoJet}, p = -1, - algorithm::JetAlgorithm.Algorithm = JetAlgorithm.AntiKt, - 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 diff --git a/src/TiledAlgoLL.jl b/src/TiledAlgoLL.jl index 25a16481..b197cb9b 100644 --- a/src/TiledAlgoLL.jl +++ b/src/TiledAlgoLL.jl @@ -286,15 +286,15 @@ function find_tile_neighbours!(tile_union, jetA, jetB, oldB, tiling) end """ - tiled_jet_reconstruct(particles::AbstractVector{T}; p::Union{Real, Nothing} = -1, - algorithm::Union{JetAlgorithm.Algorithm, Nothing} = nothing, - R = 1.0, recombine = addjets, preprocess = nothing) where {T} + tiled_jet_reconstruct(particles::AbstractVector{T}; + algorithm::JetAlgorithm.Algorithm, + p::Union{Real, Nothing} = nothing, R = 1.0, + recombine = addjets, preprocess = nothing) where {T} Main jet reconstruction algorithm entry point for reconstructing jets using the tiled strategy for generic jet type T. -This code will use the `k_t` algorithm types, operating in `(rapidity, φ)` -space. +This code will use the `k_t` algorithm types, operating in `(rapidity, φ)` space. It is not necessary to specify both the `algorithm` and the `p` (power) value. If both are given they must be consistent or an exception is thrown. @@ -303,12 +303,10 @@ If both are given they must be consistent or an exception is thrown. - `particles::AbstractVector{T}`: A vector of particles used as input for jet reconstruction. T must support methods px, py, pz and energy (defined in the JetReconstruction namespace) -- `p::Union{Real, Nothing} = -1`: The power parameter for the jet reconstruction - algorithm, thus switching between different algorithms. -- `algorithm::Union{JetAlgorithm.Algorithm, Nothing} = nothing`: The explicit - jet algorithm to use. -- `R::Float64 = 1.0`: The jet radius parameter for the jet reconstruction - algorithm. +- `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. +- `R = 1.0`: The jet radius parameter for the jet reconstruction algorithm. - `recombine::Function = addjets`: The recombination function used to combine particles into a new jet. - `preprocess::Function = nothing`: A function to preprocess the input particles. @@ -321,12 +319,13 @@ If both are given they must be consistent or an exception is thrown. tiled_jet_reconstruct(particles::Vector{LorentzVectorHEP}; p = -1, R = 0.4) ``` """ -function tiled_jet_reconstruct(particles::AbstractVector{T}; p::Union{Real, Nothing} = -1, - algorithm::Union{JetAlgorithm.Algorithm, Nothing} = nothing, - R = 1.0, recombine = addjets, preprocess = nothing) where {T} +function tiled_jet_reconstruct(particles::AbstractVector{T}; + algorithm::JetAlgorithm.Algorithm, + p::Union{Real, Nothing} = nothing, R = 1.0, + recombine = addjets, preprocess = nothing) where {T} - # Check for consistency between algorithm and power - (p, algorithm) = get_algorithm_power_consistency(p = p, algorithm = algorithm) + # Get consistent algorithm power + p = get_algorithm_power(p = p, algorithm = algorithm) if isnothing(preprocess) if T == PseudoJet @@ -354,14 +353,14 @@ function tiled_jet_reconstruct(particles::AbstractVector{T}; p::Union{Real, Noth end end - _tiled_jet_reconstruct(recombination_particles; p = p, R = R, algorithm = algorithm, + _tiled_jet_reconstruct(recombination_particles; algorithm = algorithm, p = p, R = R, recombine = recombine) end """ - _tiled_jet_reconstruct(particles::AbstractVector{PseudoJet}; p::Real = -1, - algorithm::JetAlgorithm.Algorithm = JetAlgorithm.AntiKt, - R = 1.0, recombine = addjets) + _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 @@ -370,12 +369,10 @@ power parameter. ## Arguments - `particles::AbstractVector{PseudoJet}`: A vector of `PseudoJet` particles used as input for jet reconstruction. -- `p::Real = -1`: The power parameter for the jet reconstruction algorithm, thus +- `algorithm::JetAlgorithm.Algorithm`: The jet reconstruction algorithm to use. +- `p::Real`: The power parameter for the jet reconstruction algorithm, thus switching between different algorithms. -- `R = 1.0`: The jet radius parameter for the jet reconstruction - algorithm. -- `algorithm::JetAlgorithm.Algorithm = JetAlgorithm.AntiKt`: The jet reconstruction - algorithm to use. +- `R = 1.0`: The jet radius parameter for the jet reconstruction algorithm. - `recombine::Function = addjets`: The recombination function used for combining pseudojets. @@ -387,9 +384,9 @@ power parameter. tiled_jet_reconstruct(particles::Vector{PseudoJet}; p = 1, R = 0.4) ``` """ -function _tiled_jet_reconstruct(particles::AbstractVector{PseudoJet}; p::Real = -1, - algorithm::JetAlgorithm.Algorithm = JetAlgorithm.AntiKt, - 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) diff --git a/test/test-algpower-consistency.jl b/test/test-algpower-consistency.jl index 18f8e6c6..fe48d71d 100644 --- a/test/test-algpower-consistency.jl +++ b/test/test-algpower-consistency.jl @@ -3,28 +3,26 @@ include("common.jl") @testset "Algorithm/power consistency" begin - @test JetReconstruction.check_algorithm_power_consistency(algorithm = JetAlgorithm.AntiKt, - p = -1) - @test JetReconstruction.check_algorithm_power_consistency(algorithm = JetAlgorithm.CA, - p = 0) - @test JetReconstruction.check_algorithm_power_consistency(algorithm = JetAlgorithm.Kt, - p = 1) + @test JetReconstruction.get_algorithm_power(algorithm = JetAlgorithm.AntiKt, + p = nothing) == -1 + @test JetReconstruction.get_algorithm_power(algorithm = JetAlgorithm.CA, + p = nothing) == 0 + @test JetReconstruction.get_algorithm_power(algorithm = JetAlgorithm.Kt, + p = nothing) == 1 + @test JetReconstruction.get_algorithm_power(algorithm = JetAlgorithm.Durham, + p = nothing) == 1 - @test JetReconstruction.check_algorithm_power_consistency(algorithm = JetAlgorithm.AntiKt, - p = nothing) - @test JetReconstruction.check_algorithm_power_consistency(algorithm = nothing, - p = -1) + @test JetReconstruction.get_algorithm_power(algorithm = JetAlgorithm.GenKt, + p = 1.5) == 1.5 + @test JetReconstruction.get_algorithm_power(algorithm = JetAlgorithm.GenKt, + p = -0.5) == -0.5 + @test JetReconstruction.get_algorithm_power(algorithm = JetAlgorithm.EEKt, + p = 1.5) == 1.5 + @test JetReconstruction.get_algorithm_power(algorithm = JetAlgorithm.EEKt, + p = -0.5) == -0.5 - @test_throws ArgumentError JetReconstruction.check_algorithm_power_consistency(algorithm = JetAlgorithm.AntiKt, - p = 0) - @test_throws ArgumentError JetReconstruction.check_algorithm_power_consistency(algorithm = JetAlgorithm.Kt, - p = 1.5) - - @test JetReconstruction.check_algorithm_power_consistency(algorithm = JetAlgorithm.GenKt, - p = 1.5) - @test JetReconstruction.check_algorithm_power_consistency(algorithm = JetAlgorithm.GenKt, - p = -0.5) - - @test_throws ArgumentError JetReconstruction.check_algorithm_power_consistency(algorithm = JetAlgorithm.GenKt, - p = nothing) + @test_throws ArgumentError JetReconstruction.get_algorithm_power(algorithm = JetAlgorithm.GenKt, + p = nothing) + @test_throws ArgumentError JetReconstruction.get_algorithm_power(algorithm = JetAlgorithm.EEKt, + p = nothing) end diff --git a/test/test-compare-types.jl b/test/test-compare-types.jl index 8ae6ebdb..a08e37e2 100644 --- a/test/test-compare-types.jl +++ b/test/test-compare-types.jl @@ -3,10 +3,10 @@ include("common.jl") function do_test_compare_types(strategy::RecoStrategy.Strategy; - algname = "Unknown", + algorithm::JetAlgorithm.Algorithm, ptmin::Float64 = 5.0, distance::Float64 = 0.4, - power::Integer = -1) + power::Union{Real, Nothing} = nothing) # Strategy if (strategy == RecoStrategy.N2Plain) @@ -25,6 +25,7 @@ function do_test_compare_types(strategy::RecoStrategy.Strategy; jet_collection = FinalJets[] for (ievt, event) in enumerate(events) finaljets = final_jets(inclusive_jets(jet_reconstruction(event, R = distance, + algorithm = algorithm, p = power), ptmin = ptmin)) sort_jets!(finaljets) push!(jet_collection, FinalJets(ievt, finaljets)) @@ -36,12 +37,13 @@ function do_test_compare_types(strategy::RecoStrategy.Strategy; jet_collection_lv = FinalJets[] for (ievt, event) in enumerate(events_lv) finaljets = final_jets(inclusive_jets(jet_reconstruction(event, R = distance, + algorithm = algorithm, p = power), ptmin = ptmin)) sort_jets!(finaljets) push!(jet_collection_lv, FinalJets(ievt, finaljets)) end - @testset "Jet Reconstruction Compare PseudoJet and LorentzVector, Strategy $strategy_name, Algorithm $algname" begin + @testset "Jet Reconstruction Compare PseudoJet and LorentzVector, Strategy $strategy_name, Algorithm $algorithm" begin # Here we test that inputting LorentzVector gave the same results as PseudoJets for (ievt, (event, event_lv)) in enumerate(zip(jet_collection, jet_collection_lv)) @testset "Event $(ievt)" begin @@ -57,5 +59,5 @@ function do_test_compare_types(strategy::RecoStrategy.Strategy; end end -do_test_compare_types(RecoStrategy.N2Plain, algname = pp_algorithms[-1], power = -1) -do_test_compare_types(RecoStrategy.N2Tiled, algname = pp_algorithms[-1], power = -1) +do_test_compare_types(RecoStrategy.N2Plain, algorithm = JetAlgorithm.AntiKt) +do_test_compare_types(RecoStrategy.N2Tiled, algorithm = JetAlgorithm.AntiKt) diff --git a/test/test-constituents.jl b/test/test-constituents.jl index 5fe53701..c64ddc5d 100644 --- a/test/test-constituents.jl +++ b/test/test-constituents.jl @@ -19,7 +19,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) diff --git a/test/test-jet_reconstruct-interface.jl b/test/test-jet_reconstruct-interface.jl index ca41b6bc..fe78b4be 100644 --- a/test/test-jet_reconstruct-interface.jl +++ b/test/test-jet_reconstruct-interface.jl @@ -4,6 +4,11 @@ include("common.jl") let inputs = JetReconstruction.read_final_state_particles(events_file_ee) @testset "jet_reconstruct() interface" begin + # No algorithm and no power will throw + @test_throws UndefKeywordError jet_reconstruct(inputs[1]) + # Power but no algorithm will throw + @test_throws UndefKeywordError jet_reconstruct(inputs[1]; p = 0.5) + # EE Algorithms @test typeof(jet_reconstruct(inputs[1]; algorithm = JetAlgorithm.Durham)) == ClusterSequence{EEJet} @@ -27,28 +32,7 @@ let inputs = JetReconstruction.read_final_state_particles(events_file_ee) ClusterSequence{PseudoJet} @test typeof(jet_reconstruct(inputs[1]; algorithm = JetAlgorithm.GenKt, p = 1.0, R = 0.4)) == ClusterSequence{PseudoJet} - - @test_throws ArgumentError jet_reconstruct(inputs[1]; - algorithm = JetAlgorithm.AntiKt, - p = 0) - @test_throws ArgumentError jet_reconstruct(inputs[1]; - algorithm = JetAlgorithm.CA, - p = 1) @test_throws ArgumentError jet_reconstruct(inputs[1]; - algorithm = JetAlgorithm.Kt, - p = -1) - @test_throws ArgumentError jet_reconstruct(inputs[1]; - algorithm = JetAlgorithm.GenKt, - R = 0.4) - - # Supported for now, but will deprecate this calling mode, where only - # the power is given, at the next major release - @test typeof(jet_reconstruct(inputs[1]; p = -1)) == ClusterSequence{PseudoJet} - @test typeof(jet_reconstruct(inputs[1]; p = 0)) == ClusterSequence{PseudoJet} - @test typeof(jet_reconstruct(inputs[1]; p = 1)) == ClusterSequence{PseudoJet} - @test_throws KeyError jet_reconstruct(inputs[1]; p = 0.5) - - # No algorithm or power will throw - @test_throws ArgumentError jet_reconstruct(inputs[1]) + algorithm = JetAlgorithm.GenKt, R = 0.4) end end diff --git a/test/test-substructure.jl b/test/test-substructure.jl index b61687ab..dff67537 100644 --- a/test/test-substructure.jl +++ b/test/test-substructure.jl @@ -35,7 +35,7 @@ for m in keys(methods) @testset "$testname" begin for (ievt, evt) in enumerate(events) - cluster_seq = jet_reconstruct(evt, p = 0, R = 1.0) + cluster_seq = jet_reconstruct(evt; algorithm = JetAlgorithm.CA, R = 1.0) jets = inclusive_jets(cluster_seq; ptmin = 5.0, T = PseudoJet) groomed = sort!([methods[m](jet, cluster_seq; groomer...) for jet in jets], by = JetReconstruction.pt2, rev = true)