Durham algorithm (#73)

Implementation of the e+e- k_T algorithm, aka Durham algorithm.

This is implemented in a new specific implementation as the distance metrics are rather different from those used for pp events. In particular, the reconstruction happens in (theta, phi) space.

A new EEjet structure is used, that stores the 4-momentum, but with different internal cache variables for optimising the distance calculations. (In particular the normalised momentum is cached, as the dot product is used to calculate the angle between two jets). Various accessors are added for components of the struct, as well as derived values, such as phi and rapidity.

initial_history() is moved to the ClusterSequence.jl file, which is where it should be (used by all algorithms/strategies).
ClusterSequence structs are updated to take jets as vectors of the abstract type FourMomentum.

The algorithm stores the vector of EEjets in the ClusterSequence as before (original and all intermediate jets) and optimised compact arrays for finding the closest jets are used.

Currently the implementation is about x3 slower than FastJet (which doesn't use any compact arrays). Timing on M2Pro for the ee->H input file per event is:

- JetReconstruction.jl Durham 59us
- FastJet Durham 21us
- JetReconstruction.jl AntiKt/N2Plain 33us

The consistency handling of p and algorithm is updated to take this new algorithm into account and new functions are added that map algorithms to either pp or e+e- cases (is_pp() and is_ee()).

Documentation is updated to favour passing the algorithm to the reconstruction instead of implying it from the p value (though that mode is still supported and will imply pp reconstruction).

There is also a major refactoring of tests to make testing Julia reconstruction against FastJet much easier and more independent of assumptions about the p <-> alg mapping. A new test file with e+e- -> H events is added and the 13TeV pp event file is renamed.

Author: graeme-a-stewart

README.md

@@ -12,41 +12,50 @@ Algorithms used are based on the C++ FastJet package (<https://fastjet.fr>,
 [arXiv:1111.6097](https://arxiv.org/abs/1111.6097)), reimplemented natively in Julia.
 
 The algorithms include anti-$`{k}_\text{T}`$, Cambridge/Aachen, inclusive
-$`k_\text{T}`$ and generalised $`k_\text{T}`$.
+$`k_\text{T}`$, generalised $`k_\text{T}`$ for pp events; and the Durham algorithm for e+e-.
 
 ### Interface
 
 The simplest interface is to call:
 
 ```julia
-cs = jet_reconstruct(particles::Vector{T}; p = -1, R = 1.0, recombine = +, strategy = RecoStrategy.Best)
+cs = jet_reconstruct(particles::Vector{T}; algorithm = JetAlgorithm.AntiKt, R = 1.0, [p = -1,] [recombine = +,] [strategy = RecoStrategy.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}$. Different values of $p$ then give different reconstruction algorithms:
-  - `-1` gives anti-$`{k}_\text{T}`$ clustering (default)
-  - `0` gives Cambridge/Aachen
-  - `1` gives inclusive $k_\text{T}$
-- `R` - the cone size parameter; no particles more geometrically distance than `R` will be merged (default 1.0)
+- `algorithm` is the name of the jet algorithm to be used (from the `JetAlgorithm` enum)
+  - `JetAlgorithm.AntiKt` anti-$`{k}_\text{T}`$ clustering (default)
+  - `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
+- `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)
 - `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 `RecoStrategy.Best`)
 
 The object returned is a `ClusterSequence`, which internally tracks all merge steps.
 
-Alternatively one can swap the `p=...` parameter for
-`algorithm=JetReconstruction.{AntiKt,CA,Kt}` for explicit algorithm selection. (Generalised $`{k}_\text{T}`$ requires `algorithm=JetReconstruction.GenKt` *and* `p=N`.)
+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}$
+
+(for the `GenKt` algorithm the `p` value *must* also be given to specify the algorithm fully).
 
 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}`.
-
 
+Only jets passing the cut $p_T > p_{Tmin}$ will be returned. The result is returned as a `Vector{LorentzVectorHEP}`.
 
 #### Sorting
 
@@ -72,7 +81,7 @@ Another option, if one wishes to use a specific strategy, is to call that strate
 
 ```julia
 # For N2Plain strategy called directly
-plain_jet_reconstruct(particles::Vector{T}; p = -1, R = 1.0, recombine = +)
+plain_jet_reconstruct(particles::Vector{T}; algorithm = JetAlgorithm.AntiKt, R = 1.0, recombine = +)
 ```
 
 Note that there is no `strategy` option in these interfaces.
@@ -84,9 +93,11 @@ In the examples directory there are a number of example scripts.
 See the `jetreco.jl` script for an example of how to call jet reconstruction.
 
 ```sh
-julia --project=. examples/jetreco.jl --maxevents=100 --nsamples=1 --strategy=N2Plain test/data/events.hepmc3
+julia --project=examples examples/jetreco.jl --algorithm=AntiKt test/data/events.pp13TeV.hepmc3.gz
+...
+julia --project=examples examples/jetreco.jl --algorithm=Durham test/data/events.eeH.hepmc3.gz
 ...
-julia --project=. examples/jetreco.jl --maxevents=100 --nsamples=1 --strategy=N2Tiled test/data/events.hepmc3
+julia --project=examples examples/jetreco.jl --maxevents=10 --strategy=N2Plain --algorithm=Kt --exclusive-njets=3 test/data/events.pp13TeV.hepmc3.gz
 ...
 ```
 

docs/src/examples.md

@@ -13,9 +13,11 @@ perform a jet reconstruction, with different algorithms and (optionally)
 strategy, producing exclusive and inclusive jet selections.
 
 ```sh
-julia --project=. examples/jetreco.jl --maxevents=100 --strategy=N2Plain test/data/events.hepmc3
+julia --project=examples examples/jetreco.jl --algorithm=AntiKt test/data/events.pp13TeV.hepmc3.gz
 ...
-julia --project=. examples/jetreco.jl --maxevents=100 --strategy=N2Tiled test/data/events.hepmc3
+julia --project=examples examples/jetreco.jl --algorithm=Durham test/data/events.eeH.hepmc3.gz
+...
+julia --project=examples examples/jetreco.jl --maxevents=10 --strategy=N2Plain --algorithm=Kt --exclusive-njets=3 test/data/events.pp13TeV.hepmc3.gz
 ...
 ```
 

docs/src/index.md

@@ -50,7 +50,9 @@ jet_reconstruct(particles; algorithm = JetAlgorithm.AntiKt, R = 1.0)
 ```
 
 Where `particles` is a collection of 4-vector objects to reconstruct and the
-algorithm is either given explicitly or implied by the power value. For the case of generalised $k_T$ both the algorithm (`JetAlgorithm.GenKt`) and `p` are needed.
+algorithm is either given explicitly or implied by the power value. For the case
+of generalised $k_T$ both the algorithm (`JetAlgorithm.GenKt`) and `p` are
+needed. For the case of the Durham algorithm the `R` value is ignored.
 
 The object returned is a [`ClusterSequence`](@ref), which internally tracks all
 merge steps.

examples/instrumented-jetreco.jl

@@ -97,7 +97,7 @@ function jet_process(events::Vector{Vector{PseudoJet}};
         for event in events
             _ = inclusive_jets(jet_reconstruct(event, R = distance, p = p,
                                                algorithm = algorithm,
-                                               strategy = strategy), ptmin)
+                                               strategy = strategy); ptmin = ptmin)
         end
     end
 
