fix GenKt not recognized as pp algorithm

src/AlgorithmStrategyEnums.jl

@@ -131,7 +131,12 @@ Check if the algorithm is a pp reconstruction algorithm.
 `true` if the algorithm is a pp reconstruction algorithm, `false` otherwise.
 """
 function is_pp(algorithm::JetAlgorithm.Algorithm)
-    return algorithm in [JetAlgorithm.AntiKt, JetAlgorithm.CA, JetAlgorithm.Kt]
+    return algorithm in [
+        JetAlgorithm.AntiKt,
+        JetAlgorithm.CA,
+        JetAlgorithm.Kt,
+        JetAlgorithm.GenKt
+    ]
 end
 
 """

src/ClusterSequence.jl

@@ -321,8 +321,12 @@
     end
 
     # Check that an algorithm was used that makes sense for exclusive jets
-    if !(clusterseq.algorithm ∈
-         (JetAlgorithm.CA, JetAlgorithm.Kt, JetAlgorithm.EEKt, JetAlgorithm.Durham))
+    if (clusterseq.algorithm ∈ (JetAlgorithm.GenKt, JetAlgorithm.EEKt)) &&
+       clusterseq.power < 0
+        throw(ArgumentError("Algorithm $(clusterseq.algorithm) requires power >= 0 for exclusive jets (power=$(clusterseq.power))"))
+    elseif clusterseq.algorithm ∉
+           (JetAlgorithm.CA, JetAlgorithm.Kt, JetAlgorithm.Durham, JetAlgorithm.GenKt,
+            JetAlgorithm.EEKt)
         throw(ArgumentError("Algorithm used is not suitable for exclusive jets ($(clusterseq.algorithm))"))
     end
 
