diff --git a/src/JetReconstruction.jl b/src/JetReconstruction.jl index e94ee70..f8a0a47 100644 --- a/src/JetReconstruction.jl +++ b/src/JetReconstruction.jl @@ -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 diff --git a/src/Pseudojet.jl b/src/Pseudojet.jl index 6e5bda7..e7c5124 100644 --- a/src/Pseudojet.jl +++ b/src/Pseudojet.jl @@ -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 diff --git a/src/Utils.jl b/src/Utils.jl index a29b3ea..daaa1ba 100644 --- a/src/Utils.jl +++ b/src/Utils.jl @@ -3,39 +3,46 @@ using CodecZlib """ - read_final_state_particles(fname; maxevents = -1, skipevents = 0) + read_final_state_particles(fname; maxevents = -1, skipevents = 0, T=PseudoJet) -Reads final state particles from a file and returns them as a vector of -PseudoJet objects. +Reads final state particles from a file and returns them as a vector of type T. # 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. +- `skipevents=0`: The number of events to skip before an event is included. +- `T=PseudoJet`: The type of object to contruct and return. # Returns -A vector of vectors of PseudoJet objects, where each inner vector represents all -the particles of a particular event. +A vector of vectors of T objects, where each inner vector represents all +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) +function read_final_state_particles(fname; maxevents = -1, skipevents = 0, T = PseudoJet) f = open(fname) if endswith(fname, ".gz") @debug "Reading gzipped file $fname" f = GzipDecompressorStream(f) end - events = Vector{PseudoJet}[] + events = Vector{T}[] ipart = 1 HepMC3.read_events(f, maxevents = maxevents, skipevents = skipevents) do parts - input_particles = PseudoJet[] + input_particles = T[] for p in parts if p.status == 1 - push!(input_particles, - PseudoJet(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) @@ -48,74 +55,6 @@ function read_final_state_particles(fname; maxevents = -1, skipevents = 0) 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 - -""" -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 - """ final_jets(jets::Vector{PseudoJet}, ptmin::AbstractFloat=0.0) diff --git a/test/runtests.jl b/test/runtests.jl index 17f77f9..1e6bf22 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -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,