diff --git a/README.md b/README.md index 83f239f1..a9d21065 100644 --- a/README.md +++ b/README.md @@ -33,7 +33,7 @@ cs = jet_reconstruct(particles::Vector{T}; algorithm = JetAlgorithm.AntiKt, R = - `particles` - a one dimensional array (vector) of input particles for the clustering - Any type that supplies the methods `pt2()`, `phi()`, `rapidity()`, `px()`, `py()`, `pz()`, `energy()` can be used - These methods have to be defined in the namespace of this package, i.e., `JetReconstruction.pt2(::T)` - - The `PseudoJet` or `EEjet` types from this package, a 4-vector from `LorentzVectorHEP`, or a `ReconstructedParticle` from [EDM4hep](https://github.com/peremato/EDM4hep.jl) are suitable (and have the appropriate definitions) + - The `PseudoJet` or `EEJet` types from this package, a 4-vector from `LorentzVectorHEP`, or a `ReconstructedParticle` from [EDM4hep](https://github.com/peremato/EDM4hep.jl) are suitable (and have the appropriate definitions) - `algorithm` is the name of the jet algorithm to be used (from the `JetAlgorithm` enum) - `JetAlgorithm.AntiKt` anti-$`{k}_\text{T}`$ clustering (default) - `JetAlgorithm.CA` Cambridge/Aachen clustering @@ -64,7 +64,7 @@ To obtain the final inclusive jets, use the `inclusive_jets` method: final_jets = inclusive_jets(cs::ClusterSequence; ptmin=0.0) ``` -Only jets passing the cut $p_T > p_{Tmin}$ will be returned. The result is returned as a `Vector{LorentzVectorHEP}`, but different return types can be specified (e.g., `T = EEjet`). +Only jets passing the cut $p_T > p_{Tmin}$ will be returned. The result is returned as a `Vector{LorentzVectorHEP}`, but different return types can be specified (e.g., `T = EEJet`). #### Sorting diff --git a/docs/src/particles.md b/docs/src/particles.md index 6dd1d294..c92a5b4a 100644 --- a/docs/src/particles.md +++ b/docs/src/particles.md @@ -11,7 +11,7 @@ to extract the 4-vector components, viz, the following are required: Currently built-in supported types are [`LorentzVectorHEP`](https://github.com/JuliaHEP/LorentzVectorHEP.jl), the -`PseudoJet` and `EEjet`s from this package, and `ReconstructedParticles` from +`PseudoJet` and `EEJet`s from this package, and `ReconstructedParticles` from [EDM4hep Inputs](@ref). If you require support for a different input collection type then ensure you diff --git a/examples/EDM4hep/EDM4hepJets.jl b/examples/EDM4hep/EDM4hepJets.jl index 439ad08f..4c508f5e 100644 --- a/examples/EDM4hep/EDM4hepJets.jl +++ b/examples/EDM4hep/EDM4hepJets.jl @@ -56,7 +56,7 @@ function main() if args[:printjets] @info begin jets = "Event $(ievt)\n" - for jet in exclusive_jets(cs; njets = args[:njets], T = EEjet) + for jet in exclusive_jets(cs; njets = args[:njets], T = EEJet) jets *= " $jet\n" end jets diff --git a/examples/EDM4hep/SimpleRecoEDM4hep.jl b/examples/EDM4hep/SimpleRecoEDM4hep.jl index 7f228446..abb6a456 100644 --- a/examples/EDM4hep/SimpleRecoEDM4hep.jl +++ b/examples/EDM4hep/SimpleRecoEDM4hep.jl @@ -14,7 +14,7 @@ recps = RootIO.get(reader, evt, "ReconstructedParticles") # Reconstruct and print the jets cs = jet_reconstruct(recps; algorithm = JetAlgorithm.Durham) -dijets = exclusive_jets(cs; njets = 2, T = EEjet) +dijets = exclusive_jets(cs; njets = 2, T = EEJet) for jet in dijets println(jet) end diff --git a/examples/instrumented-jetreco.jl b/examples/instrumented-jetreco.jl index b5c31f5e..110d4f7a 100644 --- a/examples/instrumented-jetreco.jl +++ b/examples/instrumented-jetreco.jl @@ -331,7 +331,7 @@ function main() # Try to read events into the correct type! if JetReconstruction.is_ee(args[:algorithm]) - jet_type = EEjet + jet_type = EEJet else jet_type = PseudoJet end diff --git a/examples/jetreco.jl b/examples/jetreco.jl index 82a6a115..37724a90 100644 --- a/examples/jetreco.jl +++ b/examples/jetreco.jl @@ -141,7 +141,7 @@ function main() global_logger(logger) # Try to read events into the correct type! if JetReconstruction.is_ee(args[:algorithm]) - jet_type = EEjet + jet_type = EEJet else jet_type = PseudoJet end diff --git a/ext/EDM4hepJets.jl b/ext/EDM4hepJets.jl index fb5d2bfa..a1788c2b 100644 --- a/ext/EDM4hepJets.jl +++ b/ext/EDM4hepJets.jl @@ -34,12 +34,12 @@ Return the energy component of a ReconstructedParticle's four vector. JetReconstruction.energy(recoparticle::ReconstructedParticle) = recoparticle.energy """ - JetReconstruction.EEjet(recoparticle::ReconstructedParticle) + JetReconstruction.EEJet(recoparticle::ReconstructedParticle) -Construct an EEjet from a ReconstructedParticle. +Construct an EEJet from a ReconstructedParticle. """ -function JetReconstruction.EEjet(recoparticle::ReconstructedParticle) - EEjet(JetReconstruction.px(recoparticle), JetReconstruction.py(recoparticle), +function JetReconstruction.EEJet(recoparticle::ReconstructedParticle) + EEJet(JetReconstruction.px(recoparticle), JetReconstruction.py(recoparticle), JetReconstruction.pz(recoparticle), JetReconstruction.energy(recoparticle)) end diff --git a/src/ClusterSequence.jl b/src/ClusterSequence.jl index 6c56de71..e593c50c 100644 --- a/src/ClusterSequence.jl +++ b/src/ClusterSequence.jl @@ -236,7 +236,7 @@ each parent jet. If the transverse momentum is greater than or equal to `ptmin`, the jet is added to the array of inclusive jets. Valid return types are `LorentzVectorCyl` and the jet type of the input `clusterseq` -(`U` - either `PseudoJet` or `EEjet` depending which algorithm was used) +(`U` - either `PseudoJet` or `EEJet` depending which algorithm was used) (N.B. this will evolve in the future to be any subtype of `FourMomentumBase`; currently unrecognised types will return `LorentzVectorCyl`). @@ -294,7 +294,7 @@ jets or a cut on the maximum distance parameter. - An array of `T` objects representing the exclusive jets. Valid return types are `LorentzVectorCyl` and the jet type of the input `clusterseq` -(`U` - either `PseudoJet` or `EEjet` depending which algorithm was used) +(`U` - either `PseudoJet` or `EEJet` depending which algorithm was used) (N.B. this will evolve in the future to be any subtype of `FourMomentumBase`; currently unrecognised types will return `LorentzVectorCyl`). diff --git a/src/EEAlgorithm.jl b/src/EEAlgorithm.jl index 7ba48f67..61c3bfa6 100644 --- a/src/EEAlgorithm.jl +++ b/src/EEAlgorithm.jl @@ -236,17 +236,17 @@ function ee_genkt_algorithm(particles::Vector{T}; p = 1, R = 4.0 end - if T == EEjet + if T == EEJet # recombination_particles will become part of the cluster sequence, so size it for # the starting particles and all N recombinations recombination_particles = copy(particles) sizehint!(recombination_particles, length(particles) * 2) else - recombination_particles = EEjet[] + recombination_particles = EEJet[] sizehint!(recombination_particles, length(particles) * 2) for i in eachindex(particles) push!(recombination_particles, - EEjet(px(particles[i]), py(particles[i]), pz(particles[i]), + EEJet(px(particles[i]), py(particles[i]), pz(particles[i]), energy(particles[i]))) end end @@ -258,13 +258,13 @@ function ee_genkt_algorithm(particles::Vector{T}; p = 1, end """ - _ee_genkt_algorithm(; particles::Vector{EEjet}, p = 1, R = 4.0, + _ee_genkt_algorithm(; particles::Vector{EEJet}, p = 1, R = 4.0, algorithm::JetAlgorithm.Algorithm = JetAlgorithm.Durham, recombine = +) This function is the actual implementation of the e+e- jet clustering algorithm. """ -function _ee_genkt_algorithm(; particles::Vector{EEjet}, p = 1, R = 4.0, +function _ee_genkt_algorithm(; particles::Vector{EEJet}, p = 1, R = 4.0, algorithm::JetAlgorithm.Algorithm = JetAlgorithm.Durham, recombine = +) # Bounds diff --git a/src/EEjet.jl b/src/EEJet.jl similarity index 62% rename from src/EEjet.jl rename to src/EEJet.jl index 7685502f..0d1d4f8d 100644 --- a/src/EEjet.jl +++ b/src/EEJet.jl @@ -1,7 +1,7 @@ """ - mutable struct EEjet + mutable struct EEJet -The `EEjet` struct is a 4-momentum object used for the e+e jet reconstruction routines. +The `EEJet` struct is a 4-momentum object used for the e+e jet reconstruction routines. # Fields - `px::Float64`: The x-component of the jet momentum. @@ -12,7 +12,7 @@ The `EEjet` struct is a 4-momentum object used for the e+e jet reconstruction ro - `_p2::Float64`: The squared momentum of the jet. - `_inv_p::Float64`: The inverse momentum of the jet. """ -mutable struct EEjet <: FourMomentum +mutable struct EEJet <: FourMomentum px::Float64 py::Float64 pz::Float64 @@ -22,30 +22,30 @@ mutable struct EEjet <: FourMomentum _cluster_hist_index::Int end -function EEjet(px::Real, py::Real, pz::Real, E::Real, _cluster_hist_index::Int) +function EEJet(px::Real, py::Real, pz::Real, E::Real, _cluster_hist_index::Int) @muladd p2 = px * px + py * py + pz * pz inv_p = @fastmath 1.0 / sqrt(p2) - EEjet(px, py, pz, E, p2, inv_p, _cluster_hist_index) + EEJet(px, py, pz, E, p2, inv_p, _cluster_hist_index) end -EEjet(px::Real, py::Real, pz::Real, E::Real) = EEjet(px, py, pz, E, 0) +EEJet(px::Real, py::Real, pz::Real, E::Real) = EEJet(px, py, pz, E, 0) -EEjet(pj::PseudoJet) = EEjet(px(pj), py(pj), pz(pj), energy(pj), cluster_hist_index(pj)) +EEJet(pj::PseudoJet) = EEJet(px(pj), py(pj), pz(pj), energy(pj), cluster_hist_index(pj)) -p2(eej::EEjet) = eej._p2 -pt2(eej::EEjet) = eej.px^2 + eej.py^2 +p2(eej::EEJet) = eej._p2 +pt2(eej::EEJet) = eej.px^2 + eej.py^2 const kt2 = pt2 -pt(eej::EEjet) = sqrt(pt2(eej)) -energy(eej::EEjet) = eej.E -px(eej::EEjet) = eej.px -py(eej::EEjet) = eej.py -pz(eej::EEjet) = eej.pz -nx(eej::EEjet) = eej.px * eej._inv_p -ny(eej::EEjet) = eej.py * eej._inv_p -nz(eej::EEjet) = eej.pz * eej._inv_p -cluster_hist_index(eej::EEjet) = eej._cluster_hist_index +pt(eej::EEJet) = sqrt(pt2(eej)) +energy(eej::EEJet) = eej.E +px(eej::EEJet) = eej.px +py(eej::EEJet) = eej.py +pz(eej::EEJet) = eej.pz +nx(eej::EEJet) = eej.px * eej._inv_p +ny(eej::EEJet) = eej.py * eej._inv_p +nz(eej::EEJet) = eej.pz * eej._inv_p +cluster_hist_index(eej::EEJet) = eej._cluster_hist_index -phi(eej::EEjet) = begin +phi(eej::EEJet) = begin phi = pt2(eej) == 0.0 ? 0.0 : atan(eej.py, eej.px) if phi < 0.0 phi += 2π @@ -53,10 +53,10 @@ phi(eej::EEjet) = begin phi end -m2(eej::EEjet) = energy(eej)^2 - p2(eej) -mass(eej::EEjet) = m2(eej) < 0.0 ? -sqrt(-m2(eej)) : sqrt(m2(eej)) +m2(eej::EEJet) = energy(eej)^2 - p2(eej) +mass(eej::EEJet) = m2(eej) < 0.0 ? -sqrt(-m2(eej)) : sqrt(m2(eej)) -function rapidity(eej::EEjet) +function rapidity(eej::EEJet) if energy(eej) == abs(pz(eej)) && iszero(pt2(eej)) MaxRapHere = _MaxRap + abs(pz(eej)) rap = (pz(eej) >= 0.0) ? MaxRapHere : -MaxRapHere @@ -75,13 +75,13 @@ function rapidity(eej::EEjet) end import Base.+; -function +(jet1::EEjet, jet2::EEjet) - EEjet(jet1.px + jet2.px, jet1.py + jet2.py, jet1.pz + jet2.pz, jet1.E + jet2.E) +function +(jet1::EEJet, jet2::EEJet) + EEJet(jet1.px + jet2.px, jet1.py + jet2.py, jet1.pz + jet2.pz, jet1.E + jet2.E) end import Base.show -function show(io::IO, eej::EEjet) - print(io, "EEjet(px: ", eej.px, " py: ", eej.py, " pz: ", eej.pz, " E: ", eej.E, +function show(io::IO, eej::EEJet) + print(io, "EEJet(px: ", eej.px, " py: ", eej.py, " pz: ", eej.pz, " E: ", eej.E, " cluster_hist_index: ", eej._cluster_hist_index, ")") end diff --git a/src/JetReconstruction.jl b/src/JetReconstruction.jl index 8d5a18fd..97f1e00a 100644 --- a/src/JetReconstruction.jl +++ b/src/JetReconstruction.jl @@ -36,11 +36,11 @@ py(p::LorentzVectorCyl) = LorentzVectorHEP.py(p) pz(p::LorentzVectorCyl) = LorentzVectorHEP.pz(p) energy(p::LorentzVectorCyl) = LorentzVectorHEP.energy(p) -# Pseudojet and EEjet types +# Pseudojet and EEJet types include("Pseudojet.jl") -include("EEjet.jl") +include("EEJet.jl") include("JetUtils.jl") -export PseudoJet, EEjet, pt_fraction, kt_scale, lorentzvector, lorentzvector_cyl +export PseudoJet, EEJet, pt_fraction, kt_scale, lorentzvector, lorentzvector_cyl # Jet reconstruction strategies and algorithms (enums!) include("AlgorithmStrategyEnums.jl") diff --git a/test/runtests.jl b/test/runtests.jl index 8d8c2038..9d2e9c60 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,6 +6,9 @@ function main() # Algorithm/power consistency checks include("test-algpower-consistency.jl") + # Basic tests for the Jet types + include("test-jet-types.jl") + # Jet utilities tests include("test-jet-utils.jl") diff --git a/test/test-jet-types.jl b/test/test-jet-types.jl new file mode 100644 index 00000000..a6a619f0 --- /dev/null +++ b/test/test-jet-types.jl @@ -0,0 +1,38 @@ +# A few very basic tests for the Jet types. + +include("common.jl") + +pj = PseudoJet(1.0, 2.0, 3.0, 10.0) +pj._cluster_hist_index = 7 +@testset "PseudoJet tests" begin + @test JetReconstruction.px(pj) ≈ 1.0 + @test JetReconstruction.py(pj) ≈ 2.0 + @test JetReconstruction.pz(pj) ≈ 3.0 + @test JetReconstruction.energy(pj) ≈ 10.0 + @test JetReconstruction.mag(pj) ≈ sqrt(1.0^2 + 2.0^2 + 3.0^2) + @test JetReconstruction.pt2(pj) ≈ 1.0^2 + 2.0^2 + @test JetReconstruction.phi(pj) ≈ atan(2.0, 1.0) + @test JetReconstruction.mass(pj) ≈ sqrt(10.0^2 - 1.0^2 - 2.0^2 - 3.0^2) + @test JetReconstruction.cluster_hist_index(pj) == 7 +end + +eej = EEJet(1.0, 2.0, 3.0, 10.0, 7) +@testset "EEJet tests" begin + @test JetReconstruction.px(eej) ≈ 1.0 + @test JetReconstruction.py(eej) ≈ 2.0 + @test JetReconstruction.pz(eej) ≈ 3.0 + @test JetReconstruction.energy(eej) ≈ 10.0 + @test JetReconstruction.p2(eej) ≈ 1.0^2 + 2.0^2 + 3.0^2 + @test JetReconstruction.pt2(eej) ≈ 1.0^2 + 2.0^2 + @test JetReconstruction.phi(eej) ≈ atan(2.0, 1.0) + @test JetReconstruction.mass(eej) ≈ sqrt(10.0^2 - 1.0^2 - 2.0^2 - 3.0^2) + @test JetReconstruction.cluster_hist_index(eej) == 7 +end + +eej_pseudojet = EEJet(pj) +@testset "Jet Interoperability" begin + @test JetReconstruction.px(eej_pseudojet) ≈ JetReconstruction.px(pj) + @test JetReconstruction.py(eej_pseudojet) ≈ JetReconstruction.py(pj) + @test JetReconstruction.pz(eej_pseudojet) ≈ JetReconstruction.pz(pj) + @test JetReconstruction.energy(eej_pseudojet) ≈ JetReconstruction.energy(pj) +end diff --git a/test/test-jet-utils.jl b/test/test-jet-utils.jl index bc99bc21..0bff3b9f 100644 --- a/test/test-jet-utils.jl +++ b/test/test-jet-utils.jl @@ -5,8 +5,8 @@ include("common.jl") pj1 = PseudoJet(1.0, 1.3, -2.3, 25.0) pj2 = PseudoJet(1.0, 0.3, -1.3, 17.0) -eej1 = EEjet(1.0, 1.3, -2.3, 25.0) -eej2 = EEjet(-1.0, 3.2, -1.2, 39.0) +eej1 = EEJet(1.0, 1.3, -2.3, 25.0) +eej2 = EEJet(-1.0, 3.2, -1.2, 39.0) @testset "Common jet utilities" begin @test JetReconstruction.pt_fraction(pj1, pj2) ≈ 0.3889609897118418 diff --git a/test/test-jet_reconstruct-interface.jl b/test/test-jet_reconstruct-interface.jl index b844287a..ca41b6bc 100644 --- a/test/test-jet_reconstruct-interface.jl +++ b/test/test-jet_reconstruct-interface.jl @@ -6,9 +6,9 @@ let inputs = JetReconstruction.read_final_state_particles(events_file_ee) @testset "jet_reconstruct() interface" begin # EE Algorithms @test typeof(jet_reconstruct(inputs[1]; algorithm = JetAlgorithm.Durham)) == - ClusterSequence{EEjet} + ClusterSequence{EEJet} @test typeof(jet_reconstruct(inputs[1]; algorithm = JetAlgorithm.EEKt, p = -1, - R = 1.0)) == ClusterSequence{EEjet} + R = 1.0)) == ClusterSequence{EEJet} @test_throws ArgumentError jet_reconstruct(inputs[1]; algorithm = JetAlgorithm.EEKt)