diff --git a/docs/src/examples.md b/docs/src/examples.md index 08e3688a..60f7dac8 100644 --- a/docs/src/examples.md +++ b/docs/src/examples.md @@ -26,6 +26,8 @@ julia --project jetreco.jl --algorithm=AntiKt ../test/data/events.pp13TeV.hepmc3 ... julia --project jetreco.jl --algorithm=Durham ../test/data/events.eeH.hepmc3.zst ... +julia --project jetreco.jl --algorithm=Valencia --p=1.2 --gamma=0.8 ../test/data/events.eeH.hepmc3.zst +... julia --project jetreco.jl --maxevents=10 --strategy=N2Plain --algorithm=Kt --exclusive-njets=3 ../test/data/events.pp13TeV.hepmc3.zst ... ``` diff --git a/docs/src/index.md b/docs/src/index.md index f8eb4854..ba4c4451 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -56,6 +56,7 @@ Each known algorithm is referenced using a `JetAlgorithm` scoped enum value. | generalised ``k_\text{T}`` | `JetAlgorithm.GenKt` | For $pp$, value of `p` must also be specified | | ``e^+e-`` ``k_\text{T}`` / Durham | `JetAlgorithm.Durham` | `R` value ignored and can be omitted | | generalised ``e^+e-`` ``k_\text{T}`` | `JetAlgorithm.EEKt` | For ``e^+e^-``, value of `p` must also be specified | +| Valencia | `JetAlgorithm.Valencia` | For ``e^+e^-``, values of `p` (β) and `γ` must be specified | ### Strategy diff --git a/src/AlgorithmStrategyEnums.jl b/src/AlgorithmStrategyEnums.jl index 64fe8fa7..757eea2e 100644 --- a/src/AlgorithmStrategyEnums.jl +++ b/src/AlgorithmStrategyEnums.jl @@ -26,8 +26,9 @@ Scoped enumeration (using EnumX) representing different jet algorithms used in t - `GenKt`: The Generalised Kt algorithm (with arbitrary power). - `EEKt`: The Generalised e+e- kt algorithm. - `Durham`: The e+e- kt algorithm, aka Durham. +- `Valencia`: The Valencia e+e- algorithm. """ -@enumx T=Algorithm JetAlgorithm AntiKt CA Kt GenKt EEKt Durham +@enumx T=Algorithm JetAlgorithm AntiKt CA Kt GenKt EEKt Durham Valencia const AllJetRecoAlgorithms = [String(Symbol(x)) for x in instances(JetAlgorithm.Algorithm)] """ @@ -36,7 +37,7 @@ const AllJetRecoAlgorithms = [String(Symbol(x)) for x in instances(JetAlgorithm. A constant array that contains the jet algorithms for which power is variable. """ -const varpower_algorithms = [JetAlgorithm.GenKt, JetAlgorithm.EEKt] +const varpower_algorithms = [JetAlgorithm.GenKt, JetAlgorithm.EEKt, JetAlgorithm.Valencia] """ algorithm2power @@ -109,7 +110,7 @@ Check if the algorithm is a e+e- reconstruction algorithm. `true` if the algorithm is a e+e- reconstruction algorithm, `false` otherwise. """ function is_ee(algorithm::JetAlgorithm.Algorithm) - return algorithm in [JetAlgorithm.EEKt, JetAlgorithm.Durham] + return algorithm in (JetAlgorithm.EEKt, JetAlgorithm.Durham, JetAlgorithm.Valencia) end """ diff --git a/src/ClusterSequence.jl b/src/ClusterSequence.jl index 47c4fd03..2476190e 100644 --- a/src/ClusterSequence.jl +++ b/src/ClusterSequence.jl @@ -332,7 +332,7 @@ function exclusive_jets(clusterseq::ClusterSequence{U}, throw(ArgumentError("Algorithm $(clusterseq.algorithm) requires power >= 0 for exclusive jets (power=$(clusterseq.power))")) elseif clusterseq.algorithm ∉ (JetAlgorithm.CA, JetAlgorithm.Kt, JetAlgorithm.Durham, JetAlgorithm.GenKt, - JetAlgorithm.EEKt) + JetAlgorithm.EEKt, JetAlgorithm.Valencia) throw(ArgumentError("Algorithm used is not suitable for exclusive jets ($(clusterseq.algorithm))")) end @@ -395,7 +395,7 @@ function n_exclusive_jets(clusterseq::ClusterSequence; dcut::AbstractFloat) throw(ArgumentError("Algorithm $(clusterseq.algorithm) requires power >= 0 for exclusive jets(power=$(clusterseq.power))")) elseif clusterseq.algorithm ∉ (JetAlgorithm.CA, JetAlgorithm.Kt, JetAlgorithm.Durham, JetAlgorithm.GenKt, - JetAlgorithm.EEKt) + JetAlgorithm.EEKt, JetAlgorithm.Valencia) throw(ArgumentError("Algorithm used is not suitable for exclusive jets ($(clusterseq.algorithm))")) end diff --git a/src/EEAlgorithm.jl b/src/EEAlgorithm.jl index 7c070479..abc714fd 100644 --- a/src/EEAlgorithm.jl +++ b/src/EEAlgorithm.jl @@ -1,5 +1,4 @@ -# Use these constants whenever we need to set a large value for the distance or -# metric distance +# Use these constants whenever we need to set a large value for the distance or metric distance const large_distance = 16.0 # = 4^2 const large_dij = 1.0e6 @@ -24,113 +23,220 @@ Calculate the angular distance between two jets `i` and `j` using the formula end """ - dij_dist(eereco, i, j, dij_factor) + dij_dist(eereco, i, j, dij_factor, algorithm::JetAlgorithm.Algorithm, invR2) -Calculate the dij distance between two ``e^+e^-``jets. +Compute dij distance for (Durham, EEKt, Valencia) using simple conditionals. +Beam index `j==0` returns a large sentinel. For Valencia we use the full +Valencia metric (independent of `dij_factor`). +""" +@inline function dij_dist(eereco, i, j, dij_factor, algorithm, invR2) + if !(algorithm isa JetAlgorithm.Algorithm) + throw(ArgumentError("algorithm must be a JetAlgorithm.Algorithm")) + end + j == 0 && return large_dij + @inbounds begin + if algorithm == JetAlgorithm.Valencia + return valencia_distance(eereco, i, j, invR2) + else + # Durham & EEKt share same form here (min(E2p_i,E2p_j) * dij_factor * angular_metric) + return min(eereco[i].E2p, eereco[j].E2p) * dij_factor * eereco[i].nndist + end + end +end + +""" + valencia_distance(eereco, i, j, invR2) -> Float64 + +Calculate the Valencia distance between two jets `i` and `j` as +``min(E_i^{2β}, E_j^{2β}) * 2 * (1 - cos(θ_{ij})) * invR2``. # Arguments - `eereco`: The array of `EERecoJet` objects. - `i`: The first jet. - `j`: The second jet. -- `dij_factor`: The scaling factor to multiply the dij distance by. +- `invR2`: The inverse square of the radius, i.e. ``1 / R^2``. # Returns -- The dij distance between `i` and `j`. +- `Float64`: The Valencia distance between `i` and `j`. """ -@inline function dij_dist(eereco, i, j, dij_factor) - # Calculate the dij distance for jet i from jet j - j == 0 && return large_dij - @inbounds min(eereco[i].E2p, eereco[j].E2p) * dij_factor * eereco[i].nndist +Base.@propagate_inbounds @inline function valencia_distance(eereco, i, j, invR2) + @muladd angular_dist = 1.0 - eereco[i].nx * eereco[j].nx - eereco[i].ny * eereco[j].ny - + eereco[i].nz * eereco[j].nz + return min(eereco[i].E2p, eereco[j].E2p) * 2 * angular_dist * invR2 +end + +""" + valencia_beam_distance(eereco, i, γ, β) -> Float64 + +Calculate the Valencia beam distance for jet `i` using the FastJet ValenciaPlugin +definition: ``d_iB = E_i^{2β} * (sin θ_i)^{2γ}``, where ``cos θ_i = nz`` +for unit direction cosines. Since ``sin^2 θ = 1 - nz^2``, we implement +``d_iB = E_i^{2β} * (1 - nz^2)^γ``. + +# Arguments +- `eereco`: The array of `EERecoJet` objects. +- `i`: The jet index. +- `γ`: The angular exponent parameter used in the Valencia beam distance. +- `β`: The energy exponent (same as `p` in our implementation). + +# Returns +- `Float64`: The Valencia beam distance for jet `i`. +""" +Base.@propagate_inbounds @inline function valencia_beam_distance(eereco, i, γ, β) + nz_i = eereco[i].nz + sin2 = 1 - nz_i * nz_i + E2p = eereco[i].E2p + if γ == 1.0 + return E2p * sin2 + elseif γ == 2.0 + return E2p * (sin2 * sin2) + else + return E2p * sin2^γ + end end -function get_angular_nearest_neighbours!(eereco, algorithm, dij_factor) +function get_angular_nearest_neighbours!(eereco, algorithm, dij_factor, p, invR2, γ = 1.0) # Get the initial nearest neighbours for each jet N = length(eereco) - # this_dist_vector = Vector{Float64}(undef, N) - # Nearest neighbour geometric distance + # Initialise sentinels so the first comparison always wins + @inbounds for i in 1:N + if algorithm == JetAlgorithm.Valencia + eereco.nndist[i] = Inf + else + eereco.nndist[i] = large_distance + end + eereco.nni[i] = i + end + # Nearest neighbour search @inbounds for i in 1:N - # TODO: Replace the 'j' loop with a vectorised operation over the appropriate array elements? - # this_dist_vector .= 1.0 .- eereco.nx[i:N] .* eereco[i + 1:end].nx .- - # eereco[i].ny .* eereco[i + 1:end].ny .- eereco[i].nz .* eereco[i + 1:end].nz - # The problem here will be avoiding allocations for the array outputs, which would easily - # kill performance @inbounds for j in (i + 1):N - this_nndist = angular_distance(eereco, i, j) + # Metric used to pick the nearest neighbour + if algorithm == JetAlgorithm.Valencia + # Use canonical Valencia distance (StructArray-aware) + this_metric = valencia_distance(eereco, i, j, invR2) + else + this_metric = angular_distance(eereco, i, j) + end # Using these ternary operators is faster than the if-else block - better_nndist_i = this_nndist < eereco[i].nndist - eereco.nndist[i] = better_nndist_i ? this_nndist : eereco.nndist[i] + better_nndist_i = this_metric < eereco[i].nndist + eereco.nndist[i] = better_nndist_i ? this_metric : eereco.nndist[i] eereco.nni[i] = better_nndist_i ? j : eereco.nni[i] - better_nndist_j = this_nndist < eereco[j].nndist - eereco.nndist[j] = better_nndist_j ? this_nndist : eereco.nndist[j] + better_nndist_j = this_metric < eereco[j].nndist + eereco.nndist[j] = better_nndist_j ? this_metric : eereco.nndist[j] eereco.nni[j] = better_nndist_j ? i : eereco.nni[j] end end # Nearest neighbour dij distance - for i in 1:N - eereco.dijdist[i] = dij_dist(eereco, i, eereco[i].nni, dij_factor) + @inbounds for i in 1:N + if algorithm == JetAlgorithm.Valencia + eereco.dijdist[i] = valencia_distance(eereco, i, eereco[i].nni, invR2) + else + eereco.dijdist[i] = dij_dist(eereco, i, eereco[i].nni, dij_factor, algorithm, + invR2) + end end - # For the EEKt algorithm, we need to check the beam distance as well - # (This is structured to only check for EEKt once) + # For the EEKt and Valencia algorithms, we need to check the beam distance as well + # (This is structured to check each algorithm's beam distance once) if algorithm == JetAlgorithm.EEKt @inbounds for i in 1:N beam_closer = eereco[i].E2p < eereco[i].dijdist eereco.dijdist[i] = beam_closer ? eereco[i].E2p : eereco.dijdist[i] eereco.nni[i] = beam_closer ? 0 : eereco.nni[i] end + elseif algorithm == JetAlgorithm.Valencia + @inbounds for i in 1:N + valencia_beam_dist = valencia_beam_distance(eereco, i, γ, p) + beam_closer = valencia_beam_dist < eereco[i].dijdist + eereco.dijdist[i] = beam_closer ? valencia_beam_dist : eereco.dijdist[i] + eereco.nni[i] = beam_closer ? 0 : eereco.nni[i] + end end end -# Update the nearest neighbour for jet i, w.r.t. all other active jets -function update_nn_no_cross!(eereco, i, N, algorithm, dij_factor) - eereco.nndist[i] = large_distance +@inline function update_nn_no_cross!(eereco, i, N, algorithm::JetAlgorithm.Algorithm, + dij_factor, invR2, β = 1.0, γ = 1.0) + # Valencia metric is unbounded, others use a large finite value + if algorithm == JetAlgorithm.Valencia + eereco.nndist[i] = Inf + else + eereco.nndist[i] = large_distance + end eereco.nni[i] = i + # Scan all other jets, update i, and cross-update j @inbounds for j in 1:N if j != i - this_nndist = angular_distance(eereco, i, j) - better_nndist_i = this_nndist < eereco[i].nndist - eereco.nndist[i] = better_nndist_i ? this_nndist : eereco.nndist[i] - eereco.nni[i] = better_nndist_i ? j : eereco.nni[i] + this_metric = algorithm == JetAlgorithm.Valencia ? + valencia_distance(eereco, i, j, invR2) : + angular_distance(eereco, i, j) + better_i = this_metric < eereco[i].nndist + eereco.nndist[i] = better_i ? this_metric : eereco.nndist[i] + eereco.nni[i] = better_i ? j : eereco.nni[i] end end - eereco.dijdist[i] = dij_dist(eereco, i, eereco[i].nni, dij_factor) + # Set dij for i using unified dispatcher and apply beam checks + eereco.dijdist[i] = dij_dist(eereco, i, eereco[i].nni, dij_factor, algorithm, invR2) if algorithm == JetAlgorithm.EEKt beam_close = eereco[i].E2p < eereco[i].dijdist eereco.dijdist[i] = beam_close ? eereco[i].E2p : eereco.dijdist[i] eereco.nni[i] = beam_close ? 0 : eereco.nni[i] + elseif algorithm == JetAlgorithm.Valencia + val_beam = valencia_beam_distance(eereco, i, γ, β) + beam_close = val_beam < eereco[i].dijdist + eereco.dijdist[i] = beam_close ? val_beam : eereco.dijdist[i] + eereco.nni[i] = beam_close ? 0 : eereco.nni[i] end end -function update_nn_cross!(eereco, i, N, algorithm, dij_factor) - # Update the nearest neighbour for jet i, w.r.t. all other active jets - # also doing the cross check for the other jet - eereco.nndist[i] = large_distance +@inline function update_nn_cross!(eereco, i, N, algorithm::JetAlgorithm.Algorithm, + dij_factor, invR2, β = 1.0, γ = 1.0) + # Valencia metric is unbounded, others use a large finite value + if algorithm == JetAlgorithm.Valencia + eereco.nndist[i] = Inf + else + eereco.nndist[i] = large_distance + end eereco.nni[i] = i + # Scan all other jets, update i, and cross-update j @inbounds for j in 1:N if j != i - this_nndist = angular_distance(eereco, i, j) - better_nndist_i = this_nndist < eereco[i].nndist - eereco.nndist[i] = better_nndist_i ? this_nndist : eereco.nndist[i] - eereco.nni[i] = better_nndist_i ? j : eereco.nni[i] - if this_nndist < eereco[j].nndist - eereco.nndist[j] = this_nndist + this_metric = algorithm == JetAlgorithm.Valencia ? + valencia_distance(eereco, i, j, invR2) : + angular_distance(eereco, i, j) + better_i = this_metric < eereco[i].nndist + eereco.nndist[i] = better_i ? this_metric : eereco.nndist[i] + eereco.nni[i] = better_i ? j : eereco.nni[i] + + if this_metric < eereco[j].nndist + eereco.nndist[j] = this_metric eereco.nni[j] = i - # j will not be revisited, so update metric distance here - eereco.dijdist[j] = dij_dist(eereco, j, i, dij_factor) + eereco.dijdist[j] = dij_dist(eereco, j, i, dij_factor, algorithm, invR2) if algorithm == JetAlgorithm.EEKt if eereco[j].E2p < eereco[j].dijdist eereco.dijdist[j] = eereco[j].E2p eereco.nni[j] = 0 end + elseif algorithm == JetAlgorithm.Valencia + val_beam_j = valencia_beam_distance(eereco, j, γ, β) + if val_beam_j < eereco[j].dijdist + eereco.dijdist[j] = val_beam_j + eereco.nni[j] = 0 + end end end end end - eereco.dijdist[i] = dij_dist(eereco, i, eereco[i].nni, dij_factor) + # Set dij for i using unified dispatcher and apply beam checks + eereco.dijdist[i] = dij_dist(eereco, i, eereco[i].nni, dij_factor, algorithm, invR2) if algorithm == JetAlgorithm.EEKt beam_close = eereco[i].E2p < eereco[i].dijdist eereco.dijdist[i] = beam_close ? eereco[i].E2p : eereco.dijdist[i] eereco.nni[i] = beam_close ? 0 : eereco.nni[i] + elseif algorithm == JetAlgorithm.Valencia + val_beam_i = valencia_beam_distance(eereco, i, γ, β) + beam_close = val_beam_i < eereco[i].dijdist + eereco.dijdist[i] = beam_close ? val_beam_i : eereco.dijdist[i] + eereco.nni[i] = beam_close ? 0 : eereco.nni[i] end end @@ -148,30 +254,44 @@ function ee_check_consistency(clusterseq, eereco, N) end end end - @debug "Consistency check passed at $msg" + @debug "Consistency check passed" end -function fill_reco_array!(eereco, particles, R2, p) - for i in eachindex(particles) +Base.@propagate_inbounds @inline function fill_reco_array!(eereco, particles, invR2, p) + @inbounds for i in eachindex(particles) eereco.index[i] = i eereco.nni[i] = 0 - eereco.nndist[i] = R2 + eereco.nndist[i] = inv(invR2) # R^2 as initial sentinel for angular algorithms # eereco.dijdist[i] = UNDEF # Does not need to be initialised eereco.nx[i] = nx(particles[i]) eereco.ny[i] = ny(particles[i]) eereco.nz[i] = nz(particles[i]) eereco.E2p[i] = energy(particles[i])^(2p) + # No precomputed beam factor; compute on demand to preserve previous behaviour end end -@inline function insert_new_jet!(eereco, i, newjet_k, R2, merged_jet, p) - eereco.index[i] = newjet_k - eereco.nni[i] = 0 - eereco.nndist[i] = R2 - eereco.nx[i] = nx(merged_jet) - eereco.ny[i] = ny(merged_jet) - eereco.nz[i] = nz(merged_jet) - eereco.E2p[i] = energy(merged_jet)^(2p) +Base.@propagate_inbounds @inline function insert_new_jet!(eereco, i, newjet_k, invR2, + merged_jet, p) + @inbounds begin + eereco.index[i] = newjet_k + eereco.nni[i] = 0 + eereco.nndist[i] = inv(invR2) + eereco.nx[i] = nx(merged_jet) + eereco.ny[i] = ny(merged_jet) + eereco.nz[i] = nz(merged_jet) + E = energy(merged_jet) + if p isa Int + if p == 1 + eereco.E2p[i] = E * E + else + E2 = E * E + eereco.E2p[i] = E2^p + end + else + eereco.E2p[i] = E^(2p) + end + end end """ @@ -179,21 +299,23 @@ end Copy the contents of slot `i` in the `eereco` array to slot `j`. """ -@inline function copy_to_slot!(eereco, i, j) - eereco.index[j] = eereco.index[i] - eereco.nni[j] = eereco.nni[i] - eereco.nndist[j] = eereco.nndist[i] - eereco.dijdist[j] = eereco.dijdist[i] - eereco.nx[j] = eereco.nx[i] - eereco.ny[j] = eereco.ny[i] - eereco.nz[j] = eereco.nz[i] - eereco.E2p[j] = eereco.E2p[i] +Base.@propagate_inbounds @inline function copy_to_slot!(eereco, i, j) + @inbounds begin + eereco.index[j] = eereco.index[i] + eereco.nni[j] = eereco.nni[i] + eereco.nndist[j] = eereco.nndist[i] + eereco.dijdist[j] = eereco.dijdist[i] + eereco.nx[j] = eereco.nx[i] + eereco.ny[j] = eereco.ny[i] + eereco.nz[j] = eereco.nz[i] + eereco.E2p[j] = eereco.E2p[i] + end end """ ee_genkt_algorithm(particles::AbstractVector{T}; algorithm::JetAlgorithm.Algorithm, p::Union{Real, Nothing} = nothing, R = 4.0, recombine = addjets, - preprocess = nothing) where {T} + preprocess = nothing, γ::Real = 1.0, β::Union{Real, Nothing} = nothing) where {T} Run an e+e- reconstruction algorithm on a set of initial particles. @@ -201,12 +323,16 @@ Run an e+e- reconstruction algorithm on a set of initial particles. - `particles::AbstractVector{T}`: A vector of particles to be clustered. - `algorithm::JetAlgorithm.Algorithm`: The jet algorithm to use. - `p::Union{Real, Nothing} = nothing`: The power parameter for the algorithm. - This is not required for the `Durham` algorithm, but must be specified for the - `EEKt`` algorithm. -- `R = 4.0`: The jet radius parameter. Not required and ignored for the `Durham` + Must be specified for EEKt algorithm. + For Valencia algorithm, this corresponds to the β parameter. + Other algorithms will ignore this value. +- `R = 4.0`: The jet radius parameter. Not required / ignored for the Durham algorithm. -- `recombine = addjets`: The recombination scheme to use. -- `preprocess = nothing`: Preprocessing function for input particles. +- `recombine`: The recombination scheme to use. +- `preprocess`: Preprocessing function for input particles. +- `γ::Real = 1.0`: The angular exponent parameter for Valencia algorithm. Ignored for other algorithms. +- `β::Union{Real, Nothing} = nothing`: Optional alias for the Valencia energy exponent; if provided for + Valencia it overrides `p`. # Returns - The result of the jet clustering as a `ClusterSequence` object. @@ -217,18 +343,27 @@ will check for consistency between the algorithm and the power parameter as needed. It will then prepare the internal EDM particles for the clustering itself, and call the actual reconstruction method `_ee_genkt_algorithm!`. -If the algorithm is Durham, `R` is nominally set to 4. If the algorithm is EEkt, -power `p` must also be specified. +If the algorithm is Durham, `R` is nominally set to 4. +If the algorithm is EEkt, power `p` must be specified. +If the algorithm is Valencia, you can provide `p` (β) and `γ`, or pass `β` explicitly to override `p`. """ function ee_genkt_algorithm(particles::AbstractVector{T}; algorithm::JetAlgorithm.Algorithm, p::Union{Real, Nothing} = nothing, R = 4.0, recombine = addjets, - preprocess = nothing) where {T} + preprocess = nothing, γ::Real = 1.0, + β::Union{Real, Nothing} = nothing) where {T} + + # For Valencia, if β is provided, overwrite p + if algorithm == JetAlgorithm.Valencia && β !== nothing + p = β + end # Check for consistency algorithm power p = get_algorithm_power(p = p, algorithm = algorithm) - # Integer p if possible - p = (round(p) == p) ? Int(p) : p + # Integer p if possible, i.e. if not running Valencia + if algorithm != JetAlgorithm.Valencia + p = (round(p) == p) ? Int(p) : p + end # For the Durham algorithm, p=1 and R is not used, but nominally set to 4 if algorithm == JetAlgorithm.Durham @@ -261,16 +396,19 @@ function ee_genkt_algorithm(particles::AbstractVector{T}; algorithm::JetAlgorith end end - # Now call the actual reconstruction method, tuned for our internal EDM - _ee_genkt_algorithm!(recombination_particles; p = p, R = R, - algorithm = algorithm, - recombine = recombine) + # Compute invR2 once and thread it through + invR2 = inv(R * R) + # Now call the unified implementation with conditional logic. + return _ee_genkt_algorithm(particles = recombination_particles, p = p, R = R, + invR2 = invR2, algorithm = algorithm, recombine = recombine, + γ = γ) end """ - _ee_genkt_algorithm!(particles::AbstractVector{EEJet}; - algorithm::JetAlgorithm.Algorithm, p::Real, R = 4.0, - recombine = addjets) + _ee_genkt_algorithm(particles::AbstractVector{EEJet}; + algorithm::JetAlgorithm.Algorithm, p::Real, R::Real, + invR2::Union{Real, Nothing} = nothing, recombine = addjets, γ::Real = 1.0, + beta::Union{Real, Nothing} = nothing) This function is the internal implementation of the e+e- jet clustering algorithm. It takes a vector of `EEJet` `particles` representing the input @@ -280,7 +418,7 @@ Users of the package should use the `ee_genkt_algorithm` function as their entry point to this jet reconstruction. # Arguments -- `particles::AbstractVector{EEJet}`: A vector of `EEJet` particles used + - `particles::AbstractVector{EEJet}`: A vector of `EEJet` particles used as input for jet reconstruction. This vector must supply the correct `cluster_hist_index` values and will be *mutated* as part of the returned `ClusterSequence`. @@ -289,19 +427,29 @@ entry point to this jet reconstruction. is raised. - `R = 4.0`: The jet radius parameter. - `recombine = addjets`: The recombination function used to merge two jets. + - `γ::Real = 1.0`: Angular exponent for the Valencia beam metric (ignored for other algorithms). + - `beta::Union{Real, Nothing} = nothing`: Optional alias for the Valencia energy exponent (β). + When provided with `algorithm == JetAlgorithm.Valencia`, it overrides `p`. # Returns - `clusterseq`: The resulting `ClusterSequence` object representing the reconstructed jets. """ -function _ee_genkt_algorithm!(particles::AbstractVector{EEJet}; - algorithm::JetAlgorithm.Algorithm, p::Real, R::Real = 4.0, - recombine = addjets) +function _ee_genkt_algorithm(; particles::AbstractVector{EEJet}, + algorithm::JetAlgorithm.Algorithm, p::Real, R::Real, + invR2::Union{Real, Nothing} = nothing, recombine = addjets, + γ::Real = 1.0, + beta::Union{Real, Nothing} = nothing) # Bounds N::Int = length(particles) - # R squared - R2 = R^2 + # invR2 provided by caller when available; otherwise compute from R once here + if invR2 === nothing + invR2 = inv(R * R) + end + if algorithm == JetAlgorithm.Valencia && beta !== nothing + p = beta + end # Constant factor for the dij metric and the beam distance function if algorithm == JetAlgorithm.Durham @@ -312,6 +460,8 @@ function _ee_genkt_algorithm!(particles::AbstractVector{EEJet}; else dij_factor = 1 / (3 + cos(R)) end + elseif algorithm == JetAlgorithm.Valencia + dij_factor = 1.0 else throw(ArgumentError("Algorithm $algorithm not supported for e+e-")) end @@ -320,7 +470,8 @@ function _ee_genkt_algorithm!(particles::AbstractVector{EEJet}; # jet information and populate it accordingly # We need N slots for this array eereco = StructArray{EERecoJet}(undef, N) - fill_reco_array!(eereco, particles, R2, p) + + fill_reco_array!(eereco, particles, invR2, p) # Setup the initial history and get the total energy history, Qtot = initial_history(particles) @@ -329,7 +480,7 @@ function _ee_genkt_algorithm!(particles::AbstractVector{EEJet}; Qtot) # Run over initial pairs of jets to find nearest neighbours - get_angular_nearest_neighbours!(eereco, algorithm, dij_factor) + get_angular_nearest_neighbours!(eereco, algorithm, dij_factor, p, invR2, γ) # Only for debugging purposes... # ee_check_consistency(clusterseq, clusterseq_index, N, nndist, nndij, nni, "Start") @@ -378,7 +529,7 @@ function _ee_genkt_algorithm!(particles::AbstractVector{EEJet}; newjet_k, dij_min) # Update the compact arrays, reusing the JetA slot - insert_new_jet!(eereco, ijetA, newjet_k, R2, merged_jet, p) + insert_new_jet!(eereco, ijetA, newjet_k, invR2, merged_jet, p) end # Squash step - copy the final jet's compact data into the jetB slot @@ -400,7 +551,7 @@ function _ee_genkt_algorithm!(particles::AbstractVector{EEJet}; # plus "belt and braces" check for an invalid NN (>N) if (eereco[i].nni == ijetA) || (eereco[i].nni == ijetB) || (eereco[i].nni > N) - update_nn_no_cross!(eereco, i, N, algorithm, dij_factor) + update_nn_no_cross!(eereco, i, N, algorithm, dij_factor, invR2, p, γ) end end end @@ -408,7 +559,7 @@ function _ee_genkt_algorithm!(particles::AbstractVector{EEJet}; # Finally, we need to update the nearest neighbours for the new jet, checking both ways # (But only if there was a new jet!) if ijetA != ijetB - update_nn_cross!(eereco, ijetA, N, algorithm, dij_factor) + update_nn_cross!(eereco, ijetA, N, algorithm, dij_factor, invR2, p, γ) end # Only for debugging purposes... diff --git a/src/GenericAlgo.jl b/src/GenericAlgo.jl index 8227c8ed..380e5a84 100644 --- a/src/GenericAlgo.jl +++ b/src/GenericAlgo.jl @@ -2,7 +2,8 @@ jet_reconstruct(particles::AbstractVector; algorithm::JetAlgorithm.Algorithm, p::Union{Real, Nothing} = nothing, R = 1.0, recombine = addjets, preprocess = nothing, - strategy::RecoStrategy.Strategy = RecoStrategy.Best) + strategy::RecoStrategy.Strategy = RecoStrategy.Best, + γ::Real = 1.0) Reconstructs jets from a collection of particles using a specified algorithm and strategy. @@ -22,6 +23,8 @@ strategy. - `strategy::RecoStrategy.Strategy = RecoStrategy.Best`: The jet reconstruction strategy to use. `RecoStrategy.Best` makes a dynamic decision based on the number of starting particles. +- `γ::Real = 1.0`: The angular exponent parameter for Valencia algorithm. Ignored + by other algorithms. Note that `p` must be specified for `GenKt` and `EEKt` algorithms, other algorithms will ignore its value. @@ -67,7 +70,8 @@ jet_reconstruct(particles; algorithm = JetAlgorithm.AntiKt, R = 1.0, preprocess function jet_reconstruct(particles::AbstractVector; algorithm::JetAlgorithm.Algorithm, p::Union{Real, Nothing} = nothing, R = 1.0, recombine = addjets, preprocess = nothing, - strategy::RecoStrategy.Strategy = RecoStrategy.Best) + strategy::RecoStrategy.Strategy = RecoStrategy.Best, + γ::Real = 1.0) if is_pp(algorithm) # We assume a pp reconstruction if strategy == RecoStrategy.Best @@ -87,6 +91,11 @@ function jet_reconstruct(particles::AbstractVector; algorithm::JetAlgorithm.Algo end # Now call the chosen algorithm, passing through the other parameters - alg(particles; p = p, algorithm = algorithm, R = R, recombine = recombine, - preprocess = preprocess) + if is_ee(algorithm) + alg(particles; p, algorithm, R, recombine, + preprocess, γ) + else + alg(particles; p, algorithm, R, recombine, + preprocess) + end end diff --git a/test/_common.jl b/test/_common.jl index 63dd6b5a..2a783fb7 100644 --- a/test/_common.jl +++ b/test/_common.jl @@ -10,6 +10,11 @@ using Logging using LorentzVectorHEP using JSON using Test +using StructArrays + +using JetReconstruction: EERecoJet, + fill_reco_array!, insert_new_jet!, copy_to_slot!, + dij_dist, valencia_distance logger = ConsoleLogger(stdout, Logging.Warn) global_logger(logger) @@ -72,6 +77,31 @@ function ComparisonTest(events_file, fastjet_outputs, algorithm, strategy, power selector_name, addjets, nothing) end +""" + struct ComparisonTestValencia + +Test parameters for Valencia jet algorithm comparison, with explicit gamma parameter. +""" +struct ComparisonTestValencia + events_file::AbstractString + fastjet_outputs::AbstractString + algorithm::JetAlgorithm.Algorithm + strategy::RecoStrategy.Strategy + power::Real + R::Real + γ::Real + selector::Any + selector_name::AbstractString + recombine::Any + preprocess::Any +end + +function ComparisonTestValencia(events_file, fastjet_outputs, algorithm, strategy, power, R, + γ, selector, selector_name) + ComparisonTestValencia(events_file, fastjet_outputs, algorithm, strategy, power, R, γ, + selector, selector_name, addjets, nothing) +end + """Read JSON file with fastjet jets in it""" function read_fastjet_outputs(fname) f = JetReconstruction.open_with_stream(fname) @@ -107,11 +137,24 @@ function run_reco_test(test::ComparisonTest; testname = nothing) # Run the jet reconstruction jet_collection = Vector{FinalJets}() for (ievent, event) in enumerate(events) - cluster_seq = JetReconstruction.jet_reconstruct(event; R = test.R, p = test.power, - algorithm = test.algorithm, - strategy = test.strategy, - recombine = test.recombine, - preprocess = test.preprocess) + if test.algorithm == JetAlgorithm.Valencia + # For VLC: pass both beta (power) and γ + cluster_seq = JetReconstruction.jet_reconstruct(event; R = test.R, + p = test.power, + γ = getfield(test, :γ), + algorithm = test.algorithm, + strategy = test.strategy, + recombine = test.recombine, + preprocess = test.preprocess) + else + # All other algorithms + cluster_seq = JetReconstruction.jet_reconstruct(event; R = test.R, + p = test.power, + algorithm = test.algorithm, + strategy = test.strategy, + recombine = test.recombine, + preprocess = test.preprocess) + end finaljets = final_jets(test.selector(cluster_seq)) sort_jets!(finaljets) push!(jet_collection, FinalJets(ievent, finaljets)) @@ -147,3 +190,95 @@ function run_reco_test(test::ComparisonTest; testname = nothing) end end end + +function run_reco_test(test::ComparisonTestValencia; testname = nothing) + events = JetReconstruction.read_final_state_particles(test.events_file) + fastjet_jets = read_fastjet_outputs(test.fastjet_outputs) + sort_jets!(fastjet_jets) + + jet_collection = Vector{FinalJets}() + for (ievent, event) in enumerate(events) + cluster_seq = JetReconstruction.jet_reconstruct(event; R = test.R, p = test.power, + γ = test.γ, + algorithm = test.algorithm, + strategy = test.strategy, + recombine = test.recombine, + preprocess = test.preprocess) + finaljets = final_jets(test.selector(cluster_seq)) + sort_jets!(finaljets) + push!(jet_collection, FinalJets(ievent, finaljets)) + end + + if isnothing(testname) + testname = "FastJet comparison: alg=$(test.algorithm), p=$(test.power), R=$(test.R), strategy=$(test.strategy)" + if test.selector_name != "" + testname *= ", $(test.selector_name)" + end + end + + @testset "$testname" begin + for (ievt, event) in enumerate(jet_collection) + @testset "Event $(ievt)" begin + @test size(event.jets) == size(fastjet_jets[ievt]["jets"]) + # For each reconstructed jet, find the closest reference jet in deltaR + ref_jets = fastjet_jets[ievt]["jets"] + used_refs = falses(length(ref_jets)) + function deltaR(rap1, phi1, rap2, phi2) + dphi = abs(phi1 - phi2) + if dphi > π + dphi -= 2π + end + return sqrt((rap1 - rap2)^2 + dphi^2) + end + for (ijet, jet) in enumerate(event.jets) + # Find closest reference jet in deltaR that hasn't been used + minidx = 0 + mindr = Inf + norm_phi = jet.phi < 0.0 ? jet.phi + 2π : jet.phi + for (ridx, rjet) in enumerate(ref_jets) + if used_refs[ridx] + continue + end + norm_phi_ref = rjet["phi"] < 0.0 ? rjet["phi"] + 2π : rjet["phi"] + dr = deltaR(jet.rap, norm_phi, rjet["rap"], norm_phi_ref) + if dr < mindr + mindr = dr + minidx = ridx + end + end + if minidx == 0 + println("No unused reference jet found for reconstructed jet $(ijet)") + continue + end + used_refs[minidx] = true + rjet = ref_jets[minidx] + rap_ref = rjet["rap"] + phi_ref = rjet["phi"] + pt_ref = rjet["pt"] + normalised_phi_ref = phi_ref < 0.0 ? phi_ref + 2π : phi_ref + rap_test = isapprox(jet.rap, rap_ref; atol = 1e-4) + phi_test = isapprox(norm_phi, normalised_phi_ref; atol = 1e-4) + pt_test = isapprox(jet.pt, pt_ref; rtol = 1e-5) + if !rap_test || !phi_test || !pt_test + println("Jet mismatch in Event $(ievt), Jet $(ijet):") + println(" Failing Jet: pt=$(jet.pt), rap=$(jet.rap), phi=$(norm_phi)") + println(" Reference Jet: pt=$(pt_ref), rap=$(rap_ref), phi=$(phi_ref)") + println(" Passes: pt=$(pt_test), rap=$(rap_test), phi=$(phi_test)") + println("\nAll jets in this event:") + for (jidx, j) in enumerate(event.jets) + norm_phi_j = j.phi < 0.0 ? j.phi + 2π : j.phi + println(" Jet $(jidx): pt=$(j.pt), rap=$(j.rap), phi=$(norm_phi_j)") + end + println("\nAll reference jets in this event:") + for (ridx, rjet2) in enumerate(ref_jets) + println(" Ref Jet $(ridx): pt=$(rjet2["pt"]), rap=$(rjet2["rap"]), phi=$(rjet2["phi"])") + end + end + @test rap_test + @test phi_test + @test pt_test + end + end + end + end +end diff --git a/test/data/jet-collections-fastjet-valencia-exclusive-d500-beta1.0-gamma1.0-R1.0-eeH.json.zst b/test/data/jet-collections-fastjet-valencia-exclusive-d500-beta1.0-gamma1.0-R1.0-eeH.json.zst new file mode 100644 index 00000000..e366cf13 Binary files /dev/null and b/test/data/jet-collections-fastjet-valencia-exclusive-d500-beta1.0-gamma1.0-R1.0-eeH.json.zst differ diff --git a/test/data/jet-collections-fastjet-valencia-exclusive-d500-beta1.2-gamma1.2-R0.8-eeH.json.zst b/test/data/jet-collections-fastjet-valencia-exclusive-d500-beta1.2-gamma1.2-R0.8-eeH.json.zst new file mode 100644 index 00000000..3728583f Binary files /dev/null and b/test/data/jet-collections-fastjet-valencia-exclusive-d500-beta1.2-gamma1.2-R0.8-eeH.json.zst differ diff --git a/test/data/jet-collections-fastjet-valencia-exclusive4-beta1.0-gamma1.0-R1.0-eeH.json.zst b/test/data/jet-collections-fastjet-valencia-exclusive4-beta1.0-gamma1.0-R1.0-eeH.json.zst new file mode 100644 index 00000000..05edc08b Binary files /dev/null and b/test/data/jet-collections-fastjet-valencia-exclusive4-beta1.0-gamma1.0-R1.0-eeH.json.zst differ diff --git a/test/data/jet-collections-fastjet-valencia-exclusive4-beta1.2-gamma1.2-R0.8-eeH.json.zst b/test/data/jet-collections-fastjet-valencia-exclusive4-beta1.2-gamma1.2-R0.8-eeH.json.zst new file mode 100644 index 00000000..9eb753e3 Binary files /dev/null and b/test/data/jet-collections-fastjet-valencia-exclusive4-beta1.2-gamma1.2-R0.8-eeH.json.zst differ diff --git a/test/data/jet-collections-fastjet-valencia-inclusive20-beta1.0-gamma1.0-R1.0-eeH.json.zst b/test/data/jet-collections-fastjet-valencia-inclusive20-beta1.0-gamma1.0-R1.0-eeH.json.zst new file mode 100644 index 00000000..725d85fc Binary files /dev/null and b/test/data/jet-collections-fastjet-valencia-inclusive20-beta1.0-gamma1.0-R1.0-eeH.json.zst differ diff --git a/test/data/jet-collections-fastjet-valencia-inclusive20-beta1.2-gamma1.2-R0.8-eeH.json.zst b/test/data/jet-collections-fastjet-valencia-inclusive20-beta1.2-gamma1.2-R0.8-eeH.json.zst new file mode 100644 index 00000000..871a1793 Binary files /dev/null and b/test/data/jet-collections-fastjet-valencia-inclusive20-beta1.2-gamma1.2-R0.8-eeH.json.zst differ diff --git a/test/runtests.jl b/test/runtests.jl index 4256853f..5390a039 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -19,6 +19,9 @@ function main() include("test-pp-reconstruction.jl") include("test-ee-reconstruction.jl") + # Valencia algorithm tests + include("test-valencia.jl") + # Test jets selection include("test-selection.jl") diff --git a/test/test-ee-reconstruction.jl b/test/test-ee-reconstruction.jl index 0a942eec..83a583f3 100644 --- a/test/test-ee-reconstruction.jl +++ b/test/test-ee-reconstruction.jl @@ -26,3 +26,57 @@ for r in [2.0, 4.0] (cs) -> exclusive_jets(cs; njets = 4), "exclusive njets") run_reco_test(eekt_njets) end + +# Optimization/helper coverage + +@testset "E2p integer fast paths in fill/insert (e+e-)" begin + E1, E2 = 3.0, 2.0 + # PseudoJet signature is (px, py, pz, E); give nonzero momentum to avoid 1/0 in EEJet + pj1 = PseudoJet(1.0, 0.0, 0.0, E1) + pj2 = PseudoJet(1.0, 0.0, 0.0, E2) + particles = EEJet[] + push!(particles, EEJet(pj1; cluster_hist_index = 1)) + push!(particles, EEJet(pj2; cluster_hist_index = 2)) + + eereco = StructArray{EERecoJet}(undef, 2) + R = 1.0 + R2 = R^2 + + fill_reco_array!(eereco, particles, R2, 1) + @test eereco.E2p[1] ≈ E1^2 + @test eereco.E2p[2] ≈ E2^2 + + fill_reco_array!(eereco, particles, R2, 2) + @test eereco.E2p[1] ≈ E1^4 + @test eereco.E2p[2] ≈ E2^4 + + merged = EEJet(1.0, 0.0, 0.0, 5.0; cluster_hist_index = 0) + insert_new_jet!(eereco, 1, 3, R2, merged, 1) + @test eereco.E2p[1] ≈ 5.0^2 + + insert_new_jet!(eereco, 2, 4, R2, merged, 2) + @test eereco.E2p[2] ≈ 5.0^4 +end + +@testset "copy_to_slot! copies fields" begin + eereco = StructArray{EERecoJet}(undef, 2) + eereco.index[1] = 10 + eereco.nni[1] = 1 + eereco.nndist[1] = 3.14 + eereco.dijdist[1] = 2.71 + eereco.nx[1] = 0.1 + eereco.ny[1] = 0.2 + eereco.nz[1] = 0.3 + eereco.E2p[1] = 7.0 + copy_to_slot!(eereco, 1, 2) + @test eereco.index[2] == 10 + @test eereco.dijdist[2] == 2.71 + @test eereco.E2p[2] == 7.0 +end + +@testset "dij_dist fallback error (non-Algorithm)" begin + eereco = StructArray{EERecoJet}(undef, 1) + eereco.nndist[1] = 0.0 + eereco.E2p[1] = 1.0 + @test_throws ArgumentError dij_dist(eereco, 1, 0, 1.0, :Foo, 1.0) +end diff --git a/test/test-valencia.jl b/test/test-valencia.jl new file mode 100644 index 00000000..35e19387 --- /dev/null +++ b/test/test-valencia.jl @@ -0,0 +1,103 @@ +""" + test-valencia.jl + +Test for the Valencia jet algorithm implementation. + +Validates clustering results and cluster sequence against FastJet C++ reference output +using jet matching to handle numerical precision differences. +""" + +include("common.jl") + +# Reference outputs for R=0.8, beta=1.2, gamma=1.2 +valencia_inclusive = ComparisonTestValencia(events_file_ee, + joinpath(@__DIR__, "data", + "jet-collections-fastjet-valencia-inclusive20-beta1.2-gamma1.2-R0.8-eeH.json.zst"), + JetAlgorithm.Valencia, RecoStrategy.N2Plain, + 1.2, 0.8, 1.2, + cs -> inclusive_jets(cs; ptmin = 20.0), + "inclusive beta=1.2 gamma=1.2 R=0.8") +run_reco_test(valencia_inclusive) + +valencia_exclusive4 = ComparisonTestValencia(events_file_ee, + joinpath(@__DIR__, "data", + "jet-collections-fastjet-valencia-exclusive4-beta1.2-gamma1.2-R0.8-eeH.json.zst"), + JetAlgorithm.Valencia, RecoStrategy.N2Plain, + 1.2, 0.8, 1.2, + cs -> exclusive_jets(cs; njets = 4), + "exclusive njets beta=1.2 gamma=1.2 R=0.8") +run_reco_test(valencia_exclusive4) + +valencia_exclusive_d500 = ComparisonTestValencia(events_file_ee, + joinpath(@__DIR__, "data", + "jet-collections-fastjet-valencia-exclusive-d500-beta1.2-gamma1.2-R0.8-eeH.json.zst"), + JetAlgorithm.Valencia, + RecoStrategy.N2Plain, 1.2, 0.8, 1.2, + cs -> exclusive_jets(cs; dcut = 500.0), + "exclusive dcut=500 beta=1.2 gamma=1.2 R=0.8") +run_reco_test(valencia_exclusive_d500) + +#Reference outputs for R=1.0, beta=1.0, gamma=1.0 +valencia_inclusive_b1g1 = ComparisonTestValencia(events_file_ee, + joinpath(@__DIR__, "data", + "jet-collections-fastjet-valencia-inclusive20-beta1.0-gamma1.0-R1.0-eeH.json.zst"), + JetAlgorithm.Valencia, + RecoStrategy.N2Plain, 1.0, 1.0, 1.0, + cs -> inclusive_jets(cs; ptmin = 20.0), + "inclusive beta=1.0 gamma=1.0 R=1.0") +run_reco_test(valencia_inclusive_b1g1) + +valencia_exclusive4_b1g1 = ComparisonTestValencia(events_file_ee, + joinpath(@__DIR__, "data", + "jet-collections-fastjet-valencia-exclusive4-beta1.0-gamma1.0-R1.0-eeH.json.zst"), + JetAlgorithm.Valencia, + RecoStrategy.N2Plain, 1.0, 1.0, 1.0, + cs -> exclusive_jets(cs; njets = 4), + "exclusive njets beta=1.0 gamma=1.0 R=1.0") +run_reco_test(valencia_exclusive4_b1g1) + +valencia_exclusive_d500_b1g1 = ComparisonTestValencia(events_file_ee, + joinpath(@__DIR__, "data", + "jet-collections-fastjet-valencia-exclusive-d500-beta1.0-gamma1.0-R1.0-eeH.json.zst"), + JetAlgorithm.Valencia, + RecoStrategy.N2Plain, 1.0, 1.0, 1.0, + cs -> exclusive_jets(cs; dcut = 500.0), + "exclusive dcut=500 beta=1.0 gamma=1.0 R=1.0") +run_reco_test(valencia_exclusive_d500_b1g1) + +# Test dij_dist for Valencia algorithm +@testset "dij_dist Valencia" begin + # Minimal eereco with two reco jets; only nx,ny,nz,E2p are used by valencia_distance + eereco = EERecoJet[EERecoJet(1, 0, Inf, Inf, 1.0, 0.0, 0.0, 1.0), + EERecoJet(2, 0, Inf, Inf, 0.0, 1.0, 0.0, 2.0)] + dij = dij_dist(eereco, 1, 2, 1.0, JetAlgorithm.Valencia, 0.8) + @test dij ≈ valencia_distance(eereco, 1, 2, 0.8) +end + +# Valencia distance wrapper coverage +@testset "Valencia distance wrappers" begin + # Minimal eereco with two reco jets and identical directions so angle=0 + eereco = EERecoJet[EERecoJet(1, 0, Inf, Inf, 1.0, 0.0, 0.0, 9.0), + EERecoJet(2, 0, Inf, Inf, 1.0, 0.0, 0.0, 4.0)] + R = 2.0 + @test valencia_distance(eereco, 1, 2, R) == 0.0 + @test valencia_distance(eereco, 1, 2, R) == 0.0 +end + +# Test ee_genkt_algorithm for Valencia algorithm (covers β override and dij_factor selection) +@testset "ee_genkt_algorithm Valencia" begin + particles = [PseudoJet(1.0, 0.0, 0.0, 1.0)] + cs = ee_genkt_algorithm(particles; algorithm = JetAlgorithm.Valencia, β = 1.2, γ = 1.2, + R = 0.8) + @test cs isa JetReconstruction.ClusterSequence +end + +# Test internal _ee_genkt_algorithm entry (touches StructArray init path too) +@testset "_ee_genkt_algorithm Valencia dij_factor" begin + # Ensure cluster_hist_index(i) == i for initial jets + jets = [EEJet(1.0, 0.0, 0.0, 1.0; cluster_hist_index = 1)] + cs = JetReconstruction._ee_genkt_algorithm(particles = jets, + algorithm = JetAlgorithm.Valencia, p = 1.2, + R = 0.8, γ = 1.2) + @test cs isa JetReconstruction.ClusterSequence +end