Implement active jet area calculation with explicit ghosts (fixes #124)#240
Implement active jet area calculation with explicit ghosts (fixes #124)#240HarshitNagpal29 wants to merge 4 commits into
Conversation
Codecov Report✅ All modified and coverable lines are covered by tests. Additional details and impacted files@@ Coverage Diff @@
## main #240 +/- ##
==========================================
+ Coverage 82.67% 84.23% +1.55%
==========================================
Files 21 22 +1
Lines 1403 1541 +138
==========================================
+ Hits 1160 1298 +138
Misses 243 243 ☔ View full report in Codecov by Sentry. 🚀 New features to boost your workflow:
|
|
|
||
| has_dangerous_particles = false | ||
| if max_ghost_pt2 > 0.0 | ||
| danger_ratio = eps(Float64)^2 |
There was a problem hiding this comment.
also, can you please guide on this its little bit confusing for me that is it fine or overaggressive?
m-fila
left a comment
There was a problem hiding this comment.
A few comments from me, I didn't check the "algorithmic" part yet
| n_ghosts::Int | ||
| end | ||
|
|
||
| function _ghost_grid_spec(spec::GhostedAreaSpec)::GhostGridSpec |
There was a problem hiding this comment.
| function _ghost_grid_spec(spec::GhostedAreaSpec)::GhostGridSpec | |
| function _ghost_grid_spec(spec::GhostedAreaSpec) |
Casting to GhostGridSpec isn't really needed here
Also, would it make sense to make this function a constructor of GhostGridSpec?
| function _copy_as_pseudojet(jet; cluster_hist_index::Int) | ||
| PseudoJet(px(jet), | ||
| py(jet), | ||
| pz(jet), | ||
| E(jet); | ||
| cluster_hist_index = cluster_hist_index,) | ||
| end |
There was a problem hiding this comment.
We have a generic constructor for PseudoJet
JetReconstruction.jl/src/PseudoJet.jl
Line 152 in 3e96272
I think there is no need to introduce a new function, the
PseudoJet already has a constructor accepting any jet
| if algorithm ∉ (JetAlgorithm.Kt, | ||
| JetAlgorithm.CA, | ||
| JetAlgorithm.AntiKt, | ||
| JetAlgorithm.GenKt) | ||
| throw(ArgumentError("jet_reconstruct_area currently supports pp algorithms only")) | ||
| end |
There was a problem hiding this comment.
| if algorithm ∉ (JetAlgorithm.Kt, | |
| JetAlgorithm.CA, | |
| JetAlgorithm.AntiKt, | |
| JetAlgorithm.GenKt) | |
| throw(ArgumentError("jet_reconstruct_area currently supports pp algorithms only")) | |
| end | |
| if !is_pp(algorithm) | |
| throw(ArgumentError("jet_reconstruct_area currently supports pp algorithms only")) | |
| end |
We have a function checking if algorithm is pp
| jets = filter(jet -> !is_pure_ghost(jet, csa), jets) | ||
| end | ||
|
|
||
| return [_convert_area_result_jet(jet, T) for jet in jets] |
There was a problem hiding this comment.
I think it could would be nice to avoid allocating a new array if T is already PseudoJet and there is nothing to do
| function _make_area_ghosts(spec::GhostedAreaSpec) | ||
| grid = _ghost_grid_spec(spec) | ||
| rng = isnothing(spec.seed) ? Random.default_rng() : Random.MersenneTwister(spec.seed) |
There was a problem hiding this comment.
More Julia-idiomatic way would be to pass an rng::AbstractRNG to a function rather than create it inside
Dear @graeme-a-stewart, @m-fila, @mattleblanc,
I have tried implementing an initial solution for #124. This PR adds a first version of active jet area calculation using explicit ghost particles, inspired by the FastJet implementation and written to fit the existing style of this package.
What This Adds
GhostedAreaSpecfor configuring the ghost gridClusterSequenceAreaas a wrapper aroundClusterSequencejet_reconstruct_area(...)area(jet, csa)area_4vector(jet, csa)is_pure_ghost(jet, csa)total_areaand ghost metadata helpersCurrent Scope
This first implementation supports pp algorithms with explicit ghost active areas. It intentionally does not yet implement passive areas, Voronoi areas, rho estimation, or pile-up subtraction.
For now,
jet_reconstruct_areaonly supportsRecoStrategy.N2Plain.Note About N2Tiled
While adding tests for jet areas, I hit an existing issue with the
N2Tiledstrategy when used with ghosted areas:My current understanding is that the tiling matrix is allocated as (n_tiles_eta, n_tiles_phi), but tile_index can allow iphi == n_tiles_phi, which can produce a linear index beyond the allocated matrix size. In my case the maximum valid index should be 90, but it reached 94.
Because of this, the area implementation currently rejects tiled strategies and uses only N2Plain.
Testing
I added tests for the new area functionality. On my machine, the full test suite passes
Docs
I have not updated the docs yet, because I thought it would be better to get feedback on the API and implementation approach first.
AI usage
AI was used to help understand parts of the FastJet implementation and to generate some of the tests. I reviewed the generated material and take responsibility for the code and tests included in this PR.