@@ -113,7 +113,7 @@ function jet_process(events::Vector{Vector{PseudoJet}};
         @timev for event in events
             _ = inclusive_jets(jet_reconstruct(event, R = distance, p = p,
                                                algorithm = algorithm,
-                                               strategy = strategy), ptmin)
+                                               strategy = strategy), ptmin = ptmin)
         end
         return nothing
     end
@@ -129,13 +129,12 @@ function jet_process(events::Vector{Vector{PseudoJet}};
         for (ievt, event) in enumerate(events)
             cs = jet_reconstruct(event, R = distance, p = p, algorithm = algorithm,
                                  strategy = strategy)
-            strategy = strategy
             if !isnothing(njets)
-                finaljets = exclusive_jets(cs, njets = njets)
+                finaljets = exclusive_jets(cs; njets = njets)
             elseif !isnothing(dcut)
-                finaljets = exclusive_jets(cs, dcut = dcut)
+                finaljets = exclusive_jets(cs; dcut = dcut)
             else
-                finaljets = inclusive_jets(cs, ptmin = ptmin)
+                finaljets = inclusive_jets(cs; ptmin = ptmin)
             end
             # Only print the jet content once
             if irun == 1

src/AlgorithmStrategyEnums.jl

@@ -70,7 +70,7 @@ Base.tryparse(E::Type{<:Enum}, str::String) =
     end
 
 """
-    set_algorithm_power_consistency(; p::Union{Real, Nothing}, algorithm::Union{JetAlgorithm, Nothing})
+    get_algorithm_power_consistency(; p::Union{Real, Nothing}, algorithm::Union{JetAlgorithm, Nothing})
 
 Get the algorithm and power consistency correct
 
@@ -103,6 +103,11 @@ function get_algorithm_power_consistency(; p::Union{Real, Nothing},
         return (p = p, algorithm = algorithm)
     end
 
+    # Some algorithms have fixed power values
+    if algorithm == JetAlgorithm.Durham
+        return (p = 1, algorithm = algorithm)
+    end
+
     # Otherwise we check the consistency between the algorithm and power
     if !isnothing(algorithm)
         power_from_alg = algorithm2power[algorithm]
@@ -127,3 +132,27 @@ function check_algorithm_power_consistency(; p::Union{Real, Nothing},
     get_algorithm_power_consistency(p = p, algorithm = algorithm)
     return true
 end
+
+"""
+    is_pp(algorithm::JetAlgorithm.Algorithm)
+
+Check if the algorithm is a pp reconstruction algorithm.
+
+# Returns
+`true` if the algorithm is a pp reconstruction algorithm, `false` otherwise.
+"""
+function is_pp(algorithm::JetAlgorithm.Algorithm)
+    return algorithm in [JetAlgorithm.AntiKt, JetAlgorithm.CA, JetAlgorithm.Kt]
+end
+
+"""
+    is_ee(algorithm::JetAlgorithm.Algorithm)
+
+Check if the algorithm is a e+e- reconstruction algorithm.
+
+# Returns
+`true` if the algorithm is a e+e- reconstruction algorithm, `false` otherwise.
+"""
+function is_ee(algorithm::JetAlgorithm.Algorithm)
+    return algorithm in [JetAlgorithm.EEKt, JetAlgorithm.Durham]
+end

src/ClusterSequence.jl

@@ -65,6 +65,37 @@ A `HistoryElement` object.
 HistoryElement(jetp_index) = HistoryElement(NonexistentParent, NonexistentParent, Invalid,
                                             jetp_index, 0.0, 0.0)
 
+"""
+    initial_history(particles)
+
+Create an initial history for the given particles.
+
+# Arguments
+- `particles`: The initial vector of stable particles.
+
+# Returns
+- `history`: An array of `HistoryElement` objects.
+- `Qtot`: The total energy in the event.
+"""
+function initial_history(particles)
+    # reserve sufficient space for everything
+    history = Vector{HistoryElement}(undef, length(particles))
+    sizehint!(history, 2 * length(particles))
+
+    Qtot::Float64 = 0
+
+    for i in eachindex(particles)
+        history[i] = HistoryElement(i)
+
+        # get cross-referencing right from the Jets
+        particles[i]._cluster_hist_index = i
+
+        # determine the total energy in the event
+        Qtot += particles[i].E
+    end
+    history, Qtot
+end
+
 """
     struct ClusterSequence
 
@@ -89,7 +120,7 @@ struct ClusterSequence
     algorithm::JetAlgorithm.Algorithm
     power::Float64
     strategy::RecoStrategy.Strategy
-    jets::Vector{PseudoJet}
+    jets::Vector{FourMomentum}
     n_initial_jets::Int
     history::Vector{HistoryElement}
     Qtot::Any