From d310cbde3dcf83d27638f4f6f2095e2dbd43ffbe Mon Sep 17 00:00:00 2001 From: Graeme A Stewart Date: Fri, 19 Apr 2024 16:11:45 +0200 Subject: [PATCH] Change generic method name to `jet_reconstruct` There is no need to have the prefix `generic` Update the README to show how to ruin reconstruction using the "Best" strategy (should be the default way) Add also a note on sorting (which is super trivial in Julia, so no need to support specific methods for `Vector{PseudoJet}` as fastjet does) --- README.md | 60 +++++++++++++++++++++++++++++----------- examples/jetreco.jl | 10 +++---- src/ClusterSequence.jl | 40 ++++++++++++++------------- src/GenericAlgo.jl | 2 +- src/JetReconstruction.jl | 2 +- test/runtests.jl | 3 ++ 6 files changed, 75 insertions(+), 42 deletions(-) diff --git a/README.md b/README.md index e96d51a..c27d1b5 100644 --- a/README.md +++ b/README.md @@ -10,21 +10,14 @@ Algorithms used are based on the C++ FastJet package (, [hep-ph/0512210](https://arxiv.org/abs/hep-ph/0512210), [arXiv:1111.6097](https://arxiv.org/abs/1111.6097)), reimplemented natively in Julia. -One algorithm is the plain N2 approach `N2Plain`, the other uses the FastJet tiling -approach, `N2Tiled`. +The algorithms include anti-${k}_\text{T}$, Cambridge/Achen and inclusive $k_\text{T}$. ### Interface -To use the package one should call the appropriate API interface of the desired algorithm +The simplest interface is to call: ```julia -# N2Plain -plain_jet_reconstruct(particles::Vector{T}; p = -1, R = 1.0, recombine = +, ptmin = 0.0) -``` - -```julia -# N2Tiled -tiled_jet_reconstruct(particles::Vector{T}; p = -1, R = 1.0, recombine = +, ptmin = 0.0) +cs = jet_reconstruct(particles::Vector{T}; p = -1, R = 1.0, recombine = +, strategy = JetRecoStrategy.Best) ``` - `particles` - a vector of input particles for the clustering @@ -32,16 +25,51 @@ tiled_jet_reconstruct(particles::Vector{T}; p = -1, R = 1.0, recombine = +, ptmi - These methods have to be defined in the namespace of this package, i.e., `JetReconstruction.pt2(::T)` - The `PseudoJet` type from this package, or a 4-vector from `LorentzVectorHEP` are suitable (and have the appropriate definitions) - `p` - the transverse momentum power used in the $d_{ij}$ metric for deciding on closest jets, as $k^{2p}_\text{T}$ - - `-1` gives anti-${k}_\text{T}$ clustering + - `-1` gives anti-${k}_\text{T}$ clustering (default) - `0` gives Cambridge/Achen - `1` gives inclusive $k_\text{T}$ -- `R` - the cone size parameter (no particles more geometrically distance than `R` will be merged) -- `recombine` - the function used to merge two pseudojets (by default this is simple 4-vector addition of $(E, \mathbf{p})$) -- `ptmin` - only jets of transverse momentum greater than or equal to `ptmin` will be returned +- `R` - the cone size parameter; no particles more geometrically distance than `R` will be merged (default 1.0) +- `recombine` - the function used to merge two pseudojets (default is a simple 4-vector addition of $(E, \mathbf{p})$) +- `strategy` - the algorithm strategy to adopt, as described below (default `JetRecoStrategy.Best`) + +The object returned is a `ClusterSequence`, which internally tracks all merge steps. + +To obtain the final inclusive jets, use the `inclusive_jets` method: + +```julia +final_jets = inclusive_jets(cs::ClusterSequence; ptmin=0.0) +``` + +Only jets passing the cut $p_T > p_{Tmin}$ will be returned. The result is returned as a `Vector{LorentzVectorHEP}`. + +#### Sorting + +As sorting vectors is trivial in Julia, no special sorting methods are provided. As an example, to sort exclusive jets of $>5.0$ (usually GeV, depending on your EDM) from highest energy to lowest: -Both of these functions return tuple of `(jets, seq)`, where these are a vector of `JetReconstruction.PseudoJet` objects and cluster sequencing information, respectively. +```julia +sorted_jets = sort!(inclusive_jets(cs::ClusterSequence; ptmin=5.0), by=JetReconstruction.energy, rev=true) +``` + +#### Strategy + +Three strategies are available for the different algorithms: + +| Strategy Name | Notes | Interface | +|---|---|---| +| `JetRecoStrategy.Best` | Dynamically switch strategy based on input particle density | `jet_reconstruct` | +| `JetRecoStrategy.N2Plain` | Global matching of particles at each interation (works well for low $N$) | `plain_jet_reconstruct` | +| `JetRecoStrategy.N2Tiled` | Use tiles of radius $R$ to limit search space (works well for higher $N$) | `tiled_jet_reconstruct` | + +Generally one can use the `jet_reconstruct` interface, shown above, as the *Best* strategy safely as the overhead is extremely low. That interface supports a `strategy` option to switch to a different option. + +Another option, if one wishes to use a specific strategy, is to call that strategy's interface directly, e.g., + +```julia +# N2Plain +plain_jet_reconstruct(particles::Vector{T}; p = -1, R = 1.0, recombine = +) +``` -Note that the `N2Plain` algorithm is usually faster at low densities, $\lessapprox 50$ input particles, otherwise the tiled algorithm is faster. +Note that there is no `strategy` option in these interfaces. ### Example diff --git a/examples/jetreco.jl b/examples/jetreco.jl index 874ef32..f4d4aa2 100644 --- a/examples/jetreco.jl +++ b/examples/jetreco.jl @@ -59,7 +59,7 @@ end """ Top level call funtion for demonstrating the use of jet reconstruction -This uses the "generic_jet_reconstruct" wrapper, so algorithm swutching +This uses the "jet_reconstruct" wrapper, so algorithm switching happens inside the JetReconstruction package itself. Some other ustilities are also supported here, such as profiling and @@ -88,19 +88,19 @@ function jet_process( if nsamples > 1 || !isnothing(profile) @info "Doing initial warm-up run" for event in events - _ = inclusive_jets(generic_jet_reconstruct(event, R = distance, p = power, strategy = strategy), ptmin) + _ = inclusive_jets(jet_reconstruct(event, R = distance, p = power, strategy = strategy), ptmin) end end if !isnothing(profile) - profile_code(profile, generic_jet_reconstruct, events, nsamples, R = distance, p = power, strategy = strategy) + profile_code(profile, jet_reconstruct, events, nsamples, R = distance, p = power, strategy = strategy) return nothing end if alloc println("Memory allocation statistics:") @timev for event in events - _ = inclusive_jets(generic_jet_reconstruct(event, R = distance, p = power, strategy = strategy), ptmin) + _ = inclusive_jets(jet_reconstruct(event, R = distance, p = power, strategy = strategy), ptmin) end return nothing end @@ -114,7 +114,7 @@ function jet_process( gcoff && GC.enable(false) t_start = time_ns() for (ievt, event) in enumerate(events) - finaljets = inclusive_jets(generic_jet_reconstruct(event, R = distance, p = power, strategy = strategy), ptmin) + finaljets = inclusive_jets(jet_reconstruct(event, R = distance, p = power, strategy = strategy), ptmin) # Only print the jet content once if irun == 1 @info begin diff --git a/src/ClusterSequence.jl b/src/ClusterSequence.jl index 74c9e8d..9224dec 100644 --- a/src/ClusterSequence.jl +++ b/src/ClusterSequence.jl @@ -76,25 +76,6 @@ struct ClusterSequence Qtot end -"""Return all inclusive jets of a ClusterSequence with pt > ptmin""" -function inclusive_jets(clusterseq::ClusterSequence, ptmin = 0.0) - dcut = ptmin * ptmin - jets_local = Vector{LorentzVectorCyl}(undef, 0) - # sizehint!(jets_local, length(clusterseq.jets)) - # For inclusive jets with a plugin algorithm, we make no - # assumptions about anything (relation of dij to momenta, - # ordering of the dij, etc.) - # for elt in Iterators.reverse(clusterseq.history) - for elt in clusterseq.history - elt.parent2 == BeamJet || continue - iparent_jet = clusterseq.history[elt.parent1].jetp_index - jet = clusterseq.jets[iparent_jet] - if pt2(jet) >= dcut - push!(jets_local, LorentzVectorCyl(pt(jet), rapidity(jet), phi(jet), mass(jet))) - end - end - jets_local -end """Add a new jet's history into the recombination sequence""" add_step_to_history!(clusterseq::ClusterSequence, parent1, parent2, jetp_index, dij) = begin @@ -133,3 +114,24 @@ add_step_to_history!(clusterseq::ClusterSequence, parent1, parent2, jetp_index, clusterseq.jets[jetp_index]._cluster_hist_index = local_step end end + +"""Return all inclusive jets of a ClusterSequence with pt > ptmin""" +function inclusive_jets(clusterseq::ClusterSequence, ptmin = 0.0) + dcut = ptmin * ptmin + jets_local = Vector{LorentzVectorCyl}(undef, 0) + # sizehint!(jets_local, length(clusterseq.jets)) + # For inclusive jets with a plugin algorithm, we make no + # assumptions about anything (relation of dij to momenta, + # ordering of the dij, etc.) + # for elt in Iterators.reverse(clusterseq.history) + for elt in clusterseq.history + elt.parent2 == BeamJet || continue + iparent_jet = clusterseq.history[elt.parent1].jetp_index + jet = clusterseq.jets[iparent_jet] + if pt2(jet) >= dcut + push!(jets_local, LorentzVectorCyl(pt(jet), rapidity(jet), phi(jet), mass(jet))) + end + end + jets_local +end + diff --git a/src/GenericAlgo.jl b/src/GenericAlgo.jl index bd4eec5..d55a87d 100644 --- a/src/GenericAlgo.jl +++ b/src/GenericAlgo.jl @@ -2,7 +2,7 @@ # switch based on the strategy, or based on the event density # if the "Best" strategy is to be employed -function generic_jet_reconstruct(particles; p = -1, R = 1.0, recombine = +, strategy = JetRecoStrategy.Best) +function jet_reconstruct(particles; p = -1, R = 1.0, recombine = +, strategy = JetRecoStrategy.Best) # Either map to the fixed algorithm corresponding to the strategy # or to an optimal choice based on the density of initial particles diff --git a/src/JetReconstruction.jl b/src/JetReconstruction.jl index 597aac5..ae37791 100644 --- a/src/JetReconstruction.jl +++ b/src/JetReconstruction.jl @@ -48,7 +48,7 @@ export tiled_jet_reconstruct ## Generic algorithm, which can switch strategy dynamically include("GenericAlgo.jl") -export generic_jet_reconstruct +export jet_reconstruct # Simple HepMC3 reader include("HepMC3.jl") diff --git a/test/runtests.jl b/test/runtests.jl index acbcb73..730c4a3 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -91,6 +91,9 @@ function do_test_compare_to_fastjet(strategy::JetRecoStrategy.Strategy, fastjet_ elseif (strategy == JetRecoStrategy.N2Tiled) jet_reconstruction = tiled_jet_reconstruct strategy_name = "N2Tiled" + elseif (strategy == JetRecoStrategy.Best) + jet_reconstruction = jet_reconstruct + strategy_name = "Best" else throw(ErrorException("Strategy not yet implemented")) end