diff --git a/docs/Project.toml b/docs/Project.toml index 4269e90..d5414ba 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,6 +1,8 @@ [deps] +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" JetReconstruction = "44e8cb2c-dfab-4825-9c70-d4808a591196" +Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" [compat] Documenter = "1.4" diff --git a/docs/make.jl b/docs/make.jl index 95af683..06030b9 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,7 +1,19 @@ using Documenter +using CairoMakie using JetReconstruction -makedocs(sitename = "JetReconstruction.jl") +push!(LOAD_PATH, "../ext/") + +include(joinpath(@__DIR__, "..", "ext", "JetVisualisation.jl")) + +makedocs(sitename = "JetReconstruction.jl", + pages = [ + "Home" => "index.md", + "Examples" => "examples.md", + "Reference" => Any["Public API" => "lib/public.md", + "Internal API" => "lib/internal.md", + "Visualisation API" => "lib/visualisation.md"] + ]) deploydocs(repo = "github.com/JuliaHEP/JetReconstruction.jl.git", push_preview = true) diff --git a/docs/src/assets/logo.png b/docs/src/assets/logo.png new file mode 100644 index 0000000..53e408c Binary files /dev/null and b/docs/src/assets/logo.png differ diff --git a/docs/src/assets/logo.svg b/docs/src/assets/logo.svg new file mode 100644 index 0000000..1418c40 --- /dev/null +++ b/docs/src/assets/logo.svg @@ -0,0 +1,180 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/src/examples.md b/docs/src/examples.md index b417f0d..87e23a5 100644 --- a/docs/src/examples.md +++ b/docs/src/examples.md @@ -44,3 +44,8 @@ Similar to `visualise-jets.jl` this notebook will produce a visualisation of jet reconstruction in the browser. This is a 3D plot where all the initial energy deposits are visualised, with colours that indicate in which final cluster the deposit ended up in. + +## `animate-reconstruction.jl` + +Performs jet reconstruction and then produces and animation of the process, +showing how the jets merge from their different constituents. diff --git a/docs/src/index.md b/docs/src/index.md index 17e4076..8035e1a 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -16,17 +16,16 @@ inclusive ``k_\text{T}``. ## Reconstruction Interface -The main interface for reconstruction is: +The main interface for reconstruction is [`jet_reconstruct`](@ref), called as, e.g., -```@docs -jet_reconstruct(particles; p = -1, R = 1.0, recombine = +, strategy = RecoStrategy.Best) +```julia +jet_reconstruct(particles; p = -1, R = 1.0) ``` -The object returned is a `ClusterSequence`, which internally tracks all merge steps. +Where `particles` is a collection of 4-vector objects to reconstruct. -```@docs -ClusterSequence -``` +The object returned is a [`ClusterSequence`](@ref), which internally tracks all +merge steps. ## Strategy @@ -38,9 +37,12 @@ Three strategies are available for the different algorithms: | `RecoStrategy.N2Plain` | Global matching of particles at each interation (works well for low $N$) | `plain_jet_reconstruct` | | `RecoStrategy.N2Tiled` | Use tiles of radius $R$ to limit search space (works well for higher $N$) | `tiled_jet_reconstruct` | -Generally one can use the `jet_reconstruct` interface, shown above, as the *Best* strategy safely as the overhead is extremely low. That interface supports a `strategy` option to switch to a different option. +Generally one can use the `jet_reconstruct` interface, shown above, as the +*Best* strategy safely as the overhead is extremely low. That interface supports +a `strategy` option to switch to a different option. -Another option, if one wishes to use a specific strategy, is to call that strategy's interface directly, e.g., +Another option, if one wishes to use a specific strategy, is to call that +strategy's interface directly, e.g., ```julia # For N2Plain strategy called directly @@ -54,19 +56,8 @@ Note that there is no `strategy` option in these interfaces. To obtain final jets both inclusive (``p_T`` cut) and exclusive (``n_{jets}`` or ``d_{ij}`` cut) selections are supported: -```@docs -inclusive_jets(clusterseq::ClusterSequence, ptmin = 0.0) -``` - -```@docs -exclusive_jets(clusterseq::ClusterSequence; dcut = nothing, njets = nothing) -``` - -The number of exclusive jets passing a particular `dcut` can be obtained: - -```@docs -n_exclusive_jets(clusterseq::ClusterSequence; dcut::AbstractFloat) -``` +- [inclusive_jets(clusterseq::ClusterSequence, ptmin = 0.0)](@ref) +- [exclusive_jets(clusterseq::ClusterSequence; dcut = nothing, njets = nothing)](@ref) ### Sorting @@ -79,7 +70,7 @@ sorted_jets = sort!(inclusive_jets(cs::ClusterSequence; ptmin=5.0), by=JetReconstruction.energy, rev=true) ``` -## Plotting +## Plotting and Animation ![illustration](assets/jetvis.png) @@ -91,6 +82,10 @@ directory. The plotting code is a package extension and will load if the one of the `Makie` modules is loaded in the environment. +The [`animatereco`](@ref) function will animate the reconstruction sequence, given a +`ClusterSequence` object. See the function documentation for the many options +that can be customised. + ## Serialisation The package also provides methods such as `loadjets`, `loadjets!`, and diff --git a/docs/src/lib/internal.md b/docs/src/lib/internal.md new file mode 100644 index 0000000..9d65716 --- /dev/null +++ b/docs/src/lib/internal.md @@ -0,0 +1,23 @@ +# Jet Reconstruction Internal Documentation + +Documentation for `JetReconstruction.jl`'s internal methods and types. + +N.B. no guarantee is made of stability of these interfaces or types. + +```@meta +CollapsedDocStrings = true +``` + +## Index + +```@index +Pages = ["internal.md"] +``` + +## Public Interface + +```@autodocs +Modules = [JetReconstruction] +Public = false +Order = [:function, :type] +``` diff --git a/docs/src/lib/public.md b/docs/src/lib/public.md new file mode 100644 index 0000000..2ab0411 --- /dev/null +++ b/docs/src/lib/public.md @@ -0,0 +1,17 @@ +# Jet Reconstruction Public Documentation + +Documentation for `JetReconstruction.jl`'s public interfaces. + +## Index + +```@index +Pages = ["public.md"] +``` + +## Public Methods and Types + +```@autodocs +Modules = [JetReconstruction] +Private = false +Order = [:function, :type] +``` diff --git a/docs/src/lib/visualisation.md b/docs/src/lib/visualisation.md new file mode 100644 index 0000000..16d1efd --- /dev/null +++ b/docs/src/lib/visualisation.md @@ -0,0 +1,16 @@ +# Jet Visualisation Documentation + +Documentation for visualisation interfaces extension module. + +## Index + +```@index +Pages = ["visualisation.md"] +``` + +## Jet Visualisation Public Interfaces + +```@autodocs +Modules = [JetVisualisation] +Order = [:function] +``` diff --git a/examples/.gitignore b/examples/.gitignore index 2a92ec2..3b724ea 100644 --- a/examples/.gitignore +++ b/examples/.gitignore @@ -1,3 +1,4 @@ # Ignore visualisation outputs *.pdf *.png +*.mp4 diff --git a/examples/README.md b/examples/README.md index 3f4f3e8..8dfc0dc 100644 --- a/examples/README.md +++ b/examples/README.md @@ -29,3 +29,8 @@ colours that indicate in which final cluster the deposit ended up in. This notebook will produce a visualisation of jet reconstruction in the browser. This is a 3D plot where all the initial energy deposits are visualised, with colours that indicate in which final cluster the deposit ended up in. + +## `animate-reconstruction.jl` + +Performs jet reconstruction and then produces and animation of the process, +showing how the jets merge from their different constituents. diff --git a/examples/animate-reconstruction.jl b/examples/animate-reconstruction.jl new file mode 100644 index 0000000..2642ae0 --- /dev/null +++ b/examples/animate-reconstruction.jl @@ -0,0 +1,64 @@ +#! /usr/bin/env julia +"""Use the visualisation tools to produce an animation of the jet reconstruction""" + +using ArgParse +using Logging +using CairoMakie + +using JetReconstruction + +# Parsing for algorithm and strategy enums +include(joinpath(@__DIR__, "parse-options.jl")) + +function main() + s = ArgParseSettings(autofix_names = true) + @add_arg_table s begin + "--event", "-e" + help = "Event in file to visualise" + arg_type = Int + default = 1 + + "--distance", "-R" + help = "Distance parameter for jet merging" + arg_type = Float64 + default = 0.4 + + "--algorithm", "-A" + help = """Algorithm to use for jet reconstruction: $(join(JetReconstruction.AllJetRecoAlgorithms, ", "))""" + arg_type = JetAlgorithm.Algorithm + default = JetAlgorithm.AntiKt + + "--strategy", "-S" + help = """Strategy for the algorithm, valid values: $(join(JetReconstruction.AllJetRecoStrategies, ", "))""" + arg_type = RecoStrategy.Strategy + default = RecoStrategy.Best + + "file" + help = "HepMC3 event file in HepMC3 to read" + required = true + + "output" + help = "File for output animation" + default = "jetreco.mp4" + end + args = parse_args(ARGS, s; as_symbols = true) + + logger = ConsoleLogger(stdout, Logging.Info) + global_logger(logger) + + events::Vector{Vector{PseudoJet}} = read_final_state_particles(args[:file], + maxevents = args[:event], + skipevents = args[:event]) + + power = JetReconstruction.algorithm2power[args[:algorithm]] + cs = jet_reconstruct(events[1], R = args[:distance], p = power, + strategy = args[:strategy]) + + animatereco(cs, args[:output]; azimuth = (1.8, 3.0), elevation = 0.5, + perspective = 0.5, framerate = 20, ancestors = true, + barsize_phi = 0.1, barsize_y = 0.3) + + @info "Saved jet reconstruction animation to $(args[:output])" +end + +main() diff --git a/examples/instrumented-jetreco.jl b/examples/instrumented-jetreco.jl index 0049e4f..d2ba765 100644 --- a/examples/instrumented-jetreco.jl +++ b/examples/instrumented-jetreco.jl @@ -25,7 +25,7 @@ include(joinpath(@__DIR__, "parse-options.jl")) function profile_code(profile, jet_reconstruction, events, niters; R = 0.4, p = -1, strategy = RecoStrategy.N2Tiled) Profile.init(n = 5 * 10^6, delay = 0.00001) - profile_events(events) = begin + function profile_events(events) for evt in events jet_reconstruction(evt, R = R, p = p, strategy = strategy) end @@ -76,7 +76,7 @@ function jet_process(events::Vector{Vector{PseudoJet}}; gcoff::Bool = false, profile = nothing, alloc::Bool = false, - dump::Union{String, Nothing} = nothing,) + dump::Union{String, Nothing} = nothing) @info "Will process $(size(events)[1]) events" # Map algorithm to power @@ -130,7 +130,8 @@ function jet_process(events::Vector{Vector{PseudoJet}}; dcut = dcut) else finaljets = inclusive_jets(jet_reconstruct(event, R = distance, p = power, - strategy = strategy), ptmin) + strategy = strategy), + ptmin) end # Only print the jet content once if irun == 1 diff --git a/examples/jetreco.jl b/examples/jetreco.jl index 95bc2bf..d036963 100644 --- a/examples/jetreco.jl +++ b/examples/jetreco.jl @@ -6,8 +6,6 @@ and performs standard jet reconstruction on the final state particles. using ArgParse using Profile -using Colors -using StatProfilerHTML using Logging using JSON @@ -32,7 +30,7 @@ function jet_process(events::Vector{Vector{PseudoJet}}; dcut = nothing, njets = nothing, strategy::RecoStrategy.Strategy, - dump::Union{String, Nothing} = nothing,) + dump::Union{String, Nothing} = nothing) @info "Will process $(size(events)[1]) events" # If we are dumping the results, setup the JSON structure @@ -56,10 +54,12 @@ function jet_process(events::Vector{Vector{PseudoJet}}; for (ievt, event) in enumerate(events) if !isnothing(njets) finaljets = exclusive_jets(jet_reconstruct(event, R = distance, p = power, - strategy = strategy), njets = njets) + strategy = strategy), + njets = njets) elseif !isnothing(dcut) finaljets = exclusive_jets(jet_reconstruct(event, R = distance, p = power, - strategy = strategy), dcut = dcut) + strategy = strategy), + dcut = dcut) else finaljets = inclusive_jets(jet_reconstruct(event, R = distance, p = power, strategy = strategy), ptmin) diff --git a/ext/JetVisualisation.jl b/ext/JetVisualisation.jl index 97225bf..c0c8726 100644 --- a/ext/JetVisualisation.jl +++ b/ext/JetVisualisation.jl @@ -5,17 +5,6 @@ module JetVisualisation using JetReconstruction using Makie -function get_all_ancestors(idx, cs::ClusterSequence) - if cs.history[idx].parent1 == JetReconstruction.NonexistentParent - return [cs.history[idx].jetp_index] - else - branch1 = get_all_ancestors(cs.history[idx].parent1, cs) - cs.history[idx].parent2 == JetReconstruction.BeamJet && return branch1 - branch2 = get_all_ancestors(cs.history[idx].parent2, cs) - return [branch1; branch2] - end -end - """ jetsplot(objects, cs::ClusterSequence; barsize_phi=0.1, barsize_eta=0.1, colormap=:glasbey_hv_n256, Module=Main) @@ -46,11 +35,11 @@ jetsplot(my_objects, cs, Module=Main) #default """ function JetReconstruction.jetsplot(objects, cs::ClusterSequence; barsize_phi = 0.1, barsize_eta = 0.1, colormap = :glasbey_hv_n256, - Module = CairoMakie) + Module = Makie) idx_arrays = Vector{Int}[] for elt in cs.history elt.parent2 == JetReconstruction.BeamJet || continue - push!(idx_arrays, get_all_ancestors(elt.parent1, cs)) + push!(idx_arrays, JetReconstruction.get_all_ancestors(elt.parent1, cs)) end jetsplot(objects, idx_arrays; barsize_phi, barsize_eta, colormap, Module) @@ -108,7 +97,208 @@ function JetReconstruction.jetsplot(objects, idx_arrays; barsize_phi = 0.1, xlabel = "ϕ", ylabel = "η", zlabel = "kt", limits = (nothing, nothing, nothing, nothing, 0, findmax(pts)[1] + 10)), - shading = NoShading,) + shading = NoShading) +end + +function JetReconstruction.jetsplot(cs::ClusterSequence, + reco_state::Dict{Int, + JetReconstruction.JetWithAncestors}; + barsize_phi = 0.1, + barsize_y = 0.1, colormap = :glasbey_category10_n256, + Module = Makie) + # Setup the marker as a square object + jet_plot_marker = Rect3f(Vec3f(0), Vec3f(1)) + + # Get the jet variables to plot + phis = [JetReconstruction.phi(x.self) for x in values(reco_state)] + ys = [JetReconstruction.rapidity(x.self) for x in values(reco_state)] + pts = [JetReconstruction.pt(x.self) for x in values(reco_state)] + + # The core points to plot are on the pt=0 axis, with the marker size + # scaled up to the p_T of the jet + jet_plot_points = Point3f.(phis .- (barsize_phi / 2), ys .- (barsize_y / 2), 0pts) + jet_plot_marker_size = Vec3f.(barsize_phi, barsize_y, pts) + + # Colours are defined from the rank of the ancestor jets + jet_plot_colours = [x.jet_rank for x in values(reco_state)] + + # Limits for rapidity and p_T are set to the maximum values in the data + # (For ϕ the limits are (0, 2π)) + min_rap = max_rap = max_pt = 0.0 + for jet in cs.jets + min_rap = min(min_rap, JetReconstruction.rapidity(jet)) + max_rap = max(max_rap, JetReconstruction.rapidity(jet)) + max_pt = max(max_pt, JetReconstruction.pt(jet)) + end + + fig, ax, plt_obj = Module.meshscatter(jet_plot_points; + markersize = jet_plot_marker_size, + marker = jet_plot_marker, + colormap = colormap, + color = jet_plot_colours, + colorrange = (1, 256), + figure = (size = (700, 600),), + axis = (type = Axis3, perspectiveness = 0.5, + azimuth = 2.7, + elevation = 0.5, + xlabel = L"\phi", ylabel = L"y", + zlabel = L"p_T", + limits = (0, 2π, min_rap - 0.5, + max_rap + 0.5, 0, max_pt + 10)), + shading = NoShading) + fig, ax, plt_obj +end + +""" + animatereco(cs::ClusterSequence, filename; + barsize_phi = 0.1, + barsize_y = 0.1, + colormap = :glasbey_category10_n256, + perspective = 0.5, + azimuth = 2.7, + elevation = 0.5, + framerate = 5, + ancestors = false, + Module = Makie) + +Animate the jet reconstruction process and save it as a video file. + +# Arguments +- `cs::ClusterSequence`: The cluster sequence object containing the jets. +- `filename`: The name of the output video file. + +# Optional Arguments +- `barsize_phi=0.1`: The size of the bars in the phi direction. +- `barsize_y=0.1`: The size of the bars in the y direction. +- `colormap=:glasbey_category10_n256`: The colormap to use for coloring the + jets. +- `perspective=0.5`: The perspective of the plot. +- `azimuth=2.7`: The azimuth angle of the plot. +- `elevation=0.5`: The elevation angle of the plot. +- `framerate=5`: The framerate of the output video. +- `end_frames=0`: The number of static frames to show at the end of the + animation. This can be useful to show the final state of the jets for a longer + time. +- `title=nothing`: The title to add to the plot. +- `ancestors=false`: Whether to include ancestors of the jets in the animation. + When `true` the ancestors of the jets will be plotted as well, as height zero + bars, with the same colour as the jet they are ancestors of. +- `Module`: The plotting module to use. Default is `Makie`. + +For `perspective`, `azimuth`, and `elevation`, a single value can be passed for +a fixed viewpoint, or a tuple of two values for a changing viewpoint. The +viewpoint will then change linearly between the two values over the course of +the animation. + +# Returns +- `fig`: The figure object representing the final frame. + +""" +function JetReconstruction.animatereco(cs::ClusterSequence, filename; + barsize_phi = 0.1, + barsize_y = 0.1, + colormap = :glasbey_category10_n256, + perspective::Union{Real, Tuple{Real, Real}} = 0.5, + azimuth::Union{Real, Tuple{Real, Real}} = 2.7, + elevation::Union{Real, Tuple{Real, Real}} = 0.5, + framerate = 5, end_frames = 0, ancestors = false, + title = nothing, + Module = Makie) + # Setup the marker as a square object + jet_plot_marker = Rect3f(Vec3f(0), Vec3f(1)) + + # Get the number of meaningful reconstruction steps and rank the initial + # particles by p_T + merge_steps = JetReconstruction.merge_steps(cs) + jet_ranks = JetReconstruction.jet_ranks(cs) + + # End point of the catagorical color map + # (How to get this programmatically from the CM Symbol?) + colormap_end = 256 + + # Keep plot limits constant + min_rap = max_rap = max_pt = 0.0 + for jet in cs.jets + min_rap = min(min_rap, JetReconstruction.rapidity(jet)) + max_rap = max(max_rap, JetReconstruction.rapidity(jet)) + max_pt = max(max_pt, JetReconstruction.pt(jet)) + end + + # Get the reconstruction state at each meaningful iteration + # And calculate the plot parameters + all_jet_plot_points = Vector{Vector{Point3f}}() + all_jet_plot_marker_size = Vector{Vector{Vec3f}}() + all_jet_plot_colours = Vector{Vector{Int}}() + for step in 0:merge_steps + reco_state = JetReconstruction.reco_state(cs, jet_ranks; iteration = step) + phis = [JetReconstruction.phi(x.self) for x in values(reco_state)] + ys = [JetReconstruction.rapidity(x.self) for x in values(reco_state)] + pts = [JetReconstruction.pt(x.self) for x in values(reco_state)] + push!(all_jet_plot_points, + Point3f.(phis .- (barsize_phi / 2), ys .- (barsize_y / 2), 0pts)) + push!(all_jet_plot_marker_size, Vec3f.(barsize_phi, barsize_y, pts)) + push!(all_jet_plot_colours, + [mod1(x.jet_rank, colormap_end) for x in values(reco_state)]) + if ancestors + for jet_entry in values(reco_state) + for ancestor in jet_entry.ancestors + ancestor_jet = cs.jets[ancestor] + push!(all_jet_plot_points[end], + Point3f(JetReconstruction.phi(ancestor_jet) - (barsize_phi / 2), + JetReconstruction.rapidity(ancestor_jet) - + (barsize_y / 2), + 0.0)) + push!(all_jet_plot_marker_size[end], + Vec3f(barsize_phi, barsize_y, 0.001max_pt)) + push!(all_jet_plot_colours[end], mod1(jet_entry.jet_rank, colormap_end)) + end + end + end + end + + # Keep plot limits constant + min_rap = max_rap = max_pt = 0.0 + for jet in cs.jets + min_rap = min(min_rap, JetReconstruction.rapidity(jet)) + max_rap = max(max_rap, JetReconstruction.rapidity(jet)) + max_pt = max(max_pt, JetReconstruction.pt(jet)) + end + + # Setup an observable for the iteration number + it_obs = Observable(0) + jet_plot_points_obs = @lift all_jet_plot_points[$it_obs + 1] + jet_plot_marker_size_obs = @lift all_jet_plot_marker_size[$it_obs + 1] + jet_plot_colours_obs = @lift all_jet_plot_colours[$it_obs + 1] + + # We may want to have a shifting viewpoint + azimuth_axis = typeof(azimuth) <: Tuple ? + @lift(azimuth[1]+$it_obs / merge_steps * (azimuth[2] - azimuth[1])) : + azimuth + elevation_axis = typeof(elevation) <: Tuple ? + @lift(elevation[1]+$it_obs / merge_steps * + (elevation[2] - elevation[1])) : elevation + perspective_axis = typeof(perspective) <: Tuple ? + @lift(perspective[1]+$it_obs / merge_steps * + (perspective[2] - perspective[1])) : perspective + + ax = (type = Axis3, title = isnothing(title) ? "" : title, + xlabel = L"\phi", ylabel = L"y", zlabel = L"p_T", + limits = (0, 2π, min_rap - 0.5, max_rap + 0.5, 0, max_pt + 10), + perspectiveness = perspective_axis, azimuth = azimuth_axis, + elevation = elevation_axis) + fig = Module.meshscatter(jet_plot_points_obs; + markersize = jet_plot_marker_size_obs, + marker = jet_plot_marker, + colormap = colormap, + color = jet_plot_colours_obs, + colorrange = (1, colormap_end), + figure = (size = (800, 600),), + axis = ax, + shading = NoShading) + record(fig, filename, 0:(merge_steps + end_frames); framerate = framerate) do iteration + it_obs[] = min(iteration, merge_steps) + end + fig end end diff --git a/src/ClusterSequence.jl b/src/ClusterSequence.jl index ff2277d..5257cd1 100644 --- a/src/ClusterSequence.jl +++ b/src/ClusterSequence.jl @@ -25,12 +25,21 @@ const BeamJet = -1 A struct holding a record of jet mergers and finalisations Fields: -- parent1: Index in history where first parent of this jet was created (NonexistentParent if this jet is an original particle) -- parent2: Index in history where second parent of this jet was created (NonexistentParent if this jet is an original particle); BeamJet if this history entry just labels the fact that the jet has recombined with the beam) -- child: Index in history where the current jet is recombined with another jet to form its child. It is Invalid if this jet does not further recombine. -- jetp_index: Index in the jets vector where we will find the PseudoJet object corresponding to this jet (i.e. the jet created at this entry of the history). NB: if this element of the history corresponds to a beam recombination, then jetp_index=Invalid. -- dij: The distance corresponding to the recombination at this stage of the clustering. -- max_dij_so_far: The largest recombination distance seen so far in the clustering history. +- `parent1`: Index in history where first parent of this jet was created + (NonexistentParent if this jet is an original particle) +- `parent2`: Index in history where second parent of this jet was created + (NonexistentParent if this jet is an original particle); BeamJet if this + history entry just labels the fact that the jet has recombined with the beam) +- `child`: Index in history where the current jet is recombined with another jet + to form its child. It is Invalid if this jet does not further recombine. +- `jetp_index`: Index in the jets vector where we will find the PseudoJet object + corresponding to this jet (i.e. the jet created at this entry of the history). + NB: if this element of the history corresponds to a beam recombination, then + `jetp_index=Invalid`. +- `dij`: The distance corresponding to the recombination at this stage of the + clustering. +- `max_dij_so_far`: The largest recombination distance seen so far in the + clustering history. """ struct HistoryElement parent1::Int @@ -331,3 +340,174 @@ function n_exclusive_jets(clusterseq::ClusterSequence; dcut::AbstractFloat) # The number of jets is then given by this formula length(clusterseq.history) - i_dcut end + +""" + get_all_ancestors(idx, cs::ClusterSequence) + +Recursively finds all ancestors of a given index in a `ClusterSequence` object. + +# Arguments +- `idx`: The index of the jet for which to find ancestors. +- `cs`: The `ClusterSequence` object containing the jet history. + +# Returns +An array of indices representing the ancestors of the given jet. +""" +function get_all_ancestors(idx, cs::ClusterSequence) + if cs.history[idx].parent1 == JetReconstruction.NonexistentParent + return [cs.history[idx].jetp_index] + else + branch1 = get_all_ancestors(cs.history[idx].parent1, cs) + cs.history[idx].parent2 == JetReconstruction.BeamJet && return branch1 + branch2 = get_all_ancestors(cs.history[idx].parent2, cs) + return [branch1; branch2] + end +end + +""" + merge_steps(clusterseq::ClusterSequence) + +Compute the number of jet-jet merge steps in a cluster sequence. This is useful +to give the number of meaningful recombination steps in a jet reconstruction +sequence (beam merge steps are not counted). + +# Arguments +- `clusterseq::ClusterSequence`: The cluster sequence object. + +# Returns +- `merge_steps::Int`: The number of merge steps. +""" +function merge_steps(clusterseq::ClusterSequence) + merge_steps = 0 + for step in clusterseq.history[(clusterseq.n_initial_jets + 1):end] + step.parent2 == BeamJet && continue + merge_steps += 1 + end + merge_steps +end + +""" + jet_ranks(clusterseq::ClusterSequence; compare_fn = JetReconstruction.pt) + +Compute the ranks of jets in a given `ClusterSequence` object based on a +specified comparison function. + +## Arguments +- `clusterseq::ClusterSequence`: The `ClusterSequence` object containing the + jets to rank. +- `compare_fn = JetReconstruction.pt`: The comparison function used to determine + the order of the jets. Defaults to `JetReconstruction.pt`, which compares jets + based on their transverse momentum. + +## Returns +A dictionary mapping each jet index to its rank. + +## Note +This is a utility function that can be used to rank initial clusters based on a specified +jet property. It can be used to assign a consistent "rank" to each reconstructed jet in +the cluster sequence, which is useful for stable plotting of jet outputs. +""" +function jet_ranks(clusterseq::ClusterSequence; compare_fn = JetReconstruction.pt) + initial_jet_list = collect(1:(clusterseq.n_initial_jets)) + sort!(initial_jet_list, by = i -> compare_fn(clusterseq.jets[i]), rev = true) + jet_ranks = Dict{Int, Int}() + for (rank, jetp_index) in enumerate(initial_jet_list) + jet_ranks[jetp_index] = rank + end + jet_ranks +end + +""" + struct JetWithAncestors + +A struct representing a jet with its origin ancestors. + +# Fields +- `self::PseudoJet`: The PseudoJet object for this jet. +- `jetp_index::Int`: The index of the jet in the corresponding cluster sequence. +- `ancestors::Set{Int}`: A set of indices representing the jetp_indexes of + ancestors of the jet (in the cluster sequence). +- `jet_rank::Int`: The rank of the jet based on a comparison of all of the jet's + ancestors + +# Note +This structure needs its associated cluster sequence origin to be useful. +""" +struct JetWithAncestors + self::PseudoJet + jetp_index::Int + ancestors::Set{Int} + jet_rank::Int +end + +""" + reco_state(cs::ClusterSequence, pt_ranks; iteration=0) + +This function returns the reconstruction state of a `ClusterSequence` object +based on a given iteration number in the reconstruction. + +# Arguments +- `cs::ClusterSequence`: The `ClusterSequence` object to update. +- `ranks`: The ranks of the original clusters, that are inherited by peudojets + during the reconstruction process. +- `iteration=0`: The iteration number to consider for updating the + reconstruction state (0 represents the initial state). +- `ignore_beam_merge=true`: Ignore beam merging steps in the reconstruction + (which produce no change in status). + +# Returns +A dictionary representing a snapshot of the reconstruction state. + +# Details +The function starts by initializing the reconstruction state with the initial +particles. Then, it walks over the iteration sequence and updates the +reconstruction state based on the history of recombination and finalization/beam +merger steps. +""" +function reco_state(cs::ClusterSequence, ranks; iteration = 0, ignore_beam_merge = true) + # Get the initial particles + reco_state = Dict{Int, JetWithAncestors}() + for jet_index in 1:(cs.n_initial_jets) + reco_state[jet_index] = JetWithAncestors(cs.jets[cs.history[jet_index].jetp_index], + jet_index, Set([]), ranks[jet_index]) + end + # Now update the reconstruction state by walking over the iteration sequence + iterations_done = 0 + for h_step in (cs.n_initial_jets + 1):length(cs.history) + h_entry = cs.history[h_step] + if h_entry.parent2 > 0 + # This is a recombination + iterations_done += 1 + # We store all of the original particle ancestors (but only the + # original ones, not intermediate merges) + my_ancestors = union(reco_state[cs.history[h_entry.parent1].jetp_index].ancestors, + reco_state[cs.history[h_entry.parent2].jetp_index].ancestors) + cs.history[h_entry.parent1].parent1 == JetReconstruction.NonexistentParent && + push!(my_ancestors, h_entry.parent1) + cs.history[h_entry.parent2].parent1 == JetReconstruction.NonexistentParent && + push!(my_ancestors, h_entry.parent2) + + # Now find the ancestor with the highest p_T value + pt_rank = cs.n_initial_jets + for ancestor in my_ancestors + (ranks[ancestor] < pt_rank) && (pt_rank = ranks[ancestor]) + end + + reco_state[h_entry.jetp_index] = JetWithAncestors(cs.jets[h_entry.jetp_index], + h_entry.jetp_index, + my_ancestors, pt_rank) + delete!(reco_state, cs.history[h_entry.parent1].jetp_index) + delete!(reco_state, cs.history[h_entry.parent2].jetp_index) + else + # This is a finalisation / beam merger, so here we do nothing + if ignore_beam_merge != true + iterations_done += 1 + end + end + # Abort when we have done the required number of iterations + if iterations_done == iteration + break + end + end + reco_state +end diff --git a/src/JetReconstruction.jl b/src/JetReconstruction.jl index 99164c3..e94ee70 100644 --- a/src/JetReconstruction.jl +++ b/src/JetReconstruction.jl @@ -75,7 +75,8 @@ export read_final_state_particles, read_final_state_particles_lv, final_jets # Jet visualisation is an extension, see ext/JetVisualisation.jl function jetsplot() end -export jetsplot +function animatereco() end +export jetsplot, animatereco # JSON results include("JSONresults.jl") diff --git a/src/TiledAlgoLL.jl b/src/TiledAlgoLL.jl index 1c719bb..e2190b3 100644 --- a/src/TiledAlgoLL.jl +++ b/src/TiledAlgoLL.jl @@ -101,7 +101,9 @@ tile_index(tiling_setup, eta::Float64, phi::Float64) = begin # - 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) + 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) return iphi * tiling_setup._n_tiles_eta + ieta