Valencia jets revisited 🍊 (#217)

* Rebased VLC implementation

Add the Valencia algorithm to the general Durham/EEKt implementation.
Uses internal algorithm switching at key points inside the distance
and update functions to add different Valencia calculations.

* Allow Valencia algorithm from CLI

Add --gamma to the reconstruction examples and pass it
through to the underlying algorithms correctly.
Add Valencia test to benchmark script

* Interface fixes after rebase

* Remove algorithm check

No need to check that algorithm is the correct type this deep in the
reconstruction. This is also a hot function, so unnecessary checks
are costly.

* Support --beta as an alias for -p/--power

* Updated documentation for Valencia

Add information on the Valencia algorithm to the documentation.

Update the reference to latest (CHEP 2024) paper.

Minor documentation tweaks.

* Add Valencia to benchmarks

Minor fixes to docstrings

* Update CITATION.cff

* Rewind citation changes

These will now be handled in a separate PR

* Improvements from review comments

---------

Co-authored-by: Matt LeBlanc <matt.leblanc@cern.ch>

Author: graeme-a-stewart

CITATION.cff

@@ -29,26 +29,18 @@ identifiers:
     description: Zenodo
 preferred-citation:
   type: article
-  title: Fast Jet Finding in Julia
-  doi: 10.1051/epjconf/202533701067
+  title: Polyglot Jet Finding
+  doi: 10.1051/epjconf/202429505017
   journal: EPJ Web Conf.
-  volume: 337
-  start: "01067"
-  year: 2025
+  volume: 295
+  start: "05017"
+  year: 2024
   authors:
     - given-names: Graeme Andrew
       family-names: Stewart
       email: graeme.andrew.stewart@cern.ch
       affiliation: CERN, Switzerland
       orcid: 'https://orcid.org/0000-0003-0182-7088'
-    - given-names: Sanmay
-      family-names: Ganguly
-      email: Sanmay.Ganguly@cern.ch
-      affiliation: Indian Institute of Technology, Kanpur, India
-    - given-names: Sattwamo
-      family-names: Ghosh
-      email: sg21ms204@iiserkol.ac.in
-      affiliation: Indian Institute of Science Education and Research, Kanpur, India
     - given-names: Atell
       family-names: Krasnopolski
       email: delta_atell@protonmail.com
@@ -59,6 +51,11 @@ preferred-citation:
       email: philippe.gras@cern.ch
       affiliation: IRFU, CEA, Université Paris-Saclay, France
       orcid: 'https://orcid.org/0000-0002-3932-5967'
+    - given-names: Benedikt
+      family-names: Hegner
+      email: benedikt.hegner@cern.ch
+      affiliation: CERN
+      orcid: 'https://orcid.org/0009-0009-2020-1620'
 repository-code: 'https://github.com/JuliaHEP/JetReconstruction.jl'
 abstract: Jet reconstruction (reclustering) with Julia
 license: MIT

README.md

@@ -41,6 +41,7 @@ cs = jet_reconstruct(particles::AbstractVector{T}; algorithm = JetAlgorithm.Anti
   - `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 (which also requires specification of `p`)
+  - `JetAlgorithm.Valencia` the Valencia $e^+e-$ pile-up resistant algorithm (which also requires specification of `p` (γ) and β)
 - `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)
 
 The object returned is a `ClusterSequence`, which internally tracks all merge steps.
@@ -65,15 +66,15 @@ sorted_jets = sort!(inclusive_jets(cs::ClusterSequence; ptmin=5.0), by=JetRecons
 
 #### Strategy
 
-Three strategies are available for the different algorithms:
+Three strategies are available for the different algorithms, but only make a difference for $`pp`$ events:
 
 | Strategy Name | Notes | Interface |
 |---|---|---|
 | `RecoStrategy.Best` | Dynamically switch strategy based on input particle density | `jet_reconstruct` |
 | `RecoStrategy.N2Plain` | Global matching of particles at each interaction (works well for low $N$) | `plain_jet_reconstruct` |
 | `RecoStrategy.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.
+Generally one can use the `jet_reconstruct` interface, shown above, as the *Best* strategy safely as the overhead is utterly negligible. 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.,
 
@@ -98,6 +99,7 @@ julia --project jetreco.jl --algorithm=Durham ../test/data/events.eeH.hepmc3.zst
 ...
 julia --project jetreco.jl --maxevents=10 --strategy=N2Plain --algorithm=Kt --exclusive-njets=3 ../test/data/events.pp13TeV.hepmc3.zst
 ...
+julia --project ./jetreco.jl --algorithm=Valencia --gamma=1.2 --beta=1.2 -R 0.8 ../test/data/events.eeH.hepmc3.zst
 ```
 
 There are options to explicitly set the algorithm (use `--help` to see these).
@@ -125,28 +127,25 @@ modules is loaded in the environment.
 
 ## Reference
 
-
-The current recommended reference for JetReconstruction.jl is:
+Although it has been developed further since the CHEP2023 conference, the CHEP
+conference proceedings, [arXiv:2309.17309](https://arxiv.org/abs/2309.17309),
+should be cited if you use this package:
 
 ```bibtex
-@article{refId0,
-  author = {{Stewart, Graeme Andrew} and {Ganguly, Sanmay} and {Ghosh, Sattwamo} and {Gras, Philippe} and {Krasnopolski, Atell}},
-  title = {Fast Jet Finding in Julia},
-  DOI= "10.1051/epjconf/202533701067",
-  url= "https://doi.org/10.1051/epjconf/202533701067",
-  journal = {EPJ Web Conf.},
-  year = 2025,
-  volume = 337,
-  pages = "01067",
+@article{Stewart:2023lgt,
+    author = "Stewart, Graeme Andrew and Gras, Philippe and Hegner, Benedikt and Krasnopolski, Atell",
+    title = "{Polyglot Jet Finding}",
+    eprint = "2309.17309",
+    archivePrefix = "arXiv",
+    primaryClass = "hep-ex",
+    doi = "10.1051/epjconf/202429505017",
+    journal = "EPJ Web Conf.",
+    volume = "295",
+    pages = "05017",
+    year = "2024"
 }
 ```
 
-Also available as [arXiv:2503.08146](https://arxiv.org/abs/2503.08146).
-
-### Other Articles
-
-- CHEP2023, *Polyglot Jet Finding*: [arXiv:2309.17309](https://arxiv.org/abs/2309.17309), [10.1051/epjconf/202429505017](https://doi.org/10.1051/epjconf)
-
 ## Copyright
 
 Code in this package is Copyright 2022-2025 Graeme A Stewart, Philippe Gras, Atell Krasnopolski, CERN.

benchmark/benchmarks.jl

@@ -9,12 +9,13 @@ const events_file_ee = joinpath(@__DIR__, "..", "test", "data", "events.eeH.hepm
 const pp_events = JetReconstruction.read_final_state_particles(events_file_pp, PseudoJet)
 const ee_events = JetReconstruction.read_final_state_particles(events_file_ee, EEJet)
 
-function jet_reconstruct_harness(events; algorithm, strategy, distance = nothing,
-                                 power = nothing,
+function jet_reconstruct_harness(events; algorithm, strategy = RecoStrategy.Best,
+                                 distance = nothing,
+                                 power = nothing, γ = nothing,
                                  recombine = RecombinationMethods[RecombinationScheme.EScheme],
                                  ptmin::Real = 5.0, dcut = nothing, njets = nothing,)
     for event in events
-        cs = jet_reconstruct(event; R = distance, p = power, algorithm = algorithm,
+        cs = jet_reconstruct(event; R = distance, p = power, γ = γ, algorithm = algorithm,
                              strategy = strategy, recombine...)
         if !isnothing(njets)
             finaljets = exclusive_jets(cs; njets = njets)
@@ -50,7 +51,6 @@ let ee_group = SUITE["ee"] = BenchmarkGroup(["ee"])
     for alg in [JetAlgorithm.Durham]
         ee_group["$alg"] = @benchmarkable jet_reconstruct_harness($ee_events;
                                                                   algorithm = $alg,
-                                                                  strategy = $RecoStrategy.Best,
                                                                   ptmin = 5.0) evals=1 samples=32
     end
     for alg in [JetAlgorithm.EEKt]
@@ -60,11 +60,27 @@ let ee_group = SUITE["ee"] = BenchmarkGroup(["ee"])
             for power in [-1.0, 0.0, 1.0]
                 distance_group["$power"] = @benchmarkable jet_reconstruct_harness($ee_events;
                                                                                   algorithm = $alg,
-                                                                                  strategy = $RecoStrategy.Best,
                                                                                   distance = $distance,
                                                                                   power = $power,
                                                                                   ptmin = 5.0) evals=1 samples=32
             end
         end
     end
+    for alg in [JetAlgorithm.Valencia]
+        alg_group = SUITE["ee"]["$alg"] = BenchmarkGroup(["$alg"])
+        for distance in [1.0]
+            distance_group = alg_group["$distance"] = BenchmarkGroup(["$distance"])
+            for β in [1.2]
+                power_group = distance_group["$β"] = BenchmarkGroup(["$β"])
+                for γ in [1.2]
+                    power_group["$γ"] = @benchmarkable jet_reconstruct_harness($ee_events;
+                                                                               algorithm = $alg,
+                                                                               distance = $distance,
+                                                                               power = $β,
+                                                                               γ = $γ,
+                                                                               ptmin = 5.0) evals=1 samples=32
+                end
+            end
+        end
+    end
 end

docs/src/examples.md

@@ -26,6 +26,8 @@ julia --project jetreco.jl --algorithm=AntiKt ../test/data/events.pp13TeV.hepmc3
 ...
 julia --project jetreco.jl --algorithm=Durham ../test/data/events.eeH.hepmc3.zst
 ...
+julia --project jetreco.jl --algorithm=Valencia --beta=1.2 --gamma=0.8 ../test/data/events.eeH.hepmc3.zst
+...
 julia --project jetreco.jl --maxevents=10 --strategy=N2Plain --algorithm=Kt --exclusive-njets=3 ../test/data/events.pp13TeV.hepmc3.zst
 ...
 ```

docs/src/index.md

@@ -13,7 +13,9 @@ Julia.
 
 The algorithms include anti-``{k}_\text{T}``, Cambridge/Aachen, inclusive
 ``k_\text{T}``, generalised ``k_\text{T}`` for ``pp`` events; and the Durham
-algorithm and generalised ``k_\text{T}`` for ``e^+e^-``.
+algorithm and generalised ``k_\text{T}`` for ``e^+e^-``. The Valencia algorithm
+for ``e^+e^-``, based on [1404.4294](https://arxiv.org/abs/1404.4294) is also
+supported.
 
 ## Reconstruction Interface
 
@@ -36,6 +38,10 @@ 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.
 
+For the Valencia algorithm, as well as `R`, the β (equivevent to the existing
+power p of the algorithm) and γ (angular exponent parameter used in the beam
+distance) are needed.
+
 The `R` value determines the cone size; in the case of the Durham algorithm the
 `R` value is ignored.
 
@@ -56,6 +62,7 @@ Each known algorithm is referenced using a `JetAlgorithm` scoped enum value.
 | generalised ``k_\text{T}`` | `JetAlgorithm.GenKt` | For $pp$, value of `p` must also be specified |
 | ``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 |
+| Valencia | `JetAlgorithm.Valencia` | For ``e^+e^-``, values of `p` (β) and `γ` must be specified |
 
 ### Strategy
 
@@ -112,27 +119,39 @@ undergoes a *beam merge* step).
 
 ## References
 
-
-The current recommended reference for JetReconstruction.jl is:
+Although it has been developed further since the CHEP2023 conference, the CHEP
+conference proceedings,
+[10.1051/epjconf/202429505017](https://doi.org/10.1051/epjconf/202429505017),
+should be cited if you use this package:
 
 ```bibtex
 @article{refId0,
-  author = {{Stewart, Graeme Andrew} and {Ganguly, Sanmay} and {Ghosh, Sattwamo} and {Gras, Philippe} and {Krasnopolski, Atell}},
-  title = {Fast Jet Finding in Julia},
-  DOI= "10.1051/epjconf/202533701067",
-  url= "https://doi.org/10.1051/epjconf/202533701067",
-  journal = {EPJ Web Conf.},
-  year = 2025,
-  volume = 337,
-  pages = "01067",
+    author = {{Stewart, Graeme Andrew} and {Gras, Philippe} and {Hegner, Benedikt} and {Krasnopolski, Atell}},
+    doi = {10.1051/epjconf/202429505017},
+    journal = {EPJ Web of Conf.},
+    pages = {05017},
+    title = {Polyglot Jet Finding},
+    url = {https://doi.org/10.1051/epjconf/202429505017},
+    volume = 295,
+    year = 2024,
+    eprint={2309.17309},
+    archivePrefix={arXiv},
+    primaryClass={hep-ex}
 }
 ```
 
-Also available as [arXiv:2503.08146](https://arxiv.org/abs/2503.08146).
-
-### Other Articles
+The original paper on [arXiv](https://arxiv.org/abs/2309.17309) is:
 
-- CHEP2023, *Polyglot Jet Finding*: [arXiv:2309.17309](https://arxiv.org/abs/2309.17309), [10.1051/epjconf/202429505017](https://doi.org/10.1051/epjconf)
+```bibtex
+@misc{stewart2023polyglot,
+      title={Polyglot Jet Finding}, 
+      author={Graeme Andrew Stewart and Philippe Gras and Benedikt Hegner and Atell Krasnopolski},
+      year={2023},
+      eprint={2309.17309},
+      archivePrefix={arXiv},
+      primaryClass={hep-ex}
+}
+```
 
 ## Community
 

docs/src/strategy.md

@@ -14,5 +14,6 @@ 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.
 
-For ``e^+e^-`` algorithms particle densities are low, so the only implementation
-is of the same type as `N2Plain`.
+For ``e^+e^-`` algorithms particle densities are low, so the only
+implementation for these algorithms is effectively of the same type as
+`N2Plain`.

examples/benchmark.sh

@@ -11,3 +11,6 @@ julia --project instrumented-jetreco.jl --algorithm=AntiKt -R 0.4 ../test/data/e
 
 echo "ee H Durham"
 julia --project instrumented-jetreco.jl --algorithm=Durham ../test/data/events.eeH.hepmc3.zst -m $trials
+
+echo "ee Valencia"
+julia --project ./instrumented-jetreco.jl --algorithm=Valencia --gamma=1.2 --power=1.2 -R 0.8 ../test/data/events.eeH.hepmc3.zst -m $trials

examples/instrumented-jetreco.jl

@@ -31,14 +31,15 @@ Profile the jet reconstruction code using the `@profile` macro and generate a
 flamegraph which is saved to the `profile/profile_subdir` directory.
 """
 function profile_code(events::Vector{Vector{T}}, profile, nsamples; R = 0.4, p = -1,
+                      γ::Union{Real, Nothing} = nothing,
                       algorithm::JetAlgorithm.Algorithm = JetAlgorithm.AntiKt,
                       strategy = RecoStrategy.N2Tiled,
                       recombine = RecombinationMethods[RecombinationScheme.ESchemeRaw]) where {T <:
                                                                                                JetReconstruction.FourMomentum}
     Profile.init(n = 5 * 10^6, delay = 0.00001)
     function profile_events(events)
         for evt in events
-            jet_reconstruct(evt; R = R, p = p, algorithm = algorithm,
+            jet_reconstruct(evt; R = R, p = p, γ = γ, algorithm = algorithm,
                             strategy = strategy, recombine...)
         end
     end
@@ -58,6 +59,7 @@ end
 """
     allocation_stats(events::Vector{Vector{T}}; algorithm::JetAlgorithm.Algorithm,
                           distance::Real = 0.4, p::Union{Real, Nothing} = nothing,
+                          γ::Union{Real, Nothing} = nothing,
                           strategy::RecoStrategy.Strategy,
                           recombine = RecombinationMethods[RecombinationScheme.ESchemeRaw],
                           ptmin::Real = 5.0) where {T <: JetReconstruction.FourMomentum}
@@ -66,12 +68,13 @@ output.
 """
 function allocation_stats(events::Vector{Vector{T}}; algorithm::JetAlgorithm.Algorithm,
                           distance::Real = 0.4, p::Union{Real, Nothing} = nothing,
+                          γ::Union{Real, Nothing} = nothing,
                           strategy::RecoStrategy.Strategy,
                           recombine = RecombinationMethods[RecombinationScheme.ESchemeRaw],
                           ptmin::Real = 5.0) where {T <: JetReconstruction.FourMomentum}
     println("Memory allocation statistics:")
     @timev for event in events
-        _ = inclusive_jets(jet_reconstruct(event; R = distance, p = p,
+        _ = inclusive_jets(jet_reconstruct(event; R = distance, p = p, γ = γ,
                                            algorithm = algorithm,
                                            strategy = strategy, recombine...),
                            ptmin = ptmin)
@@ -84,6 +87,7 @@ end
                             algorithm::JetAlgorithm.Algorithm,
                             distance::Real = 0.4,
                             p::Union{Real, Nothing} = nothing,
+                            γ::Union{Real, Nothing} = nothing,
                             ptmin::Real = 5.0,
                             dcut = nothing,
                             njets = nothing,
@@ -106,6 +110,7 @@ function benchmark_jet_reco(events::Vector{Vector{T}};
                             algorithm::JetAlgorithm.Algorithm,
                             distance::Real = 0.4,
                             p::Union{Real, Nothing} = nothing,
+                            γ::Union{Real, Nothing} = nothing,
                             strategy::RecoStrategy.Strategy,
                             recombine = RecombinationMethods[RecombinationScheme.ESchemeRaw],
                             ptmin::Real = 5.0,
@@ -136,7 +141,7 @@ function benchmark_jet_reco(events::Vector{Vector{T}};
         gcoff && GC.enable(false)
         t_start = time_ns()
         for (ievt, event) in enumerate(events)
-            cs = jet_reconstruct(event; R = distance, p = p, algorithm = algorithm,
+            cs = jet_reconstruct(event; R = distance, p = p, γ = γ, algorithm = algorithm,
                                  strategy = strategy, recombine...)
             if !isnothing(njets)
                 selectedjets = exclusive_jets(cs; njets = njets)
@@ -245,8 +250,12 @@ function parse_command_line(args)
         help = """Algorithm to use for jet reconstruction: $(join(JetReconstruction.AllJetRecoAlgorithms, ", "))"""
         arg_type = JetAlgorithm.Algorithm
 
-        "--power", "-p"
-        help = """Power value for jet reconstruction"""
+        "--power", "--beta", "-p"
+        help = """Power value for jet reconstruction (a.k.a. β for Valencia algorithm)"""
+        arg_type = Float64
+
+        "--gamma"
+        help = """Gamma, γ, only for Valencia algorithm)"""
         arg_type = Float64
 
         "--strategy", "-S"
@@ -328,18 +337,18 @@ function main()
     # Major switch between modes of running 
     if args[:alloc]
         allocation_stats(events; distance = args[:distance],
-                         p = args[:power], algorithm = args[:algorithm],
+                         p = args[:power], γ = args[:gamma], algorithm = args[:algorithm],
                          strategy = args[:strategy],
                          recombine = JetReconstruction.RecombinationMethods[args[:recombine]],
                          ptmin = args[:ptmin])
     elseif !isnothing(args[:profile])
         profile_code(events, args[:profile], args[:nsamples];
-                     R = args[:distance], p = args[:power],
+                     R = args[:distance], p = args[:power], γ = args[:gamma],
                      algorithm = args[:algorithm], strategy = args[:strategy],
                      recombine = JetReconstruction.RecombinationMethods[args[:recombine]])
     else
         benchmark_jet_reco(events, distance = args[:distance], algorithm = args[:algorithm],
-                           p = args[:power],
+                           p = args[:power], γ = args[:gamma],
                            strategy = args[:strategy],
                            recombine = JetReconstruction.RecombinationMethods[args[:recombine]],
                            ptmin = args[:ptmin], dcut = args[:exclusive_dcut],