Skip to content

Commit

Permalink
Documentation improvements and clean up
Browse files Browse the repository at this point in the history
Document Utils.jl. Remove unused pseudojets2vectors() method.
  • Loading branch information
graeme-a-stewart committed Jun 7, 2024
1 parent f58578d commit 5c79f93
Show file tree
Hide file tree
Showing 2 changed files with 94 additions and 45 deletions.
2 changes: 1 addition & 1 deletion src/JetReconstruction.jl
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,7 @@ export savejets, loadjets!, loadjets

# utility functions, useful for different primary scripts
include("Utils.jl")
export read_final_state_particles, read_final_state_particles_lv, pseudojets2vectors, final_jets
export read_final_state_particles, read_final_state_particles_lv, final_jets

# Jet visualisation is an extension, see ext/JetVisualisation.jl
function jetsplot() end
Expand Down
137 changes: 93 additions & 44 deletions src/Utils.jl
Original file line number Diff line number Diff line change
@@ -1,6 +1,23 @@
# Utility functions, which can be used by different top level scripts

"""Read HepMC3 event and keep final state particles (return PseudoJets)"""

"""
read_final_state_particles(fname; maxevents = -1, skipevents = 0)
Reads final state particles from a file and returns them as a vector of
PseudoJet objects.
# Arguments
- `fname`: The name of the file to read particles from.
- `maxevents=-1`: The maximum number of events to read. -1 means all events will
be read.
- `skipevents`: The number of events to skip before an event is included.
Default is 0.
# Returns
A vector of vectors of PseudoJet objects, where each inner vector represents all
the particles of a particular event.
"""
function read_final_state_particles(fname; maxevents = -1, skipevents = 0)
f = open(fname)

Expand All @@ -27,7 +44,23 @@ function read_final_state_particles(fname; maxevents = -1, skipevents = 0)
events
end

"""Read HepMC3 event and keep final state particles (return LorentzVectors)"""
"""
read_final_state_particles_lv(fname; maxevents = -1, skipevents = 0)
Reads final state particles from a file and returns them as a vector of
LorentzVector objects.
# Arguments
- `fname`: The name of the file to read particles from.
- `maxevents=-1`: The maximum number of events to read. -1 means all events will
be read.
- `skipevents`: The number of events to skip before an event is included.
Default is 0.
# Returns
A vector of vectors of LorentzVector objects, where each inner vector represents
all the particles of a particular event.
"""
function read_final_state_particles_lv(fname; maxevents = -1, skipevents = 0)
f = open(fname)

Expand All @@ -54,45 +87,46 @@ function read_final_state_particles_lv(fname; maxevents = -1, skipevents = 0)
events
end

function pseudojets2vectors(events::Vector{Vector{PseudoJet}})
"""Convert events in PseudoJets into deep vectors"""
event_vector = Vector{Vector{Vector{Float64}}}(undef, size(events)[1])
for (ievent, event) in enumerate(events)
jet_struct = Vector{Vector{Float64}}(undef, size(event)[1])
for (ipart, initial_particle) in enumerate(event)
jet_struct[ipart] = [
initial_particle.px,
initial_particle.py,
initial_particle.pz,
initial_particle.E,
]
end
event_vector[ievent] = jet_struct
end
event_vector
end


"""
Return the list of jets passing a pt cut
The ptmin cut in these functions is slightly legacy as often the
input jets were already filtered on pt
"""
function final_jets(jets::Vector{Vector{Float64}}, ptmin::AbstractFloat=0.0)
count = 0
final_jets = Vector{FinalJet}()
sizehint!(final_jets, 6)
for jet in jets
dcut = ptmin^2
p = PseudoJet(jet[1], jet[2], jet[3], jet[4])
if p._pt2 > dcut
count += 1
push!(final_jets, FinalJet(rap(p), phi(p), sqrt(pt2(p))))
end
end
final_jets
end

# function final_jets(jets::Vector{Vector{Float64}}, ptmin::AbstractFloat=0.0)
# count = 0
# final_jets = Vector{FinalJet}()
# sizehint!(final_jets, 6)
# for jet in jets
# dcut = ptmin^2
# p = PseudoJet(jet[1], jet[2], jet[3], jet[4])
# if p._pt2 > dcut
# count += 1
# push!(final_jets, FinalJet(rap(p), phi(p), sqrt(pt2(p))))
# end
# end
# final_jets
# end

"""
final_jets(jets::Vector{PseudoJet}, ptmin::AbstractFloat=0.0)
This function takes a vector of `PseudoJet` objects and a minimum transverse
momentum `ptmin` as input. It returns a vector of `FinalJet` objects that
satisfy the transverse momentum condition.
# Arguments
- `jets::Vector{PseudoJet}`: A vector of `PseudoJet` objects representing the
input jets.
- `ptmin::AbstractFloat=0.0`: The minimum transverse momentum required for a jet
to be included in the final jets vector.
# Returns
A vector of `FinalJet` objects that satisfy the transverse momentum condition.
"""
function final_jets(jets::Vector{PseudoJet}, ptmin::AbstractFloat=0.0)
count = 0
final_jets = Vector{FinalJet}()
Expand All @@ -107,6 +141,7 @@ function final_jets(jets::Vector{PseudoJet}, ptmin::AbstractFloat=0.0)
final_jets
end

"""Specialisation for final jets from LorentzVectors (TODO: merge into more general function)"""
function final_jets(jets::Vector{LorentzVector}, ptmin::AbstractFloat=0.0)
count = 0
final_jets = Vector{FinalJet}()
Expand All @@ -121,6 +156,7 @@ function final_jets(jets::Vector{LorentzVector}, ptmin::AbstractFloat=0.0)
final_jets
end

"""Specialisation for final jets from LorentzVectorCyl (TODO: merge into more general function)"""
function final_jets(jets::Vector{LorentzVectorCyl}, ptmin::AbstractFloat=0.0)
count = 0
final_jets = Vector{FinalJet}()
Expand All @@ -135,16 +171,29 @@ function final_jets(jets::Vector{LorentzVectorCyl}, ptmin::AbstractFloat=0.0)
final_jets
end

## Faster findmin function
"""Find the lowest value in the array, returning the value and the index"""

"""
fast_findmin(dij, n)
Find the minimum value and its index in the first `n` elements of the `dij`
array. The use of `@turbo` macro gives a significiant performance boost.
# Arguments
- `dij`: An array of values.
- `n`: The number of elements to consider in the `dij` array.
# Returns
- `dij_min`: The minimum value in the first `n` elements of the `dij` array.
- `best`: The index of the minimum value in the `dij` array.
"""
fast_findmin(dij, n) = begin
# findmin(@inbounds @view dij[1:n])
best = 1
@inbounds dij_min = dij[1]
@turbo for here in 2:n
newmin = dij[here] < dij_min
best = newmin ? here : best
dij_min = newmin ? dij[here] : dij_min
end
dij_min, best
# findmin(@inbounds @view dij[1:n])
best = 1
@inbounds dij_min = dij[1]
@turbo for here in 2:n
newmin = dij[here] < dij_min
best = newmin ? here : best
dij_min = newmin ? dij[here] : dij_min
end
dij_min, best
end

0 comments on commit 5c79f93

Please sign in to comment.