diff --git a/.gitignore b/.gitignore index e6c70b7..cd76dd6 100644 --- a/.gitignore +++ b/.gitignore @@ -2,7 +2,6 @@ /jetsavetest*.dat /demo.jl .vscode -*.json /notebooks/* /statprof/* /debug/* diff --git a/Project.toml b/Project.toml index e31697b..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,6 +18,7 @@ ProfileSVG = "132c30aa-f267-4189-9183-c8a63c7e05e6" StatProfilerHTML = "a8a75453-ed82-57c9-9e16-4cd1196ecbf5" [compat] +LorentzVectorHEP = "0.1.6" julia = "1.8" [extras] 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 85% rename from chep.jl rename to examples/jetreco.jl index f35a7bd..0f092ad 100644 --- a/chep.jl +++ b/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() 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/Algo.jl b/src/Algo.jl deleted file mode 100644 index 622f406..0000000 --- a/src/Algo.jl +++ /dev/null @@ -1,334 +0,0 @@ - -Base.@propagate_inbounds function _dist(i, j, _eta, _phi) - deta = _eta[i] - _eta[j] - dphi = abs(_phi[i] - _phi[j]) - dphi = ifelse(dphi > pi, 2pi - dphi, dphi) - muladd(deta, deta, dphi*dphi) -end - -# d_{ij} distance with i's NN (times R^2) -Base.@propagate_inbounds function _dij(i, _kt2, _nn, _nndist) - j = _nn[i] - d = _nndist[i] - d*min(_kt2[i], _kt2[j]) -end - -# finds new nn for i and checks everyone additionally -Base.@propagate_inbounds function _upd_nn_crosscheck!(i::Int, from::Int, to::Int, _eta, _phi, _R2, _nndist, _nn) - nndist = _R2 - nn = i - @inbounds @simd for j in from:to - Δ2 = _dist(i, j, _eta, _phi) - if Δ2 < nndist - nn = j - nndist = Δ2 - end - if Δ2 < _nndist[j] - _nndist[j] = Δ2 - _nn[j] = i - end - end - _nndist[i] = nndist - _nn[i] = nn; - #nothing -end - -# finds new nn for i -Base.@propagate_inbounds function _upd_nn_nocross!(i::Int, from::Int, to::Int, _eta, _phi, _R2, _nndist, _nn) - nndist = _R2 - nn = i - @inbounds @simd for j in from:(i-1) - Δ2 = _dist(i, j, _eta, _phi) - if Δ2 <= nndist - nn = j - nndist = Δ2 - end - end - @inbounds @simd for j in (i+1):to - Δ2 = _dist(i, j, _eta, _phi) - f = Δ2 <= nndist - nn = ifelse(f, j, nn) - nndist = ifelse(f, Δ2, nndist) - end - _nndist[i] = nndist - _nn[i] = nn; - #nothing -end - -# entire NN update step -Base.@propagate_inbounds function _upd_nn_step!(i, j, k, N, Nn, _kt2, _eta, _phi, _R2, _nndist, _nn, _nndij) - nnk = _nn[k] - if nnk == i || nnk == j - _upd_nn_nocross!(k, 1, N, _eta, _phi, _R2, _nndist, _nn) # update dist and nn - _nndij[k] = _dij(k, _kt2, _nn, _nndist) - nnk = _nn[k] - end - - if j != i && k != i - Δ2 = _dist(i, k, _eta, _phi) - if Δ2 < _nndist[k] - _nndist[k] = Δ2 - nnk = _nn[k] = i - _nndij[k] = _dij(k, _kt2, _nn, _nndist) - end - - cond = Δ2 < _nndist[i] - _nndist[i], _nn[i] = ifelse(cond, (Δ2,k), (_nndist[i],_nn[i])) - end - - nnk == Nn && (_nn[k] = j); - #nothing -end - -function sequential_jet_reconstruct(objects::AbstractArray{T}; p=-1.0, R=1., recombine=+) where T - # bounds - N::Int = length(objects) - - # returned values - jets = T[] # result - sequences = Vector{Int}[] # recombination sequences, WARNING: first index in the sequence is not necessarily the seed - - # params - _R2 = R*R - _p = (round(p) == p) ? Int(p) : p # integer p if possible - - # data - _objects = copy(objects) - _kt2 = (JetReconstruction.pt.(_objects) .^ 2) .^ _p - # _kt2 = 1.0 ./ (JetReconstruction.pt.(_objects) .^ 2) <- Demo code for antikt talks (i.e., _p = -1) - _phi = JetReconstruction.phi.(_objects) - _eta = JetReconstruction.eta.(_objects) - _nn = Vector(1:N) # nearest neighbours - _nndist = fill(float(_R2), N) # distances to the nearest neighbour - _sequences = Vector{Int}[[x] for x in 1:N] - - # initialize _nn - @simd for i in 1:N - _upd_nn_crosscheck!(i, 1, i-1, _eta, _phi, _R2, _nndist, _nn) - end - - # diJ table *_R2 - _nndij = zeros(N) - @inbounds @simd for i in 1:N - _nndij[i] = _dij(i, _kt2, _nn, _nndist) - end - - iteration::Int = 1 - while N != 0 - # findmin - i = 1 - dij_min = _nndij[1] - @inbounds @simd for k in 2:N - cond = _nndij[k] < dij_min - dij_min, i = ifelse(cond, (_nndij[k], k), (dij_min, i)) - end - - j::Int = _nn[i] - - ## Needed for certain tricky debugging situations - # if iteration==1 - # debug_jets(_nn, _nndist, _nndij) - # end - - if i != j - # swap if needed - if j < i - i, j = j, i - end - - # update ith jet, replacing it with the new one - _objects[i] = recombine(_objects[i], _objects[j]) - _phi[i] = JetReconstruction.phi(_objects[i]) - _eta[i] = JetReconstruction.eta(_objects[i]) - _kt2[i] = (JetReconstruction.pt(_objects[i])^2)^_p - - _nndist[i] = _R2 - _nn[i] = i - - @inbounds for x in _sequences[j] # WARNING: first index in the sequence is not necessarily the seed - push!(_sequences[i], x) - end - else # i == j - push!(jets, _objects[i]) - push!(sequences, _sequences[i]) - end - - # copy jet N to j - _objects[j] = _objects[N] - - _phi[j] = _phi[N] - _eta[j] = _eta[N] - _kt2[j] = _kt2[N] - _nndist[j] = _nndist[N] - _nn[j] = _nn[N] - _nndij[j] = _nndij[N] - - _sequences[j] = _sequences[N] - - Nn::Int = N - N -= 1 - iteration += 1 - - # update nearest neighbours step - @inbounds @simd for k in 1:N - _upd_nn_step!(i, j, k, N, Nn, _kt2, _eta, _phi, _R2, _nndist, _nn, _nndij) - end - - _nndij[i] = _dij(i, _kt2, _nn, _nndist) - end - - 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., 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., 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., recombine=+) - sequential_jet_reconstruct(objects, p=0, R=R, recombine=recombine) -end - -## EXPERIMENTAL ZONE -# typically sequential_jet_reconstruct_alt is used when developing an alternative (possibly better) way of reclustureing to keep the previous working version intact - -function sequential_jet_reconstruct_alt(objects::AbstractArray{T}; p=-1, R=1, recombine=+) where T - # bounds - N::Int = length(objects) - - # returned values - jets = T[] # result - sequences = Vector{Int}[] # recombination sequences, WARNING: first index in the sequence is not necessarily the seed - - # params - _R2::Float64 = R*R - _p = (round(p) == p) ? Int(p) : p # integer p if possible - ap = abs(_p); # absolute p - - # data - _objects = copy(objects) - _kt2 = (JetReconstruction.pt.(_objects) .^ 2) .^ _p - _phi = JetReconstruction.phi.(_objects) - _eta = JetReconstruction.eta.(_objects) - _nn = Vector(1:N) # nearest neighbours - _nndist = fill(float(_R2), N) # distances to the nearest neighbour - _sequences = Vector{Int}[[x] for x in 1:N] - - # initialize _nn - for i::Int in 1:N - _upd_nn_crosscheck!(i, 1, i-1, _eta, _phi, _R2, _nndist, _nn) - end - - # diJ table *_R2 - _nndij = zeros(N) - for i::Int in 1:N - _nndij[i] = _dij(i, _kt2, _nn, _nndist) - end - - while N != 0 - # findmin - i::Int = 1 - dij_min = _nndij[1] - for k::Int in 2:N - if _nndij[k] < dij_min - dij_min = _nndij[k] - i = k - end - end - - j::Int = _nn[i] - - if i != j - # swap if needed - if j < i - i, j = j, i - end - - # update ith jet, replacing it with the new one - _objects[i] = recombine(_objects[i], _objects[j]) - _phi[i] = JetReconstruction.phi(_objects[i]) - _eta[i] = JetReconstruction.eta(_objects[i]) - _kt2[i] = (JetReconstruction.pt(_objects[i])^2)^_p - - _nndist[i] = _R2 - _nn[i] = i - - for x in _sequences[j] # WARNING: first index in the sequence is not necessarily the seed - push!(_sequences[i], x) - end - else # i == j - push!(jets, _objects[i]) - push!(sequences, _sequences[i]) - end - - # copy jet N to j - _objects[j] = _objects[N] - - _phi[j] = _phi[N] - _eta[j] = _eta[N] - _kt2[j] = _kt2[N] - _nndist[j] = _nndist[N] - _nn[j] = _nn[N] - _nndij[j] = _nndij[N] - - _sequences[j] = _sequences[N] - - Nn::Int = N - N -= 1 - - # update nearest neighbours step - for k::Int in 1:N - _upd_nn_step!(i, j, k, N, Nn, _kt2, _eta, _phi, _R2, _nndist, _nn, _nndij) - end - - _nndij[i] = _dij(i, _kt2, _nn, _nndist) - end - - jets, sequences -end - -""" -Not for usage. Use `anti_kt_algo` instead. Correctness is not guaranteed. -""" -function anti_kt_algo_alt(objects; p=-1, R=1, recombine=+) - sequential_jet_reconstruct_alt(objects, p=p, R=R, recombine=recombine) -end - -""" -Jet state debugger -""" -function debug_jets(nn, nndist, dijdist) - for ijet ∈ eachindex(nn) - println("$ijet: $(nn[ijet]); $(nndist[ijet]); $(dijdist[ijet])") - end -end diff --git a/src/JetReconstruction.jl b/src/JetReconstruction.jl index 456c63b..e3e424d 100644 --- a/src/JetReconstruction.jl +++ b/src/JetReconstruction.jl @@ -3,38 +3,44 @@ Jet reconstruction (reclustering) in Julia. """ module JetReconstruction -# particle type definition -include("Particle.jl") -export energy, px, py, pz, pt, phi, mass, eta, kt, ϕ, η +using LorentzVectorHEP + +# Import from LorentzVectorHEP methods for those 4-vector types +pt2(p::LorentzVector) = LorentzVectorHEP.pt2(p) +phi(p::LorentzVector) = LorentzVectorHEP.phi(p) +rapidity(p::LorentzVector) = LorentzVectorHEP.rapidity(p) +px(p::LorentzVector) = LorentzVectorHEP.px(p) +py(p::LorentzVector) = LorentzVectorHEP.py(p) +pz(p::LorentzVector) = LorentzVectorHEP.pz(p) +energy(p::LorentzVector) = LorentzVectorHEP.energy(p) + +pt2(p::LorentzVectorCyl) = LorentzVectorHEP.pt2(p) +phi(p::LorentzVectorCyl) = LorentzVectorHEP.phi(p) +rapidity(p::LorentzVectorCyl) = LorentzVectorHEP.rapidity(p) +px(p::LorentzVectorCyl) = LorentzVectorHEP.px(p) +py(p::LorentzVectorCyl) = LorentzVectorHEP.py(p) +pz(p::LorentzVectorCyl) = LorentzVectorHEP.pz(p) +energy(p::LorentzVectorCyl) = LorentzVectorHEP.energy(p) # Philipp's pseudojet include("Pseudojet.jl") -export PseudoJet, rap, phi, pt2 +export PseudoJet # 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 SoA per tile -include("TiledAlgoSoATile.jl") -export tiled_jet_reconstruct_soa_tile - -# Algorithmic part, tiled reconstruction strategy with global SoA -include("TiledAlgoSoAGlobal.jl") -export tiled_jet_reconstruct_soa_global - # 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") @@ -42,7 +48,7 @@ export savejets, loadjets!, loadjets # utility functions, useful for different primary scripts include("Utils.jl") -export read_final_state_particles, pseudojets2vectors, final_jets +export read_final_state_particles, read_final_state_particles_lv, pseudojets2vectors, final_jets # jet visualisation include("JetVis.jl") @@ -53,7 +59,8 @@ include("JSONresults.jl") export FinalJet, FinalJets, JSON3 # Strategy to be used -@enum JetRecoStrategy Best N2Plain N2Tiled N2TiledSoAGlobal N2TiledSoATile -export JetRecoStrategy, Best, N2Plain, N2Tiled, N2TiledSoAGlobal, N2TiledSoATile +## Maybe an enum is not the best idea, use type dispatch instead? +@enum JetRecoStrategy Best N2Plain N2Tiled +export JetRecoStrategy, Best, N2Plain, N2Tiled end diff --git a/src/Particle.jl b/src/Particle.jl deleted file mode 100644 index 1d4d2a8..0000000 --- a/src/Particle.jl +++ /dev/null @@ -1,68 +0,0 @@ -""" -This module defines the basic approach to particle and pseudojet representation. That is, in the form of an indexable (either mutable or immutable) object `obj` where -``` -obj[1] # px -obj[2] # py -obj[3] # pz -obj[4] # energy -``` -Therefore, a simple vector [x, y, z, E] or a [static array](https://github.com/JuliaArrays/StaticArrays.jl) will do fine. - -You can define your own specific data type to represent physical objects. For the anti-kt algorithm to work with a data type, for instance, it only needs to have `pt`, `phi`, `eta`, and `+` operations defined. Be sure, however, to define those functions as extentions to the already existing `JetReconstruction.pt`, `JetReconstruction.phi`, `JetReconstruction.eta` and `Base.:+` respectively: -``` -import JetReconstruction - -JetReconstruction.eta(x::YourType) = # your eta definition -JetReconstruction.phi(x::YourType) = # your phi definition -JetReconstruction.pt(x::YourType) = # your pt definition - -Base.:+(x::YourType, y::YourType) = # your + definition (only used as a recombination function) -``` -Since `+` is only used in recombination, you can leave it undefined, if you use a custom recombination routine. -""" - -""" -""" -@inline energy(p) = p[4] - -@inline px(p) = p[1] - -@inline py(p) = p[2] - -@inline pz(p) = p[3] - -@inline pt2(p) = @fastmath px(p)^2 + py(p)^2 -@inline pt(p) = @fastmath sqrt(pt2(p)) -const kt = pt - -# @inline phi(p) = @fastmath atan(py(p), px(p)) -# Fix to return in range [0, 2π) -phi(p) = begin - if pt2(p) == 0.0 - phi = 0.0 - else - phi = atan(py(p), px(p)) - end - if phi < 0.0 - phi += 2π - end - phi -end -const ϕ = phi - -@inline mass(p) = @fastmath sqrt(max(energy(p)^2 - px(p)^2 - py(p)^2 - pz(p)^2, 0)) - -#@inline pseudorap(p) = asinh(pz(p)/pt(p)) # pseudorapidity - -# Rapidity -eta(p) = begin - _pt2 = pt2(p) - _abspz = abs(pz(p)) - if (energy(p) == _abspz) && (_pt2 == 0.0) - return (-1)^(pz(p) < 0)*(1e5 + _abspz) # a very large value that depends on pz - end - _m2 = max((energy(p) + pz(p))*(energy(p) - pz(p)) - _pt2, 0.0) # mass^2 - E_plus_z = energy(p) + _abspz - return (-1)^(pz(p) > 0) * 0.5*log((_pt2 + _m2)/(E_plus_z^2)) -end -const η = eta diff --git a/src/PlainAlgo.jl b/src/PlainAlgo.jl new file mode 100644 index 0000000..33744b4 --- /dev/null +++ b/src/PlainAlgo.jl @@ -0,0 +1,210 @@ +Base.@propagate_inbounds function dist(i, j, rapidity_array, phi_array) + drapidity = rapidity_array[i] - rapidity_array[j] + dphi = abs(phi_array[i] - phi_array[j]) + dphi = ifelse(dphi > pi, 2pi - dphi, dphi) + muladd(drapidity, drapidity, dphi * dphi) +end + +# d_{ij} distance with i's NN (times R^2) +Base.@propagate_inbounds function dij(i, kt2_array, nn, nndist) + j = nn[i] + d = nndist[i] + d * min(kt2_array[i], kt2_array[j]) +end + +# finds new nn for i and checks everyone additionally +Base.@propagate_inbounds function upd_nn_crosscheck!(i::Int, from::Int, to::Int, rapidity_array, phi_array, R2, nndist, nn) + nndist_min = R2 + nn_min = i + @inbounds @simd for j in from:to + Δ2 = dist(i, j, rapidity_array, phi_array) + if Δ2 < nndist_min + nn_min = j + nndist_min = Δ2 + end + if Δ2 < nndist[j] + nndist[j] = Δ2 + nn[j] = i + end + end + nndist[i] = nndist_min + nn[i] = nn_min +end + +# finds new nn for i +Base.@propagate_inbounds function upd_nn_nocross!(i::Int, from::Int, to::Int, rapidity_array, phi_array, R2, nndist, nn) + nndist_min = R2 + nn_min = i + @inbounds @simd for j in from:(i-1) + Δ2 = dist(i, j, rapidity_array, phi_array) + if Δ2 <= nndist_min + nn_min = j + nndist_min = Δ2 + end + end + @inbounds @simd for j in (i+1):to + Δ2 = dist(i, j, rapidity_array, phi_array) + f = Δ2 <= nndist_min + nn_min = ifelse(f, j, nn_min) + nndist_min = ifelse(f, Δ2, nndist_min) + end + nndist[i] = nndist_min + nn[i] = nn_min +end + +# entire NN update step +Base.@propagate_inbounds function upd_nn_step!(i, j, k, N, Nn, kt2_array, rapidity_array, phi_array, R2, nndist, nn, nndij) + nnk = nn[k] + if nnk == i || nnk == j + upd_nn_nocross!(k, 1, N, rapidity_array, phi_array, R2, nndist, nn) # update dist and nn + nndij[k] = dij(k, kt2_array, nn, nndist) + nnk = nn[k] + end + + if j != i && k != i + Δ2 = dist(i, k, rapidity_array, phi_array) + if Δ2 < nndist[k] + nndist[k] = Δ2 + nnk = nn[k] = i + nndij[k] = dij(k, kt2_array, nn, nndist) + end + + cond = Δ2 < nndist[i] + nndist[i], nn[i] = ifelse(cond, (Δ2, k), (nndist[i], nn[i])) + end + + nnk == Nn && (nn[k] = j) +end + +""" +This is the N2Plain jet reconstruction algorithm interface, called with an arbitrary array +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 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.(particles) .^ p + phi_array::Vector{Float64} = phi.(particles) + rapidity_array::Vector{Float64} = rapidity.(particles) + + objects_array = copy(particles) + + # Now call the actual reconstruction method, tuned for our internal EDM + 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 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) + + # Returned values + jets = PseudoJet[] + sequences = Vector{Int}[] # recombination sequences, WARNING: first index in the sequence is not necessarily the seed + + # Parameters + R2 = R^2 + ptmin2 = ptmin^2 + + # Data + nn = Vector(1:N) # nearest neighbours + nndist = fill(float(R2), N) # distances to the nearest neighbour + sequences = Vector{Int}[[x] for x in 1:N] + + # initialize _nn + @simd for i in 1:N + upd_nn_crosscheck!(i, 1, i - 1, rapidity_array, phi_array, R2, nndist, nn) + end + + # diJ table *_R2 + nndij::Vector{Float64} = zeros(N) + @inbounds @simd for i in 1:N + nndij[i] = dij(i, kt2_array, nn, nndist) + end + + iteration::Int = 1 + while N != 0 + # findmin + i = 1 + dij_min = nndij[1] + @inbounds @simd for k in 2:N + cond = nndij[k] < dij_min + dij_min, i = ifelse(cond, (nndij[k], k), (dij_min, i)) + end + + j::Int = nn[i] + + ## Needed for certain tricky debugging situations + # if iteration==1 + # debug_jets(_nn, _nndist, _nndij) + # end + + if i != j + # swap if needed + if j < i + i, j = j, i + end + + # update ith jet, replacing it with the new one + objects_array[i] = recombine(objects_array[i], objects_array[j]) + phi_array[i] = phi(objects_array[i]) + rapidity_array[i] = rapidity(objects_array[i]) + kt2_array[i] = pt2(objects_array[i]) ^ p + + nndist[i] = R2 + nn[i] = i + + @inbounds for x in sequences[j] # WARNING: first index in the sequence is not necessarily the seed + push!(sequences[i], x) + end + else # i == j + if (pt2(objects_array[i]) >= ptmin2) + # We return PseudoJets, so if we were not passed these then we need to convert (N.B. this is costly!) + if J == PseudoJet + push!(jets, objects_array[i]) + else + push!(jets, PseudoJet(px(objects_array[i]), py(objects_array[i]), pz(objects_array[i]), energy(objects_array[i]))) + end + end + push!(sequences, sequences[i]) + end + + # copy jet N to j + objects_array[j] = objects_array[N] + + phi_array[j] = phi_array[N] + rapidity_array[j] = rapidity_array[N] + kt2_array[j] = kt2_array[N] + nndist[j] = nndist[N] + nn[j] = nn[N] + nndij[j] = nndij[N] + + sequences[j] = sequences[N] + + Nn::Int = N + N -= 1 + iteration += 1 + + # update nearest neighbours step + @inbounds @simd for k in 1:N + upd_nn_step!(i, j, k, N, Nn, kt2_array, rapidity_array, phi_array, R2, nndist, nn, nndij) + end + # @infiltrate + + nndij[i] = dij(i, kt2_array, nn, nndist) + end + + jets, sequences +end diff --git a/src/Pseudojet.jl b/src/Pseudojet.jl index d2a9298..d8a712e 100644 --- a/src/Pseudojet.jl +++ b/src/Pseudojet.jl @@ -1,4 +1,4 @@ -# Copyright (c) 2022, Philippe Gras +# Copyright (c) 2022-2023, Philippe Gras and CERN # # Adapted from PseudoJet class of c++ code from the Fastjet # software (https://fastjet.fr, hep-ph/0512210, arXiv:1111.6097) @@ -6,30 +6,7 @@ # Copyright (c) 2005-2020, Matteo Cacciari, Gavin P. Salam and Gregory # Soyez # -# -#---------------------------------------------------------------------- -# This file is part of AntiKt.jl. -# -# AntiKt.jl is free software; you can redistribute it and/or modify -# it under the terms of the GNU General Public License as published by -# the Free Software Foundation; either version 2 of the License, or -# (at your option) any later version. -# -# The algorithms that underlie FastJet have required considerable -# development. They are described in the original FastJet paper, -# hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use -# FastJet as part of work towards a scientific publication, please -# quote the version you use and include a citation to the manual and -# optionally also to hep-ph/0512210. -# -# AntiKet.jl is distributed in the hope that it will be useful, -# but WITHOUT ANY WARRANTY; without even the implied warranty of -# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -# GNU General Public License for more details. -# -# You should have received a copy of the GNU General Public License -# along with FastJet. If not, see . -#---------------------------------------------------------------------- +# Some implementation is taken from LorentzVectorHEP.jl, (c) Jerry Ling """Interface for composite types that includes fields px, py, py, and E @@ -77,7 +54,7 @@ show(io::IO, jet::PseudoJet) = begin end -set_momentum(j::PseudoJet, px, py, pz, E) = begin +set_momentum!(j::PseudoJet, px, py, pz, E) = begin j.px = px j.py = py j.pz = pz @@ -125,7 +102,7 @@ phi_02pi(p::PseudoJet) = begin return p._phi end -rap(p::PseudoJet) = begin +rapidity(p::PseudoJet) = begin _ensure_valid_rap_phi(p) return p._rap end @@ -138,13 +115,42 @@ pt(p::PseudoJet) = sqrt(p._pt2) "Returns the squared invariant mass" m2(p::PseudoJet) = (p.E + p.pz)*(p.E-p.pz) - p._pt2 +"Returns the magnitude of the momentum, |p|" +mag(p::PseudoJet) = sqrt(muladd(p.px, p.px, p.py^2) + p.pz^2) + +@inline function CosTheta(p::PseudoJet) + fZ = p.pz + ptot = mag(p) + return ifelse(ptot == 0.0, 1.0, fZ / ptot) +end + +"Returns pseudorapidity, η" +function eta(p::PseudoJet) + cosTheta = CosTheta(p) + (cosTheta^2 < 1.0) && return -0.5 * log((1.0 - cosTheta) / (1.0 + cosTheta)) + fZ = p.pz + iszero(fZ) && return 0.0 + # Warning("PseudoRapidity","transverse momentum = 0! return +/- 10e10"); + fZ > 0.0 && return 10e10 + return -10e10 +end +const η = eta + # returns the invariant mass -# (If m2() is negative then -sqrt(-m2()) is returned, as in CLHEP) +# (If m2() is negative then -sqrt(-m2()) is returned, as in CLHEP, SciKitHEP Particle and ROOT) m(p::PseudoJet) = begin x = m2(p) x < 0. ? -sqrt(-x) : sqrt(x) end +# Ensure we have accessors for jet parameters +px(p::PseudoJet) = p.px +py(p::PseudoJet) = p.py +pz(p::PseudoJet) = p.pz +mass(p::PseudoJet) = m(p) +const mass2 = m2 +energy(p::PseudoJet) = p.E + import Base.+; +(j1::PseudoJet, j2::PseudoJet) = begin diff --git a/src/TiledAlgoLL.jl b/src/TiledAlgoLL.jl index ae17696..29871a6 100644 --- a/src/TiledAlgoLL.jl +++ b/src/TiledAlgoLL.jl @@ -5,6 +5,7 @@ using Logging using Accessors +using LoopVectorization # Include struct definitions and basic operations include("TiledAlglLLStructs.jl") @@ -17,7 +18,7 @@ Returns the history and the total event energy. function initial_history(particles) # reserve sufficient space for everything history = Vector{HistoryElement}(undef, length(particles)) - sizehint!(history, 2*length(particles)) + sizehint!(history, 2 * length(particles)) Qtot::Float64 = 0 @@ -38,7 +39,7 @@ between two jets.""" _tj_dist(jetA, jetB) = begin dphi = π - abs(π - abs(jetA.phi - jetB.phi)) deta = jetA.eta - jetB.eta - return dphi*dphi + deta*deta + return dphi * dphi + deta * deta end _tj_diJ(jet) = begin @@ -56,17 +57,17 @@ tile_index(tiling_setup, eta::Float64, phi::Float64) = begin # - eta can be out of range by construction (open end bins) # - phi is protection against bad rounding ieta = clamp(1 + unsafe_trunc(Int, (eta - tiling_setup._tiles_eta_min) / tiling_setup._tile_size_eta), 1, tiling_setup._n_tiles_eta) - iphi = clamp(unsafe_trunc(Int, phi / tiling_setup._tile_size_phi), 0, tiling_setup._n_tiles_phi) + iphi = clamp(unsafe_trunc(Int, phi / tiling_setup._tile_size_phi), 0, tiling_setup._n_tiles_phi) return iphi * tiling_setup._n_tiles_eta + ieta end """Initialise a tiled jet from a PseudoJet (using an index into our ClusterSequence)""" tiledjet_set_jetinfo!(jet::TiledJet, clusterseq::ClusterSequence, jets_index, R2) = begin - @inbounds jet.eta = rap(clusterseq.jets[jets_index]) - @inbounds jet.phi = phi_02pi(clusterseq.jets[jets_index]) - @inbounds jet.kt2 = pt2(clusterseq.jets[jets_index]) > 1.e-300 ? 1. / pt2(clusterseq.jets[jets_index]) : 1.e300 - jet.jets_index = jets_index + @inbounds jet.eta = rapidity(clusterseq.jets[jets_index]) + @inbounds jet.phi = phi_02pi(clusterseq.jets[jets_index]) + @inbounds jet.kt2 = pt2(clusterseq.jets[jets_index]) > 1.e-300 ? 1.0 / pt2(clusterseq.jets[jets_index]) : 1.e300 + jet.jets_index = jets_index # Initialise NN info as well jet.NN_dist = R2 jet.NN = noTiledJet @@ -77,7 +78,9 @@ tiledjet_set_jetinfo!(jet::TiledJet, clusterseq::ClusterSequence, jets_index, R2 # Insert it into the tile's linked list of jets (at the beginning) jet.previous = noTiledJet @inbounds jet.next = clusterseq.tiling.tiles[jet.tile_index] - if isvalid(jet.next) jet.next.previous = jet; end + if isvalid(jet.next) + jet.next.previous = jet + end @inbounds clusterseq.tiling.tiles[jet.tile_index] = jet nothing end @@ -90,7 +93,9 @@ function set_nearest_neighbours!(clusterseq::ClusterSequence, tiledjets::Vector{ isvalid(tile) || continue for jetA in tile for jetB in tile - if jetB == jetA break; end + if jetB == jetA + break + end dist = _tj_dist(jetA, jetB) if (dist < jetA.NN_dist) jetA.NN_dist = dist @@ -133,7 +138,7 @@ function set_nearest_neighbours!(clusterseq::ClusterSequence, tiledjets::Vector{ diJ[i] = _tj_diJ(jetA) # kt distance * R^2 # our compact diJ table will not be in one-to-one corresp. with non-compact jets, # so set up bi-directional correspondence here. - @inbounds NNs[i] = jetA + @inbounds NNs[i] = jetA jetA.dij_posn = i end NNs, diJ @@ -142,10 +147,10 @@ end """Carries out the bookkeeping associated with the step of recombining jet_i and jet_j (assuming a distance dij) and returns the index of the recombined jet, newjet_k.""" -do_ij_recombination_step!(clusterseq::ClusterSequence, jet_i, jet_j, dij) = begin +do_ij_recombination_step!(clusterseq::ClusterSequence, jet_i, jet_j, dij, recombine=+) = begin # Create the new jet by recombining the first two with # the E-scheme - push!(clusterseq.jets, clusterseq.jets[jet_i] + clusterseq.jets[jet_j]) + push!(clusterseq.jets, recombine(clusterseq.jets[jet_i], clusterseq.jets[jet_j])) # Get its index and the history index newjet_k = length(clusterseq.jets) @@ -159,7 +164,7 @@ do_ij_recombination_step!(clusterseq::ClusterSequence, jet_i, jet_j, dij) = begi hist_j = clusterseq.jets[jet_j]._cluster_hist_index add_step_to_history!(clusterseq, minmax(hist_i, hist_j)..., - newjet_k, dij) + newjet_k, dij) newjet_k end @@ -169,14 +174,14 @@ jet_i with the beam (i.e., finalising the jet)""" do_iB_recombination_step!(clusterseq::ClusterSequence, jet_i, diB) = begin # Recombine the jet with the beam add_step_to_history!(clusterseq, clusterseq.jets[jet_i]._cluster_hist_index, BeamJet, - Invalid, diB) + Invalid, diB) end """Add a new jet's history into the recombination sequence""" add_step_to_history!(clusterseq::ClusterSequence, parent1, parent2, jetp_index, dij) = begin max_dij_so_far = max(dij, clusterseq.history[end].max_dij_so_far) push!(clusterseq.history, HistoryElement(parent1, parent2, Invalid, - jetp_index, dij, max_dij_so_far)) + jetp_index, dij, max_dij_so_far)) local_step = length(clusterseq.history) @@ -195,7 +200,10 @@ add_step_to_history!(clusterseq::ClusterSequence, parent1, parent2, jetp_index, clusterseq.history[parent1] = @set hist_elem.child = local_step if parent2 >= 1 - clusterseq.history[parent2].child == Invalid || error("Internal error. Trying to recombine an object that has previsously been recombined. Parent " * string(parent2) * "'s child index " * string(clusterseq.history[parent1].child) * ". Parent jet index: " * string(clusterseq.history[parent2].jetp_index) * ".") + clusterseq.history[parent2].child == Invalid || error( + "Internal error. Trying to recombine an object that has previsously been recombined. Parent " * string(parent2) * "'s child index " * string(clusterseq.history[parent1].child) * ". Parent jet index: " * + string(clusterseq.history[parent2].jetp_index) * ".", + ) hist_elem = clusterseq.history[parent2] clusterseq.history[parent2] = @set hist_elem.child = local_step end @@ -214,8 +222,8 @@ you go along. When a neighbour is added its tagged status is set to true. Returns the updated number of near_tiles.""" function add_untagged_neighbours_to_tile_union(center_index, - tile_union, n_near_tiles, - tiling) + tile_union, n_near_tiles, + tiling) for tile_index in surrounding(center_index, tiling) @inbounds if !tiling.tags[tile_index] n_near_tiles += 1 @@ -237,15 +245,15 @@ Updates tile_union and returns n_near_tiles """ function find_tile_neighbours!(tile_union, jetA, jetB, oldB, tiling) n_near_tiles = add_untagged_neighbours_to_tile_union(jetA.tile_index, - tile_union, 0, tiling) + tile_union, 0, tiling) if isvalid(jetB) if jetB.tile_index != jetA.tile_index n_near_tiles = add_untagged_neighbours_to_tile_union(jetB.tile_index, - tile_union, n_near_tiles, tiling) + tile_union, n_near_tiles, tiling) end if oldB.tile_index != jetA.tile_index && oldB.tile_index != jetB.tile_index n_near_tiles = add_untagged_neighbours_to_tile_union(oldB.tile_index, - tile_union, n_near_tiles, tiling) + tile_union, n_near_tiles, tiling) end end n_near_tiles @@ -267,9 +275,9 @@ end """Return all inclusive jets of a ClusterSequence with pt > ptmin""" -function inclusive_jets(clusterseq::ClusterSequence, ptmin = 0.) - dcut = ptmin*ptmin - jets_local = Vector{PseudoJet}(undef, 0) +function inclusive_jets(clusterseq::ClusterSequence, ptmin = 0.0) + dcut = ptmin * ptmin + jets_local = Vector{LorentzVectorCyl}(undef, 0) # sizehint!(jets_local, length(clusterseq.jets)) # For inclusive jets with a plugin algorithm, we make no # assumptions about anything (relation of dij to momenta, @@ -280,32 +288,55 @@ function inclusive_jets(clusterseq::ClusterSequence, ptmin = 0.) iparent_jet = clusterseq.history[elt.parent1].jetp_index jet = clusterseq.jets[iparent_jet] if pt2(jet) >= dcut - push!(jets_local, jet) + push!(jets_local, LorentzVectorCyl(pt(jet), rapidity(jet), phi(jet), mass(jet))) + # push!(jets_local, jet) end end jets_local end + +""" +Main jet reconstruction algorithm entry point for generic data types + +`particles` must support methods px, py, pz and energy (N.B. these must be in the +JetReconstruction namespace). In particular Cartesian LorentzVector structs can +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(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(pseudojets, p = p, R = R, recombine = recombine, ptmin = ptmin) +end + """ -Main jet reconstruction algorithm +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)" + N::Int = length(particles) + # @debug "Initial particles: $(N)" - # Algorithm parameters - R2::Float64 = R * R - p = (round(p) == p) ? Int(p) : p # integer p if possible + # Algorithm parameters + R2::Float64 = R * R + p = (round(p) == p) ? Int(p) : p # integer p if possible # This will be used quite deep inside loops, but declare it here so that # memory (de)allocation gets done only once - tile_union = Vector{Int}(undef, 3*_n_tile_neighbours) + tile_union = Vector{Int}(undef, 3 * _n_tile_neighbours) # Container for pseudojets, sized for all initial particles, plus all of the # merged jets that can be created during reconstruction jets = PseudoJet[] - sizehint!(jets, N*2) + sizehint!(jets, N * 2) resize!(jets, N) # Copy input data into the jets container @@ -317,7 +348,7 @@ function tiled_jet_reconstruct_ll(particles::Vector{PseudoJet}; p = -1, R = 1.0, # Now get the tiling setup _eta = Vector{Float64}(undef, length(particles)) for ijet in 1:length(particles) - _eta[ijet] = rap(particles[ijet]) + _eta[ijet] = rapidity(particles[ijet]) end tiling = Tiling(setup_tiling(_eta, R)) @@ -342,7 +373,7 @@ function tiled_jet_reconstruct_ll(particles::Vector{PseudoJet}; p = -1, R = 1.0, for iteration in 1:N # Last slot holds the index of the final valid entry in the # compact NNs and dij arrays - ilast = N - (iteration-1) + ilast = N - (iteration - 1) # Search for the lowest value of min_dij_ijet dij_min, ibest = find_lowest(dij, ilast) @inbounds jetA = NNs[ibest] @@ -360,17 +391,17 @@ function tiled_jet_reconstruct_ll(particles::Vector{PseudoJet}; p = -1, R = 1.0, # the new jet that is added as jetB is inserted in a position that # has a future! if jetA.id < jetB.id - jetA, jetB = jetB, jetA; + jetA, jetB = jetB, jetA end # Recombine jetA and jetB and retrieves the new index, nn - nn = do_ij_recombination_step!(clusterseq, jetA.jets_index, jetB.jets_index, dij_min) + nn = do_ij_recombination_step!(clusterseq, jetA.jets_index, jetB.jets_index, dij_min, recombine) tiledjet_remove_from_tiles!(clusterseq.tiling, jetA) oldB = copy(jetB) # take a copy because we will need it... tiledjet_remove_from_tiles!(clusterseq.tiling, jetB) tiledjet_set_jetinfo!(jetB, clusterseq, nn, R2) # cause jetB to become _jets[nn] - # (in addition, registers the jet in the tiling) + # (in addition, registers the jet in the tiling) else # Jet-beam recombination do_iB_recombination_step!(clusterseq, jetA.jets_index, dij_min) @@ -391,7 +422,7 @@ function tiled_jet_reconstruct_ll(particles::Vector{PseudoJet}; p = -1, R = 1.0, # other particles. # Run over all tiles in our union for itile in 1:n_near_tiles - @inbounds tile = tiling.tiles[ @inbounds tile_union[itile]] #TAKES 5μs + @inbounds tile = tiling.tiles[@inbounds tile_union[itile]] #TAKES 5μs @inbounds tiling.tags[tile_union[itile]] = false # reset tag, since we're done with unions isvalid(tile) || continue #Probably not required @@ -415,14 +446,14 @@ function tiled_jet_reconstruct_ll(particles::Vector{PseudoJet}; p = -1, R = 1.0, end end # next jetJ end # next near_tile - dij[jetI.dij_posn] = _tj_diJ(jetI) # update diJ kt-dist + dij[jetI.dij_posn] = _tj_diJ(jetI) # update diJ kt-dist end #jetI.NN == jetA || (jetI.NN == jetB && !isnothing(jetB)) # check whether new jetB is closer than jetI's current NN and # if jetI is closer than jetB's current (evolving) nearest # neighbour. Where relevant update things. if isvalid(jetB) - dist = _tj_dist(jetI,jetB) + dist = _tj_dist(jetI, jetB) if dist < jetI.NN_dist if jetI != jetB jetI.NN_dist = dist diff --git a/src/TiledAlgoSoAGlobal.jl b/src/TiledAlgoSoAGlobal.jl deleted file mode 100644 index 14aceaa..0000000 --- a/src/TiledAlgoSoAGlobal.jl +++ /dev/null @@ -1,578 +0,0 @@ -# Tiled jet reconstruction, linked list data structure approach - -using Logging -using LoopVectorization - -""" -Structure holding the flat jets for a tiled reconstruction -""" -struct FlatJets - # Physics quantities - kt2::Vector{Float64} # p_t^(-2*power) - eta::Vector{Float64} # Rapidity - phi::Vector{Float64} # Phi coordinate - - # Mapping to original jets in (px, py, pz, E) - jet_index::Vector{Int} - - # Tiles and linked list - tile_index::Vector{Int} # My tile index (this is the linear index) - next_jet::Vector{Int} # This is the linked list index of the next jet for this tile (0→end) - prev_jet::Vector{Int} # This is the linked list index of the previous jet for this tile (0→first) - - # Reconstruction parameters - nearest_neighbour::Vector{Int} - nn_distance::Vector{Float64} - dij_distance::Vector{Float64} -end - -# Accessor functions for the jet array -kt2(jets::FlatJets, n::Int) = jets.kt2[n] -eta(jets::FlatJets, n::Int) = jets.eta[n] -phi(jets::FlatJets, n::Int) = jets.phi[n] -jet_index(jets::FlatJets, n::Int) = jets.jet_index[n] -tile_index(jets::FlatJets, n::Int) = jets.tile_index[n] -next_jet(jets::FlatJets, n::Int) = jets.next_jet[n] -prev_jet(jets::FlatJets, n::Int) = jets.prev_jet[n] -nearest_neighbour(jets::FlatJets, n::Int) = jets.nearest_neighbour[n] -nn_distance(jets::FlatJets, n::Int) = jets.nn_distance[n] -dij_distance(jets::FlatJets, n::Int) = jets.dij_distance[n] - -# Setters -set_kt2!(jets::FlatJets, n::Int, v) = jets.kt2[n] = v -set_eta!(jets::FlatJets, n::Int, v) = jets.eta[n] = v -set_phi!(jets::FlatJets, n::Int, v) = jets.phi[n] = v - -set_jet_index!(jets::FlatJets, n::Int, i) = jets.jet_index[n] = i -set_tile_index!(jets::FlatJets, n::Int, i) = jets.tile_index[n] = i -set_next_jet!(jets::FlatJets, n::Int, next::Int) = jets.next_jet[n] = next -set_prev_jet!(jets::FlatJets, n::Int, prev::Int) = jets.prev_jet[n] = prev - -set_nearest_neighbour!(jets::FlatJets, n::Int, v::Int) = jets.nearest_neighbour[n] = v -set_nn_distance!(jets::FlatJets, n::Int, v::Float64) = jets.nn_distance[n] = v -set_dij_distance!(jets::FlatJets, n::Int, v::Float64) = jets.dij_distance[n] = v - -"""Add jet from (px, py, pz, E) values""" -function insert_flatjet!(jets::FlatJets, tiling_setup, p, n::Int, jet_index::Int, newjet::Vector{Float64}) - set_kt2!(jets, n, JetReconstruction.pt2(newjet) ^ p) - set_eta!(jets, n, JetReconstruction.eta(newjet)) - set_phi!(jets, n, JetReconstruction.phi(newjet)) - set_jet_index!(jets, n, jet_index) - tile_η, tile_ϕ = get_tile(tiling_setup, eta(jets, n), phi(jets, n)) - set_tile_index!(jets, n, get_tile_linear_index(tiling_setup, tile_η, tile_ϕ)) - set_next_jet!(jets, n, 0) - set_prev_jet!(jets, n, 0) - set_nearest_neighbour!(jets, n, 0) - set_nn_distance!(jets, n, 1e9) - set_dij_distance!(jets, n, 1e9) -end - -"""Suppress jet at index, copying in the last jet is needed""" -function suppress_flatjet!(jets::FlatJets, n::Int) - # Is the jet we want to get rid of the final jet? In this case the job is trivial - ilast = lastindex(jets.kt2) - tainted_index::Int = 0 - if n != ilast - # Not the last jet - need to shuffle... - set_kt2!(jets, n, kt2(jets, ilast)) - set_eta!(jets, n, eta(jets, ilast)) - set_phi!(jets, n, phi(jets, ilast)) - set_jet_index!(jets, n, jet_index(jets, ilast)) - set_tile_index!(jets, n, tile_index(jets, ilast)) - set_next_jet!(jets, n, next_jet(jets, ilast)) - set_prev_jet!(jets, n, prev_jet(jets, ilast)) - set_nearest_neighbour!(jets, n, nearest_neighbour(jets, ilast)) - set_nn_distance!(jets, n, nn_distance(jets, ilast)) - set_dij_distance!(jets, n, dij_distance(jets, ilast)) - tainted_index = ilast - end - pop!(jets.kt2) - pop!(jets.eta) - pop!(jets.phi) - pop!(jets.jet_index) - pop!(jets.tile_index) - pop!(jets.next_jet) - pop!(jets.prev_jet) - pop!(jets.nearest_neighbour) - pop!(jets.nn_distance) - pop!(jets.dij_distance) - tainted_index -end - -""" -Push a jet onto a tile's linked list -""" -function push_to_tile!(tiles, itile, flatjets::FlatJets, ijet::Int) - if tiles[itile] != 0 - set_prev_jet!(flatjets, tiles[itile], ijet) - set_next_jet!(flatjets, ijet, tiles[itile]) - end - tiles[itile] = ijet -end - -""" -Delete a jet from a tile's linked list -""" -function delete_from_tile!(tiles, itile, flatjets::FlatJets, ijet::Int) - # Are we the first jet for this tile? - if tiles[itile] == ijet - tiles[itile] = next_jet(flatjets, ijet) - if tiles[itile] != 0 - set_prev_jet!(flatjets, tiles[itile], 0) - end - return - end - # If not, then join our prev ⇄ next - if next_jet(flatjets, ijet) != 0 - set_prev_jet!(flatjets, next_jet(flatjets, ijet), prev_jet(flatjets, ijet)) - end - set_next_jet!(flatjets, prev_jet(flatjets, ijet), next_jet(flatjets, ijet)) -end - -""" -Populate the tiles with nearest neighbour information -""" -function populate_tile_nn!(tiles, tiling_setup) - # To help with later iterations, we now find and cache neighbour tile indexes - @inbounds for iη in 1:tiling_setup._n_tiles_eta - @inbounds for iϕ in 1:tiling_setup._n_tiles_phi - # Clamping ensures we don't go beyond the limits of the eta tiling (which do not wrap) - @inbounds for jη in clamp(iη - 1, 1, tiling_setup._n_tiles_eta):clamp(iη + 1, 1, tiling_setup._n_tiles_eta) - δη = jη - iη - @inbounds for jϕ in iϕ-1:iϕ+1 - if (jη == iη && jϕ == iϕ) - continue - end - # Phi tiles wrap around to meet each other - δϕ = jϕ - iϕ # Hold this unwrapped value for rightmost comparison - if (jϕ == 0) - jϕ = tiling_setup._n_tiles_phi - elseif (jϕ == tiling_setup._n_tiles_phi + 1) - jϕ = 1 - end - # Tile is a neighbour - tile_index = tiling_setup._tile_linear_indexes[jη, jϕ] - push!(tiles[iη, iϕ].neighbour_tiles, tile_index) - # Only the tile directly above or to the right are _righttiles - if (((δη == -1) && (δϕ == 0)) || (δϕ == 1)) - push!(tiles[iη, iϕ].right_neighbour_tiles, tile_index) - end - end - end - end - end -end - -""" -Populate the tiles with their first jets and construct the linked lists in the flat jet array - -Note that we *push* new jets into the tile, which avoids walking through the list -""" -function populate_tile_lists!(tiles, flatjets::FlatJets, tiling_setup::TilingDef) - for ijet in eachindex(flatjets.eta) - tile_η, tile_ϕ = get_tile(tiling_setup, eta(flatjets, ijet), phi(flatjets, ijet)) - set_tile_index!(flatjets, ijet, get_tile_linear_index(tiling_setup, tile_η, tile_ϕ)) - if tiles[tile_η, tile_ϕ] != 0 - set_next_jet!(flatjets, ijet, tiles[tile_η, tile_ϕ]) - set_prev_jet!(flatjets, tiles[tile_η, tile_ϕ], ijet) - end - tiles[tile_η, tile_ϕ] = ijet - end -end - - -""" -Check that the state of tiling is consistent -- Every live pseudojet should be in one and only one tile's LL -- Each tile should have the correct start point for its pseudojets -""" -function debug_tiles(tiles, flatjets::FlatJets) - msg = "Testing tile structures\n" - jet_tile_index = similar(flatjets.kt2, Int) - active_jets = length(flatjets.next_jet) - fill!(jet_tile_index, 0) - for itile in eachindex(tiles) - ijet = tiles[itile] - while ijet != 0 - if ijet > active_jets - msg *= "Error: Found out of bounds jet $ijet in tile $itile\n" - continue - end - if jet_tile_index[ijet] != 0 - msg *= "Error: Found jet $ijet already in tile $(jet_tile_index[ijet])\n" - else - jet_tile_index[ijet] = itile - end - ijet = next_jet(flatjets, ijet) - end - end - for ijet in eachindex(jet_tile_index) - if jet_tile_index[ijet] == 0 - msg *= "Error: Jet $ijet is not associated with a tile\n" - end - end - msg -end - -""" -Debug jet state -""" -function debug_jets(flatjets::FlatJets) - for ijet ∈ eachindex(flatjets.dij_distance) - println("$ijet: $(nearest_neighbour(flatjets, ijet)); $(nn_distance(flatjets, ijet)); $(dij_distance(flatjets, ijet))") - end -end - -""" -Complete scan over all tiles to setup the nearest neighbour mappings at the beginning -""" -function find_all_tiled_nearest_neighbours!(tiles, flatjets::FlatJets, tiling_setup::TilingDef, R2) - # Iterate tile by tile... - for itile ∈ eachindex(tiles) - itile_cartesian = get_tile_cartesian_indices(tiling_setup, itile) - - ijet = tiles[itile] - while ijet!= 0 - jjet = next_jet(flatjets, ijet) - while jjet!=0 - nn_dist = geometric_distance(eta(flatjets, ijet), phi(flatjets, ijet), - eta(flatjets, jjet), phi(flatjets, jjet)) - if nn_dist < nn_distance(flatjets, ijet) - set_nn_distance!(flatjets, ijet, nn_dist) - set_nearest_neighbour!(flatjets, ijet, jjet) - end - if nn_dist < nn_distance(flatjets, jjet) - set_nn_distance!(flatjets, jjet, nn_dist) - set_nearest_neighbour!(flatjets, jjet, ijet) - end - jjet = next_jet(flatjets, jjet) - end - - # Now scan over rightmost neighbour tiles - # This means moving in (η,ϕ) indexes, so now we need to map back from the 1D index - # to the 2D one... - for jtile_cartesian in rightmost_tiles(tiling_setup._n_tiles_eta, tiling_setup._n_tiles_phi, itile_cartesian[1], itile_cartesian[2]) - jtile = get_tile_linear_index(tiling_setup, jtile_cartesian[1], jtile_cartesian[2]) - ## Debug for checking that my index calculations are correct - # @assert jtile == tiling_setup._tile_linear_indexes[jtile_cartesian[1], jtile_cartesian[2]] - jjet = tiles[jtile] - while jjet!=0 - nn_dist = geometric_distance(eta(flatjets, ijet), phi(flatjets, ijet), - eta(flatjets, jjet), phi(flatjets, jjet)) - if nn_dist < nn_distance(flatjets, ijet) - set_nn_distance!(flatjets, ijet, nn_dist) - set_nearest_neighbour!(flatjets, ijet, jjet) - end - if nn_dist < nn_distance(flatjets, jjet) - set_nn_distance!(flatjets, jjet, nn_dist) - set_nearest_neighbour!(flatjets, jjet, ijet) - end - jjet = next_jet(flatjets, jjet) - end - end - ijet = next_jet(flatjets, ijet) - end - end - # Done with nearest neighbours, now calculate dij_distances - for ijet in 1:length(flatjets.kt2) - set_dij_distance!(flatjets, ijet, - get_dij_dist( - nn_distance(flatjets, ijet), - kt2(flatjets, ijet), - nearest_neighbour(flatjets, ijet) == 0 ? 0.0 : kt2(flatjets, nearest_neighbour(flatjets, ijet)), - R2 - ) - ) - - end -end - -""" -Find neighbour for a particular jet index -""" -const updated_pair_jets = Int[] # If any reverese mappings change their NN distance, need to record this to update dij -function find_jet_nearest_neighbour!(tiles, flatjets::FlatJets, tiling_setup::TilingDef, ijet::Int, R2) - @debug "Finding neighbours of $ijet" - set_nearest_neighbour!(flatjets, ijet, 0) - set_nn_distance!(flatjets, ijet, R2) - empty!(updated_pair_jets) - # Jet's in own tile (I may not be the first jet here) - itile_cartesian = get_tile_cartesian_indices(tiling_setup, tile_index(flatjets, ijet)) - jjet = tiles[tile_index(flatjets, ijet)] - while jjet != 0 - if ijet == jjet - jjet = next_jet(flatjets, jjet) - continue - end - nn_dist = geometric_distance(eta(flatjets, ijet), phi(flatjets, ijet), - eta(flatjets, jjet), phi(flatjets, jjet)) - if nn_dist < nn_distance(flatjets, ijet) - set_nn_distance!(flatjets, ijet, nn_dist) - set_nearest_neighbour!(flatjets, ijet, jjet) - end - if nn_dist < nn_distance(flatjets, jjet) - push!(updated_pair_jets, jjet) - set_nn_distance!(flatjets, jjet, nn_dist) - set_nearest_neighbour!(flatjets, jjet, ijet) - end - jjet = next_jet(flatjets, jjet) - end - for jtile_cartesian in neighbour_tiles(tiling_setup._n_tiles_eta, tiling_setup._n_tiles_phi, itile_cartesian[1], itile_cartesian[2]) - jtile = get_tile_linear_index(tiling_setup, jtile_cartesian[1], jtile_cartesian[2]) - jjet = tiles[jtile] - while jjet != 0 - nn_dist = geometric_distance(eta(flatjets, ijet), phi(flatjets, ijet), - eta(flatjets, jjet), phi(flatjets, jjet)) - if nn_dist < nn_distance(flatjets, ijet) - set_nn_distance!(flatjets, ijet, nn_dist) - set_nearest_neighbour!(flatjets, ijet, jjet) - end - if nn_dist < nn_distance(flatjets, jjet) - push!(updated_pair_jets, jjet) - set_nn_distance!(flatjets, jjet, nn_dist) - set_nearest_neighbour!(flatjets, jjet, ijet) - end - jjet = next_jet(flatjets, jjet) - end - end - set_dij_distance!(flatjets, ijet, - get_dij_dist( - nn_distance(flatjets, ijet), - kt2(flatjets, ijet), - nearest_neighbour(flatjets, ijet) == 0 ? 0.0 : kt2(flatjets, nearest_neighbour(flatjets, ijet)), - R2 - ) - ) - for jjet in updated_pair_jets - set_dij_distance!(flatjets, jjet, - get_dij_dist( - nn_distance(flatjets, jjet), - kt2(flatjets, jjet), - kt2(flatjets, nearest_neighbour(flatjets, jjet)), - R2 - ) - ) - end -end - -""" -Find all of the jets with a particular nearest neighbour -""" -const neighbours = Vector{Int}() -find_neighbours_of(flatjets::FlatJets, innA::Int, innB::Int) = begin - empty!(neighbours) - @inbounds for (ijet, nn) in enumerate(flatjets.nearest_neighbour) - if ((nn == innA) || (nn == innB)) && ((ijet != innA) && (ijet != innB)) - push!(neighbours, ijet) - end - end - neighbours -end - -find_neighbours_of(flatjets::FlatJets, inn::Int) = begin - empty!(neighbours) - @inbounds for (ijet, nn) in enumerate(flatjets.nearest_neighbour) - if nn == inn - push!(neighbours, ijet) - end - end - neighbours -end - -""" -Look for any jets which had a shuffled jet as a NN or next/prev jet, which has -moved from old_index to new_index -""" -function move_indexes!(flatjets::FlatJets, old_index, new_index) - # Do this in a simple loop - it's contiguous memory, very fast and can be vectorised - @debug "Moving jet $old_index to $new_index" - @turbo for ijet ∈ eachindex(flatjets.nearest_neighbour) - moveme = flatjets.nearest_neighbour[ijet] == old_index - flatjets.nearest_neighbour[ijet] = moveme ? new_index : flatjets.nearest_neighbour[ijet] - end - @turbo for ijet ∈ eachindex(flatjets.next_jet) - moveme = flatjets.next_jet[ijet] == old_index - flatjets.next_jet[ijet] = moveme ? new_index : flatjets.next_jet[ijet] - end - @turbo for ijet ∈ eachindex(flatjets.prev_jet) - moveme = flatjets.prev_jet[ijet] == old_index - flatjets.prev_jet[ijet] = moveme ? new_index : flatjets.prev_jet[ijet] - end -end - -""" -Move any tiles which have a first_jet that has been shuffled to the new -index -""" -function move_tile_first_jet!(tiles, old_index, new_index) - @debug "Moving first_jet $old_index to $new_index" - @turbo for itile ∈ eachindex(tiles) - moveme = tiles[itile] == old_index - tiles[itile] = moveme ? new_index : tiles[itile] - end -end - -""" -Find the index in the vector which has the lowest value of dij -""" -function find_lowest_dij_index(dij_distances) - imin = 0 - vmin = typemax(typeof(dij_distances[1])) - @turbo for i ∈ eachindex(dij_distances) - newmin = dij_distances[i] < vmin - vmin = newmin ? dij_distances[i] : vmin - imin = newmin ? i : imin - end - imin -end - -""" -Tiled jet reconstruction -""" -function tiled_jet_reconstruct_soa_global(objects::AbstractArray{T}; p = -1, R = 1.0, recombine = +) where T - # bounds - N::Int = length(objects) - @debug "Initial particles: $(N)" - - # params - R2::Float64 = R * R - p = (round(p) == p) ? Int(p) : p # integer p if possible - - # Input data - jet_objects = copy(objects) # We will append merged jets, so we need a copy - _kt2 = (JetReconstruction.pt.(objects) .^ 2) .^ p - _phi = JetReconstruction.phi.(objects) - _eta = JetReconstruction.eta.(objects) - _jet_index = collect(1:N) # Initial jets are just numbered 1:N, mapping directly to jet_objects - - # Each jet stores which tile it is in, so need the usual container for that - _tile_index = zeros(Int, N) - - # Linked list: this is the index of the next/previous jet for this tile (0 = the end/beginning) - _next_jet = zeros(Int, N) - _prev_jet = zeros(Int, N) - - # Nearest neighbour parameters - _nearest_neighbour = zeros(Int, N) - _nn_distance = fill(R2, N) - _dij_distance = fill(1.0e9, N) - - # returned values - jets = T[] # result - sequences = zeros(Int, N) - - flatjets = FlatJets(_kt2, _eta, _phi, _jet_index, _tile_index, - _next_jet, _prev_jet, - _nearest_neighbour, _nn_distance, _dij_distance) - - # Tiling - tiling_setup= setup_tiling(_eta, R) - @debug("Tiles: $(tiling_setup._n_tiles_eta)x$(tiling_setup._n_tiles_phi)") - - # Setup the tiling array (which is simply the index of the first jet in the tile) - tiles = zeros(Int, (tiling_setup._n_tiles_eta, tiling_setup._n_tiles_phi)) - - # Populate tile jet lists, from the initial particles - populate_tile_lists!(tiles, flatjets, tiling_setup) - - # Setup initial nn, nndist and dij values - find_all_tiled_nearest_neighbours!(tiles, flatjets, tiling_setup, R2) - - # A few allocations outside the loop - itouched_jets = Set{Int}() - - # At each iteration we either merge two jets to one, or finalise a jet - # Thus each time we lose one jet, and it therefore takes N iterations to complete - # the algorithm - for iteration in 1:N - @debug "Iteration $(iteration) - Active Jets $(lastindex(flatjets.kt2))" - - iclosejetA = find_lowest_dij_index(flatjets.dij_distance) - iclosejetB = nearest_neighbour(flatjets, iclosejetA) - @debug "Closest jets $iclosejetA, $iclosejetB: $(dij_distance(flatjets, iclosejetA))" - - # Finalise jet or merge jets? - if iclosejetB != 0 - # Merge jets A and B - jet-jet recombination - # If necessary relabel A & B to ensure jetB < jetA, that way if - # the larger of them == newtail then that ends up being jetA and - # the new jet that is added as jetB is inserted in a position that - # has a future! - if iclosejetA < iclosejetB - iclosejetA, iclosejetB = iclosejetB, iclosejetA - end - @debug "Merge jets: jet indexes $(jet_index(flatjets, iclosejetA)) in $(get_tile_cartesian_indices(tiling_setup, tile_index(flatjets, iclosejetA))), $(jet_index(flatjets, iclosejetB)) in $(get_tile_cartesian_indices(tiling_setup, tile_index(flatjets, iclosejetB)))" - newjet = recombine(jet_objects[jet_index(flatjets, iclosejetA)], - jet_objects[jet_index(flatjets, iclosejetB)]) - push!(jet_objects, newjet) - push!(sequences, 0) - inewjet = lastindex(jet_objects) - ## Set the sequence for the two merged jets, which is the merged jet index - sequences[jet_index(flatjets, iclosejetA)] = inewjet - sequences[jet_index(flatjets, iclosejetB)] = inewjet - # Remove merged jets from their tile jet lists - delete_from_tile!(tiles, tile_index(flatjets, iclosejetA), flatjets, iclosejetA) - delete_from_tile!(tiles, tile_index(flatjets, iclosejetB), flatjets, iclosejetB) - # Now find out which jets had A or B as their nearest neighbour - they will - # need to be rescanned - empty!(itouched_jets) - union!(itouched_jets, find_neighbours_of(flatjets, iclosejetA, iclosejetB)) - @debug "Jets to update from A/B neighbours: $(itouched_jets)" - #### - # Now push the newjet into jetB's slot, add it to its tile's jet list and to - # the set of jets that need to rescan their neighbours - @debug "Adding merged jet to slot $iclosejetB" - insert_flatjet!(flatjets, tiling_setup, p, iclosejetB, inewjet, newjet) - push_to_tile!(tiles, tile_index(flatjets, iclosejetB), flatjets, iclosejetB) - push!(itouched_jets, iclosejetB) # This is the newjet! - # Now kill jetA, shuffling if needed - shuffled_jet = suppress_flatjet!(flatjets, iclosejetA) - @debug "Killed jet $iclosejetA and shuffled jet $shuffled_jet" - # If we had a shuffled jet, we need to update any jet who had this jet as their - # neighbour - if shuffled_jet != 0 - move_indexes!(flatjets, shuffled_jet, iclosejetA) - move_tile_first_jet!(tiles, shuffled_jet, iclosejetA) - if shuffled_jet in itouched_jets - delete!(itouched_jets, shuffled_jet) - push!(itouched_jets, iclosejetA) - end - end - ## Updates - for ijet in itouched_jets - find_jet_nearest_neighbour!(tiles, flatjets, tiling_setup, ijet, R2) - end - else - # Finalise jet A - @debug "Finalise jet: jet index $(jet_index(flatjets, iclosejetA))" - push!(jets, jet_objects[jet_index(flatjets, iclosejetA)]) - # Remove finalsied jet from its tile jet lists - delete_from_tile!(tiles, tile_index(flatjets, iclosejetA), flatjets, iclosejetA) - # Now find out which jets had A as their nearest neighbour - they will - # need to be rescanned - empty!(itouched_jets) - union!(itouched_jets, find_neighbours_of(flatjets, iclosejetA)) - @debug "Jets to update from A neighbours: $(itouched_jets)" - # Now kill jetA, shuffling if needed - shuffled_jet = suppress_flatjet!(flatjets, iclosejetA) - @debug "Killed jet $iclosejetA and shuffled jet $shuffled_jet" - # If we had a shuffled jet, we need to update any jet who had this jet as their - # neighbour - if shuffled_jet != 0 - move_indexes!(flatjets, shuffled_jet, iclosejetA) - move_tile_first_jet!(tiles, shuffled_jet, iclosejetA) - if shuffled_jet in itouched_jets - delete!(itouched_jets, shuffled_jet) - push!(itouched_jets, iclosejetA) - end - end - ## Updates - for ijet in itouched_jets - find_jet_nearest_neighbour!(tiles, flatjets, tiling_setup, ijet, R2) - end - end - end - - # The sequences return value is a list of all jets that merged to this one - jets, sequences -end diff --git a/src/TiledAlgoSoATile.jl b/src/TiledAlgoSoATile.jl deleted file mode 100644 index c913246..0000000 --- a/src/TiledAlgoSoATile.jl +++ /dev/null @@ -1,335 +0,0 @@ -### Tiled Jet Reconstruction, SoA - -include("TiledUtilsSoA.jl") - -using Logging - -""" -Do a complete scan over tiles to find all nearest neighbour distances -""" -function find_all_nearest_neighbours!(tile_jets::Array{TiledJetSoA, 2}, tiling_setup::TilingDef, - flat_jets::FlatJetSoA, R2::AbstractFloat) - # We march over all tiles, evaluating the distance to - # - each jet in this tile - # - each jet in the rightmost neighbour tiles - # As we compare jet-to-jet, in both directions, the rightmost tiles ensure that the - # whole space is swept - @inbounds for itile in eachindex(tile_jets) - tile = tile_jets[itile] - if (tile._size == 0) - continue - end - @inbounds for ijet in 1:tile._size - # Could we do this in a broadcast way...? Would it be faster? - @inbounds for jjet in ijet+1:tile._size - _dist = geometric_distance(tile._eta[ijet], tile._phi[ijet], - tile._eta[jjet], tile._phi[jjet]) - if (_dist < tile._nndist[ijet]) - tile._nndist[ijet] = _dist - set_nn!(tile._nn[ijet], itile, jjet) - end - if (_dist < tile._nndist[jjet]) - tile._nndist[jjet] = _dist - set_nn!(tile._nn[jjet], itile, ijet) - end - end - @inbounds for jtile in tile._righttiles - tile2 = tile_jets[jtile] - @inbounds for jjet in 1:tile2._size - _dist = geometric_distance(tile._eta[ijet], tile._phi[ijet], - tile2._eta[jjet], tile2._phi[jjet]) - if (_dist < tile._nndist[ijet]) - tile._nndist[ijet] = _dist - set_nn!(tile._nn[ijet], jtile, jjet) - end - if (_dist < tile2._nndist[jjet]) - tile2._nndist[jjet] = _dist - set_nn!(tile2._nn[jjet], itile, ijet) - end - end - end - end - end - - # Now calculate the dij distances - min_dij = 1e20 - min_dij_itile = 0 - min_dij_ijet = 0 - @inbounds for (itile, tile) in enumerate(tile_jets) - @inbounds for ijet in 1:tile._size - if valid_nn(tile._nn[ijet]) - tile._dij[ijet] = tile._nndist[ijet] * - min(tile._kt2[ijet], tile_jets[tile._nn[ijet]._itile]._kt2[tile._nn[ijet]._ijet]) - else - tile._dij[ijet] = tile._nndist[ijet] * tile._kt2[ijet] - end - # And as this is a complete scan, find the minimum dij as well - if tile._dij[ijet] < min_dij - min_dij = tile._dij[ijet] - min_dij_itile = itile - min_dij_ijet = ijet - end - end - end - min_dij_itile, min_dij_ijet, min_dij -end - -""" -Scan a tile and all of its neighbours for jet distances -""" -function scan_neighbors!(tile_jets::Array{TiledJetSoA, 2}, jet_tile_index::TiledSoACoord, R2::AbstractFloat) - itile = jet_tile_index._itile - tile = tile_jets[itile] - ijet = jet_tile_index._ijet - # println("$(tile) $(ijet) $(tile._size)") - @inbounds for jjet in 1:tile._size - if jjet == ijet - continue - end - _dist = geometric_distance(tile._eta[ijet], tile._phi[ijet], - tile._eta[jjet], tile._phi[jjet]) - if (_dist < tile._nndist[ijet]) - tile._nndist[ijet] = _dist - set_nn!(tile._nn[ijet], itile, jjet) - end - if (_dist < tile._nndist[jjet]) - tile._nndist[jjet] = _dist - set_nn!(tile._nn[jjet], itile, ijet) - tile._dij[jjet] = tile._nndist[jjet] * min(tile._kt2[ijet], tile._kt2[jjet]) - end - end - @inbounds for jtile in tile._nntiles - tile_j = tile_jets[jtile] - @inbounds for jjet in 1:tile_j._size - _dist = geometric_distance(tile._eta[ijet], tile._phi[ijet], - tile_j._eta[jjet], tile_j._phi[jjet]) - if (_dist < tile._nndist[ijet]) - tile._nndist[ijet] = _dist - set_nn!(tile._nn[ijet], jtile, jjet) - end - if (_dist < tile_j._nndist[jjet]) - tile_j._nndist[jjet] = _dist - set_nn!(tile_j._nn[jjet], itile, ijet) - tile_j._dij[jjet] = tile_j._nndist[jjet] * min(tile._kt2[ijet], tile_j._kt2[jjet]) - end - end - end - if tile._nndist[ijet] != R2 - tile._dij[ijet] = tile._nndist[ijet] * min(tile._kt2[ijet], - tile_jets[tile._nn[ijet]._itile]._kt2[tile._nn[ijet]._ijet]) - end - -end - -"""Dump NN and dij distances for the current tiled jets, for debugging""" -function get_nn_str(tile_jets::Array{TiledJetSoA, 2}) - jet_nn_strings = Vector{Tuple{Int64, String}}(undef, 0) - for itile in eachindex(tile_jets) - tile = tile_jets[itile] - for ijet in 1:tile._size - push!(jet_nn_strings, (tile._index[ijet], "Jet $(tile._index[ijet]), NN is $(nnindex(tile_jets, itile, ijet)), geo $(tile._nndist[ijet]) dij $(tile._dij[ijet])\n")) - end - end - sort!(jet_nn_strings) - nn_str = "" - for (ijet, jet_str) in jet_nn_strings - nn_str *= jet_str - end - nn_str -end - -"""Validate distance calculations, for debugging""" -function validate_distances(tile_jets::Array{TiledJetSoA, 2}, flat_jets, ijetA::Int, ijetB::Int) - itileA = itjetA = itileB = itjetB = 0 - for itile in eachindex(tile_jets) - tile = tile_jets[itile] - for ijet in 1:tile._size - if tile._index[ijet] == ijetA - itileA = itile - itjetA = ijet - elseif tile._index[ijet] == ijetB - itileB = itile - itjetB = ijet - end - end - end - dist_string = "$(ijetA) $(flat_jets._eta[ijetA]) $(flat_jets._phi[ijetA]) $(flat_jets._kt2[ijetA]) ($(itileA),$(itjetA))\n" - dist_string *= "$(ijetB) $(flat_jets._eta[ijetB]) $(flat_jets._phi[ijetB]) $(flat_jets._kt2[ijetB]) ($(itileB),$(itjetB))\n" - dist_string *= "$(tile_jets[itileA]._eta[itjetA]) $(tile_jets[itileA]._phi[itjetA])\n" - dist_string *= "$(tile_jets[itileB]._eta[itjetB]) $(tile_jets[itileB]._phi[itjetB])\n" - dist = geometric_distance(tile_jets[itileA]._eta[itjetA], tile_jets[itileA]._phi[itjetA], tile_jets[itileB]._eta[itjetB], tile_jets[itileB]._phi[itjetB]) - dist_string *= "Geo $(dist)" - dist_string -end - -""" -Tiled jet reconstruction -""" -function tiled_jet_reconstruct_soa_tile(objects::AbstractArray{T}; p = -1, R = 1.0, recombine = +) where T - # bounds - N::Int = length(objects) - @debug "Initial particles: $(N)" - - # params - _R2::Float64 = R * R - _p = (round(p) == p) ? Int(p) : p # integer p if possible - ap = abs(_p) # absolute p - - # Input data - _objects = copy(objects) - sizehint!(_objects, N * 2) - _kt2 = (JetReconstruction.pt.(_objects) .^ 2) .^ _p - sizehint!(_kt2, N * 2) - _phi = JetReconstruction.phi.(_objects) - sizehint!(_phi, N * 2) - _eta = JetReconstruction.eta.(_objects) - sizehint!(_eta, N * 2) - _index = collect(1:N) # Initial jets are just numbered 1:N - sizehint!(_index, N * 2) - - # returned values - jets = T[] # result - _sequences = Vector{Int}[[x] for x in 1:N] - - flat_jets = FlatJetSoA(N, _kt2, _eta, _phi, _index) - - # Tiling - tiling_setup = setup_tiling(_eta, R) - tile_jets = Array{TiledJetSoA, 2}(undef, tiling_setup._n_tiles_eta, tiling_setup._n_tiles_phi) - - # Populate tiles, from the initial particles - populate_tiles!(tile_jets, tiling_setup, flat_jets, _R2) - - # Setup initial nn, nndist and dij values - min_dij_itile, min_dij_ijet, min_dij = find_all_nearest_neighbours!(tile_jets, tiling_setup, flat_jets, _R2) - - # Move some variables outside the loop, to avoid per-loop allocations - itouched_tiles = Set{Int}() - sizehint!(itouched_tiles, 12) - tainted_slots = Set{TiledSoACoord}() - sizehint!(tainted_slots, 4) - - # At each iteration we either merge two jets to one, or finalise a jet - # Thus each time we lose one jet, and it therefore takes N iterations to complete - # the algorithm - for iteration in 1:N - @debug "Iteration $(iteration)" - - # For the first iteration the nearest neighbour is known - if iteration > 1 - # Now find again the new nearest dij jets - min_dij = 1.0e20 - min_dij_itile = 0 - min_dij_ijet = 0 - @inbounds for itile in eachindex(tile_jets) - @inbounds for ijet in 1:tile_jets[itile]._size - if tile_jets[itile]._dij[ijet] < min_dij - min_dij_itile = itile - min_dij_ijet = ijet - min_dij = tile_jets[itile]._dij[ijet] - end - end - end - end - - @debug "$(min_dij) at ($(min_dij_itile), $(min_dij_ijet)) $(tile_jets[min_dij_itile]._index[min_dij_ijet]) -> $(tile_jets[min_dij_itile]._nn[min_dij_ijet])" - # Is this a merger or a final jet? - if tile_jets[min_dij_itile]._nn[min_dij_ijet]._itile == 0 - # Final jet - jet_merger = false - index_tile_jetA = TiledSoACoord(min_dij_itile, min_dij_ijet) - index_jetA = tile_jets[min_dij_itile]._index[min_dij_ijet] - empty!(tainted_slots) - push!(tainted_slots, index_tile_jetA) - push!(jets, _objects[index_jetA]) - push!(_sequences[index_jetA], 0) - @debug "Finalise jet $(tile_jets[min_dij_itile]._index[min_dij_ijet]) ($(_sequences[index_jetA])) $(JetReconstruction.pt(_objects[index_jetA]))" - push!(tainted_slots, remove_jet!(tile_jets, index_tile_jetA)) - else - # Merge two jets - jet_merger = true - index_tile_jetA = TiledSoACoord(min_dij_itile, min_dij_ijet) - index_tile_jetB = tile_jets[min_dij_itile]._nn[min_dij_ijet] - index_jetA = tile_jets[min_dij_itile]._index[min_dij_ijet] - index_jetB = nnindex(tile_jets, min_dij_itile, min_dij_ijet) - @debug "Merge jets $(index_jetA) ($(index_tile_jetA)) and $(index_jetB) ($(index_tile_jetB))" - merged_jet = recombine(_objects[index_jetA], _objects[index_jetB]) - - # If A and B are in the same tile, ensure that A is the earlier slot - # so that slots are filled up correctly - if (index_tile_jetA._itile == index_tile_jetB._itile) && (index_tile_jetA._ijet > index_tile_jetB._ijet) - index_tile_jetA, index_tile_jetB = index_tile_jetB, index_tile_jetA - index_jetA, index_jetB = index_jetB, index_jetA - end - - push!(_objects, merged_jet) - push!(flat_jets._index, length(_objects)) - push!(flat_jets._phi, JetReconstruction.phi(merged_jet)) - push!(flat_jets._eta, JetReconstruction.eta(merged_jet)) - push!(flat_jets._kt2, (JetReconstruction.pt(merged_jet)^2)^_p) - merged_jet_index = lastindex(_objects) - - ieta_merged_jet, iphi_merged_jet = get_tile(tiling_setup, flat_jets._eta[merged_jet_index], - flat_jets._phi[merged_jet_index]) - itile_merged_jet = get_tile_linear_index(tiling_setup, ieta_merged_jet, iphi_merged_jet) - - # Set the _sequence for the two merged jets, which is the merged jet index - push!(_sequences, [_sequences[index_jetA]; _sequences[index_jetB]; merged_jet_index]) - push!(_sequences[index_jetA], merged_jet_index) - push!(_sequences[index_jetB], merged_jet_index) - - - # Delete jetA and jetB from their tiles - empty!(tainted_slots) - push!(tainted_slots, index_tile_jetA) - push!(tainted_slots, index_tile_jetB) - if itile_merged_jet == index_tile_jetA._itile - # Put the new jet into jetA's slot - insert_jet!(tile_jets[itile_merged_jet], index_tile_jetA._ijet, merged_jet_index, flat_jets, _R2) - index_tile_merged_jet = TiledSoACoord(itile_merged_jet, index_tile_jetA._ijet) - # Now zap jetB - push!(tainted_slots, remove_jet!(tile_jets, index_tile_jetB)) - elseif itile_merged_jet == index_tile_jetB._itile - # Use jetB's slot - insert_jet!(tile_jets[itile_merged_jet], index_tile_jetB._ijet, merged_jet_index, flat_jets, _R2) - index_tile_merged_jet = TiledSoACoord(itile_merged_jet, index_tile_jetB._ijet) - # Now zap jetA - push!(tainted_slots, remove_jet!(tile_jets, index_tile_jetA)) - else - # Merged jet is in a different tile - add_jet!(tile_jets[itile_merged_jet], merged_jet_index, flat_jets, _R2) - index_tile_merged_jet = TiledSoACoord(itile_merged_jet, tile_jets[itile_merged_jet]._size) - # Now zap both A and B - push!(tainted_slots, remove_jet!(tile_jets, index_tile_jetA)) - push!(tainted_slots, remove_jet!(tile_jets, index_tile_jetB)) - end - - # For our new merged jet, scan for nearest neighbours - # Remember, this is pair-wise, so it will update all jets in its tile and neighbours - scan_neighbors!(tile_jets, index_tile_merged_jet, _R2) - end - - # Now take care of tainted neighbours - empty!(itouched_tiles) - push!(itouched_tiles, index_tile_jetA._itile) - union!(itouched_tiles, tile_jets[index_tile_jetA._itile]._nntiles) - if (jet_merger && index_tile_jetB._itile != index_tile_jetA._itile) - push!(itouched_tiles, index_tile_jetB._itile) - union!(itouched_tiles, tile_jets[index_tile_jetB._itile]._nntiles) - end - - # Scan over the touched tiles, look for jets whose _nn is tainted - @inbounds for itouched_tile in itouched_tiles - tile = tile_jets[itouched_tile] - @inbounds for ijet in 1:tile._size - if tile._nn[ijet] in tainted_slots - tile._nn[ijet] = TiledSoACoord(0, 0) - tile._nndist[ijet] = _R2 - scan_neighbors!(tile_jets, TiledSoACoord(itouched_tile, ijet), _R2) - end - end - end - end - # The sequences return value is a list of all jets that merged to this one - jets, _sequences -end diff --git a/src/TiledUtilsSoA.jl b/src/TiledUtilsSoA.jl deleted file mode 100644 index 58ccb22..0000000 --- a/src/TiledUtilsSoA.jl +++ /dev/null @@ -1,240 +0,0 @@ -""" -Structures and utilities used by the per-tile SoA implementation -""" - -import Base.== -import Base.hash - - -"""Nearest neighbour coordinates""" -mutable struct TiledSoACoord - _itile::Int # Jet tile index (flattened) - _ijet::Int # Jet position in this tile -end - -"""Correct hash for the coordinate pair""" -hash(x::TiledSoACoord, h::UInt) = begin - my_hash = hash(hash(x._itile, UInt(0)), hash(x._ijet, UInt(0))) - hash(my_hash, h) -end - -"""Equality operator for tiled coordinates""" -==(x::TiledSoACoord,y::TiledSoACoord) = hash(x, UInt(0))==hash(y, UInt(0)) - -"""Setter for nearest neighbour""" -set_nn!(mynn::TiledSoACoord, itile, ijet) = begin - mynn._itile = itile - mynn._ijet = ijet -end - -"""Do we have a NN, or not""" -valid_nn(mynn::TiledSoACoord) = begin - return mynn._itile > 0 -end - -abstract type JetSoA end - -"""Structure of arrays for tiled jet parameters, using an SoA layout -for computational efficiency""" -mutable struct TiledJetSoA <: JetSoA - _size::Int # Active jet count (can be less than the vector length) - _kt2::Vector{Float64} # p_t^-2p - _eta::Vector{Float64} # Rapidity - _phi::Vector{Float64} # Phi coordinate - _index::Vector{Int} # My jet index - _nn::Vector{TiledSoACoord} # Nearest neighbour location (if (0,0) no nearest neighbour) - _nndist::Vector{Float64} # Distance to my nearest neighbour - _dij::Vector{Float64} # Jet metric distance to my nearest neighbour - _righttiles::Vector{Int} # Indexes of all tiles to my right - _nntiles::Vector{Int} # Indexes of all neighbour tiles - # Inner constructor to ensure the sizehint for tile caches - function TiledJetSoA(n::Int) - my_tiledjet = new(n, - Vector{Float64}(undef, n), - Vector{Float64}(undef, n), - Vector{Float64}(undef, n), - Vector{Int}(undef, n), - Vector{TiledSoACoord}(undef, n), - Vector{Float64}(undef, n), - Vector{Float64}(undef, n), - Vector{Int}(undef, 0), - Vector{Int}(undef, 0) - ) - sizehint!(my_tiledjet._righttiles, 4) - sizehint!(my_tiledjet._nntiles, 8) - my_tiledjet - end -end - -"""Return the flat jet index of the nearest neighbour tile of a jet""" -nnindex(tile_jets::Array{TiledJetSoA, 2}, itile, ijet) = begin - if tile_jets[itile]._nn[ijet]._itile == 0 - return 0 - end - return tile_jets[tile_jets[itile]._nn[ijet]._itile]._index[tile_jets[itile]._nn[ijet]._ijet] -end - -"""Structure for the flatjet SoA, which is convenient""" -mutable struct FlatJetSoA <: JetSoA - _size::Int # Number of active entries (may be less than the vector size!) - _kt2::Vector{Float64} # p_t^-2 - _eta::Vector{Float64} # Rapidity - _phi::Vector{Float64} # Phi coordinate - _index::Vector{Int} # My jet index -end - -"""Insert a jet into a tile at the given slot""" -insert_jet!(tile::TiledJetSoA, slot::Int, index::Int, flat_jets::FlatJetSoA, R2::AbstractFloat) = begin - tile._kt2[slot] = flat_jets._kt2[index] - tile._eta[slot] = flat_jets._eta[index] - tile._phi[slot] = flat_jets._phi[index] - tile._index[slot] = index - tile._nn[slot] = TiledSoACoord(0, 0) - tile._nndist[slot] = R2 - tile._dij[slot] = R2 * tile._kt2[slot] -end - -"""Add a jet to a tile beyond the active slots""" -add_jet!(tile::TiledJetSoA, index::Int, flat_jets::FlatJetSoA, R2::AbstractFloat) = begin - tile._size += 1 - push!(tile._kt2, flat_jets._kt2[index]) - push!(tile._eta, flat_jets._eta[index]) - push!(tile._phi, flat_jets._phi[index]) - push!(tile._index, index) - push!(tile._nn, TiledSoACoord(0, 0)) - push!(tile._nndist, R2) - push!(tile._dij, R2 * flat_jets._kt2[index]) -end - -"""Remove a jet from the tile at the given tile index and repack if needed""" -remove_jet!(tile_jets::Array{TiledJetSoA, 2}, tile_index::TiledSoACoord) = begin - tile = tile_jets[tile_index._itile] - if tile._size != tile_index._ijet - # Need to copy the last jet into the slot, to ensure a contiguous array - # These two slots become tainted - tile._kt2[tile_index._ijet] = tile._kt2[tile._size] - tile._eta[tile_index._ijet] = tile._eta[tile._size] - tile._phi[tile_index._ijet] = tile._phi[tile._size] - tile._index[tile_index._ijet] = tile._index[tile._size] - tile._nn[tile_index._ijet] = tile._nn[tile._size] - tile._nndist[tile_index._ijet] = tile._nndist[tile._size] - tile._dij[tile_index._ijet] = tile._dij[tile._size] - tainted = TiledSoACoord(tile_index._itile, tile._size) - else - tainted = TiledSoACoord(0, 0) - end - # Physically reduce the arrays - safer while developing to avoid accidents! - # For optimised code, consider not doing this, as we can use _size as a bound when scanning - deleteat!(tile._kt2, tile._size) - deleteat!(tile._eta, tile._size) - - deleteat!(tile._phi, tile._size) - deleteat!(tile._index, tile._size) - deleteat!(tile._nn, tile._size) - deleteat!(tile._nndist, tile._size) - deleteat!(tile._dij, tile._size) - # Removing the last entry - tile._size -= 1 - return tainted -end - -get_jet(tile_jets::Array{TiledJetSoA, 2}, tile_index::TiledSoACoord) = begin - tile = tile_jets[tile_index._itile] - ijet = tile_index._ijet - return "($(tile._eta[ijet]), $(tile._phi[ijet]) [$(tile._index[ijet])]) -> $(tile._nn[ijet])" -end - - -"""Return the nth jet in the SoA""" -get_jet(j, n::Int) = begin - return j._index[n], j._eta[n], j._phi[n], j._kt2[n] -end - -# Getters - will work for both SoAs -kt2(j::JetSoA, n::Int) = j._kt2[n] -eta(j::JetSoA, n::Int) = j._eta[n] -phi(j::JetSoA, n::Int) = j._phi[n] -index(j::JetSoA, n::Int) = j._index[n] -nn(j::JetSoA, n::Int) = j._nn[n] -nndist(j::JetSoA, n::Int) = j._nndist[n] - -# Setters -set_kt2!(j::JetSoA, n::Int, v) = j._kt2[n] = v -set_eta!(j::JetSoA, n::Int, v) = j._eta[n] = v -set_phi!(j::JetSoA, n::Int, v) = j._phi[n] = v -set_index!(j::JetSoA, n::Int, v) = j._index[n] = v -set_nn!(j::JetSoA, n::Int, v::TiledSoACoord) = j._nn[n] = v -set_nndist!(j::JetSoA, n::Int, v) = j._nndist[n] = v - - -""" -Populate tiling structure with our initial jets and setup neighbour tile caches -""" -function populate_tiles!(tile_jets::Array{TiledJetSoA, 2}, tiling_setup::TilingDef, - flat_jets::FlatJetSoA, R2::AbstractFloat) - # This is a special case, where the initial particles are all - # "linear" in the flat_jets structure, so we scan through that - # and match each jet to a tile, so that we can assign correct size - # vectors in the tiled jets structure - tile_jet_count = Array{Vector{Int}, 2}(undef, tiling_setup._n_tiles_eta, tiling_setup._n_tiles_phi) - # Using fill() doesn't work as we fill all tiles with the same vector! - @inbounds for itile in eachindex(tile_jet_count) - tile_jet_count[itile] = Int[] - end - - # Find out where each jet lives, then push its index value to the correct tile - @inbounds for ijet in 1:flat_jets._size - ieta, iphi = get_tile(tiling_setup, eta(flat_jets, ijet), phi(flat_jets, ijet)) - push!(tile_jet_count[ieta, iphi], index(flat_jets, ijet)) - end - - # Now use the cached indexes to assign and fill the tiles - @inbounds for itile in eachindex(tile_jet_count) - ijets = tile_jet_count[itile] - this_tile_jets = TiledJetSoA(length(ijets)) - @inbounds for (itilejet, ijet) in enumerate(ijets) - set_kt2!(this_tile_jets, itilejet, kt2(flat_jets, ijet)) - set_eta!(this_tile_jets, itilejet, eta(flat_jets, ijet)) - set_phi!(this_tile_jets, itilejet, phi(flat_jets, ijet)) - set_index!(this_tile_jets, itilejet, index(flat_jets, ijet)) - set_nn!(this_tile_jets, itilejet, TiledSoACoord(0,0)) - set_nndist!(this_tile_jets, itilejet, R2) - end - tile_jets[itile] = this_tile_jets - end - populate_tile_cache!(tile_jets, tiling_setup) -end - -""" -For each tile, populate a cache of the nearest tile neighbours -""" -function populate_tile_cache!(tile_jets::Array{TiledJetSoA, 2}, tiling_setup::TilingDef) - # To help with later iterations, we now find and cache neighbour tile indexes - @inbounds for ieta in 1:tiling_setup._n_tiles_eta - @inbounds for iphi in 1:tiling_setup._n_tiles_phi - # Clamping ensures we don't go beyond the limits of the eta tiling (which do not wrap) - @inbounds for jeta in clamp(ieta - 1, 1, tiling_setup._n_tiles_eta):clamp(ieta + 1, 1, tiling_setup._n_tiles_eta) - δeta = jeta - ieta - @inbounds for jphi in iphi-1:iphi+1 - if (jeta == ieta && jphi == iphi) - continue - end - # Phi tiles wrap around to meet each other - δphi = jphi - iphi # Hold this unwrapped value for rightmost comparison - if (jphi == 0) - jphi = tiling_setup._n_tiles_phi - elseif (jphi == tiling_setup._n_tiles_phi + 1) - jphi = 1 - end - # Tile is a neighbour - tile_index = get_tile_linear_index(tiling_setup, jeta, jphi) - push!(tile_jets[ieta, iphi]._nntiles, tile_index) - # Only the tile directly above or to the right are _righttiles - if (((δeta == -1) && (δphi == 0)) || (δphi == 1)) - push!(tile_jets[ieta, iphi]._righttiles, tile_index) - end - end - end - end - end -end diff --git a/src/Utils.jl b/src/Utils.jl index b48e702..c9c34a6 100644 --- a/src/Utils.jl +++ b/src/Utils.jl @@ -1,6 +1,6 @@ # Utility functions, which can be used by different top level scripts -"""Read HepMC3 event and keep final state particles""" +"""Read HepMC3 event and keep final state particles (return PseudoJets)""" function read_final_state_particles(fname; maxevents = -1, skipevents = 0) f = open(fname) @@ -22,6 +22,33 @@ function read_final_state_particles(fname; maxevents = -1, skipevents = 0) end @info "Total Events: $(length(events))" + @debug events + events +end + +"""Read HepMC3 event and keep final state particles (return LorentzVectors)""" +function read_final_state_particles_lv(fname; maxevents = -1, skipevents = 0) + f = open(fname) + + events = Vector{LorentzVector{Float64}}[] + + ipart = 1 + HepMC3.read_events(f, maxevents = maxevents, skipevents = skipevents) do parts + input_particles = LorentzVector{Float64}[] + for p in parts + if p.status == 1 + push!( + input_particles, + LorentzVector(p.momentum.t, p.momentum.x, p.momentum.y, p.momentum.z), + ) + end + end + push!(events, input_particles) + ipart += 1 + end + + @info "Total Events: $(length(events))" + @debug events events end @@ -67,7 +94,35 @@ function final_jets(jets::Vector{PseudoJet}, ptmin::AbstractFloat) dcut = ptmin^2 if pt2(jet) > dcut count += 1 - push!(final_jets, FinalJet(rap(jet), phi(jet), sqrt(pt2(jet)))) + push!(final_jets, FinalJet(rapidity(jet), phi(jet), sqrt(pt2(jet)))) + end + end + final_jets +end + +function final_jets(jets::Vector{LorentzVector}, ptmin::AbstractFloat) + count = 0 + final_jets = Vector{FinalJet}() + sizehint!(final_jets, 6) + dcut = ptmin^2 + for jet in jets + if LorentzVectorHEP.pt(jet)^2 > dcut + count += 1 + push!(final_jets, FinalJet(LorentzVectorHEP.rapidity(jet), LorentzVectorHEP.phi(jet), LorentzVectorHEP.pt(jet))) + end + end + final_jets +end + +function final_jets(jets::Vector{LorentzVectorCyl}, ptmin::AbstractFloat) + count = 0 + final_jets = Vector{FinalJet}() + sizehint!(final_jets, 6) + dcut = ptmin^2 + for jet in jets + if LorentzVectorHEP.pt(jet)^2 > dcut + count += 1 + push!(final_jets, FinalJet(LorentzVectorHEP.eta(jet), LorentzVectorHEP.phi(jet), LorentzVectorHEP.pt(jet))) end end final_jets diff --git a/test/data/jet_collections_fastjet.json b/test/data/jet_collections_fastjet.json new file mode 100644 index 0000000..86f0d34 --- /dev/null +++ b/test/data/jet_collections_fastjet.json @@ -0,0 +1,4362 @@ +[ + { + "jetid": 1, + "jets": [ + { + "rap": 3.35635648, + "phi": 2.68836383, + "pt": 24.83429978 + }, + { + "rap": -2.18503736, + "phi": 5.30438152, + "pt": 19.89959594 + }, + { + "rap": 1.36876323, + "phi": 6.15824888, + "pt": 8.95819427 + } + ] + }, + { + "jetid": 2, + "jets": [ + { + "rap": -1.29396553, + "phi": 5.49161728, + "pt": 35.68347579 + }, + { + "rap": 1.32327026, + "phi": 2.43974319, + "pt": 21.58688573 + }, + { + "rap": 1.01245908, + "phi": 1.9905492, + "pt": 18.1274308 + }, + { + "rap": -3.35185572, + "phi": 2.19540272, + "pt": 11.98498134 + }, + { + "rap": 1.77143648, + "phi": 5.06224577, + "pt": 11.67129138 + } + ] + }, + { + "jetid": 3, + "jets": [ + { + "rap": 2.52621253, + "phi": 5.55618673, + "pt": 18.48635867 + }, + { + "rap": 2.52921435, + "phi": 2.5858984, + "pt": 17.53067774 + }, + { + "rap": 4.57319341, + "phi": 4.26972618, + "pt": 8.87757818 + }, + { + "rap": 1.75010192, + "phi": 2.0279478, + "pt": 5.06765805 + } + ] + }, + { + "jetid": 4, + "jets": [ + { + "rap": -3.6000889, + "phi": 4.41465604, + "pt": 22.34660671 + }, + { + "rap": -1.06879143, + "phi": 0.58334587, + "pt": 16.37930922 + }, + { + "rap": 1.03423013, + "phi": 4.38693121, + "pt": 10.08822235 + }, + { + "rap": -2.75293618, + "phi": 4.98340453, + "pt": 9.62035452 + }, + { + "rap": -0.86106785, + "phi": 2.86081819, + "pt": 8.74462895 + }, + { + "rap": 4.15755763, + "phi": 0.93531031, + "pt": 7.40693086 + }, + { + "rap": 0.78224427, + "phi": 0.20976581, + "pt": 7.20178622 + }, + { + "rap": 0.0361121, + "phi": 1.0691108, + "pt": 7.16731332 + }, + { + "rap": 4.90928412, + "phi": 1.31763908, + "pt": 6.731507 + }, + { + "rap": -4.10121684, + "phi": 6.13769142, + "pt": 6.69228155 + }, + { + "rap": -2.34515022, + "phi": 0.46615119, + "pt": 5.46843643 + }, + { + "rap": -0.00793253, + "phi": 1.801414, + "pt": 5.26841956 + }, + { + "rap": 2.42589704, + "phi": 2.43385383, + "pt": 5.26840475 + }, + { + "rap": -4.07518127, + "phi": 3.17324965, + "pt": 5.23440031 + } + ] + }, + { + "jetid": 5, + "jets": [ + { + "rap": 1.51746024, + "phi": 6.10728928, + "pt": 29.2915043 + }, + { + "rap": 4.00989773, + "phi": 3.35607135, + "pt": 25.68535118 + }, + { + "rap": -0.20033213, + "phi": 2.46449272, + "pt": 12.32853543 + }, + { + "rap": -1.21413834, + "phi": 4.67988556, + "pt": 12.11037323 + }, + { + "rap": 3.59652224, + "phi": 5.95274508, + "pt": 9.23877412 + }, + { + "rap": 1.92154289, + "phi": 2.97898007, + "pt": 9.1967848 + }, + { + "rap": 2.64235234, + "phi": 5.59566664, + "pt": 9.12489762 + }, + { + "rap": -2.12735144, + "phi": 2.43214047, + "pt": 7.43036564 + }, + { + "rap": 2.4175818, + "phi": 3.67759239, + "pt": 7.31584383 + }, + { + "rap": -2.42047257, + "phi": 1.20195366, + "pt": 7.17763041 + }, + { + "rap": -1.85261632, + "phi": 3.24197682, + "pt": 6.97909616 + }, + { + "rap": 4.57492099, + "phi": 1.22681112, + "pt": 6.46286469 + }, + { + "rap": 0.56547878, + "phi": 2.25956615, + "pt": 6.44094567 + }, + { + "rap": -1.48846441, + "phi": 5.19698648, + "pt": 6.24309049 + }, + { + "rap": -2.52037432, + "phi": 2.20600879, + "pt": 6.01960524 + }, + { + "rap": -1.52136689, + "phi": 4.22641631, + "pt": 5.45736783 + } + ] + }, + { + "jetid": 6, + "jets": [ + { + "rap": -1.17413991, + "phi": 1.09505304, + "pt": 23.50089909 + }, + { + "rap": 1.84431528, + "phi": 3.2278215, + "pt": 15.35349284 + }, + { + "rap": -0.81249159, + "phi": 5.89571528, + "pt": 13.30928934 + }, + { + "rap": -0.94731344, + "phi": 2.6408007, + "pt": 8.67400488 + }, + { + "rap": 0.79032299, + "phi": 4.77835786, + "pt": 7.92850192 + }, + { + "rap": 1.93519206, + "phi": 2.30272459, + "pt": 6.18020988 + }, + { + "rap": -0.69522031, + "phi": 0.92155634, + "pt": 5.12728461 + } + ] + }, + { + "jetid": 7, + "jets": [ + { + "rap": 3.78093834, + "phi": 4.72933236, + "pt": 31.58144663 + }, + { + "rap": 1.78627745, + "phi": 1.83187236, + "pt": 27.1442145 + }, + { + "rap": -2.63721673, + "phi": 1.30448676, + "pt": 6.34748385 + } + ] + }, + { + "jetid": 8, + "jets": [ + { + "rap": 4.04329648, + "phi": 1.49198388, + "pt": 31.50568192 + }, + { + "rap": 2.66043875, + "phi": 4.37723279, + "pt": 21.32095741 + }, + { + "rap": 3.10408462, + "phi": 4.57804853, + "pt": 10.98017471 + }, + { + "rap": -3.12696594, + "phi": 0.6368243, + "pt": 8.12284964 + }, + { + "rap": 2.43707449, + "phi": 5.15613495, + "pt": 7.44645611 + }, + { + "rap": 1.59542154, + "phi": 2.08646564, + "pt": 6.77863691 + }, + { + "rap": 2.65640256, + "phi": 5.84409582, + "pt": 6.20807881 + }, + { + "rap": -2.43686775, + "phi": 6.18810947, + "pt": 5.71334783 + }, + { + "rap": 0.25861072, + "phi": 5.26967131, + "pt": 5.43953451 + }, + { + "rap": -1.56188444, + "phi": 3.08380551, + "pt": 5.28543681 + } + ] + }, + { + "jetid": 9, + "jets": [ + { + "rap": -0.33820031, + "phi": 3.60418377, + "pt": 22.89306266 + }, + { + "rap": 0.04483094, + "phi": 0.23396081, + "pt": 9.43407103 + }, + { + "rap": -3.64255933, + "phi": 0.59048654, + "pt": 6.90464536 + }, + { + "rap": -3.62167164, + "phi": 5.59277433, + "pt": 6.33654109 + }, + { + "rap": 0.63890587, + "phi": 0.26769197, + "pt": 6.2570091 + } + ] + }, + { + "jetid": 10, + "jets": [ + { + "rap": 1.1795804, + "phi": 2.17021382, + "pt": 36.5507518 + }, + { + "rap": -0.39693619, + "phi": 6.20943497, + "pt": 32.85482814 + }, + { + "rap": 2.52831789, + "phi": 5.01691758, + "pt": 27.3831182 + }, + { + "rap": -2.22039708, + "phi": 2.53402339, + "pt": 20.24838086 + }, + { + "rap": 2.96112835, + "phi": 4.23065675, + "pt": 15.76016901 + }, + { + "rap": 1.7374121, + "phi": 1.79886129, + "pt": 5.01854043 + } + ] + }, + { + "jetid": 11, + "jets": [ + { + "rap": 1.664624, + "phi": 0.7755939, + "pt": 17.14605685 + }, + { + "rap": 2.70489795, + "phi": 3.94732797, + "pt": 8.82339015 + }, + { + "rap": 4.89506769, + "phi": 2.73660932, + "pt": 6.1018262 + }, + { + "rap": 2.18013493, + "phi": 2.66823489, + "pt": 5.10759484 + } + ] + }, + { + "jetid": 12, + "jets": [ + { + "rap": -3.14962744, + "phi": 0.5214671, + "pt": 56.3254634 + }, + { + "rap": -0.24165705, + "phi": 3.57899081, + "pt": 32.9362266 + }, + { + "rap": -3.87369308, + "phi": 3.78932748, + "pt": 20.43134274 + }, + { + "rap": -3.41431675, + "phi": 0.90273946, + "pt": 9.25199302 + }, + { + "rap": -3.47684312, + "phi": 4.7766388, + "pt": 6.17510367 + } + ] + }, + { + "jetid": 13, + "jets": [ + { + "rap": 0.1374073, + "phi": 1.87320506, + "pt": 11.53466838 + }, + { + "rap": -1.36497747, + "phi": 0.13186385, + "pt": 9.9417517 + }, + { + "rap": -2.91939246, + "phi": 1.62156565, + "pt": 9.02007682 + }, + { + "rap": -2.0795162, + "phi": 4.54409845, + "pt": 8.61481286 + }, + { + "rap": -1.3026033, + "phi": 5.29487551, + "pt": 7.36161522 + }, + { + "rap": -3.28388267, + "phi": 2.78877243, + "pt": 7.19858138 + }, + { + "rap": -2.45485764, + "phi": 4.78343435, + "pt": 6.75181354 + }, + { + "rap": -0.39568218, + "phi": 1.88951263, + "pt": 5.47255466 + }, + { + "rap": -1.79125906, + "phi": 0.4112483, + "pt": 5.1989312 + } + ] + }, + { + "jetid": 14, + "jets": [ + { + "rap": -0.47688425, + "phi": 1.3094454, + "pt": 39.59518833 + }, + { + "rap": 1.56994294, + "phi": 2.58322313, + "pt": 16.22571243 + }, + { + "rap": -0.79471851, + "phi": 0.30290435, + "pt": 15.51036651 + }, + { + "rap": 5.39781894, + "phi": 3.79072736, + "pt": 15.49635858 + }, + { + "rap": -0.24543333, + "phi": 5.43352191, + "pt": 14.93208354 + }, + { + "rap": -1.51240299, + "phi": 0.50957812, + "pt": 14.77758354 + }, + { + "rap": 1.64725745, + "phi": 4.71766372, + "pt": 12.69470443 + }, + { + "rap": -2.22341966, + "phi": 5.66617838, + "pt": 11.77801131 + }, + { + "rap": -0.25326752, + "phi": 3.86875156, + "pt": 9.19985078 + }, + { + "rap": -3.6430709, + "phi": 3.36606275, + "pt": 9.15242887 + }, + { + "rap": 2.13951606, + "phi": 3.41991361, + "pt": 7.63068382 + }, + { + "rap": -1.5846547, + "phi": 4.45151124, + "pt": 7.34152512 + }, + { + "rap": -2.77804142, + "phi": 5.27729622, + "pt": 7.17426973 + }, + { + "rap": 1.08144286, + "phi": 4.15756839, + "pt": 6.91748004 + }, + { + "rap": 2.48425822, + "phi": 4.34390269, + "pt": 6.79104516 + }, + { + "rap": -1.53684177, + "phi": 5.84561369, + "pt": 6.65588683 + }, + { + "rap": -0.29665688, + "phi": 1.81411897, + "pt": 6.6321586 + }, + { + "rap": 1.0373415, + "phi": 2.19975137, + "pt": 6.60352738 + }, + { + "rap": -2.54023384, + "phi": 2.7605771, + "pt": 6.39099513 + }, + { + "rap": -4.44039707, + "phi": 0.15972957, + "pt": 6.20537604 + }, + { + "rap": -2.01065082, + "phi": 1.8325315, + "pt": 6.10110275 + }, + { + "rap": 0.8893975, + "phi": 0.87211139, + "pt": 5.51201193 + }, + { + "rap": 2.95521251, + "phi": 2.748836, + "pt": 5.25143272 + }, + { + "rap": 0.44771018, + "phi": 5.87537805, + "pt": 5.07186595 + } + ] + }, + { + "jetid": 15, + "jets": [ + { + "rap": 0.87239604, + "phi": 5.45041439, + "pt": 36.58676106 + }, + { + "rap": -1.06422928, + "phi": 2.33090777, + "pt": 11.29175446 + }, + { + "rap": 3.33526431, + "phi": 2.37276887, + "pt": 8.32223005 + }, + { + "rap": -4.63224457, + "phi": 6.05890784, + "pt": 7.59987102 + }, + { + "rap": -1.8781512, + "phi": 5.47575228, + "pt": 7.53702604 + }, + { + "rap": -4.88029608, + "phi": 2.31237405, + "pt": 6.41566541 + }, + { + "rap": 2.70160945, + "phi": 3.71330845, + "pt": 6.38376342 + }, + { + "rap": -1.60718561, + "phi": 2.32493851, + "pt": 6.36918325 + }, + { + "rap": 1.37938293, + "phi": 1.87758848, + "pt": 5.55588597 + }, + { + "rap": -0.47995291, + "phi": 3.03236152, + "pt": 5.3207285 + } + ] + }, + { + "jetid": 16, + "jets": [ + { + "rap": -0.63153165, + "phi": 0.84341198, + "pt": 21.14094711 + }, + { + "rap": 0.79215635, + "phi": 3.34987295, + "pt": 19.9078355 + }, + { + "rap": 3.25339733, + "phi": 1.03864084, + "pt": 18.61709566 + }, + { + "rap": -3.59035004, + "phi": 5.36701912, + "pt": 11.69501706 + }, + { + "rap": 1.95420095, + "phi": 3.60683591, + "pt": 10.54419295 + }, + { + "rap": 1.32516489, + "phi": 3.83109345, + "pt": 10.41300781 + }, + { + "rap": -1.94719758, + "phi": 4.12325465, + "pt": 7.40752426 + }, + { + "rap": 2.15985747, + "phi": 4.04354489, + "pt": 7.22413907 + }, + { + "rap": 3.71509017, + "phi": 0.60786413, + "pt": 5.35303273 + }, + { + "rap": 0.40228037, + "phi": 1.73595104, + "pt": 5.34972673 + }, + { + "rap": -0.96259256, + "phi": 3.92349533, + "pt": 5.0296159 + } + ] + }, + { + "jetid": 17, + "jets": [ + { + "rap": 0.84807318, + "phi": 2.98516228, + "pt": 29.18359101 + }, + { + "rap": 3.73672205, + "phi": 6.1492938, + "pt": 15.79482524 + }, + { + "rap": -2.39378808, + "phi": 0.26015374, + "pt": 10.07934654 + }, + { + "rap": -2.23777394, + "phi": 6.07084393, + "pt": 9.55249132 + } + ] + }, + { + "jetid": 18, + "jets": [ + { + "rap": -2.39391497, + "phi": 5.99387906, + "pt": 21.78813277 + }, + { + "rap": -1.74490796, + "phi": 2.56246273, + "pt": 18.32064021 + }, + { + "rap": -3.87159821, + "phi": 2.14712334, + "pt": 8.88999564 + }, + { + "rap": 1.12338326, + "phi": 6.08296975, + "pt": 5.60471139 + }, + { + "rap": 1.196868, + "phi": 3.93238459, + "pt": 5.50034218 + }, + { + "rap": 3.09167658, + "phi": 3.33951272, + "pt": 5.47586256 + }, + { + "rap": -2.21196333, + "phi": 1.35655249, + "pt": 5.33386356 + } + ] + }, + { + "jetid": 19, + "jets": [ + { + "rap": -1.92721086, + "phi": 2.51010234, + "pt": 29.52029045 + }, + { + "rap": 0.39448, + "phi": 5.59479588, + "pt": 23.42478757 + }, + { + "rap": -1.18938206, + "phi": 1.53451494, + "pt": 5.77973735 + } + ] + }, + { + "jetid": 20, + "jets": [ + { + "rap": -3.82016853, + "phi": 5.79978032, + "pt": 14.48225854 + }, + { + "rap": -3.96121231, + "phi": 0.15208152, + "pt": 13.55785452 + }, + { + "rap": -3.94049556, + "phi": 1.38028187, + "pt": 8.23788142 + }, + { + "rap": -1.9836822, + "phi": 2.78622261, + "pt": 6.05901817 + }, + { + "rap": 1.59859754, + "phi": 1.49082161, + "pt": 5.75424023 + }, + { + "rap": -2.78024846, + "phi": 3.3650565, + "pt": 5.64603349 + }, + { + "rap": -3.44204681, + "phi": 2.85397428, + "pt": 5.52271648 + }, + { + "rap": -4.02908312, + "phi": 0.64438159, + "pt": 5.25295326 + } + ] + }, + { + "jetid": 21, + "jets": [ + { + "rap": 1.83811983, + "phi": 0.88966879, + "pt": 27.03914544 + }, + { + "rap": 3.51903519, + "phi": 3.87753782, + "pt": 17.04595072 + }, + { + "rap": -1.36780549, + "phi": 4.46083451, + "pt": 12.19889738 + } + ] + }, + { + "jetid": 22, + "jets": [ + { + "rap": -3.73368334, + "phi": 4.68085841, + "pt": 26.08180488 + }, + { + "rap": -2.81084787, + "phi": 2.12756484, + "pt": 17.44326063 + }, + { + "rap": -3.83116815, + "phi": 1.31716677, + "pt": 11.61887415 + }, + { + "rap": -1.86613782, + "phi": 6.26216283, + "pt": 11.38888156 + }, + { + "rap": -1.98242425, + "phi": 1.79559166, + "pt": 9.47057162 + }, + { + "rap": -1.21582963, + "phi": 5.55815706, + "pt": 6.66419303 + }, + { + "rap": 4.23760661, + "phi": 3.63790724, + "pt": 5.58246889 + }, + { + "rap": -1.90860599, + "phi": 1.40091439, + "pt": 5.089263 + } + ] + }, + { + "jetid": 23, + "jets": [ + { + "rap": 2.90234363, + "phi": 2.58135218, + "pt": 21.41762504 + }, + { + "rap": -1.74431734, + "phi": 0.17026693, + "pt": 11.86706235 + }, + { + "rap": 2.11923121, + "phi": 3.72465919, + "pt": 10.02977198 + }, + { + "rap": -2.01600967, + "phi": 5.77900884, + "pt": 9.07214436 + }, + { + "rap": 3.09512109, + "phi": 4.30970119, + "pt": 8.91173841 + }, + { + "rap": -5.48129275, + "phi": 1.03183521, + "pt": 7.16706255 + }, + { + "rap": 3.83169534, + "phi": 3.98604398, + "pt": 6.49372897 + } + ] + }, + { + "jetid": 24, + "jets": [ + { + "rap": -0.57051877, + "phi": 2.13290237, + "pt": 20.72314149 + }, + { + "rap": 0.08115969, + "phi": 5.24683426, + "pt": 18.8030857 + } + ] + }, + { + "jetid": 25, + "jets": [ + { + "rap": 0.51668219, + "phi": 2.97232316, + "pt": 30.03658562 + }, + { + "rap": -3.74458779, + "phi": 5.03805258, + "pt": 19.92017897 + }, + { + "rap": -2.21775184, + "phi": 0.61250408, + "pt": 15.8096051 + }, + { + "rap": -1.28582133, + "phi": 1.05151496, + "pt": 7.38704859 + }, + { + "rap": 0.46636776, + "phi": 3.72000128, + "pt": 6.36944074 + }, + { + "rap": -0.91573074, + "phi": 0.42891385, + "pt": 5.4937362 + }, + { + "rap": 0.08801557, + "phi": 4.36890549, + "pt": 5.29979113 + }, + { + "rap": 0.33487218, + "phi": 2.43213647, + "pt": 5.02786122 + } + ] + }, + { + "jetid": 26, + "jets": [ + { + "rap": 2.45960516, + "phi": 3.88142326, + "pt": 42.77793814 + }, + { + "rap": 3.86992729, + "phi": 1.42992821, + "pt": 18.39773156 + }, + { + "rap": 1.85517971, + "phi": 0.73863771, + "pt": 15.52035935 + }, + { + "rap": 2.94904558, + "phi": 6.17423036, + "pt": 11.35315763 + } + ] + }, + { + "jetid": 27, + "jets": [ + { + "rap": 0.77131642, + "phi": 3.82705972, + "pt": 35.40628248 + }, + { + "rap": -2.30284287, + "phi": 5.55441357, + "pt": 17.41682074 + }, + { + "rap": 1.60568491, + "phi": 0.89038949, + "pt": 16.60658105 + }, + { + "rap": 0.4429229, + "phi": 1.89389961, + "pt": 11.9219618 + }, + { + "rap": -1.79646225, + "phi": 3.38526895, + "pt": 10.43231234 + }, + { + "rap": -0.35565198, + "phi": 0.30971223, + "pt": 9.21983223 + }, + { + "rap": -2.45247311, + "phi": 0.0516068, + "pt": 8.99210991 + }, + { + "rap": 4.60204888, + "phi": 3.73833896, + "pt": 8.88260906 + }, + { + "rap": -0.03136169, + "phi": 2.90580482, + "pt": 8.67297574 + }, + { + "rap": 0.97476096, + "phi": 0.77461309, + "pt": 7.84177375 + }, + { + "rap": 1.01681213, + "phi": 1.97497602, + "pt": 7.00798304 + }, + { + "rap": 1.64822399, + "phi": 5.03502731, + "pt": 6.1622655 + }, + { + "rap": 3.37398049, + "phi": 1.80814064, + "pt": 5.86162932 + }, + { + "rap": -3.18718446, + "phi": 4.31494815, + "pt": 5.75137174 + }, + { + "rap": -4.40159767, + "phi": 0.59099845, + "pt": 5.41783526 + }, + { + "rap": -1.72540479, + "phi": 6.19995813, + "pt": 5.27371495 + } + ] + }, + { + "jetid": 28, + "jets": [ + { + "rap": -3.16670119, + "phi": 2.89621113, + "pt": 20.67691172 + }, + { + "rap": -4.61156429, + "phi": 5.93009143, + "pt": 20.29633517 + }, + { + "rap": -0.42387519, + "phi": 6.20284475, + "pt": 13.34944706 + }, + { + "rap": -2.29257205, + "phi": 3.06547666, + "pt": 12.19823739 + }, + { + "rap": -0.05940605, + "phi": 0.30734166, + "pt": 5.62275548 + } + ] + }, + { + "jetid": 29, + "jets": [ + { + "rap": 3.39979791, + "phi": 4.96066435, + "pt": 10.97783759 + }, + { + "rap": 2.89103054, + "phi": 0.05938771, + "pt": 6.84793428 + }, + { + "rap": 3.1467668, + "phi": 1.96194274, + "pt": 6.57017415 + } + ] + }, + { + "jetid": 30, + "jets": [ + { + "rap": 0.88639522, + "phi": 2.49324638, + "pt": 56.72544863 + }, + { + "rap": 1.20342743, + "phi": 5.37423351, + "pt": 31.05368026 + }, + { + "rap": -1.11689034, + "phi": 5.72939699, + "pt": 19.66646727 + }, + { + "rap": -0.43369899, + "phi": 2.83151956, + "pt": 13.63770202 + }, + { + "rap": 2.72697135, + "phi": 5.28052275, + "pt": 12.81533378 + }, + { + "rap": -0.93168849, + "phi": 0.08490773, + "pt": 10.68147754 + }, + { + "rap": 3.03949964, + "phi": 1.55181026, + "pt": 8.94473255 + }, + { + "rap": -2.71709288, + "phi": 6.04400089, + "pt": 8.93608142 + }, + { + "rap": 2.35732261, + "phi": 5.78856525, + "pt": 7.784274 + }, + { + "rap": 4.31716503, + "phi": 1.63165571, + "pt": 7.06754033 + }, + { + "rap": -2.17789145, + "phi": 3.5170153, + "pt": 6.70266396 + }, + { + "rap": -0.11866952, + "phi": 2.16158359, + "pt": 6.33582519 + } + ] + }, + { + "jetid": 31, + "jets": [ + { + "rap": 0.49773401, + "phi": 1.1442432, + "pt": 22.86750303 + }, + { + "rap": 0.44723607, + "phi": 3.67145839, + "pt": 20.0397718 + }, + { + "rap": 0.86538834, + "phi": 5.16620833, + "pt": 11.43269605 + } + ] + }, + { + "jetid": 32, + "jets": [ + { + "rap": 0.51036746, + "phi": 3.68487252, + "pt": 24.83743363 + }, + { + "rap": 1.5859647, + "phi": 2.82007718, + "pt": 18.24654642 + }, + { + "rap": 2.54438399, + "phi": 0.56549233, + "pt": 18.24200057 + }, + { + "rap": 0.05362132, + "phi": 4.08040715, + "pt": 13.13862998 + }, + { + "rap": 0.14550491, + "phi": 6.0374313, + "pt": 11.8287509 + }, + { + "rap": 0.52997843, + "phi": 1.43101166, + "pt": 9.35451277 + }, + { + "rap": 3.33199475, + "phi": 5.92164695, + "pt": 9.1614528 + }, + { + "rap": 1.02121642, + "phi": 6.00327249, + "pt": 8.40910455 + }, + { + "rap": 3.10764407, + "phi": 0.5088296, + "pt": 8.30894232 + }, + { + "rap": -0.55073287, + "phi": 0.16601996, + "pt": 8.26803666 + }, + { + "rap": -3.15844579, + "phi": 1.36204941, + "pt": 7.5578621 + }, + { + "rap": -2.80210726, + "phi": 3.6249838, + "pt": 7.43724241 + }, + { + "rap": 1.08656054, + "phi": 3.89942435, + "pt": 6.84052394 + }, + { + "rap": 5.72935121, + "phi": 4.20697281, + "pt": 6.83235922 + }, + { + "rap": -1.33367469, + "phi": 4.23388764, + "pt": 6.43973876 + }, + { + "rap": -5.30510096, + "phi": 4.36541285, + "pt": 6.12065483 + }, + { + "rap": -2.56284802, + "phi": 5.73023999, + "pt": 6.0631507 + }, + { + "rap": -2.41082527, + "phi": 1.82827956, + "pt": 5.87571416 + }, + { + "rap": -2.96686986, + "phi": 5.22095929, + "pt": 5.57502596 + }, + { + "rap": -3.60553747, + "phi": 2.41909382, + "pt": 5.46858846 + }, + { + "rap": 1.45859935, + "phi": 1.34446168, + "pt": 5.41285118 + }, + { + "rap": 2.26697694, + "phi": 1.09667571, + "pt": 5.37432271 + }, + { + "rap": -5.45405941, + "phi": 1.72160264, + "pt": 5.351616 + }, + { + "rap": 2.08196204, + "phi": 0.26274384, + "pt": 5.32906741 + }, + { + "rap": 1.47389698, + "phi": 0.34920989, + "pt": 5.29002908 + }, + { + "rap": -0.67742893, + "phi": 3.44326635, + "pt": 5.26111873 + }, + { + "rap": -0.69086218, + "phi": 4.51629164, + "pt": 5.0218307 + } + ] + }, + { + "jetid": 33, + "jets": [ + { + "rap": 1.66430731, + "phi": 4.73928565, + "pt": 28.45117737 + }, + { + "rap": -1.43293955, + "phi": 2.08269443, + "pt": 8.11675385 + }, + { + "rap": -2.15252347, + "phi": 1.33059069, + "pt": 7.6583696 + }, + { + "rap": -4.76420487, + "phi": 1.31276342, + "pt": 6.9808209 + }, + { + "rap": 1.84865364, + "phi": 1.68299047, + "pt": 6.57733749 + }, + { + "rap": -0.64503515, + "phi": 1.0780381, + "pt": 6.55982626 + }, + { + "rap": -0.68736899, + "phi": 2.72466085, + "pt": 5.54875926 + }, + { + "rap": 0.74676526, + "phi": 3.17084634, + "pt": 5.4505977 + }, + { + "rap": -4.15381903, + "phi": 0.71616185, + "pt": 5.31770637 + } + ] + }, + { + "jetid": 34, + "jets": [ + { + "rap": 2.03490181, + "phi": 1.79491832, + "pt": 25.81534548 + }, + { + "rap": 1.83343754, + "phi": 0.93174776, + "pt": 14.47791791 + }, + { + "rap": -0.0253806, + "phi": 5.56790301, + "pt": 14.412503 + }, + { + "rap": -1.12717839, + "phi": 4.97443268, + "pt": 12.79744532 + }, + { + "rap": 1.29848795, + "phi": 4.81957233, + "pt": 10.90106866 + }, + { + "rap": 0.30170933, + "phi": 0.64065005, + "pt": 7.50952703 + }, + { + "rap": 2.32043336, + "phi": 2.79165492, + "pt": 6.78467806 + }, + { + "rap": -5.81090373, + "phi": 3.14512546, + "pt": 6.10710491 + }, + { + "rap": -1.97985578, + "phi": 5.80277102, + "pt": 5.85082422 + }, + { + "rap": -0.07495323, + "phi": 3.67585111, + "pt": 5.73727029 + } + ] + }, + { + "jetid": 35, + "jets": [ + { + "rap": -2.58569454, + "phi": 5.00156339, + "pt": 17.33592925 + }, + { + "rap": 3.76050564, + "phi": 0.87822525, + "pt": 14.95224798 + }, + { + "rap": -2.14356653, + "phi": 4.80050443, + "pt": 8.1924767 + } + ] + }, + { + "jetid": 36, + "jets": [ + { + "rap": -1.56994311, + "phi": 1.50843711, + "pt": 16.61790238 + }, + { + "rap": -1.33026985, + "phi": 3.88721542, + "pt": 12.92005267 + }, + { + "rap": -2.60764562, + "phi": 2.26532948, + "pt": 12.58575578 + }, + { + "rap": -2.77874325, + "phi": 0.16047235, + "pt": 8.35479442 + }, + { + "rap": -2.13765048, + "phi": 1.32008407, + "pt": 7.76053424 + }, + { + "rap": -3.00070299, + "phi": 3.33346135, + "pt": 6.288284 + }, + { + "rap": 5.13248497, + "phi": 5.29107121, + "pt": 5.47978873 + }, + { + "rap": -2.50000059, + "phi": 5.64188276, + "pt": 5.44314359 + }, + { + "rap": -1.75950115, + "phi": 0.11376169, + "pt": 5.42662751 + }, + { + "rap": -3.43659578, + "phi": 6.21215179, + "pt": 5.40521047 + }, + { + "rap": 3.97965504, + "phi": 4.07499249, + "pt": 5.25350165 + }, + { + "rap": -1.88529164, + "phi": 5.30610997, + "pt": 5.04156919 + }, + { + "rap": -0.87398583, + "phi": 1.4700037, + "pt": 5.00894952 + } + ] + }, + { + "jetid": 37, + "jets": [ + { + "rap": -0.13546515, + "phi": 1.9495252, + "pt": 22.83851067 + }, + { + "rap": -2.44352517, + "phi": 5.39619386, + "pt": 12.36617919 + }, + { + "rap": -1.54934592, + "phi": 5.26265593, + "pt": 7.4803184 + }, + { + "rap": -1.04517734, + "phi": 4.54464209, + "pt": 6.66659542 + }, + { + "rap": -2.40934265, + "phi": 5.9105635, + "pt": 5.3576065 + } + ] + }, + { + "jetid": 38, + "jets": [ + { + "rap": 4.74673296, + "phi": 3.66890201, + "pt": 18.59047591 + }, + { + "rap": 4.974745, + "phi": 0.95283077, + "pt": 12.18384151 + }, + { + "rap": -0.73703643, + "phi": 1.43499018, + "pt": 8.68146797 + }, + { + "rap": 2.78351921, + "phi": 6.01644277, + "pt": 8.47890531 + } + ] + }, + { + "jetid": 39, + "jets": [ + { + "rap": 0.89828866, + "phi": 3.04658767, + "pt": 20.98614198 + }, + { + "rap": -2.44777364, + "phi": 2.45473259, + "pt": 13.2147192 + }, + { + "rap": 2.07097153, + "phi": 5.06481913, + "pt": 13.11395901 + }, + { + "rap": -1.76348252, + "phi": 6.16215674, + "pt": 11.5058729 + }, + { + "rap": -1.65547441, + "phi": 0.33379478, + "pt": 10.72289175 + }, + { + "rap": -0.44288895, + "phi": 4.7781182, + "pt": 8.01711185 + }, + { + "rap": 1.11292115, + "phi": 4.58679437, + "pt": 6.92874226 + }, + { + "rap": -0.12507023, + "phi": 2.39706328, + "pt": 5.99806407 + }, + { + "rap": -2.56055249, + "phi": 0.62456236, + "pt": 5.91700331 + }, + { + "rap": -2.0535143, + "phi": 2.05055505, + "pt": 5.57190138 + }, + { + "rap": 1.36287543, + "phi": 3.42672344, + "pt": 5.27656043 + } + ] + }, + { + "jetid": 40, + "jets": [ + { + "rap": 3.62108225, + "phi": 0.71965088, + "pt": 31.06375699 + }, + { + "rap": -1.544335, + "phi": 4.12685476, + "pt": 7.39951742 + }, + { + "rap": 0.76608775, + "phi": 4.25002184, + "pt": 5.96919204 + } + ] + }, + { + "jetid": 41, + "jets": [ + { + "rap": -4.99421586, + "phi": 5.98255257, + "pt": 18.97582593 + }, + { + "rap": -2.21873405, + "phi": 2.90565437, + "pt": 14.9775239 + }, + { + "rap": 1.57845244, + "phi": 3.41295713, + "pt": 10.33563658 + }, + { + "rap": -1.9442346, + "phi": 2.48571061, + "pt": 6.77981942 + }, + { + "rap": -2.43986922, + "phi": 2.54053038, + "pt": 6.06103345 + }, + { + "rap": -2.42548297, + "phi": 3.45094499, + "pt": 5.36418537 + }, + { + "rap": 5.65189522, + "phi": 4.69202859, + "pt": 5.03478313 + } + ] + }, + { + "jetid": 42, + "jets": [ + { + "rap": 1.80273941, + "phi": 5.17091788, + "pt": 24.04791266 + }, + { + "rap": 3.58446644, + "phi": 1.8132847, + "pt": 22.46519562 + }, + { + "rap": 2.70064486, + "phi": 2.61543446, + "pt": 9.41353751 + }, + { + "rap": 4.61428882, + "phi": 5.18670875, + "pt": 6.35704344 + }, + { + "rap": -3.11937292, + "phi": 3.03495483, + "pt": 5.10336593 + } + ] + }, + { + "jetid": 43, + "jets": [ + { + "rap": 5.17500976, + "phi": 5.4989582, + "pt": 28.84646362 + }, + { + "rap": 5.08833118, + "phi": 2.19851288, + "pt": 20.21301587 + } + ] + }, + { + "jetid": 44, + "jets": [ + { + "rap": -3.49016614, + "phi": 2.92534278, + "pt": 31.43910703 + }, + { + "rap": -2.11619961, + "phi": 0.20000564, + "pt": 22.78395715 + }, + { + "rap": -2.4511975, + "phi": 0.72249281, + "pt": 16.5512101 + }, + { + "rap": 2.99696345, + "phi": 5.22849495, + "pt": 14.51002376 + }, + { + "rap": -4.50092665, + "phi": 2.90728866, + "pt": 14.16174111 + }, + { + "rap": 3.04861805, + "phi": 1.31769873, + "pt": 10.41417392 + }, + { + "rap": 0.26510373, + "phi": 3.31777131, + "pt": 9.23363224 + }, + { + "rap": -4.82729828, + "phi": 0.92109355, + "pt": 9.17640949 + }, + { + "rap": -2.11243445, + "phi": 5.98691646, + "pt": 7.45151133 + }, + { + "rap": 2.62895802, + "phi": 4.22858422, + "pt": 7.35231344 + }, + { + "rap": -2.86603696, + "phi": 0.23295109, + "pt": 7.15102193 + }, + { + "rap": 0.74464063, + "phi": 2.93338285, + "pt": 5.90204627 + }, + { + "rap": -0.01932336, + "phi": 4.05045303, + "pt": 5.44207316 + }, + { + "rap": -2.30659038, + "phi": 1.53909373, + "pt": 5.26975889 + }, + { + "rap": 0.6534349, + "phi": 3.73533287, + "pt": 5.19233188 + }, + { + "rap": 2.44753885, + "phi": 1.76426394, + "pt": 5.06300207 + } + ] + }, + { + "jetid": 45, + "jets": [ + { + "rap": -0.84450174, + "phi": 3.10497745, + "pt": 38.21407069 + }, + { + "rap": 2.66564854, + "phi": 5.89703968, + "pt": 16.91785372 + }, + { + "rap": -3.67354145, + "phi": 0.30154748, + "pt": 10.05339211 + }, + { + "rap": 1.8981449, + "phi": 2.88540114, + "pt": 9.38748353 + }, + { + "rap": -4.19774211, + "phi": 1.00894894, + "pt": 8.61695014 + }, + { + "rap": 2.19393679, + "phi": 2.07602432, + "pt": 6.90812695 + }, + { + "rap": -4.06476757, + "phi": 4.67437213, + "pt": 6.69618445 + }, + { + "rap": 0.27940337, + "phi": 0.34196289, + "pt": 6.48333129 + }, + { + "rap": -0.59869769, + "phi": 0.5051657, + "pt": 6.3857297 + }, + { + "rap": 1.66809815, + "phi": 0.93301802, + "pt": 5.61871696 + } + ] + }, + { + "jetid": 46, + "jets": [ + { + "rap": -0.29174278, + "phi": 5.15233629, + "pt": 21.69417966 + }, + { + "rap": 0.27562083, + "phi": 5.91538083, + "pt": 12.70145295 + }, + { + "rap": 2.91217184, + "phi": 3.44082046, + "pt": 10.64502176 + }, + { + "rap": -0.11753471, + "phi": 2.5368116, + "pt": 10.49085215 + }, + { + "rap": 2.65477053, + "phi": 1.10724752, + "pt": 9.32132341 + }, + { + "rap": 2.94072156, + "phi": 1.95507108, + "pt": 7.90848016 + }, + { + "rap": -0.14195953, + "phi": 1.85948773, + "pt": 7.8987001 + }, + { + "rap": -0.41783496, + "phi": 0.05750623, + "pt": 6.40065666 + }, + { + "rap": 1.93554243, + "phi": 3.04287668, + "pt": 5.79829227 + }, + { + "rap": 1.56833451, + "phi": 4.61760168, + "pt": 5.40264161 + }, + { + "rap": -1.85454097, + "phi": 5.04060847, + "pt": 5.30128847 + }, + { + "rap": 1.39236023, + "phi": 3.025858, + "pt": 5.27617843 + }, + { + "rap": 3.55086487, + "phi": 0.60512744, + "pt": 5.1791226 + } + ] + }, + { + "jetid": 47, + "jets": [ + { + "rap": -1.11922336, + "phi": 3.82634863, + "pt": 19.03560955 + }, + { + "rap": -0.21128955, + "phi": 0.88073639, + "pt": 10.79803033 + } + ] + }, + { + "jetid": 48, + "jets": [ + { + "rap": 0.44869347, + "phi": 2.64258807, + "pt": 16.98495478 + }, + { + "rap": 3.03896878, + "phi": 2.75147085, + "pt": 13.92362033 + }, + { + "rap": 1.70591755, + "phi": 0.22958573, + "pt": 10.06987181 + }, + { + "rap": -2.62528738, + "phi": 1.42650912, + "pt": 7.54969312 + }, + { + "rap": -1.59336587, + "phi": 5.59339413, + "pt": 7.47863427 + }, + { + "rap": -0.11627512, + "phi": 1.2292313, + "pt": 7.4633832 + }, + { + "rap": 3.26115003, + "phi": 2.23517513, + "pt": 6.88838611 + }, + { + "rap": -1.11727123, + "phi": 4.19924981, + "pt": 6.87532856 + }, + { + "rap": 0.35036854, + "phi": 0.40601207, + "pt": 5.71034299 + }, + { + "rap": 1.65353783, + "phi": 4.77388987, + "pt": 5.34748487 + } + ] + }, + { + "jetid": 49, + "jets": [ + { + "rap": -0.98484561, + "phi": 1.33584918, + "pt": 51.94172514 + }, + { + "rap": -0.84148244, + "phi": 4.84139855, + "pt": 49.5963835 + }, + { + "rap": -1.13156109, + "phi": 4.26418938, + "pt": 9.26656312 + } + ] + }, + { + "jetid": 50, + "jets": [ + { + "rap": 2.85211211, + "phi": 3.07527979, + "pt": 17.0668719 + }, + { + "rap": 3.38220423, + "phi": 0.18039629, + "pt": 16.72930797 + }, + { + "rap": 2.45452399, + "phi": 3.90607493, + "pt": 5.70676246 + }, + { + "rap": 4.06284835, + "phi": 5.83161284, + "pt": 5.43751824 + } + ] + }, + { + "jetid": 51, + "jets": [ + { + "rap": -4.12071362, + "phi": 3.37711062, + "pt": 22.68186246 + }, + { + "rap": -4.09845373, + "phi": 1.14331499, + "pt": 18.12353909 + }, + { + "rap": -3.51526931, + "phi": 5.52503261, + "pt": 10.01720543 + }, + { + "rap": 0.93968293, + "phi": 5.4438079, + "pt": 7.64128356 + }, + { + "rap": 1.57591912, + "phi": 2.56640306, + "pt": 7.03327944 + }, + { + "rap": 3.36011947, + "phi": 0.18570879, + "pt": 6.56795322 + }, + { + "rap": -3.10095974, + "phi": 1.12785355, + "pt": 5.38942046 + }, + { + "rap": 5.21531101, + "phi": 3.97884661, + "pt": 5.31908901 + }, + { + "rap": -3.89514003, + "phi": 5.15025731, + "pt": 5.23402934 + } + ] + }, + { + "jetid": 52, + "jets": [ + { + "rap": 0.52990759, + "phi": 0.76543626, + "pt": 27.54351413 + }, + { + "rap": -1.7416513, + "phi": 4.04040869, + "pt": 22.05839225 + }, + { + "rap": -3.0784054, + "phi": 4.39735382, + "pt": 9.93650986 + }, + { + "rap": -0.04527219, + "phi": 0.92620056, + "pt": 6.92486187 + }, + { + "rap": -0.62322574, + "phi": 1.97930699, + "pt": 6.70572406 + }, + { + "rap": 0.10821429, + "phi": 1.6977548, + "pt": 6.66116083 + }, + { + "rap": 0.60099739, + "phi": 1.41967361, + "pt": 5.62495498 + } + ] + }, + { + "jetid": 53, + "jets": [ + { + "rap": -0.48703241, + "phi": 0.68391087, + "pt": 28.75099429 + }, + { + "rap": 3.14189985, + "phi": 4.56042508, + "pt": 13.32810081 + }, + { + "rap": 5.27680869, + "phi": 1.98776389, + "pt": 11.45905881 + }, + { + "rap": 2.75620654, + "phi": 3.29882401, + "pt": 10.35585886 + }, + { + "rap": 2.37163836, + "phi": 0.81861679, + "pt": 10.11293059 + }, + { + "rap": 1.92887993, + "phi": 4.18973327, + "pt": 8.4875459 + }, + { + "rap": -1.27043261, + "phi": 3.85176272, + "pt": 7.88154462 + }, + { + "rap": 0.7350985, + "phi": 5.98359276, + "pt": 6.74073373 + }, + { + "rap": 2.16750575, + "phi": 4.89634451, + "pt": 6.48924176 + }, + { + "rap": 0.11942529, + "phi": 0.16603529, + "pt": 6.1762165 + }, + { + "rap": -2.88940129, + "phi": 1.60248258, + "pt": 5.55701166 + }, + { + "rap": -1.22587675, + "phi": 4.2868427, + "pt": 5.52679742 + } + ] + }, + { + "jetid": 54, + "jets": [ + { + "rap": 0.03877767, + "phi": 5.05389401, + "pt": 20.8196713 + }, + { + "rap": -1.25422288, + "phi": 2.78801414, + "pt": 12.70429969 + }, + { + "rap": -4.51970076, + "phi": 1.42630985, + "pt": 5.39521276 + } + ] + }, + { + "jetid": 55, + "jets": [ + { + "rap": 2.83546981, + "phi": 2.58969232, + "pt": 20.36715562 + }, + { + "rap": 2.76778413, + "phi": 5.39202269, + "pt": 17.77371977 + }, + { + "rap": 2.25606961, + "phi": 0.55785543, + "pt": 5.45395047 + } + ] + }, + { + "jetid": 56, + "jets": [ + { + "rap": -1.6058401, + "phi": 5.32486983, + "pt": 22.783447 + }, + { + "rap": 1.90188526, + "phi": 2.19497599, + "pt": 21.69123455 + }, + { + "rap": -0.34655377, + "phi": 0.35218659, + "pt": 10.73899398 + }, + { + "rap": -0.58167905, + "phi": 5.23035221, + "pt": 10.50026601 + } + ] + }, + { + "jetid": 57, + "jets": [ + { + "rap": -2.07054339, + "phi": 6.0073307, + "pt": 24.35334165 + }, + { + "rap": -1.69562613, + "phi": 2.73898951, + "pt": 11.57144793 + }, + { + "rap": -1.22861657, + "phi": 2.31642363, + "pt": 6.97839143 + } + ] + }, + { + "jetid": 58, + "jets": [ + { + "rap": -1.45886822, + "phi": 0.10870344, + "pt": 22.82473754 + }, + { + "rap": -1.16573181, + "phi": 2.93725253, + "pt": 15.22924094 + }, + { + "rap": -3.02872857, + "phi": 3.73182454, + "pt": 11.26602029 + } + ] + }, + { + "jetid": 59, + "jets": [ + { + "rap": -1.44621807, + "phi": 3.88795181, + "pt": 23.97727633 + }, + { + "rap": 0.92442504, + "phi": 0.6632397, + "pt": 19.58838573 + }, + { + "rap": -0.27290976, + "phi": 1.89676568, + "pt": 8.82229736 + }, + { + "rap": 2.43562458, + "phi": 1.3573952, + "pt": 7.24928406 + }, + { + "rap": 3.04475866, + "phi": 2.96505854, + "pt": 6.90598212 + }, + { + "rap": -0.92902223, + "phi": 5.11249516, + "pt": 5.8038351 + }, + { + "rap": -4.67171076, + "phi": 1.93356567, + "pt": 5.42155748 + }, + { + "rap": -2.05798142, + "phi": 3.93996048, + "pt": 5.30540181 + } + ] + }, + { + "jetid": 60, + "jets": [ + { + "rap": 0.60973356, + "phi": 4.1837463, + "pt": 31.73637158 + }, + { + "rap": 0.85096608, + "phi": 0.91159551, + "pt": 27.23080175 + } + ] + }, + { + "jetid": 61, + "jets": [ + { + "rap": 1.12893587, + "phi": 1.23686802, + "pt": 29.95013126 + }, + { + "rap": 2.80411458, + "phi": 3.55243048, + "pt": 12.56115138 + }, + { + "rap": 3.72791224, + "phi": 3.80188028, + "pt": 7.9250092 + }, + { + "rap": 5.20958376, + "phi": 6.19762015, + "pt": 7.61934619 + }, + { + "rap": -1.0131503, + "phi": 5.15369363, + "pt": 6.5810922 + }, + { + "rap": 2.59419556, + "phi": 4.56102655, + "pt": 6.43316224 + } + ] + }, + { + "jetid": 62, + "jets": [ + { + "rap": -1.51444928, + "phi": 1.51555336, + "pt": 17.53638083 + }, + { + "rap": -1.46774834, + "phi": 0.79804757, + "pt": 16.02325604 + }, + { + "rap": -3.6691301, + "phi": 5.09792877, + "pt": 12.36910821 + }, + { + "rap": -2.91347442, + "phi": 3.03035573, + "pt": 10.22547952 + }, + { + "rap": -0.9130224, + "phi": 4.24943498, + "pt": 9.15051926 + }, + { + "rap": -5.07686284, + "phi": 4.75367654, + "pt": 7.59710373 + }, + { + "rap": 2.53714375, + "phi": 0.21841523, + "pt": 6.8532384 + }, + { + "rap": -1.6705541, + "phi": 5.15177251, + "pt": 6.30624737 + }, + { + "rap": 4.92550835, + "phi": 0.92484368, + "pt": 5.67212552 + }, + { + "rap": 2.64814149, + "phi": 5.86822711, + "pt": 5.46404386 + }, + { + "rap": 1.13371152, + "phi": 0.40777635, + "pt": 5.3163966 + }, + { + "rap": -1.3035186, + "phi": 0.15957488, + "pt": 5.24787737 + }, + { + "rap": -5.62157201, + "phi": 2.2793931, + "pt": 5.14671486 + } + ] + }, + { + "jetid": 63, + "jets": [ + { + "rap": -0.61464715, + "phi": 5.98179196, + "pt": 27.63800158 + }, + { + "rap": -2.3299799, + "phi": 2.76536262, + "pt": 25.64663752 + }, + { + "rap": -0.48412237, + "phi": 2.40874591, + "pt": 9.72312331 + } + ] + }, + { + "jetid": 64, + "jets": [ + { + "rap": -4.31620325, + "phi": 2.57045971, + "pt": 18.33202094 + }, + { + "rap": 3.55720812, + "phi": 0.4536455, + "pt": 12.90612095 + }, + { + "rap": -0.01713831, + "phi": 5.53805367, + "pt": 10.16764206 + }, + { + "rap": 1.30922987, + "phi": 5.89654538, + "pt": 6.49381399 + }, + { + "rap": 0.14841271, + "phi": 2.29816441, + "pt": 6.16124857 + }, + { + "rap": -0.6545446, + "phi": 5.09682118, + "pt": 5.00787838 + } + ] + }, + { + "jetid": 65, + "jets": [ + { + "rap": 1.61437102, + "phi": 0.10564807, + "pt": 26.60798494 + }, + { + "rap": 1.01051932, + "phi": 3.37907586, + "pt": 15.85780664 + }, + { + "rap": 1.00983095, + "phi": 4.26813829, + "pt": 10.56302789 + }, + { + "rap": -2.77726056, + "phi": 3.37968796, + "pt": 6.34981022 + }, + { + "rap": 2.92617753, + "phi": 1.31364452, + "pt": 5.77478261 + }, + { + "rap": 0.02272524, + "phi": 5.05224034, + "pt": 5.75346713 + } + ] + }, + { + "jetid": 66, + "jets": [ + { + "rap": -2.13754906, + "phi": 5.94131867, + "pt": 18.11350687 + }, + { + "rap": -0.53941413, + "phi": 3.1595856, + "pt": 16.75609778 + }, + { + "rap": 4.17649134, + "phi": 1.75155139, + "pt": 8.56033487 + }, + { + "rap": -0.76331085, + "phi": 3.54843894, + "pt": 6.21454458 + }, + { + "rap": 2.98630952, + "phi": 3.21316895, + "pt": 6.11950658 + }, + { + "rap": 2.40264638, + "phi": 5.09090004, + "pt": 5.24280907 + }, + { + "rap": 3.04198687, + "phi": 0.63230325, + "pt": 5.14451928 + }, + { + "rap": 1.73621211, + "phi": 2.96956367, + "pt": 5.14231997 + } + ] + }, + { + "jetid": 67, + "jets": [ + { + "rap": 4.16657983, + "phi": 2.35010871, + "pt": 25.81761573 + }, + { + "rap": 1.08243428, + "phi": 5.56285266, + "pt": 9.74933685 + }, + { + "rap": -3.67724082, + "phi": 4.42326616, + "pt": 9.39552291 + }, + { + "rap": -0.70387734, + "phi": 6.27976207, + "pt": 6.34978231 + }, + { + "rap": 3.24136089, + "phi": 4.55592752, + "pt": 5.99395028 + }, + { + "rap": -0.27026742, + "phi": 4.10675774, + "pt": 5.83983568 + }, + { + "rap": -1.08867464, + "phi": 2.2761876, + "pt": 5.73209852 + }, + { + "rap": 6.28215524, + "phi": 5.73702121, + "pt": 5.16282579 + }, + { + "rap": -0.78299542, + "phi": 4.35831353, + "pt": 5.14867458 + }, + { + "rap": 0.56922577, + "phi": 5.02463516, + "pt": 5.05592132 + } + ] + }, + { + "jetid": 68, + "jets": [ + { + "rap": -0.15066877, + "phi": 0.20684908, + "pt": 23.26776193 + }, + { + "rap": -1.75497233, + "phi": 3.18831426, + "pt": 21.51822056 + }, + { + "rap": -1.76250848, + "phi": 3.75675716, + "pt": 8.83835089 + } + ] + }, + { + "jetid": 69, + "jets": [ + { + "rap": 1.86172798, + "phi": 5.54471468, + "pt": 24.1878514 + }, + { + "rap": -1.66270525, + "phi": 3.54129451, + "pt": 14.37826989 + }, + { + "rap": -1.71293395, + "phi": 3.03443062, + "pt": 13.01315446 + }, + { + "rap": -0.12465186, + "phi": 1.51460706, + "pt": 10.16346863 + }, + { + "rap": -0.9892985, + "phi": 2.95210397, + "pt": 10.06382264 + }, + { + "rap": -1.33937668, + "phi": 0.45053166, + "pt": 9.1616915 + }, + { + "rap": -4.65020552, + "phi": 1.02620366, + "pt": 8.30518489 + }, + { + "rap": 0.12186764, + "phi": 0.23368328, + "pt": 7.53058488 + }, + { + "rap": 4.70461528, + "phi": 3.19244572, + "pt": 6.91127975 + }, + { + "rap": 2.77370398, + "phi": 2.08198679, + "pt": 6.08138833 + }, + { + "rap": 4.11927396, + "phi": 6.00822648, + "pt": 5.87100231 + }, + { + "rap": 5.72022883, + "phi": 0.20790395, + "pt": 5.81541639 + }, + { + "rap": -1.50977948, + "phi": 5.24268679, + "pt": 5.65814598 + }, + { + "rap": 3.26392468, + "phi": 4.53904432, + "pt": 5.51322879 + }, + { + "rap": 5.51802476, + "phi": 2.86522774, + "pt": 5.41039087 + }, + { + "rap": 1.83785959, + "phi": 3.88978008, + "pt": 5.27407913 + }, + { + "rap": -4.88776399, + "phi": 5.42409357, + "pt": 5.16161046 + }, + { + "rap": -2.75901021, + "phi": 0.59786967, + "pt": 5.13029712 + }, + { + "rap": -4.25709126, + "phi": 2.86765818, + "pt": 5.04019773 + } + ] + }, + { + "jetid": 70, + "jets": [ + { + "rap": -0.46982661, + "phi": 4.28531916, + "pt": 22.8991046 + }, + { + "rap": -1.2592364, + "phi": 1.32114744, + "pt": 9.77912996 + }, + { + "rap": -3.57004522, + "phi": 1.19476003, + "pt": 5.37080453 + } + ] + }, + { + "jetid": 71, + "jets": [ + { + "rap": 2.05981358, + "phi": 4.50569096, + "pt": 21.35134558 + }, + { + "rap": 1.66211429, + "phi": 2.44810138, + "pt": 15.46021883 + }, + { + "rap": 1.67717561, + "phi": 1.46451984, + "pt": 15.41910803 + }, + { + "rap": 3.60425352, + "phi": 5.17877356, + "pt": 6.97183596 + }, + { + "rap": 4.36606787, + "phi": 6.1630021, + "pt": 5.19176911 + } + ] + }, + { + "jetid": 72, + "jets": [ + { + "rap": -1.75442841, + "phi": 5.82592889, + "pt": 24.05771447 + }, + { + "rap": -1.00062383, + "phi": 2.85023607, + "pt": 12.78380297 + }, + { + "rap": -0.53248434, + "phi": 1.65729956, + "pt": 9.69695855 + }, + { + "rap": 3.79593172, + "phi": 2.82545105, + "pt": 5.0633652 + } + ] + }, + { + "jetid": 73, + "jets": [ + { + "rap": 0.12298645, + "phi": 0.47014221, + "pt": 19.29034251 + }, + { + "rap": 0.73376168, + "phi": 3.39250374, + "pt": 9.35662502 + }, + { + "rap": 0.1156221, + "phi": 3.87302638, + "pt": 6.34984744 + }, + { + "rap": -2.79168561, + "phi": 4.26718772, + "pt": 5.26317185 + } + ] + }, + { + "jetid": 74, + "jets": [ + { + "rap": -2.7719803, + "phi": 2.57018006, + "pt": 32.31670914 + }, + { + "rap": -2.51402731, + "phi": 5.25640537, + "pt": 26.20629216 + }, + { + "rap": -3.11057576, + "phi": 0.05244812, + "pt": 5.51161601 + } + ] + }, + { + "jetid": 75, + "jets": [ + { + "rap": 1.56505609, + "phi": 2.99472397, + "pt": 36.11547782 + }, + { + "rap": -1.55894028, + "phi": 0.66055711, + "pt": 10.69433621 + }, + { + "rap": -5.00652819, + "phi": 5.12290183, + "pt": 10.49363539 + }, + { + "rap": -0.26724342, + "phi": 0.37293061, + "pt": 8.63622835 + }, + { + "rap": 1.57513631, + "phi": 5.42864926, + "pt": 8.08114735 + }, + { + "rap": -0.36506162, + "phi": 5.99048358, + "pt": 6.71897449 + }, + { + "rap": -4.45445858, + "phi": 3.5454506, + "pt": 6.57009869 + }, + { + "rap": -0.74610756, + "phi": 2.76300852, + "pt": 6.34934044 + }, + { + "rap": 2.20305295, + "phi": 1.4553232, + "pt": 5.60816092 + }, + { + "rap": -1.3250387, + "phi": 1.74093731, + "pt": 5.50127828 + }, + { + "rap": 1.16291014, + "phi": 3.24139039, + "pt": 5.35970296 + }, + { + "rap": 1.06563577, + "phi": 5.76661433, + "pt": 5.30660593 + }, + { + "rap": 2.61084349, + "phi": 3.83911627, + "pt": 5.29850845 + }, + { + "rap": -1.43913319, + "phi": 3.55894277, + "pt": 5.20666209 + } + ] + }, + { + "jetid": 76, + "jets": [ + { + "rap": -0.69493967, + "phi": 4.17120081, + "pt": 25.05346336 + }, + { + "rap": 2.12366137, + "phi": 0.43746161, + "pt": 15.03475605 + }, + { + "rap": -1.74693558, + "phi": 1.10194692, + "pt": 11.43470599 + }, + { + "rap": -0.94644565, + "phi": 3.61461305, + "pt": 11.20522251 + } + ] + }, + { + "jetid": 77, + "jets": [ + { + "rap": -1.14500448, + "phi": 5.06562644, + "pt": 16.24723524 + }, + { + "rap": -3.07222809, + "phi": 1.99451543, + "pt": 16.08875886 + }, + { + "rap": -1.7275963, + "phi": 3.81811752, + "pt": 7.86219412 + }, + { + "rap": -2.2418064, + "phi": 1.09344353, + "pt": 5.69293328 + }, + { + "rap": 1.50231288, + "phi": 0.47059978, + "pt": 5.6843573 + }, + { + "rap": -4.96836663, + "phi": 2.2606741, + "pt": 5.00910919 + } + ] + }, + { + "jetid": 78, + "jets": [ + { + "rap": -4.78590484, + "phi": 4.98634272, + "pt": 27.53362053 + }, + { + "rap": -0.52549806, + "phi": 1.2446939, + "pt": 18.77306536 + }, + { + "rap": -1.75878681, + "phi": 2.03380469, + "pt": 10.39355656 + }, + { + "rap": -1.10189738, + "phi": 1.09779514, + "pt": 9.28656838 + }, + { + "rap": 0.3921684, + "phi": 2.06651499, + "pt": 8.65028521 + }, + { + "rap": -3.47684571, + "phi": 2.59563701, + "pt": 8.36642243 + }, + { + "rap": -3.43271212, + "phi": 5.84021085, + "pt": 8.01097524 + }, + { + "rap": -2.86897555, + "phi": 2.98567191, + "pt": 7.52576784 + }, + { + "rap": -3.08890434, + "phi": 2.18665567, + "pt": 7.31410368 + }, + { + "rap": 0.6081653, + "phi": 4.90120462, + "pt": 6.81568441 + }, + { + "rap": 1.87012446, + "phi": 5.18081774, + "pt": 6.59760342 + }, + { + "rap": -0.96931738, + "phi": 4.18200653, + "pt": 6.34013752 + }, + { + "rap": 1.05225939, + "phi": 2.37440966, + "pt": 6.21824529 + }, + { + "rap": -4.77888843, + "phi": 5.62077527, + "pt": 5.13606472 + }, + { + "rap": -2.64541133, + "phi": 1.80766928, + "pt": 5.07537871 + } + ] + }, + { + "jetid": 79, + "jets": [ + { + "rap": -0.05815909, + "phi": 0.24031823, + "pt": 34.26353952 + }, + { + "rap": -1.80255713, + "phi": 3.56179851, + "pt": 32.27306895 + }, + { + "rap": 1.27406208, + "phi": 0.01087646, + "pt": 6.82996946 + }, + { + "rap": 1.94445493, + "phi": 3.62130484, + "pt": 5.44294098 + } + ] + }, + { + "jetid": 80, + "jets": [ + { + "rap": -0.05084612, + "phi": 5.45294447, + "pt": 24.35326702 + }, + { + "rap": -0.54278545, + "phi": 1.60488832, + "pt": 10.7348801 + }, + { + "rap": -1.60790614, + "phi": 2.16429529, + "pt": 8.56290251 + }, + { + "rap": 1.94441302, + "phi": 3.59553891, + "pt": 6.43420735 + }, + { + "rap": -2.15928686, + "phi": 3.8030523, + "pt": 6.07491351 + }, + { + "rap": -0.81810431, + "phi": 0.69455898, + "pt": 5.09166486 + }, + { + "rap": -3.59993722, + "phi": 0.80364981, + "pt": 5.01127274 + } + ] + }, + { + "jetid": 81, + "jets": [ + { + "rap": -2.61166492, + "phi": 5.47221968, + "pt": 22.46480128 + }, + { + "rap": -2.79984666, + "phi": 2.4424784, + "pt": 20.29775365 + }, + { + "rap": 4.41258907, + "phi": 5.56391204, + "pt": 5.27698362 + } + ] + }, + { + "jetid": 82, + "jets": [ + { + "rap": -0.57475777, + "phi": 2.17553175, + "pt": 54.39512081 + }, + { + "rap": -1.8638639, + "phi": 4.68857483, + "pt": 21.82765901 + }, + { + "rap": 0.81346465, + "phi": 5.66842022, + "pt": 10.73023433 + }, + { + "rap": -0.3150661, + "phi": 5.57224692, + "pt": 7.87824805 + }, + { + "rap": -0.71787753, + "phi": 5.64702622, + "pt": 5.30875166 + } + ] + }, + { + "jetid": 83, + "jets": [ + { + "rap": -1.35021236, + "phi": 0.86631546, + "pt": 26.36382179 + }, + { + "rap": 1.67233389, + "phi": 2.6715205, + "pt": 23.5886375 + }, + { + "rap": 0.32089721, + "phi": 1.76636413, + "pt": 13.32835644 + }, + { + "rap": -4.38355951, + "phi": 4.62734557, + "pt": 11.92584968 + }, + { + "rap": -2.54480934, + "phi": 4.52108408, + "pt": 10.15575033 + }, + { + "rap": 0.30085617, + "phi": 5.22900527, + "pt": 7.02366103 + }, + { + "rap": -5.59372071, + "phi": 4.88864955, + "pt": 5.29384108 + }, + { + "rap": 1.28178311, + "phi": 2.2460548, + "pt": 5.0765227 + } + ] + }, + { + "jetid": 84, + "jets": [ + { + "rap": -0.62440189, + "phi": 0.92978201, + "pt": 20.14416569 + }, + { + "rap": -2.32019008, + "phi": 4.78634424, + "pt": 17.54422847 + }, + { + "rap": -2.73936612, + "phi": 1.66264792, + "pt": 12.42145881 + }, + { + "rap": -1.73264384, + "phi": 3.52839467, + "pt": 7.57641221 + } + ] + }, + { + "jetid": 85, + "jets": [ + { + "rap": -1.65598606, + "phi": 6.24402928, + "pt": 19.51729754 + }, + { + "rap": -1.53649648, + "phi": 3.32006967, + "pt": 12.9710666 + }, + { + "rap": -2.63642875, + "phi": 3.82734463, + "pt": 10.10859142 + }, + { + "rap": 1.35395698, + "phi": 3.7607081, + "pt": 8.51142227 + }, + { + "rap": -5.66459093, + "phi": 1.41705621, + "pt": 7.08248093 + }, + { + "rap": -2.02375507, + "phi": 5.7686963, + "pt": 5.49190704 + }, + { + "rap": 0.11353261, + "phi": 0.23936849, + "pt": 5.42017905 + } + ] + }, + { + "jetid": 86, + "jets": [ + { + "rap": -2.93277518, + "phi": 4.51589997, + "pt": 46.44394715 + }, + { + "rap": -2.90379873, + "phi": 1.25551491, + "pt": 39.48836268 + }, + { + "rap": -2.42406577, + "phi": 4.55808281, + "pt": 7.88615876 + }, + { + "rap": 0.57013796, + "phi": 1.61062353, + "pt": 6.50429323 + }, + { + "rap": 0.66927196, + "phi": 2.90675865, + "pt": 6.39523436 + }, + { + "rap": 4.21377218, + "phi": 5.44338945, + "pt": 6.34941646 + }, + { + "rap": 1.71621707, + "phi": 3.31372713, + "pt": 5.70973238 + } + ] + }, + { + "jetid": 87, + "jets": [ + { + "rap": 2.00910066, + "phi": 0.32757274, + "pt": 41.47723011 + }, + { + "rap": -0.9772689, + "phi": 3.42987667, + "pt": 35.19997036 + }, + { + "rap": -1.5706189, + "phi": 2.25231286, + "pt": 5.53356679 + }, + { + "rap": -1.46964894, + "phi": 3.40368977, + "pt": 5.18413927 + } + ] + }, + { + "jetid": 88, + "jets": [ + { + "rap": -1.20416812, + "phi": 3.97249666, + "pt": 34.30385457 + }, + { + "rap": 1.57909705, + "phi": 6.02646529, + "pt": 19.35396758 + }, + { + "rap": 0.89520971, + "phi": 1.56955899, + "pt": 15.78294139 + }, + { + "rap": -2.09038742, + "phi": 1.10082856, + "pt": 11.82670463 + }, + { + "rap": -2.6498172, + "phi": 4.36218595, + "pt": 9.33467874 + }, + { + "rap": 3.56951709, + "phi": 0.53493114, + "pt": 7.24366614 + }, + { + "rap": 2.03858231, + "phi": 3.19512514, + "pt": 6.60619622 + }, + { + "rap": 0.76248185, + "phi": 3.23556127, + "pt": 6.33454579 + }, + { + "rap": -1.50732388, + "phi": 1.05719491, + "pt": 5.94109734 + }, + { + "rap": -5.77981264, + "phi": 6.0174733, + "pt": 5.58054981 + }, + { + "rap": -0.96696939, + "phi": 2.61335883, + "pt": 5.37769821 + }, + { + "rap": -0.12050175, + "phi": 1.64224527, + "pt": 5.21265076 + }, + { + "rap": 0.61177165, + "phi": 2.41986836, + "pt": 5.13349928 + } + ] + }, + { + "jetid": 89, + "jets": [ + { + "rap": 0.6099434, + "phi": 5.84376337, + "pt": 32.47571275 + }, + { + "rap": 1.49897356, + "phi": 2.62542888, + "pt": 24.01653258 + }, + { + "rap": 0.46682098, + "phi": 2.09664603, + "pt": 22.4674145 + }, + { + "rap": 1.51144527, + "phi": 0.55305408, + "pt": 13.37244758 + }, + { + "rap": 0.36514927, + "phi": 4.62244346, + "pt": 8.45211163 + }, + { + "rap": -0.21302278, + "phi": 5.38661058, + "pt": 8.24594213 + }, + { + "rap": -2.25190902, + "phi": 3.19365785, + "pt": 7.96326299 + }, + { + "rap": -2.58688438, + "phi": 1.24956361, + "pt": 5.86690374 + }, + { + "rap": -1.82844778, + "phi": 4.62151586, + "pt": 5.67094599 + }, + { + "rap": 1.94386946, + "phi": 0.16324486, + "pt": 5.64525576 + }, + { + "rap": 0.90656121, + "phi": 4.80948952, + "pt": 5.04133688 + } + ] + }, + { + "jetid": 90, + "jets": [ + { + "rap": 1.65652237, + "phi": 5.62175264, + "pt": 35.61534779 + }, + { + "rap": -0.28590594, + "phi": 2.62804705, + "pt": 25.03972348 + }, + { + "rap": 0.95165092, + "phi": 0.35402968, + "pt": 9.53168081 + }, + { + "rap": 1.35959312, + "phi": 0.11137391, + "pt": 8.30243334 + }, + { + "rap": -0.30763609, + "phi": 4.08076607, + "pt": 7.24431012 + }, + { + "rap": 1.03337474, + "phi": 5.75077915, + "pt": 7.21438996 + }, + { + "rap": 2.47115873, + "phi": 1.71463831, + "pt": 6.94079446 + }, + { + "rap": 3.08530764, + "phi": 5.66608526, + "pt": 6.93107373 + }, + { + "rap": 2.83311884, + "phi": 2.73639926, + "pt": 6.43631882 + }, + { + "rap": -2.32756975, + "phi": 2.27869505, + "pt": 6.42859668 + }, + { + "rap": 0.31742948, + "phi": 3.25493522, + "pt": 5.90365269 + }, + { + "rap": -1.7551644, + "phi": 3.07572508, + "pt": 5.84924081 + }, + { + "rap": 2.11141937, + "phi": 2.43625928, + "pt": 5.65180212 + } + ] + }, + { + "jetid": 91, + "jets": [ + { + "rap": 1.76222679, + "phi": 0.87410257, + "pt": 33.79601729 + }, + { + "rap": 1.0587973, + "phi": 4.40140166, + "pt": 32.82482432 + }, + { + "rap": 1.76672236, + "phi": 3.84906806, + "pt": 13.19158504 + }, + { + "rap": 3.2774358, + "phi": 6.09117064, + "pt": 9.4619856 + }, + { + "rap": -1.36215403, + "phi": 4.14607783, + "pt": 8.34910462 + }, + { + "rap": 0.37483886, + "phi": 4.88917573, + "pt": 7.91848609 + }, + { + "rap": -0.48195678, + "phi": 4.96363275, + "pt": 7.56938281 + }, + { + "rap": -0.28287109, + "phi": 4.02474303, + "pt": 7.24892274 + }, + { + "rap": 3.35290474, + "phi": 1.00635694, + "pt": 6.76549049 + }, + { + "rap": 1.51114836, + "phi": 3.02905199, + "pt": 6.53380851 + }, + { + "rap": -3.38764669, + "phi": 4.17405448, + "pt": 6.30770573 + }, + { + "rap": -0.86319972, + "phi": 5.69529015, + "pt": 6.21795371 + }, + { + "rap": 0.55025625, + "phi": 2.68252832, + "pt": 6.10451 + }, + { + "rap": 1.93581131, + "phi": 1.60685199, + "pt": 5.92587915 + }, + { + "rap": -2.73813892, + "phi": 0.95189082, + "pt": 5.8978053 + }, + { + "rap": 2.22888757, + "phi": 1.13634374, + "pt": 5.58867562 + }, + { + "rap": 0.75167335, + "phi": 0.09149569, + "pt": 5.56698118 + }, + { + "rap": -1.06702545, + "phi": 3.35000329, + "pt": 5.55144637 + }, + { + "rap": -3.2704607, + "phi": 1.1350432, + "pt": 5.52125423 + }, + { + "rap": -0.92671562, + "phi": 1.39757494, + "pt": 5.30637881 + }, + { + "rap": -1.88979165, + "phi": 3.96255411, + "pt": 5.02519413 + } + ] + }, + { + "jetid": 92, + "jets": [ + { + "rap": -0.23214577, + "phi": 1.5530677, + "pt": 24.07133124 + }, + { + "rap": 4.97891567, + "phi": 5.45226022, + "pt": 20.95387513 + }, + { + "rap": -3.6947431, + "phi": 5.51765668, + "pt": 10.89440573 + }, + { + "rap": -0.82711917, + "phi": 4.71637931, + "pt": 7.25153698 + }, + { + "rap": -2.9060308, + "phi": 4.16451714, + "pt": 6.32674035 + }, + { + "rap": 0.23165143, + "phi": 1.89490557, + "pt": 6.21003556 + }, + { + "rap": 2.97881309, + "phi": 1.5239026, + "pt": 5.83612948 + }, + { + "rap": -1.50068942, + "phi": 5.86374739, + "pt": 5.72164464 + }, + { + "rap": 1.69313299, + "phi": 5.05519143, + "pt": 5.56577513 + }, + { + "rap": -2.88395652, + "phi": 0.48004369, + "pt": 5.3352474 + }, + { + "rap": -3.53109, + "phi": 2.16251747, + "pt": 5.26592671 + }, + { + "rap": -2.97599626, + "phi": 1.61331552, + "pt": 5.25351024 + }, + { + "rap": -2.15120124, + "phi": 1.75227953, + "pt": 5.16380563 + } + ] + }, + { + "jetid": 93, + "jets": [ + { + "rap": -5.48442158, + "phi": 2.66300837, + "pt": 16.66496985 + }, + { + "rap": -4.31869505, + "phi": 5.11248151, + "pt": 10.23382828 + }, + { + "rap": 0.568638, + "phi": 5.40641458, + "pt": 9.07621482 + }, + { + "rap": 4.17860349, + "phi": 1.75276139, + "pt": 9.06239309 + }, + { + "rap": -0.01748142, + "phi": 1.07115688, + "pt": 8.98236902 + }, + { + "rap": 0.05247943, + "phi": 5.07530149, + "pt": 7.93733944 + }, + { + "rap": 1.11427199, + "phi": 6.06463718, + "pt": 7.22731973 + }, + { + "rap": 2.26441327, + "phi": 3.80886133, + "pt": 6.82027505 + }, + { + "rap": -0.85495475, + "phi": 0.33410243, + "pt": 6.20510721 + }, + { + "rap": -3.00400815, + "phi": 1.22399868, + "pt": 5.36392278 + }, + { + "rap": -3.3610046, + "phi": 0.96664676, + "pt": 5.12989605 + }, + { + "rap": 3.14341268, + "phi": 4.84867272, + "pt": 5.08998152 + } + ] + }, + { + "jetid": 94, + "jets": [ + { + "rap": -0.58557583, + "phi": 3.81252153, + "pt": 24.4613419 + }, + { + "rap": 2.83486786, + "phi": 1.65235337, + "pt": 17.23777744 + }, + { + "rap": 5.01365726, + "phi": 0.39564819, + "pt": 16.92673402 + }, + { + "rap": 3.13644559, + "phi": 2.01369685, + "pt": 14.60281354 + }, + { + "rap": 0.50234715, + "phi": 4.0074401, + "pt": 12.60431345 + }, + { + "rap": 0.59707161, + "phi": 2.20968761, + "pt": 10.58081394 + }, + { + "rap": 4.93363493, + "phi": 3.73896901, + "pt": 10.27737674 + }, + { + "rap": -2.05844964, + "phi": 3.80415836, + "pt": 9.82987965 + }, + { + "rap": 4.21063991, + "phi": 0.3679943, + "pt": 9.60541985 + }, + { + "rap": 0.01701359, + "phi": 0.8366188, + "pt": 8.93847696 + }, + { + "rap": -0.12167286, + "phi": 3.72814299, + "pt": 8.52998029 + }, + { + "rap": 0.99529612, + "phi": 2.63574605, + "pt": 8.34217213 + }, + { + "rap": 2.19669502, + "phi": 0.15951296, + "pt": 8.32410313 + }, + { + "rap": 3.36574559, + "phi": 0.36773986, + "pt": 8.18827873 + }, + { + "rap": -1.29489326, + "phi": 0.28600148, + "pt": 8.03929003 + }, + { + "rap": 0.80750315, + "phi": 1.03148915, + "pt": 7.76370663 + }, + { + "rap": -0.10430129, + "phi": 4.29048863, + "pt": 7.62447831 + }, + { + "rap": 3.82825949, + "phi": 4.89245827, + "pt": 7.56986904 + }, + { + "rap": 0.27977537, + "phi": 4.85572345, + "pt": 7.56386661 + }, + { + "rap": 1.7052144, + "phi": 1.48859073, + "pt": 7.54758312 + }, + { + "rap": 0.75521198, + "phi": 4.51268131, + "pt": 7.34752951 + }, + { + "rap": -0.5881971, + "phi": 5.67957786, + "pt": 7.33279628 + }, + { + "rap": -0.24519851, + "phi": 2.84343907, + "pt": 6.36974973 + }, + { + "rap": 2.80928831, + "phi": 0.18724231, + "pt": 6.2393796 + }, + { + "rap": 0.57855184, + "phi": 5.57229134, + "pt": 6.2389455 + }, + { + "rap": -1.56681945, + "phi": 3.28134897, + "pt": 5.95444807 + }, + { + "rap": -0.99287493, + "phi": 0.69015883, + "pt": 5.58326197 + }, + { + "rap": -1.4787118, + "phi": 5.76490146, + "pt": 5.35286386 + }, + { + "rap": 1.08088956, + "phi": 2.02983158, + "pt": 5.15233704 + }, + { + "rap": -1.13644465, + "phi": 2.83303969, + "pt": 5.13726778 + } + ] + }, + { + "jetid": 95, + "jets": [ + { + "rap": -1.31767191, + "phi": 3.11053963, + "pt": 29.72785205 + }, + { + "rap": -0.08886018, + "phi": 0.45634152, + "pt": 24.06849223 + }, + { + "rap": 2.92238925, + "phi": 5.20108267, + "pt": 11.24975421 + }, + { + "rap": -0.89794779, + "phi": 3.79712289, + "pt": 7.92451285 + }, + { + "rap": -5.76215518, + "phi": 6.2429974, + "pt": 7.81819248 + }, + { + "rap": -2.58083387, + "phi": 1.10369184, + "pt": 6.70884023 + } + ] + }, + { + "jetid": 96, + "jets": [ + { + "rap": -4.42627153, + "phi": 1.44275131, + "pt": 22.90169969 + }, + { + "rap": -0.43132041, + "phi": 4.07279612, + "pt": 20.89434516 + }, + { + "rap": -4.13708274, + "phi": 6.2241099, + "pt": 11.03194265 + }, + { + "rap": 3.84708688, + "phi": 4.92812625, + "pt": 9.98124544 + } + ] + }, + { + "jetid": 97, + "jets": [ + { + "rap": -0.58016363, + "phi": 0.39531652, + "pt": 62.96613366 + }, + { + "rap": -2.05543358, + "phi": 4.0185368, + "pt": 40.22893167 + }, + { + "rap": -0.97181225, + "phi": 2.99028235, + "pt": 26.16807412 + }, + { + "rap": 0.14656718, + "phi": 1.66084646, + "pt": 7.6323789 + }, + { + "rap": 4.37233876, + "phi": 4.78475181, + "pt": 5.27702676 + } + ] + }, + { + "jetid": 98, + "jets": [ + { + "rap": 0.36623023, + "phi": 2.59996628, + "pt": 22.9775586 + }, + { + "rap": 1.47651468, + "phi": 6.14462163, + "pt": 13.182296 + }, + { + "rap": 1.63379294, + "phi": 5.60419037, + "pt": 6.38754417 + } + ] + }, + { + "jetid": 99, + "jets": [ + { + "rap": -1.67548151, + "phi": 2.92674075, + "pt": 31.63867634 + }, + { + "rap": -1.62828031, + "phi": 5.38144848, + "pt": 22.25137197 + }, + { + "rap": 2.49950673, + "phi": 0.69203587, + "pt": 18.77465183 + }, + { + "rap": 0.43085178, + "phi": 2.37968129, + "pt": 14.23410076 + }, + { + "rap": 0.43839415, + "phi": 3.9532408, + "pt": 9.98833187 + }, + { + "rap": -0.52051885, + "phi": 3.23942862, + "pt": 9.25009352 + }, + { + "rap": -4.44198614, + "phi": 5.91804361, + "pt": 9.08973099 + }, + { + "rap": -2.22176303, + "phi": 2.21759371, + "pt": 6.90958168 + }, + { + "rap": 2.97813895, + "phi": 1.86057601, + "pt": 5.91342959 + }, + { + "rap": 0.27306655, + "phi": 3.37390185, + "pt": 5.25913207 + }, + { + "rap": 1.33965552, + "phi": 4.61473586, + "pt": 5.20977775 + }, + { + "rap": 3.22340199, + "phi": 3.60283432, + "pt": 5.14644876 + } + ] + }, + { + "jetid": 100, + "jets": [ + { + "rap": -0.74746481, + "phi": 6.12278123, + "pt": 65.45917884 + }, + { + "rap": 1.21432316, + "phi": 2.60725593, + "pt": 32.70404135 + }, + { + "rap": -0.09913624, + "phi": 2.77254168, + "pt": 19.9110653 + }, + { + "rap": 2.8753562, + "phi": 4.7193044, + "pt": 17.30978393 + }, + { + "rap": 4.95717688, + "phi": 0.27136724, + "pt": 7.68169839 + }, + { + "rap": 5.7067729, + "phi": 2.97176675, + "pt": 7.15267437 + }, + { + "rap": -5.57103521, + "phi": 2.43610538, + "pt": 5.67666524 + }, + { + "rap": 0.48552972, + "phi": 3.32968196, + "pt": 5.11498298 + } + ] + } +] \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index 24ad069..8ea0bea 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,6 +3,7 @@ using JetReconstruction using Test using JSON +using LorentzVectorHEP """Read JSON file with fastjet jets in it""" function read_fastjet_outputs(;fname="test/data/jet_collections_fastjet.json") @@ -30,12 +31,7 @@ function main() # Test each stratgy... do_jet_test(N2Plain, fastjet_jets) - do_jet_test(N2TiledSoAGlobal, fastjet_jets) - do_jet_test(N2TiledSoATile, fastjet_jets) do_jet_test(N2Tiled, fastjet_jets) - - # Atell's original test - original_tests() end function do_jet_test(strategy::JetRecoStrategy, fastjet_jets; @@ -45,40 +41,37 @@ 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 == N2TiledSoAGlobal) - jet_reconstruction = tiled_jet_reconstruct_soa_global - strategy_name = "N2TiledSoAGlobal" - elseif (strategy == N2TiledSoATile) - jet_reconstruction = tiled_jet_reconstruct_soa_tile - strategy_name = "N2TiledSoATile" elseif (strategy == N2Tiled) - jet_reconstruction = tiled_jet_reconstruct_ll + jet_reconstruction = tiled_jet_reconstruct strategy_name = "N2Tiled" else throw(ErrorException("Strategy not yet implemented")) end # Now run our jet reconstruction... + # From PseudoJets events::Vector{Vector{PseudoJet}} = read_final_state_particles("test/data/events.hepmc3") - if strategy == N2Tiled - event_vector = events - else - # First, convert all events into the Vector of Vectors that Atell's - # code likes - event_vector = pseudojets2vectors(events) - end - # event_vector = pseudojets2vectors(events) jet_collection = FinalJets[] - for (ievt, event) in enumerate(event_vector) - finaljets, _ = jet_reconstruction(event, R=distance, p=power) + for (ievt, event) in enumerate(events) + finaljets, _ = jet_reconstruction(event, R=distance, p=power, ptmin=ptmin) fj = final_jets(finaljets, ptmin) sort_jets!(fj) push!(jet_collection, FinalJets(ievt, fj)) end - @testset "Jet Reconstruction, $strategy_name" begin + # From LorentzVector + events_lv::Vector{Vector{LorentzVector}} = read_final_state_particles_lv("test/data/events.hepmc3") + jet_collection_lv = FinalJets[] + for (ievt, event) in enumerate(events_lv) + finaljets, _ = jet_reconstruction(event, R=distance, p=power, ptmin=ptmin) + fj = final_jets(finaljets, ptmin) + sort_jets!(fj) + push!(jet_collection_lv, FinalJets(ievt, fj)) + end + + @testset "Jet Reconstruction PseudoJet, $strategy_name" begin # Test each event in turn... for (ievt, event) in enumerate(jet_collection) @testset "Event $(ievt)" begin @@ -90,61 +83,33 @@ function do_jet_test(strategy::JetRecoStrategy, fastjet_jets; # tolerance into the isapprox() function # Use atol for position as this is absolute, but rtol for # the momentum + # Sometimes phi could be in the range [-π, π], but FastJet always is [0, 2π] + normalised_phi = jet.phi < 0.0 ? jet.phi + 2π : jet.phi @test jet.rap ≈ fastjet_jets[ievt]["jets"][ijet]["rap"] atol=1e-7 - @test jet.phi ≈ fastjet_jets[ievt]["jets"][ijet]["phi"] atol=1e-7 + @test normalised_phi ≈ fastjet_jets[ievt]["jets"][ijet]["phi"] atol=1e-7 @test jet.pt ≈ fastjet_jets[ievt]["jets"][ijet]["pt"] rtol=1e-6 end end end end end -end -"""Original test implementation""" -function original_tests() - """ - A function for test comparison - """ - function arrcompare(y, yt; eps=0.1) - for i in eachindex(y) - if !(sum(abs.(y[i] .- yt[i]) .< eps) == length(y[i])) - return false + @testset "Jet Reconstruction LorentzVector, $strategy_name" begin + # Here we test that inputing LorentzVector gave the same results as PseudoJets + for (ievt, (event, event_lv)) in enumerate(zip(jet_collection, jet_collection_lv)) + @testset "Event $(ievt)" begin + @test size(event.jets) == size(event_lv.jets) + # Test each jet in turn + for (jet, jet_lv) in zip(event.jets, event_lv.jets) + @test jet.rap ≈ jet_lv.rap atol=1e-7 + @test jet.phi ≈ jet_lv.phi atol=1e-7 + @test jet.pt ≈ jet_lv.pt rtol=1e-6 + end end end - true end - NUMBER_OF_TESTS = 12 # number of test files in the data folder - - @testset "N2Plain Original Tests" begin - # @test anti_kt(X) == Y - - # additional simple test - @test arrcompare( - anti_kt_algo([ - [99.0, 0.1, 0, 100], - [4.0, -0.1, 0, 5.0], - [-99., 0., 0.0, 99] - ], R=0.7)[1], - - [[103.0, 0.0, 0.0, 105.0], [-99.0, 0.0, 0.0, 99.0]] - ) - # actual testing - for i in 1:NUMBER_OF_TESTS - istr = string(i) - - @test arrcompare( - sort( - anti_kt_algo( - loadjets("test/data/"*istr*".dat"), R=1 - )[1], lt=(a,b)->(pt(a)>pt(b)) - ), - - loadjets("test/data/"*istr*"-fj-result.dat") - ) - end - end end if abspath(PROGRAM_FILE) == @__FILE__