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

@@ -10,38 +10,66 @@ Algorithms used are based on the C++ FastJet package (<https://fastjet.fr>,
 [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
   - 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` 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
 

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

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
+

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
 

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")

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