From 7b1b3df830c1ae352ec4066fb5a89f52a0dd537b Mon Sep 17 00:00:00 2001 From: Graeme A Stewart Date: Wed, 25 Oct 2023 14:38:34 +0200 Subject: [PATCH] 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 | 3 +- README.md | 93 +++++++++++++++++++++++---- anti-kt.pseudo => docs/anti-kt.pseudo | 0 chep.jl => examples/jetreco.jl | 6 +- main.jl | 79 ----------------------- src/JetReconstruction.jl | 8 +-- src/{Algo.jl => PlainAlgo.jl} | 61 ++++-------------- src/TiledAlgoLL.jl | 6 +- test/runtests.jl | 4 +- 9 files changed, 107 insertions(+), 153 deletions(-) rename anti-kt.pseudo => docs/anti-kt.pseudo (100%) rename chep.jl => examples/jetreco.jl (96%) delete mode 100644 main.jl rename src/{Algo.jl => PlainAlgo.jl} (71%) diff --git a/Project.toml b/Project.toml index 661b475..d1cd28e 100644 --- a/Project.toml +++ b/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" diff --git a/README.md b/README.md index 8d9fa2a..e2c46df 100644 --- a/README.md +++ b/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 (, +[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 +- Graeme A Stewart +- Philippe Gras + +and is Copyright 2022-2023 The Authors, CERN. + +The code is under the MIT License. diff --git a/anti-kt.pseudo b/docs/anti-kt.pseudo similarity index 100% rename from anti-kt.pseudo rename to docs/anti-kt.pseudo diff --git a/chep.jl b/examples/jetreco.jl similarity index 96% rename from chep.jl rename to examples/jetreco.jl index a57de14..0f092ad 100644 --- a/chep.jl +++ b/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], diff --git a/main.jl b/main.jl deleted file mode 100644 index 35ad608..0000000 --- a/main.jl +++ /dev/null @@ -1,79 +0,0 @@ -# choose one: -# this -#= -include("src/JetReconstruction.jl") -using .JetReconstruction -=# -# or this -using Revise; import Pkg -Pkg.activate(".") -using JetReconstruction - -using CairoMakie - -## array constructor -makejet(pt, ϕ, y) = [pt*cos(ϕ), pt*sin(ϕ), pt*sinh(y), sqrt((pt*cos(ϕ))^2 + (pt*sin(ϕ))^2 + (pt*sinh(y))^2)] - -## choose a data sample and save it -smalldata = [ - makejet(25, 1.2, 0), - makejet(0.1, 1, 0), - makejet(20, 0, 0) -] -smalldata = [ - makejet(25, 1.2, 0), - makejet(0.1, 0.6, 0), - makejet(20, 0, 0) -] -smalldata = [ - makejet(5, 0.8, 0), - makejet(2, 0, 0.8), -] -smalldata = [ - makejet(30, -0.5, -2.5), - makejet(30, -0.5, -3.5) -] -smalldata = [ - makejet(30, 0.5, -2.5), - makejet(30, 0.5, -3.5), - makejet(30, -0.5, -2.5), - makejet(30, -0.5, -3.5), - makejet(0.1, 0, -3), -] -smalldata = [ - [0.2839991726, -0.4668202293, -49.6978846131, 49.709744138], - [0.1434555273, -0.0312880942, -6.9382351613, 6.9411919273], -] -smalldata = [ - makejet(0.8, 0.8, 0.1), - makejet(0.1, 0.1, 0.8), -] - -smalldata = [ - [-0.8807412236, -1.2331262152, -157.431315651, 157.4393822839], - [-0.0051712611, 0.23815508, -9.7396045662, 9.7435168946], # 2 - [0.0362280943, 0.2694752057, -6.9243427525, 6.9310844534], # 3 - [-0.2206628664, -0.1438198985, -0.6838608666, 0.7460038429], - [1.2716787521, 1.0422298083, -6.1740167274, 6.3907254797], # 5 - [-0.5695590845, -0.3627761836, -58.5430479911, 58.5544811606], - [0.2839991726, -0.4668202293, -49.6978846131, 49.709744138], - [0.6510530003, 1.3970949413, -62.7226079598, 62.7485783532], # 8 - [0.1434555273, -0.0312880942, -6.9382351613, 6.9411919273], - [0.4931562547, 2.1627817414, -14.8865871635, 15.0516040711], # 10 - [0.2396813608, -0.0786236784, -1.9340954697, 1.9554625817], - [0.3355486441, 0.0516402769, -0.8346540063, 0.9118040941], - [-0.7853865645, -0.7810520475, -1.5367790662, 1.8994852039], -] -smalldata = [ - [1.2716787521, 1.0422298083, -6.1740167274, 6.3907254797], - [0.4931562547, 2.1627817414, -14.8865871635, 15.0516040711], - [0.3355486441, 0.0516402769, -0.8346540063, 0.9118040941], -] - -savejets("small.dat", smalldata, format="px py pz E") -## run -smalljets, smallind = anti_kt_algo(smalldata, R=1); img = jetsplot(smalldata, smallind, Module=CairoMakie) - -smalljets, smallind = anti_kt_algo_alt(smalldata, R=1); img = jetsplot(smalldata, smallind, Module=CairoMakie) - -smalljets, smallind = kt_algo(smalldata, R=1); img = jetsplot(smalldata, smallind, Module=CairoMakie) diff --git a/src/JetReconstruction.jl b/src/JetReconstruction.jl index c916bb7..e3e424d 100644 --- a/src/JetReconstruction.jl +++ b/src/JetReconstruction.jl @@ -25,16 +25,14 @@ 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 @@ -42,7 +40,7 @@ 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") diff --git a/src/Algo.jl b/src/PlainAlgo.jl similarity index 71% rename from src/Algo.jl rename to src/PlainAlgo.jl index d189188..33744b4 100644 --- a/src/Algo.jl +++ b/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 diff --git a/src/TiledAlgoLL.jl b/src/TiledAlgoLL.jl index 2dcf379..29871a6 100644 --- a/src/TiledAlgoLL.jl +++ b/src/TiledAlgoLL.jl @@ -306,7 +306,7 @@ 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)) @@ -314,13 +314,13 @@ function tiled_jet_reconstruct_ll(particles::Vector{T}; p = -1, R = 1.0, recombi 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)" diff --git a/test/runtests.jl b/test/runtests.jl index 5bb3846..8ea0bea 100644 --- a/test/runtests.jl +++ b/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"))