From ff41739cff72ee00fce9af6a9609921247b54f38 Mon Sep 17 00:00:00 2001 From: Graeme A Stewart Date: Fri, 15 Sep 2023 15:28:15 +0200 Subject: [PATCH 01/21] Remove sub-optimal algorithms Delete the less performant tiled algorithms --- chep.jl | 34 +-- src/JetReconstruction.jl | 12 +- src/TiledAlgoLL.jl | 1 + src/TiledAlgoSoAGlobal.jl | 578 -------------------------------------- src/TiledAlgoSoATile.jl | 335 ---------------------- src/TiledUtilsSoA.jl | 240 ---------------- test/runtests.jl | 8 - 7 files changed, 20 insertions(+), 1188 deletions(-) delete mode 100644 src/TiledAlgoSoAGlobal.jl delete mode 100644 src/TiledAlgoSoATile.jl delete mode 100644 src/TiledUtilsSoA.jl diff --git a/chep.jl b/chep.jl index f35a7bd..0607a9e 100644 --- a/chep.jl +++ b/chep.jl @@ -71,11 +71,7 @@ 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) + elseif (strategy == N2Tiled || stragegy == Best) jet_reconstruction = tiled_jet_reconstruct_ll else throw(ErrorException("Strategy not yet implemented")) @@ -121,6 +117,7 @@ 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() @@ -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/src/JetReconstruction.jl b/src/JetReconstruction.jl index 456c63b..9e9bae2 100644 --- a/src/JetReconstruction.jl +++ b/src/JetReconstruction.jl @@ -24,14 +24,6 @@ export sequential_jet_reconstruct, kt_algo, anti_kt_algo, anti_kt_algo_alt, camb # 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 @@ -53,7 +45,7 @@ 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 +@enum JetRecoStrategy Best N2Plain N2Tiled +export JetRecoStrategy, Best, N2Plain, N2Tiled end diff --git a/src/TiledAlgoLL.jl b/src/TiledAlgoLL.jl index ae17696..437d005 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") 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/test/runtests.jl b/test/runtests.jl index 24ad069..2619bbb 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -30,8 +30,6 @@ 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 @@ -47,12 +45,6 @@ function do_jet_test(strategy::JetRecoStrategy, fastjet_jets; if (strategy == N2Plain) jet_reconstruction = sequential_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 strategy_name = "N2Tiled" From 6a90da921ca7c4da64d82e7e4ddbad10a59b99df Mon Sep 17 00:00:00 2001 From: Graeme A Stewart Date: Fri, 6 Oct 2023 16:47:43 +0200 Subject: [PATCH 02/21] Switch to LorentzVector inputs Using LV inputs gives a standard interface to the reconstruction algorithms Implement HepMC reader that can return LVs Adapt N2Plain algorithm to digest these inputs (there is a slight performance regression, 780->833 us on my machine) Clean up all of the "dev' code from the N2Plain algorithm --- chep.jl | 30 ++++---- src/Algo.jl | 145 +++++---------------------------------- src/JetReconstruction.jl | 4 +- src/Utils.jl | 43 +++++++++++- 4 files changed, 79 insertions(+), 143 deletions(-) diff --git a/chep.jl b/chep.jl index 0607a9e..53e6236 100644 --- a/chep.jl +++ b/chep.jl @@ -14,6 +14,7 @@ using StatProfilerHTML using Logging using JSON +using LorentzVectorHEP using JetReconstruction function profile_code(jet_reconstruction, events, niters) @@ -55,7 +56,7 @@ function profile_code(jet_reconstruction, events, niters) end function jet_process( - events::Vector{Vector{PseudoJet}}; + events::Vector{Vector{LorentzVector}}; ptmin::Float64 = 5.0, distance::Float64 = 0.4, power::Integer = -1, @@ -77,13 +78,14 @@ function jet_process( throw(ErrorException("Strategy not yet implemented")) end + # Change of tactic - we adapt all interfaces to swallow LorentzVectors directly # The N2Tiled algorithm uses PseudoJets so pass these directly - if strategy == N2Tiled - event_vector = events - else - # The other algorithms swallow 4-vectors instead - event_vector = pseudojets2vectors(events) - end + # if strategy == N2Tiled + # 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 if !isnothing(dump) @@ -93,20 +95,20 @@ function jet_process( # Warmup code if we are doing a multi-sample timing run if nsamples > 1 || profile @info "Doing initial warm-up run" - for event in event_vector + for event in events finaljets, _ = jet_reconstruction(event, R = distance, p = power) final_jets(finaljets, ptmin) end end if profile - profile_code(jet_reconstruction, event_vector, nsamples) + profile_code(jet_reconstruction, events, nsamples) return nothing end if alloc println("Memory allocation statistics:") - @timev for event in event_vector + @timev for event in events finaljets, _ = jet_reconstruction(event, R = distance, p = power) final_jets(finaljets, ptmin) end @@ -121,7 +123,7 @@ function jet_process( for irun ∈ 1:nsamples gcoff && GC.enable(false) t_start = time_ns() - for (ievt, event) in enumerate(event_vector) + for (ievt, event) in enumerate(events) finaljets, _ = jet_reconstruction(event, R = distance, p = power) fj = final_jets(finaljets, ptmin) # Only print the jet content once @@ -268,8 +270,10 @@ main() = begin logger = ConsoleLogger(stdout, Logging.Warn) end global_logger(logger) - events::Vector{Vector{PseudoJet}} = - read_final_state_particles(args[:file], maxevents = args[:maxevents], skipevents = args[:skip]) + # events::Vector{Vector{PseudoJet}} = + # read_final_state_particles(args[:file], maxevents = args[:maxevents], skipevents = args[:skip]) + events::Vector{Vector{LorentzVector}} = + read_final_state_particles_lv(args[:file], maxevents = args[:maxevents], skipevents = args[:skip]) jet_process(events, ptmin = args[:ptmin], distance = args[:distance], power = args[:power], strategy = args[:strategy], nsamples = args[:nsamples], gcoff = args[:gcoff], profile = args[:profile], diff --git a/src/Algo.jl b/src/Algo.jl index 622f406..97cce90 100644 --- a/src/Algo.jl +++ b/src/Algo.jl @@ -1,4 +1,3 @@ - Base.@propagate_inbounds function _dist(i, j, _eta, _phi) deta = _eta[i] - _eta[j] dphi = abs(_phi[i] - _phi[j]) @@ -56,8 +55,10 @@ Base.@propagate_inbounds function _upd_nn_nocross!(i::Int, from::Int, to::Int, _ end # entire NN update step +# Base.@propagate_inbounds function _upd_nn_step!(i::Int, j::Int, k::Int, N::Int, Nn::Int, +# _kt2::Vector{Float64}, _eta::Vector{Float64}, _phi::Vector{Float64}, _R2::Float64, _nndist::Vector{Float64}, _nn::Vector{Int}, _nndij::Vector{Float64}) Base.@propagate_inbounds function _upd_nn_step!(i, j, k, N, Nn, _kt2, _eta, _phi, _R2, _nndist, _nn, _nndij) - nnk = _nn[k] + 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) @@ -80,24 +81,25 @@ Base.@propagate_inbounds function _upd_nn_step!(i, j, k, N, Nn, _kt2, _eta, _phi #nothing end + function sequential_jet_reconstruct(objects::AbstractArray{T}; p=-1.0, R=1., recombine=+) where T - # bounds + # Bounds N::Int = length(objects) - # returned values + # Returned values jets = T[] # result sequences = Vector{Int}[] # recombination sequences, WARNING: first index in the sequence is not necessarily the seed - # params + # Parameters _R2 = R*R _p = (round(p) == p) ? Int(p) : p # integer p if possible - # data + # 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) + ## When working with LorentzVectorHEP we make sure these arrays are type stable + _kt2::Vector{Float64} = (LorentzVectorHEP.pt.(_objects) .^ 2) .^ _p + _phi::Vector{Float64} = LorentzVectorHEP.phi.(_objects) + _eta::Vector{Float64} = LorentzVectorHEP.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] @@ -108,7 +110,7 @@ function sequential_jet_reconstruct(objects::AbstractArray{T}; p=-1.0, R=1., rec end # diJ table *_R2 - _nndij = zeros(N) + _nndij::Vector{Float64} = zeros(N) @inbounds @simd for i in 1:N _nndij[i] = _dij(i, _kt2, _nn, _nndist) end @@ -138,9 +140,9 @@ function sequential_jet_reconstruct(objects::AbstractArray{T}; p=-1.0, R=1., rec # 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 + _phi[i] = LorentzVectorHEP.phi(_objects[i]) + _eta[i] = LorentzVectorHEP.eta(_objects[i]) + _kt2[i] = (LorentzVectorHEP.pt(_objects[i])^2)^_p _nndist[i] = _R2 _nn[i] = i @@ -173,6 +175,7 @@ function sequential_jet_reconstruct(objects::AbstractArray{T}; p=-1.0, R=1., rec @inbounds @simd for k in 1:N _upd_nn_step!(i, j, k, N, Nn, _kt2, _eta, _phi, _R2, _nndist, _nn, _nndij) end + # @infiltrate _nndij[i] = _dij(i, _kt2, _nn, _nndist) end @@ -218,117 +221,3 @@ Returns: 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 9e9bae2..72f0716 100644 --- a/src/JetReconstruction.jl +++ b/src/JetReconstruction.jl @@ -3,6 +3,8 @@ Jet reconstruction (reclustering) in Julia. """ module JetReconstruction +using LorentzVectorHEP + # particle type definition include("Particle.jl") export energy, px, py, pz, pt, phi, mass, eta, kt, ϕ, η @@ -34,7 +36,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") diff --git a/src/Utils.jl b/src/Utils.jl index b48e702..46d83c5 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 @@ -72,3 +99,17 @@ function final_jets(jets::Vector{PseudoJet}, ptmin::AbstractFloat) 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.eta(jet), LorentzVectorHEP.phi(jet), LorentzVectorHEP.pt(jet))) + end + end + final_jets +end From f1fc60d930fe29906e2097c8f22be99d5bb6995e Mon Sep 17 00:00:00 2001 From: Graeme A Stewart Date: Tue, 10 Oct 2023 22:17:01 +0200 Subject: [PATCH 03/21] Change to input LorentzVectors These are much more standard and supported objects from the LorentzVectorHEP package This drops use of the internal "Particle.jl" struct. Do not export PseudoJet functions. Some standard reformatting applied --- src/Algo.jl | 45 ++++++++++---------- src/JetReconstruction.jl | 8 ++-- src/Pseudojet.jl | 7 ++- src/TiledAlgoLL.jl | 92 ++++++++++++++++++++++++---------------- src/Utils.jl | 14 ++++++ test/runtests.jl | 21 ++++----- 6 files changed, 114 insertions(+), 73 deletions(-) diff --git a/src/Algo.jl b/src/Algo.jl index 97cce90..949b3d0 100644 --- a/src/Algo.jl +++ b/src/Algo.jl @@ -2,14 +2,14 @@ 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) + 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]) + d * min(_kt2[i], _kt2[j]) end # finds new nn for i and checks everyone additionally @@ -28,8 +28,7 @@ Base.@propagate_inbounds function _upd_nn_crosscheck!(i::Int, from::Int, to::Int end end _nndist[i] = nndist - _nn[i] = nn; - #nothing + _nn[i] = nn end # finds new nn for i @@ -50,15 +49,14 @@ Base.@propagate_inbounds function _upd_nn_nocross!(i::Int, from::Int, to::Int, _ nndist = ifelse(f, Δ2, nndist) end _nndist[i] = nndist - _nn[i] = nn; - #nothing + _nn[i] = nn end # entire NN update step # Base.@propagate_inbounds function _upd_nn_step!(i::Int, j::Int, k::Int, N::Int, Nn::Int, # _kt2::Vector{Float64}, _eta::Vector{Float64}, _phi::Vector{Float64}, _R2::Float64, _nndist::Vector{Float64}, _nn::Vector{Int}, _nndij::Vector{Float64}) Base.@propagate_inbounds function _upd_nn_step!(i, j, k, N, Nn, _kt2, _eta, _phi, _R2, _nndist, _nn, _nndij) - nnk = _nn[k] + 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) @@ -74,24 +72,23 @@ Base.@propagate_inbounds function _upd_nn_step!(i, j, k, N, Nn, _kt2, _eta, _phi end cond = Δ2 < _nndist[i] - _nndist[i], _nn[i] = ifelse(cond, (Δ2,k), (_nndist[i],_nn[i])) + _nndist[i], _nn[i] = ifelse(cond, (Δ2, k), (_nndist[i], _nn[i])) end - nnk == Nn && (_nn[k] = j); - #nothing + nnk == Nn && (_nn[k] = j) end -function sequential_jet_reconstruct(objects::AbstractArray{T}; p=-1.0, R=1., recombine=+) where T +function sequential_jet_reconstruct(objects::AbstractArray{T}; p = -1.0, R = 1.0, recombine = +) where T # Bounds N::Int = length(objects) # Returned values - jets = T[] # result + jets = T[] sequences = Vector{Int}[] # recombination sequences, WARNING: first index in the sequence is not necessarily the seed # Parameters - _R2 = R*R + _R2 = R * R _p = (round(p) == p) ? Int(p) : p # integer p if possible # Data @@ -99,14 +96,14 @@ function sequential_jet_reconstruct(objects::AbstractArray{T}; p=-1.0, R=1., rec ## When working with LorentzVectorHEP we make sure these arrays are type stable _kt2::Vector{Float64} = (LorentzVectorHEP.pt.(_objects) .^ 2) .^ _p _phi::Vector{Float64} = LorentzVectorHEP.phi.(_objects) - _eta::Vector{Float64} = LorentzVectorHEP.eta.(_objects) + _eta::Vector{Float64} = LorentzVectorHEP.rap.(_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) + _upd_nn_crosscheck!(i, 1, i - 1, _eta, _phi, _R2, _nndist, _nn) end # diJ table *_R2 @@ -121,7 +118,7 @@ function sequential_jet_reconstruct(objects::AbstractArray{T}; p=-1.0, R=1., rec i = 1 dij_min = _nndij[1] @inbounds @simd for k in 2:N - cond = _nndij[k] < dij_min + cond = _nndij[k] < dij_min dij_min, i = ifelse(cond, (_nndij[k], k), (dij_min, i)) end @@ -141,7 +138,7 @@ function sequential_jet_reconstruct(objects::AbstractArray{T}; p=-1.0, R=1., rec # update ith jet, replacing it with the new one _objects[i] = recombine(_objects[i], _objects[j]) _phi[i] = LorentzVectorHEP.phi(_objects[i]) - _eta[i] = LorentzVectorHEP.eta(_objects[i]) + _eta[i] = LorentzVectorHEP.rap(_objects[i]) _kt2[i] = (LorentzVectorHEP.pt(_objects[i])^2)^_p _nndist[i] = _R2 @@ -151,6 +148,8 @@ function sequential_jet_reconstruct(objects::AbstractArray{T}; p=-1.0, R=1., rec push!(_sequences[i], x) end else # i == j + # push!(jets, LorentzVectorCyl(LorentzVectorHEP.pt(_objects[i]), LorentzVectorHEP.eta(_objects[i]), + # LorentzVectorHEP.phi(_objects[i]), LorentzVectorHEP.mt(_objects[i]))) push!(jets, _objects[i]) push!(sequences, _sequences[i]) end @@ -192,8 +191,8 @@ 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) +function anti_kt_algo(objects; R = 1.0, recombine = +) + sequential_jet_reconstruct(objects, p = -1, R = R, recombine = recombine) end """ @@ -205,8 +204,8 @@ 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) +function kt_algo(objects; R = 1.0, recombine = +) + sequential_jet_reconstruct(objects, p = 1, R = R, recombine = recombine) end """ @@ -218,6 +217,6 @@ 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) +function cambridge_aachen_algo(objects; R = 1.0, recombine = +) + sequential_jet_reconstruct(objects, p = 0, R = R, recombine = recombine) end diff --git a/src/JetReconstruction.jl b/src/JetReconstruction.jl index 72f0716..56c769e 100644 --- a/src/JetReconstruction.jl +++ b/src/JetReconstruction.jl @@ -6,12 +6,13 @@ module JetReconstruction using LorentzVectorHEP # particle type definition -include("Particle.jl") -export energy, px, py, pz, pt, phi, mass, eta, kt, ϕ, η +# include("Particle.jl") +# export energy, px, py, pz, pt, phi, mass, eta, kt, ϕ, η # Philipp's pseudojet include("Pseudojet.jl") -export PseudoJet, rap, phi, pt2 +## As this is an internal EDM class, we perhaps shouldn't export this stuff... +# export PseudoJet, rap, phi, pt, pt2, px, py, pz, pt, phi, mass, eta # Simple HepMC3 reader include("HepMC3.jl") @@ -47,6 +48,7 @@ include("JSONresults.jl") export FinalJet, FinalJets, JSON3 # Strategy to be used +## Maybe an enum is not the best idea, use type dispatch instead? @enum JetRecoStrategy Best N2Plain N2Tiled export JetRecoStrategy, Best, N2Plain, N2Tiled diff --git a/src/Pseudojet.jl b/src/Pseudojet.jl index d2a9298..d90a0f0 100644 --- a/src/Pseudojet.jl +++ b/src/Pseudojet.jl @@ -77,7 +77,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 @@ -145,6 +145,11 @@ m(p::PseudoJet) = begin x < 0. ? -sqrt(-x) : sqrt(x) end +px(p::PseudoJet) = p.px +py(p::PseudoJet) = p.py +pz(p::PseudoJet) = p.pz +mass(p::PseudoJet) = m(p) + import Base.+; +(j1::PseudoJet, j2::PseudoJet) = begin diff --git a/src/TiledAlgoLL.jl b/src/TiledAlgoLL.jl index 437d005..37f1bec 100644 --- a/src/TiledAlgoLL.jl +++ b/src/TiledAlgoLL.jl @@ -18,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 @@ -39,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 @@ -57,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 = 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.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 @@ -78,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 @@ -91,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 @@ -134,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 @@ -160,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 @@ -170,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) @@ -196,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 @@ -215,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 @@ -238,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 @@ -268,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, @@ -281,7 +288,8 @@ 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), rap(jet), phi(jet), mass(jet))) + # push!(jets_local, jet) end end jets_local @@ -290,23 +298,35 @@ end """ Main jet reconstruction algorithm """ +function tiled_jet_reconstruct_ll(particles::Vector{LorentzVector}; p = -1, R = 1.0, recombine = +, ptmin = 0.0) + # 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(LorentzVectorHEP.px(particle), LorentzVectorHEP.py(particle), + LorentzVectorHEP.pz(particle), LorentzVectorHEP.energy(particle)) + end + tiled_jet_reconstruct_ll(pseudojets, p = p, R = R, recombine = recombine, ptmin = ptmin) +end + + function tiled_jet_reconstruct_ll(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 @@ -343,7 +363,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] @@ -361,7 +381,7 @@ 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 @@ -371,7 +391,7 @@ function tiled_jet_reconstruct_ll(particles::Vector{PseudoJet}; p = -1, R = 1.0, 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) @@ -392,7 +412,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 @@ -416,14 +436,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/Utils.jl b/src/Utils.jl index 46d83c5..febecac 100644 --- a/src/Utils.jl +++ b/src/Utils.jl @@ -101,6 +101,20 @@ function final_jets(jets::Vector{PseudoJet}, ptmin::AbstractFloat) 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.rap(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) diff --git a/test/runtests.jl b/test/runtests.jl index 2619bbb..f32f6ac 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") @@ -33,7 +34,7 @@ function main() do_jet_test(N2Tiled, fastjet_jets) # Atell's original test - original_tests() + # original_tests() end function do_jet_test(strategy::JetRecoStrategy, fastjet_jets; @@ -53,17 +54,17 @@ function do_jet_test(strategy::JetRecoStrategy, fastjet_jets; end # Now run our jet reconstruction... - 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 + events::Vector{Vector{LorentzVector}} = read_final_state_particles_lv("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) + for (ievt, event) in enumerate(events) finaljets, _ = jet_reconstruction(event, R=distance, p=power) fj = final_jets(finaljets, ptmin) sort_jets!(fj) From 8025d4ca70cbd3619bfb4ad0ab43d855febd7f63 Mon Sep 17 00:00:00 2001 From: Graeme A Stewart Date: Wed, 11 Oct 2023 15:43:28 +0200 Subject: [PATCH 04/21] Update to use pt2 (pre-squared) for kt2 This is used in N2Plain Adapt tests to allow for algorithms returning phi in the [-\pi, \pi] range --- src/Algo.jl | 4 ++-- test/runtests.jl | 4 +++- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/src/Algo.jl b/src/Algo.jl index 949b3d0..34fd832 100644 --- a/src/Algo.jl +++ b/src/Algo.jl @@ -94,7 +94,7 @@ function sequential_jet_reconstruct(objects::AbstractArray{T}; p = -1.0, R = 1.0 # Data _objects = copy(objects) ## When working with LorentzVectorHEP we make sure these arrays are type stable - _kt2::Vector{Float64} = (LorentzVectorHEP.pt.(_objects) .^ 2) .^ _p + _kt2::Vector{Float64} = LorentzVectorHEP.pt2.(_objects) .^ _p _phi::Vector{Float64} = LorentzVectorHEP.phi.(_objects) _eta::Vector{Float64} = LorentzVectorHEP.rap.(_objects) _nn = Vector(1:N) # nearest neighbours @@ -139,7 +139,7 @@ function sequential_jet_reconstruct(objects::AbstractArray{T}; p = -1.0, R = 1.0 _objects[i] = recombine(_objects[i], _objects[j]) _phi[i] = LorentzVectorHEP.phi(_objects[i]) _eta[i] = LorentzVectorHEP.rap(_objects[i]) - _kt2[i] = (LorentzVectorHEP.pt(_objects[i])^2)^_p + _kt2[i] = LorentzVectorHEP.pt2(_objects[i])^_p _nndist[i] = _R2 _nn[i] = i diff --git a/test/runtests.jl b/test/runtests.jl index f32f6ac..93a4f1c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -83,8 +83,10 @@ 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 From e201573283dd4d245fae9c8b18fd858b03957e1e Mon Sep 17 00:00:00 2001 From: Graeme A Stewart Date: Fri, 20 Oct 2023 11:45:22 +0200 Subject: [PATCH 05/21] Update LorentzVectorHEP New version supports rapidity() method, as well as pt2() and more consistent handling of mass() --- Project.toml | 1 + src/Algo.jl | 4 ++-- src/Utils.jl | 2 +- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index e31697b..661b475 100644 --- a/Project.toml +++ b/Project.toml @@ -20,6 +20,7 @@ StatProfilerHTML = "a8a75453-ed82-57c9-9e16-4cd1196ecbf5" [compat] julia = "1.8" +LorentzVectorHEP = "0.1.6" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/src/Algo.jl b/src/Algo.jl index 34fd832..ad824e4 100644 --- a/src/Algo.jl +++ b/src/Algo.jl @@ -96,7 +96,7 @@ function sequential_jet_reconstruct(objects::AbstractArray{T}; p = -1.0, R = 1.0 ## When working with LorentzVectorHEP we make sure these arrays are type stable _kt2::Vector{Float64} = LorentzVectorHEP.pt2.(_objects) .^ _p _phi::Vector{Float64} = LorentzVectorHEP.phi.(_objects) - _eta::Vector{Float64} = LorentzVectorHEP.rap.(_objects) + _eta::Vector{Float64} = LorentzVectorHEP.rapidity.(_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] @@ -138,7 +138,7 @@ function sequential_jet_reconstruct(objects::AbstractArray{T}; p = -1.0, R = 1.0 # update ith jet, replacing it with the new one _objects[i] = recombine(_objects[i], _objects[j]) _phi[i] = LorentzVectorHEP.phi(_objects[i]) - _eta[i] = LorentzVectorHEP.rap(_objects[i]) + _eta[i] = LorentzVectorHEP.rapidity(_objects[i]) _kt2[i] = LorentzVectorHEP.pt2(_objects[i])^_p _nndist[i] = _R2 diff --git a/src/Utils.jl b/src/Utils.jl index febecac..8f209a6 100644 --- a/src/Utils.jl +++ b/src/Utils.jl @@ -108,7 +108,7 @@ function final_jets(jets::Vector{LorentzVector}, ptmin::AbstractFloat) for jet in jets if LorentzVectorHEP.pt(jet)^2 > dcut count += 1 - push!(final_jets, FinalJet(LorentzVectorHEP.rap(jet), LorentzVectorHEP.phi(jet), LorentzVectorHEP.pt(jet))) + push!(final_jets, FinalJet(LorentzVectorHEP.rapidity(jet), LorentzVectorHEP.phi(jet), LorentzVectorHEP.pt(jet))) end end final_jets From 304759982e69e835c421259a4dfa69a249257ddf Mon Sep 17 00:00:00 2001 From: Graeme A Stewart Date: Sun, 22 Oct 2023 21:12:52 +0200 Subject: [PATCH 06/21] Add missing FastJet JSON output file This was mistakenly masked in .gitignore --- .gitignore | 1 - test/data/jet_collections_fastjet.json | 4362 ++++++++++++++++++++++++ 2 files changed, 4362 insertions(+), 1 deletion(-) create mode 100644 test/data/jet_collections_fastjet.json 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/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 From 6afac594a9e265e1124eff960ecc7958d9b37642 Mon Sep 17 00:00:00 2001 From: Graeme A Stewart Date: Sun, 22 Oct 2023 21:14:30 +0200 Subject: [PATCH 07/21] Refactor exports Do not export "internal" pieces from JetReconstruction Add "import" of relevant methods from LorentzVectorHEP for use in the algorithms --- src/JetReconstruction.jl | 14 +++++++++----- 1 file changed, 9 insertions(+), 5 deletions(-) diff --git a/src/JetReconstruction.jl b/src/JetReconstruction.jl index 56c769e..5ee693e 100644 --- a/src/JetReconstruction.jl +++ b/src/JetReconstruction.jl @@ -5,14 +5,18 @@ module JetReconstruction using LorentzVectorHEP -# particle type definition -# include("Particle.jl") -# export energy, px, py, pz, pt, phi, mass, eta, kt, ϕ, η +# 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) + +pt2(p::LorentzVectorCyl) = LorentzVectorHEP.pt2(p) +phi(p::LorentzVectorCyl) = LorentzVectorHEP.phi(p) +rapidity(p::LorentzVectorCyl) = LorentzVectorHEP.rapidity(p) # Philipp's pseudojet include("Pseudojet.jl") -## As this is an internal EDM class, we perhaps shouldn't export this stuff... -# export PseudoJet, rap, phi, pt, pt2, px, py, pz, pt, phi, mass, eta +## As this is an internal EDM class, we don't export anything # Simple HepMC3 reader include("HepMC3.jl") From 1da44bbe2ee823c40db59073be888b4f624072ae Mon Sep 17 00:00:00 2001 From: Graeme A Stewart Date: Sun, 22 Oct 2023 21:17:08 +0200 Subject: [PATCH 08/21] Add algorithm wrapper and refactor variables Add an outer "wrapper" that writes the internal EDM objects and allow for a direct call when the internal EDM is already constructed (this is fairer for timings) Refactor many of the variable names to get rid of the underscore prefixes and to make it clearer when variables are arrays, rather than single values --- src/Algo.jl | 193 ++++++++++++++++++++++++++++------------------------ 1 file changed, 103 insertions(+), 90 deletions(-) diff --git a/src/Algo.jl b/src/Algo.jl index ad824e4..34096a1 100644 --- a/src/Algo.jl +++ b/src/Algo.jl @@ -1,128 +1,143 @@ -Base.@propagate_inbounds function _dist(i, j, _eta, _phi) - deta = _eta[i] - _eta[j] - dphi = abs(_phi[i] - _phi[j]) +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(deta, deta, 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, _nn, _nndist) - j = _nn[i] - d = _nndist[i] - d * min(_kt2[i], _kt2[j]) +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, _eta, _phi, _R2, _nndist, _nn) - nndist = _R2 - nn = i +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, _eta, _phi) - if Δ2 < nndist - nn = j - nndist = Δ2 + Δ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 + if Δ2 < nndist[j] + nndist[j] = Δ2 + nn[j] = i end end - _nndist[i] = nndist - _nn[i] = nn + 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, _eta, _phi, _R2, _nndist, _nn) - nndist = _R2 - nn = 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, _eta, _phi) - if Δ2 <= nndist - nn = j - nndist = Δ2 + Δ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, _eta, _phi) - f = Δ2 <= nndist - nn = ifelse(f, j, nn) - nndist = ifelse(f, Δ2, nndist) + Δ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 - _nn[i] = nn + nndist[i] = nndist_min + nn[i] = nn_min end # entire NN update step -# Base.@propagate_inbounds function _upd_nn_step!(i::Int, j::Int, k::Int, N::Int, Nn::Int, -# _kt2::Vector{Float64}, _eta::Vector{Float64}, _phi::Vector{Float64}, _R2::Float64, _nndist::Vector{Float64}, _nn::Vector{Int}, _nndij::Vector{Float64}) -Base.@propagate_inbounds function _upd_nn_step!(i, j, k, N, Nn, _kt2, _eta, _phi, _R2, _nndist, _nn, _nndij) - nnk = _nn[k] +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, _eta, _phi, _R2, _nndist, _nn) # update dist and nn - _nndij[k] = _dij(k, _kt2, _nn, _nndist) - nnk = _nn[k] + 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, _eta, _phi) - if Δ2 < _nndist[k] - _nndist[k] = Δ2 - nnk = _nn[k] = i - _nndij[k] = _dij(k, _kt2, _nn, _nndist) + Δ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])) + cond = Δ2 < nndist[i] + nndist[i], nn[i] = ifelse(cond, (Δ2, k), (nndist[i], nn[i])) end - nnk == Nn && (_nn[k] = j) + nnk == Nn && (nn[k] = j) end +""" +This is the N2Plain jet reconstruction algorithm interface, called with an arbitrary array +of objects, which supports the methods pt2(), phi(), rapidity() for each element. +""" +function sequential_jet_reconstruct(objects::AbstractArray{T}; p = -1, R = 1.0, recombine = +) where T + # Integer p if possible + p = (round(p) == p) ? Int(p) : p + + # We make sure these arrays are type stable - have seen issues where, depending on the values + # returned by the methods, they can become unstable and performance degrades + kt2_array::Vector{Float64} = pt2.(objects) .^ p + phi_array::Vector{Float64} = phi.(objects) + rapidity_array::Vector{Float64} = rapidity.(objects) + + objects_array = copy(objects) + + # Now call the actual reconstruction method, tuned for our internal EDM + sequential_jet_reconstruct(objects_array=objects_array, kt2_array=kt2_array, phi_array=phi_array, + rapidity_array=rapidity_array, p=p, R=R, recombine=recombine) +end + + -function sequential_jet_reconstruct(objects::AbstractArray{T}; p = -1.0, R = 1.0, recombine = +) where T +function sequential_jet_reconstruct(;objects_array::AbstractArray{J}, kt2_array::Vector{F}, + phi_array::Vector{F}, rapidity_array::Vector{F}, p = -1, R = 1.0, recombine = +) where {J, F<:AbstractFloat} # Bounds - N::Int = length(objects) + N::Int = length(objects_array) # Returned values - jets = T[] + jets = J[] sequences = Vector{Int}[] # recombination sequences, WARNING: first index in the sequence is not necessarily the seed # Parameters - _R2 = R * R - _p = (round(p) == p) ? Int(p) : p # integer p if possible + R2 = R * R # Data - _objects = copy(objects) - ## When working with LorentzVectorHEP we make sure these arrays are type stable - _kt2::Vector{Float64} = LorentzVectorHEP.pt2.(_objects) .^ _p - _phi::Vector{Float64} = LorentzVectorHEP.phi.(_objects) - _eta::Vector{Float64} = LorentzVectorHEP.rapidity.(_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] + 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) + upd_nn_crosscheck!(i, 1, i - 1, rapidity_array, phi_array, R2, nndist, nn) end # diJ table *_R2 - _nndij::Vector{Float64} = zeros(N) + nndij::Vector{Float64} = zeros(N) @inbounds @simd for i in 1:N - _nndij[i] = _dij(i, _kt2, _nn, _nndist) + nndij[i] = dij(i, kt2_array, nn, nndist) end iteration::Int = 1 while N != 0 # findmin i = 1 - dij_min = _nndij[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)) + cond = nndij[k] < dij_min + dij_min, i = ifelse(cond, (nndij[k], k), (dij_min, i)) end - j::Int = _nn[i] + j::Int = nn[i] ## Needed for certain tricky debugging situations # if iteration==1 @@ -136,35 +151,33 @@ function sequential_jet_reconstruct(objects::AbstractArray{T}; p = -1.0, R = 1.0 end # update ith jet, replacing it with the new one - _objects[i] = recombine(_objects[i], _objects[j]) - _phi[i] = LorentzVectorHEP.phi(_objects[i]) - _eta[i] = LorentzVectorHEP.rapidity(_objects[i]) - _kt2[i] = LorentzVectorHEP.pt2(_objects[i])^_p + objects_array[i] = recombine(objects_array[i], objects_array[j]) + phi_array[i] = LorentzVectorHEP.phi(objects_array[i]) + rapidity_array[i] = LorentzVectorHEP.rapidity(objects_array[i]) + kt2_array[i] = LorentzVectorHEP.pt2(objects_array[i]) ^ p - _nndist[i] = _R2 - _nn[i] = i + 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) + @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, LorentzVectorCyl(LorentzVectorHEP.pt(_objects[i]), LorentzVectorHEP.eta(_objects[i]), - # LorentzVectorHEP.phi(_objects[i]), LorentzVectorHEP.mt(_objects[i]))) - push!(jets, _objects[i]) - push!(sequences, _sequences[i]) + push!(jets, objects_array[i]) + push!(sequences, sequences[i]) end # copy jet N to j - _objects[j] = _objects[N] + objects_array[j] = objects_array[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] + 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] + sequences[j] = sequences[N] Nn::Int = N N -= 1 @@ -172,11 +185,11 @@ function sequential_jet_reconstruct(objects::AbstractArray{T}; p = -1.0, R = 1.0 # 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) + 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, _nn, _nndist) + nndij[i] = dij(i, kt2_array, nn, nndist) end jets, sequences From 09e0e463c3b84cc993abde56037218f3b8d8be67 Mon Sep 17 00:00:00 2001 From: Graeme A Stewart Date: Tue, 24 Oct 2023 16:57:30 +0200 Subject: [PATCH 09/21] Benchmark with internal EDM for N2Tiled The conversion from LorentzVectorHEP to PseudoJet is rather expensive and should be factored out of the benchmark for a fair comparison of the actual algorithm --- chep.jl | 30 ++++++++++++++++++------------ src/TiledAlgoLL.jl | 9 +++++++-- 2 files changed, 25 insertions(+), 14 deletions(-) diff --git a/chep.jl b/chep.jl index 53e6236..868cd34 100644 --- a/chep.jl +++ b/chep.jl @@ -78,14 +78,20 @@ function jet_process( throw(ErrorException("Strategy not yet implemented")) end - # Change of tactic - we adapt all interfaces to swallow LorentzVectors directly - # The N2Tiled algorithm uses PseudoJets so pass these directly - # if strategy == N2Tiled - # event_vector = events - # else - # # The other algorithms swallow 4-vectors instead - # event_vector = pseudojets2vectors(events) - # end + # Build internal EDM structures for timing measurements + if strategy == N2Tiled + event_vector = Vector{Vector{JetReconstruction.PseudoJet}}(undef, 0) + for event in events + pseudojet_event = Vector{JetReconstruction.PseudoJet}(undef, length(event)) + for (i, particle) in enumerate(event) + pseudojet_event[i] = JetReconstruction.PseudoJet(LorentzVectorHEP.px(particle), LorentzVectorHEP.py(particle), + LorentzVectorHEP.pz(particle), LorentzVectorHEP.energy(particle)) + end + push!(event_vector, pseudojet_event) + end + elseif strategy == N2Plain + event_vector = events + end # If we are dumping the results, setup the JSON structure if !isnothing(dump) @@ -95,20 +101,20 @@ function jet_process( # Warmup code if we are doing a multi-sample timing run if nsamples > 1 || profile @info "Doing initial warm-up run" - for event in events + for event in event_vector finaljets, _ = jet_reconstruction(event, R = distance, p = power) final_jets(finaljets, ptmin) end end if profile - profile_code(jet_reconstruction, events, nsamples) + profile_code(jet_reconstruction, event_vector, nsamples) return nothing end if alloc println("Memory allocation statistics:") - @timev for event in events + @timev for event in event_vector finaljets, _ = jet_reconstruction(event, R = distance, p = power) final_jets(finaljets, ptmin) end @@ -123,7 +129,7 @@ function jet_process( for irun ∈ 1:nsamples gcoff && GC.enable(false) t_start = time_ns() - for (ievt, event) in enumerate(events) + for (ievt, event) in enumerate(event_vector) finaljets, _ = jet_reconstruction(event, R = distance, p = power) fj = final_jets(finaljets, ptmin) # Only print the jet content once diff --git a/src/TiledAlgoLL.jl b/src/TiledAlgoLL.jl index 37f1bec..01a394d 100644 --- a/src/TiledAlgoLL.jl +++ b/src/TiledAlgoLL.jl @@ -295,8 +295,11 @@ function inclusive_jets(clusterseq::ClusterSequence, ptmin = 0.0) jets_local end + """ -Main jet reconstruction algorithm +Main jet reconstruction algorithm entry point + +This function will construct the internal EDM of PseudoJets from LorentzVectors """ function tiled_jet_reconstruct_ll(particles::Vector{LorentzVector}; p = -1, R = 1.0, recombine = +, ptmin = 0.0) # Here we need to populate the vector of PseudoJets that are the internal @@ -309,7 +312,9 @@ function tiled_jet_reconstruct_ll(particles::Vector{LorentzVector}; p = -1, R = tiled_jet_reconstruct_ll(pseudojets, p = p, R = R, recombine = recombine, ptmin = ptmin) end - +""" +Main jet reconstruction algorithm, using PseudoJets +""" function tiled_jet_reconstruct_ll(particles::Vector{PseudoJet}; p = -1, R = 1.0, recombine = +, ptmin = 0.0) # Bounds N::Int = length(particles) From 24069d8f7060acad24f09b4eecea20c6ef4f4758 Mon Sep 17 00:00:00 2001 From: Graeme A Stewart Date: Tue, 24 Oct 2023 21:43:53 +0200 Subject: [PATCH 10/21] Import additional methods from LorentzVectorHEP We want to have the full suite of relevant methods available in our namespace --- src/JetReconstruction.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/src/JetReconstruction.jl b/src/JetReconstruction.jl index 5ee693e..09071b5 100644 --- a/src/JetReconstruction.jl +++ b/src/JetReconstruction.jl @@ -9,10 +9,18 @@ using LorentzVectorHEP 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") From d7f1da1f409f8298e95131d3cd2076181731e20c Mon Sep 17 00:00:00 2001 From: Graeme A Stewart Date: Tue, 24 Oct 2023 21:44:24 +0200 Subject: [PATCH 11/21] Remove erroneous namespace from methods --- src/Algo.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Algo.jl b/src/Algo.jl index 34096a1..4685b19 100644 --- a/src/Algo.jl +++ b/src/Algo.jl @@ -152,9 +152,9 @@ function sequential_jet_reconstruct(;objects_array::AbstractArray{J}, kt2_array: # update ith jet, replacing it with the new one objects_array[i] = recombine(objects_array[i], objects_array[j]) - phi_array[i] = LorentzVectorHEP.phi(objects_array[i]) - rapidity_array[i] = LorentzVectorHEP.rapidity(objects_array[i]) - kt2_array[i] = LorentzVectorHEP.pt2(objects_array[i]) ^ p + 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 From 1e2df3bce49cab7cd2f54be24fec7b40392277ed Mon Sep 17 00:00:00 2001 From: Graeme A Stewart Date: Wed, 25 Oct 2023 09:04:40 +0200 Subject: [PATCH 12/21] Update PseudoJet method name Change rap() to rapidity() for consistency --- src/Pseudojet.jl | 2 +- src/TiledAlgoLL.jl | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Pseudojet.jl b/src/Pseudojet.jl index d90a0f0..4a90948 100644 --- a/src/Pseudojet.jl +++ b/src/Pseudojet.jl @@ -125,7 +125,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 diff --git a/src/TiledAlgoLL.jl b/src/TiledAlgoLL.jl index 01a394d..04dea8e 100644 --- a/src/TiledAlgoLL.jl +++ b/src/TiledAlgoLL.jl @@ -64,7 +64,7 @@ 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.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 @@ -288,7 +288,7 @@ function inclusive_jets(clusterseq::ClusterSequence, ptmin = 0.0) iparent_jet = clusterseq.history[elt.parent1].jetp_index jet = clusterseq.jets[iparent_jet] if pt2(jet) >= dcut - push!(jets_local, LorentzVectorCyl(pt(jet), rap(jet), phi(jet), mass(jet))) + push!(jets_local, LorentzVectorCyl(pt(jet), rapidity(jet), phi(jet), mass(jet))) # push!(jets_local, jet) end end @@ -343,7 +343,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)) From 098307dc5f88c2f17c9c0abd7493420a699cd232 Mon Sep 17 00:00:00 2001 From: Graeme A Stewart Date: Wed, 25 Oct 2023 10:06:17 +0200 Subject: [PATCH 13/21] Add pseudorapidity and accessors to PseudoJet This makes the struct more generic and avoids directly accessing member variables --- src/Pseudojet.jl | 53 ++++++++++++++++++++++++------------------------ 1 file changed, 27 insertions(+), 26 deletions(-) diff --git a/src/Pseudojet.jl b/src/Pseudojet.jl index 4a90948..bd88e7a 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 @@ -138,17 +115,41 @@ 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) +mass2 = m2 +energy(p::PseudoJet) = p.E import Base.+; From 645ac095be68912e5793a504f4d1085a538ad30a Mon Sep 17 00:00:00 2001 From: Graeme A Stewart Date: Wed, 25 Oct 2023 10:13:31 +0200 Subject: [PATCH 14/21] Make N2Tiled type agnostic Changing back to PseudoJets, make the generic method for N2Tiled type agnostic (anything that supports energy(), px(), py(), pz() is fine) --- src/TiledAlgoLL.jl | 14 ++++++++------ 1 file changed, 8 insertions(+), 6 deletions(-) diff --git a/src/TiledAlgoLL.jl b/src/TiledAlgoLL.jl index 04dea8e..1fadaad 100644 --- a/src/TiledAlgoLL.jl +++ b/src/TiledAlgoLL.jl @@ -297,23 +297,25 @@ end """ -Main jet reconstruction algorithm entry point +Main jet reconstruction algorithm entry point for generic data types -This function will construct the internal EDM of PseudoJets from LorentzVectors +`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. """ -function tiled_jet_reconstruct_ll(particles::Vector{LorentzVector}; p = -1, R = 1.0, recombine = +, ptmin = 0.0) +function tiled_jet_reconstruct_ll(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(LorentzVectorHEP.px(particle), LorentzVectorHEP.py(particle), - LorentzVectorHEP.pz(particle), LorentzVectorHEP.energy(particle)) + pseudojets[i] = PseudoJet(px(particle), py(particle), + pz(particle), energy(particle)) end tiled_jet_reconstruct_ll(pseudojets, p = p, R = R, recombine = recombine, ptmin = ptmin) end """ -Main jet reconstruction algorithm, using PseudoJets +Main jet reconstruction algorithm, using PseudoJet objects """ function tiled_jet_reconstruct_ll(particles::Vector{PseudoJet}; p = -1, R = 1.0, recombine = +, ptmin = 0.0) # Bounds From d8bcc4424cec5db4edefc02e7356f2f5e908dac6 Mon Sep 17 00:00:00 2001 From: Graeme A Stewart Date: Wed, 25 Oct 2023 10:25:32 +0200 Subject: [PATCH 15/21] Move back to using PseudoJet structs As N2Plain is written, it's completely happy to use PseudoJets directly --- chep.jl | 26 +++++++++----------------- src/JetReconstruction.jl | 1 + src/Utils.jl | 2 +- 3 files changed, 11 insertions(+), 18 deletions(-) diff --git a/chep.jl b/chep.jl index 868cd34..c6d83c2 100644 --- a/chep.jl +++ b/chep.jl @@ -56,7 +56,7 @@ function profile_code(jet_reconstruction, events, niters) end function jet_process( - events::Vector{Vector{LorentzVector}}; + events::Vector{Vector{PseudoJet}}; ptmin::Float64 = 5.0, distance::Float64 = 0.4, power::Integer = -1, @@ -78,18 +78,10 @@ function jet_process( throw(ErrorException("Strategy not yet implemented")) end - # Build internal EDM structures for timing measurements - if strategy == N2Tiled - event_vector = Vector{Vector{JetReconstruction.PseudoJet}}(undef, 0) - for event in events - pseudojet_event = Vector{JetReconstruction.PseudoJet}(undef, length(event)) - for (i, particle) in enumerate(event) - pseudojet_event[i] = JetReconstruction.PseudoJet(LorentzVectorHEP.px(particle), LorentzVectorHEP.py(particle), - LorentzVectorHEP.pz(particle), LorentzVectorHEP.energy(particle)) - end - push!(event_vector, pseudojet_event) - end - elseif strategy == N2Plain + # 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 end @@ -276,10 +268,10 @@ main() = begin logger = ConsoleLogger(stdout, Logging.Warn) end global_logger(logger) - # events::Vector{Vector{PseudoJet}} = - # read_final_state_particles(args[:file], maxevents = args[:maxevents], skipevents = args[:skip]) - events::Vector{Vector{LorentzVector}} = - read_final_state_particles_lv(args[:file], maxevents = args[:maxevents], skipevents = args[:skip]) + events::Vector{Vector{PseudoJet}} = + read_final_state_particles(args[:file], maxevents = args[:maxevents], skipevents = args[:skip]) + # events::Vector{Vector{LorentzVector}} = + # read_final_state_particles_lv(args[:file], maxevents = args[:maxevents], skipevents = args[:skip]) jet_process(events, ptmin = args[:ptmin], distance = args[:distance], power = args[:power], strategy = args[:strategy], nsamples = args[:nsamples], gcoff = args[:gcoff], profile = args[:profile], diff --git a/src/JetReconstruction.jl b/src/JetReconstruction.jl index 09071b5..c916bb7 100644 --- a/src/JetReconstruction.jl +++ b/src/JetReconstruction.jl @@ -24,6 +24,7 @@ energy(p::LorentzVectorCyl) = LorentzVectorHEP.energy(p) # Philipp's pseudojet include("Pseudojet.jl") +export PseudoJet ## As this is an internal EDM class, we don't export anything # Simple HepMC3 reader diff --git a/src/Utils.jl b/src/Utils.jl index 8f209a6..c9c34a6 100644 --- a/src/Utils.jl +++ b/src/Utils.jl @@ -94,7 +94,7 @@ 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 From dd02e8d06d2841715f179ad8197a15792c73b974 Mon Sep 17 00:00:00 2001 From: Graeme A Stewart Date: Wed, 25 Oct 2023 10:40:37 +0200 Subject: [PATCH 16/21] Correctly pass recombine function for N2Tiled Anything non-standard will require recombine(::PseudoJet, ::PseudoJet) to be defined. --- src/TiledAlgoLL.jl | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/TiledAlgoLL.jl b/src/TiledAlgoLL.jl index 1fadaad..2dcf379 100644 --- a/src/TiledAlgoLL.jl +++ b/src/TiledAlgoLL.jl @@ -147,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) @@ -302,6 +302,9 @@ 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_ll(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 @@ -392,7 +395,7 @@ function tiled_jet_reconstruct_ll(particles::Vector{PseudoJet}; p = -1, R = 1.0, 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... From 71c8bec5950cd65e165ddb6ae002a17db57e53b6 Mon Sep 17 00:00:00 2001 From: Graeme A Stewart Date: Wed, 25 Oct 2023 10:53:24 +0200 Subject: [PATCH 17/21] Add ptmin filtering to N2Plain To make the interface the same as N2Tiled support directly filtering on the final jets' ptmin values (which actually speeds up the code, slightly!) --- chep.jl | 2 +- src/Algo.jl | 13 ++++++++----- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/chep.jl b/chep.jl index c6d83c2..a57de14 100644 --- a/chep.jl +++ b/chep.jl @@ -122,7 +122,7 @@ function jet_process( 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 diff --git a/src/Algo.jl b/src/Algo.jl index 4685b19..050281f 100644 --- a/src/Algo.jl +++ b/src/Algo.jl @@ -80,7 +80,7 @@ end This is the N2Plain jet reconstruction algorithm interface, called with an arbitrary array of objects, which supports the methods pt2(), phi(), rapidity() for each element. """ -function sequential_jet_reconstruct(objects::AbstractArray{T}; p = -1, R = 1.0, recombine = +) where T +function sequential_jet_reconstruct(objects::AbstractArray{T}; p = -1, R = 1.0, recombine = +, ptmin = 0.0) where T # Integer p if possible p = (round(p) == p) ? Int(p) : p @@ -94,13 +94,13 @@ function sequential_jet_reconstruct(objects::AbstractArray{T}; p = -1, R = 1.0, # Now call the actual reconstruction method, tuned for our internal EDM sequential_jet_reconstruct(objects_array=objects_array, kt2_array=kt2_array, phi_array=phi_array, - rapidity_array=rapidity_array, p=p, R=R, recombine=recombine) + rapidity_array=rapidity_array, p=p, R=R, recombine=recombine, ptmin=ptmin) end function sequential_jet_reconstruct(;objects_array::AbstractArray{J}, kt2_array::Vector{F}, - phi_array::Vector{F}, rapidity_array::Vector{F}, p = -1, R = 1.0, recombine = +) where {J, F<:AbstractFloat} + 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) @@ -109,7 +109,8 @@ function sequential_jet_reconstruct(;objects_array::AbstractArray{J}, kt2_array: sequences = Vector{Int}[] # recombination sequences, WARNING: first index in the sequence is not necessarily the seed # Parameters - R2 = R * R + R2 = R^2 + ptmin2 = ptmin^2 # Data nn = Vector(1:N) # nearest neighbours @@ -163,7 +164,9 @@ function sequential_jet_reconstruct(;objects_array::AbstractArray{J}, kt2_array: push!(sequences[i], x) end else # i == j - push!(jets, objects_array[i]) + if (pt2(objects_array[i]) >= ptmin2) + push!(jets, objects_array[i]) + end push!(sequences, sequences[i]) end From 547107cf750690c84818f8169484b88eaad75671 Mon Sep 17 00:00:00 2001 From: Graeme A Stewart Date: Wed, 25 Oct 2023 11:08:01 +0200 Subject: [PATCH 18/21] Return PseudoJets always Even from N2Plain we return PseudoJet types (but avoid construction of new objects if possible) --- src/Algo.jl | 9 +++++++-- test/runtests.jl | 12 ++---------- 2 files changed, 9 insertions(+), 12 deletions(-) diff --git a/src/Algo.jl b/src/Algo.jl index 050281f..d189188 100644 --- a/src/Algo.jl +++ b/src/Algo.jl @@ -105,7 +105,7 @@ function sequential_jet_reconstruct(;objects_array::AbstractArray{J}, kt2_array: N::Int = length(objects_array) # Returned values - jets = J[] + jets = PseudoJet[] sequences = Vector{Int}[] # recombination sequences, WARNING: first index in the sequence is not necessarily the seed # Parameters @@ -165,7 +165,12 @@ function sequential_jet_reconstruct(;objects_array::AbstractArray{J}, kt2_array: end else # i == j if (pt2(objects_array[i]) >= ptmin2) - push!(jets, objects_array[i]) + # 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 diff --git a/test/runtests.jl b/test/runtests.jl index 93a4f1c..9be7d4e 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -54,18 +54,10 @@ function do_jet_test(strategy::JetRecoStrategy, fastjet_jets; end # Now run our jet reconstruction... - events::Vector{Vector{LorentzVector}} = read_final_state_particles_lv("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) + events::Vector{Vector{PseudoJet}} = read_final_state_particles("test/data/events.hepmc3") jet_collection = FinalJets[] for (ievt, event) in enumerate(events) - finaljets, _ = jet_reconstruction(event, R=distance, p=power) + finaljets, _ = jet_reconstruction(event, R=distance, p=power, ptmin=ptmin) fj = final_jets(finaljets, ptmin) sort_jets!(fj) push!(jet_collection, FinalJets(ievt, fj)) From 9ce9346dddcef4abb70750a871aa42c5c6e2390e Mon Sep 17 00:00:00 2001 From: Graeme A Stewart Date: Wed, 25 Oct 2023 13:24:20 +0200 Subject: [PATCH 19/21] Add tests for LorentzVector input Check that reconstruction from LorentzVector gives the same as PseudoJet --- src/Particle.jl | 68 ------------------------------------------------ test/runtests.jl | 68 ++++++++++++++++-------------------------------- 2 files changed, 23 insertions(+), 113 deletions(-) delete mode 100644 src/Particle.jl 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/test/runtests.jl b/test/runtests.jl index 9be7d4e..5bb3846 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -32,9 +32,6 @@ function main() # Test each stratgy... do_jet_test(N2Plain, fastjet_jets) do_jet_test(N2Tiled, fastjet_jets) - - # Atell's original test - # original_tests() end function do_jet_test(strategy::JetRecoStrategy, fastjet_jets; @@ -54,6 +51,7 @@ function do_jet_test(strategy::JetRecoStrategy, fastjet_jets; end # Now run our jet reconstruction... + # From PseudoJets events::Vector{Vector{PseudoJet}} = read_final_state_particles("test/data/events.hepmc3") jet_collection = FinalJets[] for (ievt, event) in enumerate(events) @@ -63,7 +61,17 @@ function do_jet_test(strategy::JetRecoStrategy, fastjet_jets; 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 @@ -85,53 +93,23 @@ function do_jet_test(strategy::JetRecoStrategy, fastjet_jets; 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__ From 7b1b3df830c1ae352ec4066fb5a89f52a0dd537b Mon Sep 17 00:00:00 2001 From: Graeme A Stewart Date: Wed, 25 Oct 2023 14:38:34 +0200 Subject: [PATCH 20/21] Pre-release cleanup Write a proper README file that describes how to use the package (real Julia documentation to come, but it's a start) Rename main algorithms to be more consistent, viz. sequential_jet_reconstruct and tiled_jet_reconstruct Also remove all of the pre-cooked versions with specific values of p (like cambridge_aachen_algo) as these add little value and would need to exist for all algorithms. Rename the chep.jl script to jetreco.jl and put it into an examples directory --- Project.toml | 3 +- README.md | 93 +++++++++++++++++++++++---- anti-kt.pseudo => docs/anti-kt.pseudo | 0 chep.jl => examples/jetreco.jl | 6 +- main.jl | 79 ----------------------- src/JetReconstruction.jl | 8 +-- src/{Algo.jl => PlainAlgo.jl} | 61 ++++-------------- src/TiledAlgoLL.jl | 6 +- test/runtests.jl | 4 +- 9 files changed, 107 insertions(+), 153 deletions(-) rename anti-kt.pseudo => docs/anti-kt.pseudo (100%) rename chep.jl => examples/jetreco.jl (96%) delete mode 100644 main.jl rename src/{Algo.jl => PlainAlgo.jl} (71%) diff --git a/Project.toml b/Project.toml index 661b475..d1cd28e 100644 --- a/Project.toml +++ b/Project.toml @@ -8,7 +8,6 @@ Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63" Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" FlameGraphs = "08572546-2f56-4bcf-ba4e-bab62c3a3f89" -HepMC3_jll = "b85c3e40-22db-5268-bacb-02bd65cb4e01" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" Logging = "56ddb016-857b-54e1-b83d-db4d58db5568" LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890" @@ -19,8 +18,8 @@ ProfileSVG = "132c30aa-f267-4189-9183-c8a63c7e05e6" StatProfilerHTML = "a8a75453-ed82-57c9-9e16-4cd1196ecbf5" [compat] -julia = "1.8" LorentzVectorHEP = "0.1.6" +julia = "1.8" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/README.md b/README.md index 8d9fa2a..e2c46df 100644 --- a/README.md +++ b/README.md @@ -1,21 +1,92 @@ -# JetReconstruction +# JetReconstruction.jl -[![Build Status](https://github.com/gojakuch/JetReconstruction.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/gojakuch/JetReconstruction.jl/actions/workflows/CI.yml?query=branch%3Amain) +[![Build Status](https://github.com/JuliaHEP/JetReconstruction.jl/actions/workflows/CI.yml/badge.svg?branch=main)](https://github.com/JuliaHEP/JetReconstruction.jl/actions/workflows/CI.yml?query=branch%3Amain) + +## This package implements sequential Jet Reconstruction (clustering) + +### Algorithms + +Algorithms used are based on the C++ FastJet package (, +[hep-ph/0512210](https://arxiv.org/abs/hep-ph/0512210), +[arXiv:1111.6097](https://arxiv.org/abs/1111.6097)), reimplemented natively in Julia. + +One algorithm is the plain N2 approach `N2Plain`, the other uses the FastJet tiling +approach, `N2Tiled`. + +### Interface + +To use the package one should call the appropriate API interface of the desired algorithm -This code is not a registered Julia package yet. If you want to use it, the valid option would be to change the `main.jl` file how you want and run it. That file should either start with ```julia -using Revise; import Pkg # you need to have Revise installed -Pkg.activate(".") -using JetReconstruction +# N2Plain +plain_jet_reconstruct(particles::Vector{T}; p = -1, R = 1.0, recombine = +, ptmin = 0.0) ``` -or with + ```julia -include("src/JetReconstruction.jl") -using .JetReconstruction +# N2Tiled +tiled_jet_reconstruct(particles::Vector{T}; p = -1, R = 1.0, recombine = +, ptmin = 0.0) ``` -There is some documentation provided for functions and submodules. Once everything is working properly and efficiently, this `README.md` will contain more details on usage and the module might get registered. -## Plotting +- `particles` - a vector of input particles for the clustering + - Any type that supplies the methods `pt2()`, `phi()`, `rapidity()`, `px()`, `py()`, `pz()`, `energy()` can be used + - These methods have to be defined in the namespace of this package, i.e., `JetReconstruction.pt2(::T)` + - The `PseudoJet` type from this package, or a 4-vector from `LorentzVectorHEP` are suitable (and have the appropriate definitions) +- `p` - the transverse momentum power used in the $d_{ij}$ metric for deciding on closest jets, as $k^{2p}_\text{T}$ + - `-1` gives anti-${k}_\text{T}$ clustering + - `0` gives Cambridge/Achen + - `1` gives inclusive $k_\text{T}$ +- `R` - the cone size parameter (no particles more geometrically distance than `R` will be merged) +- `recombine` - the function used to merge two pseudojets (by default this is simple 4-vector addition of $(E, \mathbf{p})$) +- `ptmin` - only jets of transverse momentum greater than or equal to `ptmin` will be returned + +Both of these functions return tuple of `(jets, seq)`, where these are a vector of `JetReconstruction.PseudoJet` objects and cluster sequencing information, respectively. + +Note that the `N2Plain` algorithm is usually faster at low densities, $\lessapprox 50$ input particles, otherwise the tiled algorithm is faster. + +### Example + +See the `examples/jetreco.jl` script for a full example of how to call jet reconstruction. + +```julia +julia --project=. examples/jetreco.jl --maxevents=100 --nsamples=1 --strategy=N2Plain test/data/events.hepmc3 +... +julia --project=. examples/jetreco.jl --maxevents=100 --nsamples=1 --strategy=N2Tiled test/data/events.hepmc3 +... +``` + +The example also shows how to use `JetReconstruction.HepMC3` to read HepMC3 ASCII files (via the `read_final_state_particles()` wrapper). + +### Plotting + +**TO BE FIXED** + ![illustration](img/illustration.jpeg) To visualise the clustered jets as a 3d bar plot (see illustration above) we now use `Makie.jl`. See the `jetsplot` function and its documentation for more. + +## Reference + +Although it has been developed further since the CHEP2023 conference, the CHEP conference proceedings, [arXiv:2309.17309](https://arxiv.org/abs/2309.17309), should be cited if you use this package: + +```bibtex +@misc{stewart2023polyglot, + title={Polyglot Jet Finding}, + author={Graeme Andrew Stewart and Philippe Gras and Benedikt Hegner and Atell Krasnopolski}, + year={2023}, + eprint={2309.17309}, + archivePrefix={arXiv}, + primaryClass={hep-ex} +} +``` + +## Authors and Copyright + +Code in this package is authored by: + +- Atell Krasnopolski +- Graeme A Stewart +- Philippe Gras + +and is Copyright 2022-2023 The Authors, CERN. + +The code is under the MIT License. diff --git a/anti-kt.pseudo b/docs/anti-kt.pseudo similarity index 100% rename from anti-kt.pseudo rename to docs/anti-kt.pseudo diff --git a/chep.jl b/examples/jetreco.jl similarity index 96% rename from chep.jl rename to examples/jetreco.jl index a57de14..0f092ad 100644 --- a/chep.jl +++ b/examples/jetreco.jl @@ -71,9 +71,9 @@ function jet_process( # Strategy if (strategy == N2Plain) - jet_reconstruction = sequential_jet_reconstruct + jet_reconstruction = plain_jet_reconstruct elseif (strategy == N2Tiled || stragegy == Best) - jet_reconstruction = tiled_jet_reconstruct_ll + jet_reconstruction = tiled_jet_reconstruct else throw(ErrorException("Strategy not yet implemented")) end @@ -270,8 +270,6 @@ main() = begin global_logger(logger) events::Vector{Vector{PseudoJet}} = read_final_state_particles(args[:file], maxevents = args[:maxevents], skipevents = args[:skip]) - # events::Vector{Vector{LorentzVector}} = - # read_final_state_particles_lv(args[:file], maxevents = args[:maxevents], skipevents = args[:skip]) jet_process(events, ptmin = args[:ptmin], distance = args[:distance], power = args[:power], strategy = args[:strategy], nsamples = args[:nsamples], gcoff = args[:gcoff], profile = args[:profile], diff --git a/main.jl b/main.jl deleted file mode 100644 index 35ad608..0000000 --- a/main.jl +++ /dev/null @@ -1,79 +0,0 @@ -# choose one: -# this -#= -include("src/JetReconstruction.jl") -using .JetReconstruction -=# -# or this -using Revise; import Pkg -Pkg.activate(".") -using JetReconstruction - -using CairoMakie - -## array constructor -makejet(pt, ϕ, y) = [pt*cos(ϕ), pt*sin(ϕ), pt*sinh(y), sqrt((pt*cos(ϕ))^2 + (pt*sin(ϕ))^2 + (pt*sinh(y))^2)] - -## choose a data sample and save it -smalldata = [ - makejet(25, 1.2, 0), - makejet(0.1, 1, 0), - makejet(20, 0, 0) -] -smalldata = [ - makejet(25, 1.2, 0), - makejet(0.1, 0.6, 0), - makejet(20, 0, 0) -] -smalldata = [ - makejet(5, 0.8, 0), - makejet(2, 0, 0.8), -] -smalldata = [ - makejet(30, -0.5, -2.5), - makejet(30, -0.5, -3.5) -] -smalldata = [ - makejet(30, 0.5, -2.5), - makejet(30, 0.5, -3.5), - makejet(30, -0.5, -2.5), - makejet(30, -0.5, -3.5), - makejet(0.1, 0, -3), -] -smalldata = [ - [0.2839991726, -0.4668202293, -49.6978846131, 49.709744138], - [0.1434555273, -0.0312880942, -6.9382351613, 6.9411919273], -] -smalldata = [ - makejet(0.8, 0.8, 0.1), - makejet(0.1, 0.1, 0.8), -] - -smalldata = [ - [-0.8807412236, -1.2331262152, -157.431315651, 157.4393822839], - [-0.0051712611, 0.23815508, -9.7396045662, 9.7435168946], # 2 - [0.0362280943, 0.2694752057, -6.9243427525, 6.9310844534], # 3 - [-0.2206628664, -0.1438198985, -0.6838608666, 0.7460038429], - [1.2716787521, 1.0422298083, -6.1740167274, 6.3907254797], # 5 - [-0.5695590845, -0.3627761836, -58.5430479911, 58.5544811606], - [0.2839991726, -0.4668202293, -49.6978846131, 49.709744138], - [0.6510530003, 1.3970949413, -62.7226079598, 62.7485783532], # 8 - [0.1434555273, -0.0312880942, -6.9382351613, 6.9411919273], - [0.4931562547, 2.1627817414, -14.8865871635, 15.0516040711], # 10 - [0.2396813608, -0.0786236784, -1.9340954697, 1.9554625817], - [0.3355486441, 0.0516402769, -0.8346540063, 0.9118040941], - [-0.7853865645, -0.7810520475, -1.5367790662, 1.8994852039], -] -smalldata = [ - [1.2716787521, 1.0422298083, -6.1740167274, 6.3907254797], - [0.4931562547, 2.1627817414, -14.8865871635, 15.0516040711], - [0.3355486441, 0.0516402769, -0.8346540063, 0.9118040941], -] - -savejets("small.dat", smalldata, format="px py pz E") -## run -smalljets, smallind = anti_kt_algo(smalldata, R=1); img = jetsplot(smalldata, smallind, Module=CairoMakie) - -smalljets, smallind = anti_kt_algo_alt(smalldata, R=1); img = jetsplot(smalldata, smallind, Module=CairoMakie) - -smalljets, smallind = kt_algo(smalldata, R=1); img = jetsplot(smalldata, smallind, Module=CairoMakie) diff --git a/src/JetReconstruction.jl b/src/JetReconstruction.jl index c916bb7..e3e424d 100644 --- a/src/JetReconstruction.jl +++ b/src/JetReconstruction.jl @@ -25,16 +25,14 @@ energy(p::LorentzVectorCyl) = LorentzVectorHEP.energy(p) # Philipp's pseudojet include("Pseudojet.jl") export PseudoJet -## As this is an internal EDM class, we don't export anything # Simple HepMC3 reader include("HepMC3.jl") -export HepMC3 ## N2Plain algorithm # Algorithmic part for simple sequential implementation -include("Algo.jl") -export sequential_jet_reconstruct, kt_algo, anti_kt_algo, anti_kt_algo_alt, cambridge_aachen_algo +include("PlainAlgo.jl") +export plain_jet_reconstruct ## Tiled algorithms # Common pieces @@ -42,7 +40,7 @@ include("TiledAlgoUtils.jl") # Algorithmic part, tiled reconstruction strategy with linked list jet objects include("TiledAlgoLL.jl") -export tiled_jet_reconstruct_ll +export tiled_jet_reconstruct # jet serialisation (saving to file) include("Serialize.jl") diff --git a/src/Algo.jl b/src/PlainAlgo.jl similarity index 71% rename from src/Algo.jl rename to src/PlainAlgo.jl index d189188..33744b4 100644 --- a/src/Algo.jl +++ b/src/PlainAlgo.jl @@ -78,28 +78,34 @@ end """ This is the N2Plain jet reconstruction algorithm interface, called with an arbitrary array -of objects, which supports the methods pt2(), phi(), rapidity() for each element. +of particles, which supports suitable 4-vector methods, viz. + - pt2(), phi(), rapidity(), px(), py(), pz(), energy() +for each element. + +Examples of suitable types are JetReconstruction.PseudoJet and LorentzVectorHEP. +N.B. these methods need to exist in the namespace of this package, i.e. JetReconstruction.pt2(::T), +which is already done for the two types above. """ -function sequential_jet_reconstruct(objects::AbstractArray{T}; p = -1, R = 1.0, recombine = +, ptmin = 0.0) where T +function plain_jet_reconstruct(particles::Vector{T}; p = -1, R = 1.0, recombine = +, ptmin = 0.0) where T # Integer p if possible p = (round(p) == p) ? Int(p) : p # We make sure these arrays are type stable - have seen issues where, depending on the values # returned by the methods, they can become unstable and performance degrades - kt2_array::Vector{Float64} = pt2.(objects) .^ p - phi_array::Vector{Float64} = phi.(objects) - rapidity_array::Vector{Float64} = rapidity.(objects) + kt2_array::Vector{Float64} = pt2.(particles) .^ p + phi_array::Vector{Float64} = phi.(particles) + rapidity_array::Vector{Float64} = rapidity.(particles) - objects_array = copy(objects) + objects_array = copy(particles) # Now call the actual reconstruction method, tuned for our internal EDM - sequential_jet_reconstruct(objects_array=objects_array, kt2_array=kt2_array, phi_array=phi_array, + plain_jet_reconstruct(objects_array=objects_array, kt2_array=kt2_array, phi_array=phi_array, rapidity_array=rapidity_array, p=p, R=R, recombine=recombine, ptmin=ptmin) end -function sequential_jet_reconstruct(;objects_array::AbstractArray{J}, kt2_array::Vector{F}, +function plain_jet_reconstruct(;objects_array::Vector{J}, kt2_array::Vector{F}, phi_array::Vector{F}, rapidity_array::Vector{F}, p = -1, R = 1.0, recombine = +, ptmin = 0.0) where {J, F<:AbstractFloat} # Bounds N::Int = length(objects_array) @@ -202,42 +208,3 @@ function sequential_jet_reconstruct(;objects_array::AbstractArray{J}, kt2_array: jets, sequences end - -""" -`anti_kt_algo(objects; R=1, recombine=(x, y)->(x + y)) -> Vector, Vector{Vector{Int}}` - -Runs the anti-kt jet reconstruction algorithm. `objects` can be any collection of *unique* elements. - -Returns: - `jets` - a vector of jets. Each jet is of the same type as elements in `objects`. - `sequences` - a vector of vectors of indices in `objects`. For all `i`, `sequences[i]` gives a sequence of indices of objects that have been combined into the i-th jet (`jets[i]`). -""" -function anti_kt_algo(objects; R = 1.0, recombine = +) - sequential_jet_reconstruct(objects, p = -1, R = R, recombine = recombine) -end - -""" -`kt_algo(objects; R=1, recombine=(x, y)->(x + y)) -> Vector, Vector{Vector{Int}}` - -Runs the kt jet reconstruction algorithm. `objects` can be any collection of *unique* elements. - -Returns: - `jets` - a vector of jets. Each jet is of the same type as elements in `objects`. - `sequences` - a vector of vectors of indices in `objects`. For all `i`, `sequences[i]` gives a sequence of indices of objects that have been combined into the i-th jet (`jets[i]`). -""" -function kt_algo(objects; R = 1.0, recombine = +) - sequential_jet_reconstruct(objects, p = 1, R = R, recombine = recombine) -end - -""" -`cambridge_aachen_algo(objects; R=1, recombine=(x, y)->(x + y)) -> Vector, Vector{Vector{Int}}` - -Runs the Cambridge/Aachen jet reconstruction algorithm. `objects` can be any collection of *unique* elements. - -Returns: - `jets` - a vector of jets. Each jet is of the same type as elements in `objects`. - `sequences` - a vector of vectors of indices in `objects`. For all `i`, `sequences[i]` gives a sequence of indices of objects that have been combined into the i-th jet (`jets[i]`). -""" -function cambridge_aachen_algo(objects; R = 1.0, recombine = +) - sequential_jet_reconstruct(objects, p = 0, R = R, recombine = recombine) -end diff --git a/src/TiledAlgoLL.jl b/src/TiledAlgoLL.jl index 2dcf379..29871a6 100644 --- a/src/TiledAlgoLL.jl +++ b/src/TiledAlgoLL.jl @@ -306,7 +306,7 @@ be used. If a non-standard recombination is used, it must be defined for JetReconstruction.PseudoJet, as this struct is used internally. """ -function tiled_jet_reconstruct_ll(particles::Vector{T}; p = -1, R = 1.0, recombine = +, ptmin = 0.0) where {T} +function tiled_jet_reconstruct(particles::Vector{T}; p = -1, R = 1.0, recombine = +, ptmin = 0.0) where {T} # Here we need to populate the vector of PseudoJets that are the internal # EDM for the main algorithm, then we call the reconstruction pseudojets = Vector{PseudoJet}(undef, length(particles)) @@ -314,13 +314,13 @@ function tiled_jet_reconstruct_ll(particles::Vector{T}; p = -1, R = 1.0, recombi pseudojets[i] = PseudoJet(px(particle), py(particle), pz(particle), energy(particle)) end - tiled_jet_reconstruct_ll(pseudojets, p = p, R = R, recombine = recombine, ptmin = ptmin) + tiled_jet_reconstruct(pseudojets, p = p, R = R, recombine = recombine, ptmin = ptmin) end """ Main jet reconstruction algorithm, using PseudoJet objects """ -function tiled_jet_reconstruct_ll(particles::Vector{PseudoJet}; p = -1, R = 1.0, recombine = +, ptmin = 0.0) +function tiled_jet_reconstruct(particles::Vector{PseudoJet}; p = -1, R = 1.0, recombine = +, ptmin = 0.0) # Bounds N::Int = length(particles) # @debug "Initial particles: $(N)" diff --git a/test/runtests.jl b/test/runtests.jl index 5bb3846..8ea0bea 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -41,10 +41,10 @@ function do_jet_test(strategy::JetRecoStrategy, fastjet_jets; # Strategy if (strategy == N2Plain) - jet_reconstruction = sequential_jet_reconstruct + jet_reconstruction = plain_jet_reconstruct strategy_name = "N2Plain" elseif (strategy == N2Tiled) - jet_reconstruction = tiled_jet_reconstruct_ll + jet_reconstruction = tiled_jet_reconstruct strategy_name = "N2Tiled" else throw(ErrorException("Strategy not yet implemented")) From e910711fe2e93382827f3b4ecef67a51b65abce5 Mon Sep 17 00:00:00 2001 From: Graeme A Stewart Date: Wed, 25 Oct 2023 15:21:18 +0200 Subject: [PATCH 21/21] Make mass2 alias const Co-authored-by: Jerry Ling --- src/Pseudojet.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Pseudojet.jl b/src/Pseudojet.jl index bd88e7a..d8a712e 100644 --- a/src/Pseudojet.jl +++ b/src/Pseudojet.jl @@ -148,7 +148,7 @@ px(p::PseudoJet) = p.px py(p::PseudoJet) = p.py pz(p::PseudoJet) = p.pz mass(p::PseudoJet) = m(p) -mass2 = m2 +const mass2 = m2 energy(p::PseudoJet) = p.E import Base.+;