Pre-release cleanup

Write a proper README file that describes how to use the package
(real Julia documentation to come, but it's a start)

Rename main algorithms to be more consistent, viz.
sequential_jet_reconstruct and tiled_jet_reconstruct

Also remove all of the pre-cooked versions with specific values of p
(like cambridge_aachen_algo) as these add little value and would need to
exist for all algorithms.

Rename the chep.jl script to jetreco.jl and put it into an examples
directory

Project.toml

@@ -8,7 +8,6 @@ Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
 ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 FlameGraphs = "08572546-2f56-4bcf-ba4e-bab62c3a3f89"
-HepMC3_jll = "b85c3e40-22db-5268-bacb-02bd65cb4e01"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
@@ -19,8 +18,8 @@ ProfileSVG = "132c30aa-f267-4189-9183-c8a63c7e05e6"
 StatProfilerHTML = "a8a75453-ed82-57c9-9e16-4cd1196ecbf5"
 
 [compat]
-julia = "1.8"
 LorentzVectorHEP = "0.1.6"
+julia = "1.8"
 
 [extras]
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"

README.md

@@ -1,21 +1,92 @@
-# JetReconstruction
+# JetReconstruction.jl
 
-[![Build Status](https://github.com/gojakuch/JetReconstruction.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/gojakuch/JetReconstruction.jl/actions/workflows/CI.yml?query=branch%3Amain)
+[![Build Status](https://github.com/JuliaHEP/JetReconstruction.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/JuliaHEP/JetReconstruction.jl/actions/workflows/CI.yml?query=branch%3Amain)
+
+## This package implements sequential Jet Reconstruction (clustering)
+
+### Algorithms
+
+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`.
+
+### Interface
+
+To use the package one should call the appropriate API interface of the desired algorithm
 
-This code is not a registered Julia package yet. If you want to use it, the valid option would be to change the `main.jl` file how you want and run it. That file should either start with
 ```julia
-using Revise; import Pkg # you need to have Revise installed
-Pkg.activate(".")
-using JetReconstruction
+# N2Plain
+plain_jet_reconstruct(particles::Vector{T}; p = -1, R = 1.0, recombine = +, ptmin = 0.0)
 ```
-or with
+
 ```julia
-include("src/JetReconstruction.jl")
-using .JetReconstruction
+# N2Tiled
+tiled_jet_reconstruct(particles::Vector{T}; p = -1, R = 1.0, recombine = +, ptmin = 0.0)
 ```
-There is some documentation provided for functions and submodules. Once everything is working properly and efficiently, this `README.md` will contain more details on usage and the module might get registered.
 
-## Plotting
+- `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
+  - `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
+
+Both of these functions return tuple of `(jets, seq)`, where these are a vector of `JetReconstruction.PseudoJet` objects and cluster sequencing information, respectively.
+
+Note that the `N2Plain` algorithm is usually faster at low densities, $\lessapprox 50$ input particles, otherwise the tiled algorithm is faster.
+
+### Example
+
+See the `examples/jetreco.jl` script for a full example of how to call jet reconstruction.
+
+```julia
+julia --project=. examples/jetreco.jl --maxevents=100 --nsamples=1 --strategy=N2Plain test/data/events.hepmc3
+...
+julia --project=. examples/jetreco.jl --maxevents=100 --nsamples=1 --strategy=N2Tiled test/data/events.hepmc3
+...
+```
+
+The example also shows how to use `JetReconstruction.HepMC3` to read HepMC3 ASCII files (via the `read_final_state_particles()` wrapper).
+
+### Plotting
+
+**TO BE FIXED**
+
 ![illustration](img/illustration.jpeg)
 
 To visualise the clustered jets as a 3d bar plot (see illustration above) we now use `Makie.jl`. See the `jetsplot` function and its documentation for more. 
+
+## Reference
+
+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
+@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}
+}
+```
+
+## Authors and Copyright
+
+Code in this package is authored by:
+
+- Atell Krasnopolski <[email protected]>
+- Graeme A Stewart <[email protected]>
+- Philippe Gras <[email protected]>
+
+and is Copyright 2022-2023 The Authors, CERN.
+
+The code is under the MIT License.

chep.jl → examples/jetreco.jl

@@ -71,9 +71,9 @@ function jet_process(
 
 	# Strategy
 	if (strategy == N2Plain)
-		jet_reconstruction = sequential_jet_reconstruct
+		jet_reconstruction = plain_jet_reconstruct
 	elseif (strategy == N2Tiled || stragegy == Best)
-		jet_reconstruction = tiled_jet_reconstruct_ll
+		jet_reconstruction = tiled_jet_reconstruct
 	else
 		throw(ErrorException("Strategy not yet implemented"))
 	end
@@ -270,8 +270,6 @@ main() = begin
 	global_logger(logger)
 	events::Vector{Vector{PseudoJet}} =
 		read_final_state_particles(args[:file], maxevents = args[:maxevents], skipevents = args[:skip])
-	# events::Vector{Vector{LorentzVector}} =
-	# 	read_final_state_particles_lv(args[:file], maxevents = args[:maxevents], skipevents = args[:skip])
 	jet_process(events, ptmin = args[:ptmin], distance = args[:distance], 
         power = args[:power], strategy = args[:strategy],
 		nsamples = args[:nsamples], gcoff = args[:gcoff], profile = args[:profile],

src/JetReconstruction.jl

@@ -25,24 +25,22 @@ energy(p::LorentzVectorCyl) = LorentzVectorHEP.energy(p)
 # Philipp's pseudojet
 include("Pseudojet.jl")
 export PseudoJet
-## As this is an internal EDM class, we don't export anything
 
 # Simple HepMC3 reader
 include("HepMC3.jl")
-export HepMC3
 
 ## N2Plain algorithm
 # Algorithmic part for simple sequential implementation
-include("Algo.jl")
-export sequential_jet_reconstruct, kt_algo, anti_kt_algo, anti_kt_algo_alt, cambridge_aachen_algo
+include("PlainAlgo.jl")
+export plain_jet_reconstruct
 
 ## Tiled algorithms
 # Common pieces
 include("TiledAlgoUtils.jl")
 
 # Algorithmic part, tiled reconstruction strategy with linked list jet objects
 include("TiledAlgoLL.jl")
-export tiled_jet_reconstruct_ll
+export tiled_jet_reconstruct
 
 # jet serialisation (saving to file)
 include("Serialize.jl")

src/Algo.jl → src/PlainAlgo.jl

@@ -78,28 +78,34 @@ end
 
 """
 This is the N2Plain jet reconstruction algorithm interface, called with an arbitrary array
-of objects, which supports the methods pt2(), phi(), rapidity() for each element.
+of particles, which supports suitable 4-vector methods, viz.
+ - pt2(), phi(), rapidity(), px(), py(), pz(), energy()
+for each element.
+
+Examples of suitable types are JetReconstruction.PseudoJet and LorentzVectorHEP.
+N.B. these methods need to exist in the namespace of this package, i.e. JetReconstruction.pt2(::T),
+which is already done for the two types above.
 """
-function sequential_jet_reconstruct(objects::AbstractArray{T}; p = -1, R = 1.0, recombine = +, ptmin = 0.0) where T
+function plain_jet_reconstruct(particles::Vector{T}; p = -1, R = 1.0, recombine = +, ptmin = 0.0) where T
     # Integer p if possible
     p = (round(p) == p) ? Int(p) : p 
 
     # We make sure these arrays are type stable - have seen issues where, depending on the values
     # returned by the methods, they can become unstable and performance degrades
-    kt2_array::Vector{Float64} = pt2.(objects) .^ p
-    phi_array::Vector{Float64} = phi.(objects)
-    rapidity_array::Vector{Float64} = rapidity.(objects)
+    kt2_array::Vector{Float64} = pt2.(particles) .^ p
+    phi_array::Vector{Float64} = phi.(particles)
+    rapidity_array::Vector{Float64} = rapidity.(particles)
 
-    objects_array = copy(objects)
+    objects_array = copy(particles)
 
     # Now call the actual reconstruction method, tuned for our internal EDM
-    sequential_jet_reconstruct(objects_array=objects_array, kt2_array=kt2_array, phi_array=phi_array, 
+    plain_jet_reconstruct(objects_array=objects_array, kt2_array=kt2_array, phi_array=phi_array, 
         rapidity_array=rapidity_array, p=p, R=R, recombine=recombine, ptmin=ptmin)
 end
 
 
 
-function sequential_jet_reconstruct(;objects_array::AbstractArray{J}, kt2_array::Vector{F}, 
+function plain_jet_reconstruct(;objects_array::Vector{J}, kt2_array::Vector{F}, 
         phi_array::Vector{F}, rapidity_array::Vector{F}, p = -1, R = 1.0, recombine = +, ptmin = 0.0) where {J, F<:AbstractFloat}
     # Bounds
     N::Int = length(objects_array)
@@ -202,42 +208,3 @@ function sequential_jet_reconstruct(;objects_array::AbstractArray{J}, kt2_array:
 
     jets, sequences
 end
-
-"""
-`anti_kt_algo(objects; R=1, recombine=(x, y)->(x + y)) -> Vector, Vector{Vector{Int}}`
-
-Runs the anti-kt jet reconstruction algorithm. `objects` can be any collection of *unique* elements.
-
-Returns:
-    `jets` - a vector of jets. Each jet is of the same type as elements in `objects`.
-    `sequences` - a vector of vectors of indices in `objects`. For all `i`, `sequences[i]` gives a sequence of indices of objects that have been combined into the i-th jet (`jets[i]`).
-"""
-function anti_kt_algo(objects; R = 1.0, recombine = +)
-    sequential_jet_reconstruct(objects, p = -1, R = R, recombine = recombine)
-end
-
-"""
-`kt_algo(objects; R=1, recombine=(x, y)->(x + y)) -> Vector, Vector{Vector{Int}}`
-
-Runs the kt jet reconstruction algorithm. `objects` can be any collection of *unique* elements.
-
-Returns:
-    `jets` - a vector of jets. Each jet is of the same type as elements in `objects`.
-    `sequences` - a vector of vectors of indices in `objects`. For all `i`, `sequences[i]` gives a sequence of indices of objects that have been combined into the i-th jet (`jets[i]`).
-"""
-function kt_algo(objects; R = 1.0, recombine = +)
-    sequential_jet_reconstruct(objects, p = 1, R = R, recombine = recombine)
-end
-
-"""
-`cambridge_aachen_algo(objects; R=1, recombine=(x, y)->(x + y)) -> Vector, Vector{Vector{Int}}`
-
-Runs the Cambridge/Aachen jet reconstruction algorithm. `objects` can be any collection of *unique* elements.
-
-Returns:
-    `jets` - a vector of jets. Each jet is of the same type as elements in `objects`.
-    `sequences` - a vector of vectors of indices in `objects`. For all `i`, `sequences[i]` gives a sequence of indices of objects that have been combined into the i-th jet (`jets[i]`).
-"""
-function cambridge_aachen_algo(objects; R = 1.0, recombine = +)
-    sequential_jet_reconstruct(objects, p = 0, R = R, recombine = recombine)
-end

src/TiledAlgoLL.jl

@@ -306,21 +306,21 @@ be used.
 If a non-standard recombination is used, it must be defined for
 JetReconstruction.PseudoJet, as this struct is used internally.
 """
-function tiled_jet_reconstruct_ll(particles::Vector{T}; p = -1, R = 1.0, recombine = +, ptmin = 0.0) where {T}
+function tiled_jet_reconstruct(particles::Vector{T}; p = -1, R = 1.0, recombine = +, ptmin = 0.0) where {T}
     # Here we need to populate the vector of PseudoJets that are the internal
     # EDM for the main algorithm, then we call the reconstruction
     pseudojets = Vector{PseudoJet}(undef, length(particles))
     for (i, particle) in enumerate(particles)
         pseudojets[i] = PseudoJet(px(particle), py(particle),
             pz(particle), energy(particle))
     end
-    tiled_jet_reconstruct_ll(pseudojets, p = p, R = R, recombine = recombine, ptmin = ptmin)
+    tiled_jet_reconstruct(pseudojets, p = p, R = R, recombine = recombine, ptmin = ptmin)
 end
 
 """
 Main jet reconstruction algorithm, using PseudoJet objects
 """
-function tiled_jet_reconstruct_ll(particles::Vector{PseudoJet}; p = -1, R = 1.0, recombine = +, ptmin = 0.0)
+function tiled_jet_reconstruct(particles::Vector{PseudoJet}; p = -1, R = 1.0, recombine = +, ptmin = 0.0)
     # Bounds
     N::Int = length(particles)
     # @debug "Initial particles: $(N)"

test/runtests.jl

@@ -41,10 +41,10 @@ function do_jet_test(strategy::JetRecoStrategy, fastjet_jets;
 
 	# Strategy
 	if (strategy == N2Plain)
-		jet_reconstruction = sequential_jet_reconstruct
+		jet_reconstruction = plain_jet_reconstruct
         strategy_name = "N2Plain"
     elseif (strategy == N2Tiled)
-		jet_reconstruction = tiled_jet_reconstruct_ll
+		jet_reconstruction = tiled_jet_reconstruct
         strategy_name = "N2Tiled"
 	else
 		throw(ErrorException("Strategy not yet implemented"))