Skip to content

Commit

Permalink
Test all algorithms
Browse files Browse the repository at this point in the history 
Ensure that tests are run against all three standard jet
finding algorithms (Akt, CA, Ikt) for both strategies

Factorise out the PseudoJet vs. LorentzVector test.

Switch runs into "Warn" level mode (avoid info messages about
number of events processed).
  • Loading branch information
@graeme-a-stewart
graeme-a-stewart committed Nov 1, 2023
1 parent bee8ab5 commit 47182d9
Showing 1 changed file with 76 additions and 40 deletions.
116 changes: 76 additions & 40 deletions test/runtests.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,9 +4,10 @@ using JetReconstruction
using Test
using JSON
using LorentzVectorHEP
using Logging

"""Read JSON file with fastjet jets in it"""
function read_fastjet_outputs(;fname="test/data/jet_collections_fastjet_akt.json")
function read_fastjet_outputs(; fname = "test/data/jet_collections_fastjet_akt.json")
f = open(fname)
JSON.parse(f)
end
Expand All @@ -15,76 +16,71 @@ end
function sort_jets!(event_jet_array)
jet_pt(jet) = jet["pt"]
for e in event_jet_array
sort!(e["jets"], by=jet_pt, rev=true)
sort!(e["jets"], by = jet_pt, rev = true)
end
end

function sort_jets!(jet_array::Vector{FinalJet})
jet_pt(jet) = jet.pt
sort!(jet_array, by=jet_pt, rev=true)
sort!(jet_array, by = jet_pt, rev = true)
end

function main()
# Read our fastjet outputs (we read for anti-kt, cambridge/achen, inclusive-kt)
algorithms = Dict(-1 => "Anti-kt",
0 => "Cambridge/Achen",
1 => "Inclusive-kt")

fastjet_alg_files = Dict(-1 => "test/data/jet_collections_fastjet_akt.json",
0 => "test/data/jet_collections_fastjet_ca.json",
1 => "test/data/jet_collections_fastjet_ikt.json")

fastjet_data = Dict{Int, Vector{Any}}()
for (k, v) in pairs(fastjet_alg_files)
fastjet_jets = read_fastjet_outputs(fname=v)
fastjet_jets = read_fastjet_outputs(fname = v)
sort_jets!(fastjet_jets)
println(typeof(fastjet_jets))
fastjet_data[k] = fastjet_jets
end

# Test each stratgy...
do_jet_test(N2Plain, fastjet_data[-1])
do_jet_test(N2Tiled, fastjet_data[-1])
for power in keys(algorithms)
do_test_compare_to_fastjet(N2Plain, fastjet_data[power], algname = algorithms[power], power = power)
do_test_compare_to_fastjet(N2Tiled, fastjet_data[power], algname = algorithms[power], power = power)
end

# Compare inputing data in PseudoJet with using a LorentzVector
do_test_compare_types(N2Plain, algname = algorithms[-1], power = -1)
do_test_compare_types(N2Tiled, algname = algorithms[-1], power = -1)
end

function do_jet_test(strategy::JetRecoStrategy, fastjet_jets;
function do_test_compare_to_fastjet(strategy::JetRecoStrategy, fastjet_jets;
algname = "Unknown",
ptmin::Float64 = 5.0,
distance::Float64 = 0.4,
power::Integer = -1)
distance::Float64 = 0.4,
power::Integer = -1)

# Strategy
if (strategy == N2Plain)
jet_reconstruction = plain_jet_reconstruct
# Strategy
if (strategy == N2Plain)
jet_reconstruction = plain_jet_reconstruct
strategy_name = "N2Plain"
elseif (strategy == N2Tiled)
jet_reconstruction = tiled_jet_reconstruct
jet_reconstruction = tiled_jet_reconstruct
strategy_name = "N2Tiled"
else
throw(ErrorException("Strategy not yet implemented"))
end
else
throw(ErrorException("Strategy not yet implemented"))
end

# Now run our jet reconstruction...
# From PseudoJets
events::Vector{Vector{PseudoJet}} = read_final_state_particles("test/data/events.hepmc3")
jet_collection = FinalJets[]
for (ievt, event) in enumerate(events)
finaljets, _ = jet_reconstruction(event, R=distance, p=power, ptmin=ptmin)
finaljets, _ = jet_reconstruction(event, R = distance, p = power, ptmin = ptmin)
fj = final_jets(finaljets, ptmin)
sort_jets!(fj)
push!(jet_collection, FinalJets(ievt, fj))
end

