Skip to content

Commit

Permalink
Allow read_final_state_particles to return different types
Browse files Browse the repository at this point in the history
A type can be passed in as a parameter, and the function will
try to construct objects of that
type to return
  • Loading branch information
graeme-a-stewart committed Jul 16, 2024
1 parent 8de4986 commit 60e3bd7
Show file tree
Hide file tree
Showing 4 changed files with 16 additions and 54 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, final_jets
export read_final_state_particles, final_jets

# Jet visualisation is an extension, see ext/JetVisualisation.jl
function jetsplot() end
Expand Down
2 changes: 1 addition & 1 deletion src/Pseudojet.jl
Original file line number Diff line number Diff line change
Expand Up @@ -111,7 +111,7 @@ Print a `PseudoJet` object to the specified IO stream.
show(io::IO, jet::PseudoJet) = begin
print(io, "Pseudojet(px: ", jet.px, " py: ", jet.py, " pz: ", jet.pz, " E: ", jet.E,
"; ",
"pt: ", sqrt(jet._pt2), " eta: ", rapidity(jet), " phi: ", phi(jet), ", m: ",
"pt: ", sqrt(jet._pt2), " rapidity: ", rapidity(jet), " phi: ", phi(jet), ", m: ",
m(jet), ")")
end

Expand Down
63 changes: 12 additions & 51 deletions src/Utils.jl
Original file line number Diff line number Diff line change
Expand Up @@ -17,9 +17,11 @@ Reads final state particles from a file and returns them as a vector of type T.
# Returns
A vector of vectors of T objects, where each inner vector represents all
the particles of a particular event.
the particles of a particular event. In particular T can be `PseudoJet` or
a `LorentzVector` type. Note, if T is not `PseudoJet`, the order of the
arguments in the constructor must be `(t, x, y, z)`.
"""
function read_final_state_particles(fname; maxevents = -1, skipevents = 0, T=PseudoJet)
function read_final_state_particles(fname; maxevents = -1, skipevents = 0, T = PseudoJet)
f = open(fname)
if endswith(fname, ".gz")
@debug "Reading gzipped file $fname"
Expand All @@ -33,8 +35,14 @@ function read_final_state_particles(fname; maxevents = -1, skipevents = 0, T=Pse
input_particles = T[]
for p in parts
if p.status == 1
push!(input_particles,
T(p.momentum.x, p.momentum.y, p.momentum.z, p.momentum.t))
# Annoyingly PseudoJet and LorentzVector constructors
# disagree on the order of arguments...
if T == PseudoJet
particle = T(p.momentum.x, p.momentum.y, p.momentum.z, p.momentum.t)
else
particle = T(p.momentum.t, p.momentum.x, p.momentum.y, p.momentum.z)
end
push!(input_particles, particle)
end
end
push!(events, input_particles)
Expand All @@ -47,53 +55,6 @@ function read_final_state_particles(fname; maxevents = -1, skipevents = 0, T=Pse
events
end

"""
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 HepMC3 ASCII file to read particles from. If the file
is gzipped, the function will automatically decompress it.
- `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)
if endswith(fname, ".gz")
@debug "Reading gzipped file $fname"
f = GzipDecompressorStream(f)
end

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
close(f)

@info "Total Events: $(length(events))"
@debug events
events
end


"""
final_jets(jets::Vector{PseudoJet}, ptmin::AbstractFloat=0.0)
Expand Down
3 changes: 2 additions & 1 deletion test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -204,7 +204,8 @@ function do_test_compare_types(strategy::RecoStrategy.Strategy;
end

# From LorentzVector
events_lv::Vector{Vector{LorentzVector}} = read_final_state_particles_lv(events_file)
events_lv::Vector{Vector{LorentzVector}} = read_final_state_particles(events_file;
T = LorentzVector)
jet_collection_lv = FinalJets[]
for (ievt, event) in enumerate(events_lv)
finaljets = final_jets(inclusive_jets(jet_reconstruction(event, R = distance,
Expand Down

0 comments on commit 60e3bd7

Please sign in to comment.