Skip to content

Commit

Permalink
Merge pull request #21 from graeme-a-stewart/graeme-chep
Browse files Browse the repository at this point in the history
Add CHEP2023 wrapper and HepMC3 reading
  • Loading branch information
gojakuch authored May 20, 2023
2 parents 91c0fa3 + 6b94212 commit 96530de
Show file tree
Hide file tree
Showing 10 changed files with 552 additions and 14 deletions.
4 changes: 4 additions & 0 deletions .gitignore
Original file line number Diff line number Diff line change
@@ -1,3 +1,7 @@
/Manifest.toml
/jetsavetest*.dat
/demo.jl
.vscode
*.json
/notebooks/*
/statprof/*
11 changes: 10 additions & 1 deletion Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,17 @@ authors = ["Atell Krasnopolski <[email protected]>"]
version = "0.1.0"

[deps]
IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
FlameGraphs = "08572546-2f56-4bcf-ba4e-bab62c3a3f89"
HepMC3_jll = "b85c3e40-22db-5268-bacb-02bd65cb4e01"
JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1"
JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899"
LorentzVectorHEP = "f612022c-142a-473f-8cfd-a09cf3793c6c"
Profile = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"
ProfileSVG = "132c30aa-f267-4189-9183-c8a63c7e05e6"
StatProfilerHTML = "a8a75453-ed82-57c9-9e16-4cd1196ecbf5"
StructTypes = "856f2bd8-1eba-4b0a-8007-ebc267875bd4"

[compat]
julia = "1.6"
Expand Down
241 changes: 241 additions & 0 deletions chep.jl
Original file line number Diff line number Diff line change
@@ -0,0 +1,241 @@
#! /usr/bin/env julia
"""
Wrapper to run jet reco code feeding in the standard set of HepMC events that
are used for benchmarking jet reco.
"""

import FlameGraphs
import ProfileSVG

using ArgParse
using Profile
using Colors
using StatProfilerHTML

using JetReconstruction

const R = 0.4
const ptmin = 5.0

function read_events(fname; maxevents=-1, skipevents=0)
f = open(fname)

events = Vector{PseudoJet}[]

ipart = 1
HepMC3.read_events(f, maxevents=maxevents, skipevents=skipevents) do parts
input_particles = PseudoJet[]
for p in parts
if p.status == 1
push!(
input_particles,
PseudoJet(p.momentum.px, p.momentum.py, p.momentum.pz, p.momentum.e),
)
end
end
push!(events, input_particles)
ipart += 1
end

println("Total Events: ", length(events))
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

function final_jets(jets::Vector{Vector{Float64}}, ptmin::AbstractFloat)
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 profile_code(events, niters)
Profile.init(n=10^6, delay=0.00001)
profile_events(events) = begin
for evt in events
anti_kt_algo(evt, R=0.4)
end
end
profile_events(events[1:1])
@profile for i = 1:niters
profile_events(events)
end
statprofilehtml()
fcolor = FlameGraphs.FlameColors(
reverse(colormap("Blues", 15))[1:5],
colorant"slategray4",
colorant"gray95",
reverse(colormap("Reds", 15))[1:5],
reverse(sequential_palette(39, 10; s=38, b=2))[1:5],#brownish pallette
)
ProfileSVG.save(
fcolor,
joinpath("statprof", "profsvg.svg");
combine=true,
timeunit=:ms,
font="Arial, Helvetica, sans-serif"
)
println(
"Flame graph from ProfileSVG.jl at file://",
abspath("statprof/profsvg.svg"),
"\n",
"""
\tRed tint: Runtime dispatch
\tBrown/yellow tint: Garbage collection
\tBlue tint: OK
""",
)
end

function jet_process(
events::Vector{Vector{PseudoJet}};
nsamples::Integer=1,
gcoff::Bool=false,
profile::Bool=false,
dump::Union{String,Nothing}=nothing
)
println("Will process $(size(events)) events")

# First, convert all events into the Vector of Vectors that Atell's
# code likes
event_vector = pseudojets2vectors(events)

# If we are dumping the results, setup the JSON structure
if !isnothing(dump)
jet_collection = FinalJets[]
end

# Warmup code if we are doing a multi-sample timing run
if nsamples > 1 || profile
println("Doing initial warm-up run")
for event in event_vector
anti_kt_algo(event, R=0.4)
end
end

GC.gc()
gcoff && GC.enable(false)

if profile
profile_code(event_vector, nsamples)
return Nothing
end

# Now setup timers and run the loop
cummulative_time = 0.0
cummulative_time2 = 0.0
for irun = 1:nsamples
print("$(irun)/$(nsamples) ")
t_start = time_ns()
for (ievt, event) in enumerate(event_vector)
finaljets, _ = anti_kt_algo(event, R=0.4)
fj = final_jets(finaljets, ptmin)
if !isnothing(dump) && irun == 1
println("Event $(ievt)")
for (ijet, jet) in enumerate(fj)
println(" $(ijet) - $(jet)")
end
push!(jet_collection, FinalJets(ievt, fj))
end
end
t_stop = time_ns()
dt_μs = convert(Float64, t_stop - t_start) * 1.e-3
println(dt_μs)
cummulative_time += dt_μs
cummulative_time2 += dt_μs^2
end

gcoff && GC.enable(true)

mean = cummulative_time / nsamples
cummulative_time2 /= nsamples
if nsamples > 1
sigma = sqrt(nsamples / (nsamples - 1) * (cummulative_time2 - mean^2))
else
sigma = 0.0
end
mean /= length(events)
sigma /= length(events)
println("Processed $(length(events)) events $(nsamples) times")
println("Time per event $(mean) ± $(sigma) μs")

if !isnothing(dump)
open(dump, "w") do io
JSON3.pretty(io, jet_collection)
end
end
end

parse_command_line(args) = begin
s = ArgParseSettings(autofix_names=true)
@add_arg_table! s begin
"--maxevents", "-n"
help = "Maximum number of events to read. -1 to read all events from the file."
arg_type = Int
default = -1

"--skip", "-s"
help = "Number of events to skip at beginning of the file."
arg_type = Int
default = 0

"--nsamples", "-m"
help = "Number of measurement points to acquire."
arg_type = Int
default = 1

"--gcoff"
help = "Turn off Julia garbage collector during each time measurement."
action = :store_true

"--profile"
help = "Profile code and generate a flame graph."
action = :store_true

"--dump"
help = "Write list of recontructed jets to a JSON formatted file and also to stdout"

"file"
help = "HepMC3 event file in HepMC3 to read."
required = true
end
return parse_args(args, s; as_symbols=true)
end

main() = begin
args = parse_command_line(ARGS)
events::Vector{Vector{PseudoJet}} =
read_events(args[:file], maxevents=args[:maxevents], skipevents=args[:skip])
jet_process(events, nsamples=args[:nsamples], gcoff=args[:gcoff], profile=args[:profile],
dump=args[:dump])
nothing
end

if abspath(PROGRAM_FILE) == @__FILE__
main()
end
27 changes: 14 additions & 13 deletions src/Algo.jl
Original file line number Diff line number Diff line change
@@ -1,9 +1,8 @@
import IfElse

Base.@propagate_inbounds function _dist(i, j, _eta, _phi)
deta = _eta[i] - _eta[j]
dphi = abs(_phi[i] - _phi[j])
dphi = IfElse.ifelse(dphi > pi, 2pi - dphi, dphi)
dphi = ifelse(dphi > pi, 2pi - dphi, dphi)
muladd(deta, deta, dphi*dphi)
end

Expand Down Expand Up @@ -40,15 +39,16 @@ Base.@propagate_inbounds function _upd_nn_nocross!(i::Int, from::Int, to::Int, _
nn = i
@inbounds @simd for j in from:(i-1)
Δ2 = _dist(i, j, _eta, _phi)
f = Δ2 <= nndist
nn = IfElse.ifelse(f, j, nn)
nndist = IfElse.ifelse(f, Δ2, nndist)
if Δ2 <= nndist
nn = j
nndist = Δ2
end
end
@inbounds @simd for j in (i+1):to
Δ2 = _dist(i, j, _eta, _phi)
f = Δ2 <= nndist
nn = IfElse.ifelse(f, j, nn)
nndist = IfElse.ifelse(f, Δ2, nndist)
nn = ifelse(f, j, nn)
nndist = ifelse(f, Δ2, nndist)
end
_nndist[i] = nndist
_nn[i] = nn;
Expand All @@ -73,14 +73,14 @@ Base.@propagate_inbounds function _upd_nn_step!(i, j, k, N, Nn, _kt2, _eta, _phi
end

cond = Δ2 < _nndist[i]
_nndist[i], _nn[i] = IfElse.ifelse(cond, (Δ2,k), (_nndist[i],_nn[i]))
_nndist[i], _nn[i] = ifelse(cond, (Δ2,k), (_nndist[i],_nn[i]))
end

nnk == Nn && (_nn[k] = j);
#nothing
end

function sequential_jet_reconstruct(objects::AbstractArray{T}; p=-1, R=1., recombine=+) where T
function sequential_jet_reconstruct(objects::AbstractArray{T}; p=-1.0, R=1., recombine=+) where T
# bounds
N::Int = length(objects)

Expand All @@ -95,6 +95,7 @@ function sequential_jet_reconstruct(objects::AbstractArray{T}; p=-1, R=1., recom
# data
_objects = copy(objects)
_kt2 = (JetReconstruction.pt.(_objects) .^ 2) .^ _p
# _kt2 = 1.0 ./ (JetReconstruction.pt.(_objects) .^ 2) <- Demo code for antikt talks (i.e., _p = -1)
_phi = JetReconstruction.phi.(_objects)
_eta = JetReconstruction.eta.(_objects)
_nn = Vector(1:N) # nearest neighbours
Expand All @@ -118,7 +119,7 @@ function sequential_jet_reconstruct(objects::AbstractArray{T}; p=-1, R=1., recom
dij_min = _nndij[1]
@inbounds @simd for k in 2:N
cond = _nndij[k] < dij_min
dij_min, i = IfElse.ifelse(cond, (_nndij[k], k), (dij_min, i))
dij_min, i = ifelse(cond, (_nndij[k], k), (dij_min, i))
end

j::Int = _nn[i]
Expand Down Expand Up @@ -179,7 +180,7 @@ Runs the anti-kt jet reconstruction algorithm. `objects` can be any collection o
Returns:
`jets` - a vector of jets. Each jet is of the same type as elements in `objects`.
`sequences` - a vector of vectors of indeces in `objects`. For all `i`, `sequences[i]` gives a sequence of indeces of objects that have been combined into the i-th jet (`jets[i]`).
`sequences` - a vector of vectors of indices in `objects`. For all `i`, `sequences[i]` gives a sequence of indices of objects that have been combined into the i-th jet (`jets[i]`).
"""
function anti_kt_algo(objects; R=1., recombine=+)
sequential_jet_reconstruct(objects, p=-1, R=R, recombine=recombine)
Expand All @@ -192,7 +193,7 @@ Runs the kt jet reconstruction algorithm. `objects` can be any collection of *un
Returns:
`jets` - a vector of jets. Each jet is of the same type as elements in `objects`.
`sequences` - a vector of vectors of indeces in `objects`. For all `i`, `sequences[i]` gives a sequence of indeces of objects that have been combined into the i-th jet (`jets[i]`).
`sequences` - a vector of vectors of indices in `objects`. For all `i`, `sequences[i]` gives a sequence of indices of objects that have been combined into the i-th jet (`jets[i]`).
"""
function kt_algo(objects; R=1., recombine=+)
sequential_jet_reconstruct(objects, p=1, R=R, recombine=recombine)
Expand All @@ -205,7 +206,7 @@ Runs the Cambridge/Aachen jet reconstruction algorithm. `objects` can be any col
Returns:
`jets` - a vector of jets. Each jet is of the same type as elements in `objects`.
`sequences` - a vector of vectors of indeces in `objects`. For all `i`, `sequences[i]` gives a sequence of indeces of objects that have been combined into the i-th jet (`jets[i]`).
`sequences` - a vector of vectors of indices in `objects`. For all `i`, `sequences[i]` gives a sequence of indices of objects that have been combined into the i-th jet (`jets[i]`).
"""
function cambridge_aachen_algo(objects; R=1., recombine=+)
sequential_jet_reconstruct(objects, p=0, R=R, recombine=recombine)
Expand Down
Loading

0 comments on commit 96530de

Please sign in to comment.