# From LorentzVector
events_lv::Vector{Vector{LorentzVector}} = read_final_state_particles_lv("test/data/events.hepmc3")
jet_collection_lv = FinalJets[]
for (ievt, event) in enumerate(events_lv)
finaljets, _ = jet_reconstruction(event, R=distance, p=power, ptmin=ptmin)
fj = final_jets(finaljets, ptmin)
sort_jets!(fj)
push!(jet_collection_lv, FinalJets(ievt, fj))
end

@testset "Jet Reconstruction PseudoJet, $strategy_name" begin
@testset "Jet Reconstruction compare to FastJet: Strategy $strategy_name, Algorithm $algname" begin
# Test each event in turn...
for (ievt, event) in enumerate(jet_collection)
@testset "Event $(ievt)" begin
Expand All @@ -98,33 +94,73 @@ function do_jet_test(strategy::JetRecoStrategy, fastjet_jets;
# the momentum
# Sometimes phi could be in the range [-π, π], but FastJet always is [0, 2π]
normalised_phi = jet.phi < 0.0 ? jet.phi + 2π : jet.phi
@test jet.rap ≈ fastjet_jets[ievt]["jets"][ijet]["rap"] atol=1e-7
@test normalised_phi ≈ fastjet_jets[ievt]["jets"][ijet]["phi"] atol=1e-7
@test jet.pt ≈ fastjet_jets[ievt]["jets"][ijet]["pt"] rtol=1e-6
@test jet.rap ≈ fastjet_jets[ievt]["jets"][ijet]["rap"] atol = 1e-7
@test normalised_phi ≈ fastjet_jets[ievt]["jets"][ijet]["phi"] atol = 1e-7
@test jet.pt ≈ fastjet_jets[ievt]["jets"][ijet]["pt"] rtol = 1e-6
end
end
end
end
end
end

@testset "Jet Reconstruction LorentzVector, $strategy_name" begin
function do_test_compare_types(strategy::JetRecoStrategy;
algname = "Unknown",
ptmin::Float64 = 5.0,
distance::Float64 = 0.4,
power::Integer = -1)

# Strategy
if (strategy == N2Plain)
jet_reconstruction = plain_jet_reconstruct
strategy_name = "N2Plain"
elseif (strategy == N2Tiled)
jet_reconstruction = tiled_jet_reconstruct
strategy_name = "N2Tiled"
else
throw(ErrorException("Strategy not yet implemented"))
end

# Now run our jet reconstruction...
# From PseudoJets
events::Vector{Vector{PseudoJet}} = read_final_state_particles("test/data/events.hepmc3")
jet_collection = FinalJets[]
for (ievt, event) in enumerate(events)
finaljets, _ = jet_reconstruction(event, R = distance, p = power, ptmin = ptmin)
fj = final_jets(finaljets, ptmin)
sort_jets!(fj)
push!(jet_collection, FinalJets(ievt, fj))
end

# From LorentzVector
events_lv::Vector{Vector{LorentzVector}} = read_final_state_particles_lv("test/data/events.hepmc3")
jet_collection_lv = FinalJets[]
for (ievt, event) in enumerate(events_lv)
finaljets, _ = jet_reconstruction(event, R = distance, p = power, ptmin = ptmin)
fj = final_jets(finaljets, ptmin)
sort_jets!(fj)
push!(jet_collection_lv, FinalJets(ievt, fj))
end

@testset "Jet Reconstruction Compare PseudoJet and LorentzVector, Strategy $strategy_name, Algorithm $algname" begin
# Here we test that inputing LorentzVector gave the same results as PseudoJets
for (ievt, (event, event_lv)) in enumerate(zip(jet_collection, jet_collection_lv))
@testset "Event $(ievt)" begin
@test size(event.jets) == size(event_lv.jets)
# Test each jet in turn
for (jet, jet_lv) in zip(event.jets, event_lv.jets)
@test jet.rap ≈ jet_lv.rap atol=1e-7
@test jet.phi ≈ jet_lv.phi atol=1e-7
@test jet.pt ≈ jet_lv.pt rtol=1e-6
@test jet.rap ≈ jet_lv.rap atol = 1e-7
@test jet.phi ≈ jet_lv.phi atol = 1e-7
@test jet.pt ≈ jet_lv.pt rtol = 1e-6
end
end
end
end


end


if abspath(PROGRAM_FILE) == @__FILE__
main()
logger = ConsoleLogger(stdout, Logging.Warn)
global_logger(logger)
main()
end

0 comments on commit 47182d9

Please sign in to comment.