@@ -378,8 +382,12 @@
 """
 function n_exclusive_jets(clusterseq::ClusterSequence; dcut::AbstractFloat)
     # Check that an algorithm was used that makes sense for exclusive jets
-    if !(clusterseq.algorithm ∈
-         (JetAlgorithm.CA, JetAlgorithm.Kt, JetAlgorithm.EEKt, JetAlgorithm.Durham))
+    if (clusterseq.algorithm ∈ (JetAlgorithm.GenKt, JetAlgorithm.EEKt)) &&
+       clusterseq.power < 0
+        throw(ArgumentError("Algorithm $(clusterseq.algorithm) requires power >= 0 for exclusive jets(power=$(clusterseq.power))"))
+    elseif clusterseq.algorithm ∉
+           (JetAlgorithm.CA, JetAlgorithm.Kt, JetAlgorithm.Durham, JetAlgorithm.GenKt,
+            JetAlgorithm.EEKt)
         throw(ArgumentError("Algorithm used is not suitable for exclusive jets ($(clusterseq.algorithm))"))
     end
 

test/runtests.jl

@@ -35,6 +35,9 @@ function main()
     include("test-pp-reconstruction.jl")
     include("test-ee-reconstruction.jl")
 
+    # Test jets selection
+    include("test-selection.jl")
+
     # Compare inputting data in PseudoJet with using a LorentzVector
     do_test_compare_types(RecoStrategy.N2Plain, algname = pp_algorithms[-1], power = -1)
     do_test_compare_types(RecoStrategy.N2Tiled, algname = pp_algorithms[-1], power = -1)
@@ -55,87 +58,6 @@ function main()
     include("test-aqua.jl")
 end
 
-"""
-Run the jet test
-
-dijmax -> apply inclusive dijmax cut
-njets -> apply inclusive njets cut
-
-If dijmax and njets are nothing, test inclusive jets with pt >= ptmin
-"""
-function do_test_compare_to_fastjet(strategy::RecoStrategy.Strategy, fastjet_jets;
-                                    algname = "Unknown",
-                                    selection = "Inclusive",
-                                    ptmin::Float64 = 5.0,
-                                    distance::Float64 = 0.4,
-                                    power::Real = -1,
-                                    dijmax = nothing,
-                                    njets = nothing)
-
-    # Strategy
-    if (strategy == RecoStrategy.N2Plain)
-        jet_reconstruction = plain_jet_reconstruct
-        strategy_name = "N2Plain"
-    elseif (strategy == RecoStrategy.N2Tiled)
-        jet_reconstruction = tiled_jet_reconstruct
-        strategy_name = "N2Tiled"
-    elseif (strategy == RecoStrategy.Best)
-        jet_reconstruction = jet_reconstruct
-        strategy_name = "Best"
-    else
-        throw(ErrorException("Strategy not yet implemented"))
-    end
-
-    # Now run our jet reconstruction...
-    # From PseudoJets
-    events::Vector{Vector{PseudoJet}} = read_final_state_particles(events_file_pp)
-    jet_collection = FinalJets[]
-    for (ievt, event) in enumerate(events)
-        # First run the reconstruction
-        if power == 1.5
-            cluster_seq = jet_reconstruction(event, R = distance,
-                                             algorithm = JetAlgorithm.GenKt, p = power)
-        else
-            cluster_seq = jet_reconstruction(event, R = distance, p = power)
-        end
-        # Now make the requested selection
-        if !isnothing(dijmax)
-            selected_jets = exclusive_jets(cluster_seq; dcut = dijmax)
-        elseif !isnothing(njets)
-            selected_jets = exclusive_jets(cluster_seq; njets = njets)
-        else
-            selected_jets = inclusive_jets(cluster_seq; ptmin = ptmin)
-        end
-        # And extract in out final_jets format
-        finaljets = final_jets(selected_jets)
-        sort_jets!(finaljets)
-        push!(jet_collection, FinalJets(ievt, finaljets))
-    end
-
-    @testset "Jet Reconstruction compare to FastJet: Strategy $strategy_name, Algorithm $algname, Selection $selection" begin
-        # Test each event in turn...
-        for (ievt, event) in enumerate(jet_collection)
-            @testset "Event $(ievt)" begin
-                @test size(event.jets) == size(fastjet_jets[ievt]["jets"])
-                # Test each jet in turn
-                for (ijet, jet) in enumerate(event.jets)
-                    if ijet <= size(fastjet_jets[ievt]["jets"])[1]
-                        # Approximate test - note that @test macro passes the 
-                        # tolerance into the isapprox() function
-                        # Use atol for position as this is absolute, but rtol for
-                        # 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
-                    end
-                end
-            end
-        end
-    end
-end
-
 function do_test_compare_types(strategy::RecoStrategy.Strategy;
                                algname = "Unknown",
                                ptmin::Float64 = 5.0,

test/test-pp-reconstruction.jl

@@ -3,15 +3,26 @@
 # Contains all common functions and necessary imports
 include("common.jl")
 
+const test_genkt_power = 1.5
+const test_cone_size = 0.4
+
 # Test inclusive jets for each algorithm and strategy
-for alg in [JetAlgorithm.AntiKt, JetAlgorithm.CA, JetAlgorithm.Kt],
+for alg in [JetAlgorithm.AntiKt, JetAlgorithm.CA, JetAlgorithm.Kt, JetAlgorithm.GenKt],
     stg in [RecoStrategy.N2Plain, RecoStrategy.N2Tiled]
 
-    fastjet_file = joinpath(@__DIR__, "data",
-                            "jet-collections-fastjet-inclusive-$(alg).json.gz")
+    if alg == JetAlgorithm.GenKt
+        power = test_genkt_power
+        fastjet_file = joinpath(@__DIR__, "data",
+                                "jet-collections-fastjet-inclusive-$(alg)-p$(power).json.gz")
+    else
+        power = JetReconstruction.algorithm2power[alg]
+        fastjet_file = joinpath(@__DIR__, "data",
+                                "jet-collections-fastjet-inclusive-$(alg).json.gz")
+    end
+
     test = ComparisonTest(events_file_pp, fastjet_file,
                           alg, stg,
-                          JetReconstruction.algorithm2power[alg], 0.4,
+                          power, test_cone_size,
                           (cs) -> inclusive_jets(cs; ptmin = 5.0), "inclusive")
     run_reco_test(test)
 end
@@ -22,7 +33,7 @@ for alg in [JetAlgorithm.CA, JetAlgorithm.Kt]
                             "jet-collections-fastjet-njets4-$(alg).json.gz")
     test = ComparisonTest(events_file_pp, fastjet_file,
                           alg, RecoStrategy.Best,
-                          JetReconstruction.algorithm2power[alg], 0.4,
+                          JetReconstruction.algorithm2power[alg], test_cone_size,
                           (cs) -> exclusive_jets(cs; njets = 4), "exclusive njets")
     run_reco_test(test)
 end
@@ -33,7 +44,7 @@ for (alg, dij_max) in zip([JetAlgorithm.CA, JetAlgorithm.Kt], ["0.99", "20.0"])
                             "jet-collections-fastjet-dij$(dij_max)-$(alg).json.gz")
     test = ComparisonTest(events_file_pp, fastjet_file,
                           alg, RecoStrategy.Best,
-                          JetReconstruction.algorithm2power[alg], 0.4,
+                          JetReconstruction.algorithm2power[alg], test_cone_size,
                           (cs) -> exclusive_jets(cs; dcut = tryparse(Float64, dij_max)),
                           "exclusive dcut")
     run_reco_test(test)

test/test-selection.jl

@@ -0,0 +1,46 @@
+# File: test/test_test-selection.jl
+
+include("common.jl")
+const test_cone_size = 0.4
+
+""" Helper function to test exclusive jets selection with given algorithm power argument."""
+function test_exclusive_jets_power(event, algo; p = nothing, should_throw::Bool)
+    test_cone_size = 0.4
+    test_strategy = RecoStrategy.Best
+    test_njets = 4
+
+    @testset "Exclusive jets selection: alg=$(algo) p=$(p)" begin
+        cluster_seq = JetReconstruction.jet_reconstruct(event; algorithm = algo, p = p,
+                                                        R = test_cone_size,
+                                                        strategy = test_strategy)
+        if should_throw
+            @test_throws ArgumentError begin
+                JetReconstruction.exclusive_jets(cluster_seq; njets = test_njets)
+            end
+        else
+            @test JetReconstruction.exclusive_jets(cluster_seq; njets = test_njets) isa Any
+        end
+    end
+end
+
+@testset "Exclusive jets selection algorithm power requirements" verbose=true begin
+    pp_event = first(JetReconstruction.read_final_state_particles(events_file_pp))
+    ee_event = first(JetReconstruction.read_final_state_particles(events_file_ee))
+
+    # Tests for pp_event
+    test_exclusive_jets_power(pp_event, JetAlgorithm.GenKt; p = -1.0,
+                              should_throw = true)
+    test_exclusive_jets_power(pp_event, JetAlgorithm.GenKt; p = 1.0,
+                              should_throw = false)
+    test_exclusive_jets_power(pp_event, JetAlgorithm.AntiKt; should_throw = true)
+    for alg in [JetAlgorithm.CA, JetAlgorithm.Kt]
+        test_exclusive_jets_power(pp_event, alg; should_throw = false)
+    end
+
+    # Tests for ee_event
+    test_exclusive_jets_power(ee_event, JetAlgorithm.EEKt; p = -1.0,
+                              should_throw = true)
+    test_exclusive_jets_power(ee_event, JetAlgorithm.EEKt; p = 1.0,
+                              should_throw = false)
+    test_exclusive_jets_power(ee_event, JetAlgorithm.Durham; should_throw = false)
+end