Generalised kT Algorithm for pp collisions (#71)

Add proper support for the generalised kT algorithm for pp collisions (i.e., operating in $(y, \phi)$ space).

Added a `JetAlgorithm` type for generalised kT, viz. `GenKt`.

Strategies are restructured to be able to be driven by only the power, or the algorithm, or (for `GenKt`) both. A new function, `get_algorithm_power_consistency` is added that will return the "partner" of a (power/algorithm) given an (algorithm/power). If inconsistent input is given an exception is thrown.

Documentation is updated explaining how to call with either `p` or `algorithm`, or both for `GenKt`.

ClusterSequence now stores the power value and the algorithm and has to be initialised with both. In addition the power is always stored as a `Float64`.

Examples are restructured to also be able to be driven by power, or algorithm, or both. However, if neither is given then the scripts will use `AntiKt` and print a warning. There is also now a clearer separation between calling the jet reconstruction itself and getting the inclusive/exclusive constituents.

A test to compare Generalised kT output with FastJet for the case `p=1.5` is added.

Tests have been added to run our code examples, to check none of these break (at least for the standalone Julia scripts), but these are currently disabled as some non-trivial problems arose:

- The dependencies are heavy and some for GLMakie actually seem to timeout
- As these are CI machines, they can't correctly setup GLMakie
- There seems to be a weird issue with setting up a different environment for the Julia nightly:

> ERROR: LoadError: ArgumentError: Package Pkg does not have Random in its dependencies

# Misc

- Notebooks are now *not* formatted by JuliaFormatter.
- For the tiled algorithm directly calling the real work interface is now discouraged (method is now `_tiled_jet_reconstruct`) to ensure that power and algorithm are handled properly in the public interface (`tiled_jet_reconstruct`). (N.B. this costs about 1us per event)

Closes #49

.JuliaFormatter.toml

@@ -1,2 +1,5 @@
 style = "sciml"
 yas_style_nesting = true
+
+# Do not format notebooks
+ignore = ["*-pluto.jl", "*-nb.jl"]

README.md

@@ -11,7 +11,8 @@ 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.
 
-The algorithms include anti-$`{k}_\text{T}`$, Cambridge/Aachen and inclusive $`k_\text{T}`$.
+The algorithms include anti-$`{k}_\text{T}`$, Cambridge/Aachen, inclusive
+$`k_\text{T}`$ and generalised $`k_\text{T}`$.
 
 ### Interface
 
@@ -35,14 +36,18 @@ cs = jet_reconstruct(particles::Vector{T}; p = -1, R = 1.0, recombine = +, strat
 
 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`.)
+
 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:

docs/src/index.md

@@ -43,7 +43,14 @@ The main interface for reconstruction is [`jet_reconstruct`](@ref), called as, e
 jet_reconstruct(particles; p = -1, R = 1.0)
 ```
 
-Where `particles` is a collection of 4-vector objects to reconstruct.
+or
+
+```julia
+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.
 
 The object returned is a [`ClusterSequence`](@ref), which internally tracks all
 merge steps.

examples/Project.toml

@@ -7,8 +7,10 @@ GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 JetReconstruction = "44e8cb2c-dfab-4825-9c70-d4808a591196"
 LorentzVectorHEP = "f612022c-142a-473f-8cfd-a09cf3793c6c"
+Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Pluto = "c3e4b0f8-55cb-11ea-2926-15256bba5781"
 PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8"
 Profile = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"
 ProfileSVG = "132c30aa-f267-4189-9183-c8a63c7e05e6"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 StatProfilerHTML = "a8a75453-ed82-57c9-9e16-4cd1196ecbf5"

examples/animate-reconstruction.jl

@@ -26,7 +26,10 @@ function main()
         "--algorithm", "-A"
         help = """Algorithm to use for jet reconstruction: $(join(JetReconstruction.AllJetRecoAlgorithms, ", "))"""
         arg_type = JetAlgorithm.Algorithm
-        default = JetAlgorithm.AntiKt
+
+        "--power", "-p"
+        help = """Power value for jet reconstruction"""
+        arg_type = Float64
 
         "--strategy", "-S"
         help = """Strategy for the algorithm, valid values: $(join(JetReconstruction.AllJetRecoStrategies, ", "))"""
@@ -49,9 +52,14 @@ function main()
     events::Vector{Vector{PseudoJet}} = read_final_state_particles(args[:file],
                                                                    maxevents = args[:event],
                                                                    skipevents = args[:event])
-
-    power = JetReconstruction.algorithm2power[args[:algorithm]]
-    cs = jet_reconstruct(events[1], R = args[:distance], p = power,
+    if isnothing(args[:algorithm]) && isnothing(args[:power])
+        @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,
                          strategy = args[:strategy])
 
     animatereco(cs, args[:output]; azimuth = (1.8, 3.0), elevation = 0.5,

examples/instrumented-jetreco.jl

@@ -23,11 +23,13 @@ using JetReconstruction
 include(joinpath(@__DIR__, "parse-options.jl"))
 
 function profile_code(profile, jet_reconstruction, events, niters; R = 0.4, p = -1,
+                      algorithm::JetAlgorithm.Algorithm = JetAlgorithm.AntiKt,
                       strategy = RecoStrategy.N2Tiled)
     Profile.init(n = 5 * 10^6, delay = 0.00001)
     function profile_events(events)
         for evt in events
-            jet_reconstruction(evt, R = R, p = p, strategy = strategy)
+            jet_reconstruction(evt, R = R, p = p, algorithm = algorithm,
+                               strategy = strategy)
         end
     end
     profile_events(events[1:1])
@@ -67,7 +69,8 @@ serialising the reconstructed jet outputs.
 """
 function jet_process(events::Vector{Vector{PseudoJet}};
                      distance::Real = 0.4,
-                     algorithm::JetAlgorithm.Algorithm = JetAlgorithm.AntiKt,
+                     algorithm::Union{JetAlgorithm.Algorithm, Nothing} = nothing,
+                     p::Union{Real, Nothing} = nothing,
                      ptmin::Real = 5.0,
                      dcut = nothing,
                      njets = nothing,
@@ -77,36 +80,40 @@ function jet_process(events::Vector{Vector{PseudoJet}};
                      profile = nothing,
                      alloc::Bool = false,
                      dump::Union{String, Nothing} = nothing)
-    @info "Will process $(size(events)[1]) events"
-
-    # Map algorithm to power
-    power = JetReconstruction.algorithm2power[algorithm]
 
     # If we are dumping the results, setup the JSON structure
     if !isnothing(dump)
         jet_collection = FinalJets[]
     end
 
+    # Set consistent algorithm and power
+    (p, algorithm) = JetReconstruction.get_algorithm_power_consistency(p = p,
+                                                                       algorithm = algorithm)
+    @info "Jet reconstruction will use $(algorithm) with power $(p)"
+
     # Warmup code if we are doing a multi-sample timing run
     if nsamples > 1 || !isnothing(profile)
         @info "Doing initial warm-up run"
         for event in events
-            _ = inclusive_jets(jet_reconstruct(event, R = distance, p = power,
-                                               strategy = strategy); ptmin = ptmin)
+            _ = inclusive_jets(jet_reconstruct(event, R = distance, p = p,
+                                               algorithm = algorithm,
+                                               strategy = strategy), ptmin)
         end
     end
 
     if !isnothing(profile)
-        profile_code(profile, jet_reconstruct, events, nsamples, R = distance, p = power,
+        profile_code(profile, jet_reconstruct, events, nsamples, algorithm = algorithm,
+                     R = distance, p = p,
                      strategy = strategy)
         return nothing
     end
 
     if alloc
         println("Memory allocation statistics:")
         @timev for event in events
-            _ = inclusive_jets(jet_reconstruct(event, R = distance, p = power,
-                                               strategy = strategy); ptmin = ptmin)
+            _ = inclusive_jets(jet_reconstruct(event, R = distance, p = p,
+                                               algorithm = algorithm,
+                                               strategy = strategy), ptmin)
         end
         return nothing
     end
@@ -120,18 +127,15 @@ function jet_process(events::Vector{Vector{PseudoJet}};
         gcoff && GC.enable(false)
         t_start = time_ns()
         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(jet_reconstruct(event, R = distance, p = power,
-                                                           strategy = strategy);
-                                           njets = njets)
+                finaljets = exclusive_jets(cs, njets = njets)
             elseif !isnothing(dcut)
-                finaljets = exclusive_jets(jet_reconstruct(event, R = distance, p = power,
-                                                           strategy = strategy);
-                                           dcut = dcut)
+                finaljets = exclusive_jets(cs, dcut = dcut)
             else
-                finaljets = inclusive_jets(jet_reconstruct(event, R = distance, p = power,
-                                                           strategy = strategy);
-                                           ptmin = ptmin)
+                finaljets = inclusive_jets(cs, ptmin = ptmin)
             end
             # Only print the jet content once
             if irun == 1
@@ -221,7 +225,10 @@ function parse_command_line(args)
         "--algorithm", "-A"
         help = """Algorithm to use for jet reconstruction: $(join(JetReconstruction.AllJetRecoAlgorithms, ", "))"""
         arg_type = JetAlgorithm.Algorithm
-        default = JetAlgorithm.AntiKt
+
+        "--power", "-p"
+        help = """Power value for jet reconstruction"""
+        arg_type = Float64
 
         "--strategy", "-S"
         help = """Strategy for the algorithm, valid values: $(join(JetReconstruction.AllJetRecoStrategies, ", "))"""
@@ -277,7 +284,12 @@ function main()
     events::Vector{Vector{PseudoJet}} = read_final_state_particles(args[:file],
                                                                    maxevents = args[:maxevents],
                                                                    skipevents = args[:skip])
+    if isnothing(args[:algorithm]) && isnothing(args[:power])
+        @warn "Neither algorithm nor power specified, defaulting to AntiKt"
+        args[:algorithm] = JetAlgorithm.AntiKt
+    end
     jet_process(events, distance = args[:distance], algorithm = args[:algorithm],
+                p = args[:power],
                 strategy = args[:strategy],
                 ptmin = args[:ptmin], dcut = args[:exclusive_dcut],
                 njets = args[:exclusive_njets],

examples/jetreco-constituents.jl

@@ -5,6 +5,10 @@
 # N.B. currently you must use the `jet-constituents` branch of `JetReconstruction`.
 using JetReconstruction
 using LorentzVectorHEP
+using Logging
+
+logger = ConsoleLogger(stdout, Logging.Info)
+global_logger(logger)
 
 input_file = joinpath(dirname(pathof(JetReconstruction)), "..", "test", "data",
                       "events.hepmc3.gz")

examples/jetreco.jl

@@ -25,35 +25,32 @@ Final jets can be serialised if the "dump" option is given
 """
 function jet_process(events::Vector{Vector{PseudoJet}};
                      distance::Real = 0.4,
-                     algorithm::JetAlgorithm.Algorithm = JetAlgorithm.AntiKt,
+                     algorithm::Union{JetAlgorithm.Algorithm, Nothing} = nothing,
+                     p::Union{Real, Nothing} = nothing,
                      ptmin::Real = 5.0,
                      dcut = nothing,
                      njets = nothing,
                      strategy::RecoStrategy.Strategy,
                      dump::Union{String, Nothing} = nothing)
-    @info "Will process $(size(events)[1]) events"
 
-    # If we are dumping the results, setup the JSON structure
-    if !isnothing(dump)
-        jet_collection = FinalJets[]
-    end
+    # Set consistent algorithm and power
+    (p, algorithm) = JetReconstruction.get_algorithm_power_consistency(p = p,
+                                                                       algorithm = algorithm)
+    @info "Jet reconstruction will use $(algorithm) with power $(p)"
 
     # A friendly label for the algorithm and final jet selection
     if !isnothing(njets)
-        @info "Running exclusive jets with n_jets = $(njets)"
+        @info "Will produce exclusive jets with n_jets = $(njets)"
     elseif !isnothing(dcut)
-        @info "Running exclusive jets with dcut = $(dcut)"
+        @info "Will produce exclusive jets with dcut = $(dcut)"
     else
-        @info "Running inclusive jets with ptmin = $(ptmin)"
+        @info "Will produce inclusive jets with ptmin = $(ptmin)"
     end
 
-    # Map algorithm to power
-    power = JetReconstruction.algorithm2power[algorithm]
-
     # Now run over each event
     for (ievt, event) in enumerate(events)
         # Run the jet reconstruction
-        cluster_seq = jet_reconstruct(event, R = distance, p = power,
+        cluster_seq = jet_reconstruct(event, R = distance, p = p, algorithm = algorithm,
                                       strategy = strategy)
         # Now select jets, with inclusive or exclusive parameters
         if !isnothing(njets)
@@ -117,7 +114,10 @@ function parse_command_line(args)
         "--algorithm", "-A"
         help = """Algorithm to use for jet reconstruction: $(join(JetReconstruction.AllJetRecoAlgorithms, ", "))"""
         arg_type = JetAlgorithm.Algorithm
-        default = JetAlgorithm.AntiKt
+
+        "--power", "-p"
+        help = """Power value for jet reconstruction"""
+        arg_type = Float64
 
         "--strategy", "-S"
         help = """Strategy for the algorithm, valid values: $(join(JetReconstruction.AllJetRecoStrategies, ", "))"""
@@ -141,7 +141,12 @@ function main()
     events::Vector{Vector{PseudoJet}} = read_final_state_particles(args[:file],
                                                                    maxevents = args[:maxevents],
                                                                    skipevents = args[:skip])
+    if isnothing(args[:algorithm]) && isnothing(args[:power])
+        @warn "Neither algorithm nor power specified, defaulting to AntiKt"
+        args[:algorithm] = JetAlgorithm.AntiKt
+    end
     jet_process(events, distance = args[:distance], algorithm = args[:algorithm],
+                p = args[:power],
                 strategy = args[:strategy],
                 ptmin = args[:ptmin], dcut = args[:exclusive_dcut],
                 njets = args[:exclusive_njets],

examples/visualise-jets.ipynb

@@ -26,8 +26,7 @@
    "outputs": [],
    "source": [
     "input_file = joinpath(dirname(pathof(JetReconstruction)), \"..\", \"test\", \"data\", \"events.hepmc3.gz\")\n",
-    "event_number = 5\n",
-    "power = JetReconstruction.algorithm2power[JetAlgorithm.AntiKt];"
+    "event_number = 5"
    ]
   },
   {
@@ -47,7 +46,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "cs = jet_reconstruct(events[1], R = 1.0, p = power);"
+    "cs = jet_reconstruct(events[1], R = 1.0, algorithm = JetAlgorithm.AntiKt);"
    ]
   },
   {

examples/visualise-jets.jl

@@ -39,7 +39,10 @@ function main()
         "--algorithm", "-A"
         help = """Algorithm to use for jet reconstruction: $(join(JetReconstruction.AllJetRecoAlgorithms, ", "))"""
         arg_type = JetAlgorithm.Algorithm
-        default = JetAlgorithm.AntiKt
+
+        "--power", "-p"
+        help = """Power value for jet reconstruction"""
+        arg_type = Float64
 
         "--strategy", "-S"
         help = """Strategy for the algorithm, valid values: $(join(JetReconstruction.AllJetRecoStrategies, ", "))"""
@@ -63,8 +66,9 @@ function main()
                                                                    maxevents = args[:event],
                                                                    skipevents = args[:event])
 
-    power = JetReconstruction.algorithm2power[args[:algorithm]]
-    cs = jet_reconstruct(events[1], R = args[:distance], p = 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,
                          strategy = args[:strategy])
 
     plt = jetsplot(events[1], cs; Module = CairoMakie)