Pre release cleanup

.gitignore

@@ -2,7 +2,6 @@
 /jetsavetest*.dat
 /demo.jl
 .vscode
-*.json
 /notebooks/*
 /statprof/*
 /debug/*

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,6 +18,7 @@ ProfileSVG = "132c30aa-f267-4189-9183-c8a63c7e05e6"
 StatProfilerHTML = "a8a75453-ed82-57c9-9e16-4cd1196ecbf5"
 
 [compat]
+LorentzVectorHEP = "0.1.6"
 julia = "1.8"
 
 [extras]

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

@@ -14,6 +14,7 @@ using StatProfilerHTML
 using Logging
 using JSON
 
+using LorentzVectorHEP
 using JetReconstruction
 
 function profile_code(jet_reconstruction, events, niters)
@@ -70,23 +71,18 @@ function jet_process(
 
  # Strategy
  if (strategy == N2Plain)
- jet_reconstruction = sequential_jet_reconstruct
- elseif (strategy == N2TiledSoAGlobal)
- jet_reconstruction = tiled_jet_reconstruct_soa_global
- elseif (strategy == N2TiledSoATile)
- jet_reconstruction = tiled_jet_reconstruct_soa_tile
- elseif (strategy == N2Tiled)
- jet_reconstruction = tiled_jet_reconstruct_ll
+ jet_reconstruction = plain_jet_reconstruct
+ elseif (strategy == N2Tiled || stragegy == Best)
+ jet_reconstruction = tiled_jet_reconstruct
  else
  throw(ErrorException("Strategy not yet implemented"))
  end
 
- # The N2Tiled algorithm uses PseudoJets so pass these directly
- if strategy == N2Tiled
+ # Build internal EDM structures for timing measurements, if needed
+ # For N2Tiled and N2Plain this is unnecessary as both these reconstruction
+ # methods can process PseudoJets directly
+ if (strategy == N2Tiled) || (strategy == N2Plain)
  event_vector = events
- else
- # The other algorithms swallow 4-vectors instead
- event_vector = pseudojets2vectors(events)
  end
 
  # If we are dumping the results, setup the JSON structure
@@ -121,11 +117,12 @@ function jet_process(
  GC.gc()
  cummulative_time = 0.0
  cummulative_time2 = 0.0
+ lowest_time = typemax(Float64)
  for irun ∈ 1:nsamples
  gcoff && GC.enable(false)
  t_start = time_ns()
  for (ievt, event) in enumerate(event_vector)
- finaljets, _ = jet_reconstruction(event, R = distance, p = power)
+ finaljets, _ = jet_reconstruction(event, R = distance, p = power, ptmin=ptmin)
  fj = final_jets(finaljets, ptmin)
  # Only print the jet content once
  if irun == 1
@@ -149,6 +146,7 @@ function jet_process(
  end
  cummulative_time += dt_μs
  cummulative_time2 += dt_μs^2
+ lowest_time = dt_μs < lowest_time ? dt_μs : lowest_time
  end
 
  mean = cummulative_time / nsamples
@@ -160,8 +158,17 @@ function jet_process(
  end
  mean /= length(events)
  sigma /= length(events)
+ lowest_time /= length(events)
+ # Why also record the lowest time? 
+ # 
+ # The argument is that on a "busy" machine, the run time of an application is
+ # always TrueRunTime+Overheads, where Overheads is a nuisance parameter that
+ # adds jitter, depending on the other things the machine is doing. Therefore
+ # the minimum value is (a) more stable and (b) reflects better the intrinsic
+ # code performance.
  println("Processed $(length(events)) events $(nsamples) times")
- println("Time per event $(mean) ± $(sigma) μs")
+ println("Average time per event $(mean) ± $(sigma) μs")
+ println("Lowest time per event $lowest_time μs")
 
  if !isnothing(dump)
  open(dump, "w") do io
@@ -246,10 +253,6 @@ function ArgParse.parse_item(::Type{JetRecoStrategy}, x::AbstractString)
  return JetRecoStrategy(1)
  elseif (x == "N2Tiled")
  return JetRecoStrategy(2)
- elseif (x == "N2TiledSoAGlobal")
- return JetRecoStrategy(3)
- elseif (x == "N2TiledSoATile")
- return JetRecoStrategy(4)
  else
  throw(ErrorException("Invalid value for strategy: $(x)"))
  end
@@ -274,12 +277,9 @@ main() = begin
  nothing
 end
 
-# The issue is that running through the debugger in VS Code actually has
+# Running through the debugger in VS Code actually has
 # ARGS[0] = "/some/path/.vscode/extensions/julialang.language-julia-1.47.2/scripts/debugger/run_debugger.jl",
 # so then the program does nothing at all if it only tests for abspath(PROGRAM_FILE) == @__FILE__
-# In addition, deep debugging with Infiltrator needs to start in an interactive session
-#
-# (Really, one starts to wonder if main() should be executed unconditionally!)
-if (abspath(PROGRAM_FILE) == @__FILE__) || (basename(PROGRAM_FILE) == "run_debugger.jl" || isinteractive())
- main()
-end
+# In addition, deep debugging with Infiltrator needs to start in an interactive session,
+# so execute main() unconditionally
+main()