From 47182d9ac470c81c549c3eb20e4b7c3fd99ce6ed Mon Sep 17 00:00:00 2001 From: Graeme A Stewart Date: Mon, 30 Oct 2023 14:40:45 +0100 Subject: [PATCH] Test all algorithms 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). --- test/runtests.jl | 116 +++++++++++++++++++++++++++++++---------------- 1 file changed, 76 insertions(+), 40 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index ef32414..ae2b9f0 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -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 @@ -15,13 +16,13 @@ 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() @@ -29,62 +30,57 @@ function main() 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 @@ -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