Improve parameter coverage in particle unit test

Project.toml

@@ -24,7 +24,8 @@ StaticArrays = "1.2.13"
 julia = "1.6"
 
 [extras]
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test"]
+test = ["SafeTestsets","Test"]

src/particles/particle_states.jl

@@ -228,8 +228,8 @@ end
 function _photon_state(pol::AllPolarization, mom::QEDbase.AbstractFourMomentum)
     cth = getCosTheta(mom)
     sth = sqrt(1 - cth^2)
+    cos_phi = getCosPhi(mom)
     sin_phi = getSinPhi(mom)
-    cos_phi = sqrt(1 - sin_phi^2)
     return SVector(
         SLorentzVector{Float64}(0.0, cth * cos_phi, cth * sin_phi, -sth),
         SLorentzVector{Float64}(0.0, -sin_phi, cos_phi, 0.0),
@@ -239,14 +239,14 @@ end
 function _photon_state(pol::PolarizationX, mom::QEDbase.AbstractFourMomentum)
     cth = getCosTheta(mom)
     sth = sqrt(1 - cth^2)
+    cos_phi = getCosPhi(mom)
     sin_phi = getSinPhi(mom)
-    cos_phi = sqrt(1 - sin_phi^2)
     return SLorentzVector{Float64}(0.0, cth * cos_phi, cth * sin_phi, -sth)
 end
 
 function _photon_state(pol::PolarizationY, mom::QEDbase.AbstractFourMomentum)
+    cos_phi = getCosPhi(mom)
     sin_phi = getSinPhi(mom)
-    cos_phi = sqrt(1 - sin_phi^2)
     return SLorentzVector{Float64}(0.0, -sin_phi, cos_phi, 0.0)
 end
 

test/dirac_tensor.jl

@@ -1,4 +1,4 @@
-
+using QEDbase
 using StaticArrays
 
 unary_methods = [-, +]

test/four_momentum.jl

@@ -1,3 +1,4 @@
+using QEDbase
 using Random
 
 const ATOL = 1e-15

test/gamma_matrices.jl

@@ -1,4 +1,4 @@
-
+using QEDbase
 using LinearAlgebra
 using Random
 using SparseArrays

test/lorentz_interface.jl

@@ -1,3 +1,4 @@
+using QEDbase
 
 lorentz_getter = [
     getT,

test/lorentz_vector.jl

@@ -1,3 +1,4 @@
+using QEDbase
 using StaticArrays
 
 unary_methods = [-, +]

test/particle_spinors.jl

@@ -1,3 +1,4 @@
+using QEDbase
 using Random
 
 const ATOL = 1e-15

test/particles.jl

@@ -1,15 +1,35 @@
 using QEDbase
+using StaticArrays
 using Random
 
+include("utils.jl")
+
 FERMION_STATES_GROUNDTRUTH_FACTORY = Dict(
     (Incoming, Electron) => IncomingFermionSpinor,
     (Outgoing, Electron) => OutgoingFermionSpinor,
     (Incoming, Positron) => IncomingAntiFermionSpinor,
     (Outgoing, Positron) => OutgoingAntiFermionSpinor,
 )
 
-rng = MersenneTwister(708583836976)
-x, y, z = rand(rng, 3)
+RNG = MersenneTwister(708583836976)
+ATOL = eps()
+RTOL = 0.0
+PHOTON_ENERGIES = (0.0, rand(RNG), rand(RNG) * 10)
+COS_THETAS = (-1.0, -rand(RNG), 0.0, rand(RNG), 1.0)
+# check every quadrant
+PHIS = (
+    0.0,
+    rand(RNG) * pi / 2,
+    pi / 2,
+    (1.0 + rand(RNG)) * pi / 2,
+    pi,
+    (2 + rand(RNG)) * pi / 2,
+    3 * pi / 2,
+    (3 + rand(RNG)) * pi / 2,
+    2 * pi,
+)
+
+X, Y, Z = rand(RNG, 3)
 
 @testset "fermion likes" begin
     @testset "fermion" begin
@@ -18,27 +38,27 @@ x, y, z = rand(rng, 3)
         @test is_particle(TestFermion())
         @test !is_anti_particle(TestFermion())
 
-        mom = SFourMomentum(sqrt(x^2 + y^2 + z^2 + mass(Electron())^2), x, y, z)
-        @testset "$P $D" for (P, D) in
+        @testset "$p $d" for (p, d) in
                              Iterators.product((Electron, Positron), (Incoming, Outgoing))
-            particle_mass = mass(P())
-            groundtruth_states = FERMION_STATES_GROUNDTRUTH_FACTORY[(D, P)](
+            mom = SFourMomentum(sqrt(mass(p()) + X^2 + Y^2 + Z^2), X, Y, Z)
+            particle_mass = mass(p())
+            groundtruth_states = FERMION_STATES_GROUNDTRUTH_FACTORY[(d, p)](
                 mom, particle_mass
             )
             groundtruth_tuple = SVector(groundtruth_states(1), groundtruth_states(2))
-            @test base_state(P(), D(), mom, AllSpin()) == groundtruth_tuple
-            @test base_state(P(), D(), mom, SpinUp()) == groundtruth_tuple[1]
-            @test base_state(P(), D(), mom, SpinDown()) == groundtruth_tuple[2]
+            @test base_state(p(), d(), mom, AllSpin()) == groundtruth_tuple
+            @test base_state(p(), d(), mom, SpinUp()) == groundtruth_tuple[1]
+            @test base_state(p(), d(), mom, SpinDown()) == groundtruth_tuple[2]
 
-            @test QEDbase._as_svec(base_state(P(), D(), mom, AllSpin())) isa SVector
-            @test QEDbase._as_svec(base_state(P(), D(), mom, SpinUp())) isa SVector
-            @test QEDbase._as_svec(base_state(P(), D(), mom, SpinDown())) isa SVector
+            @test QEDbase._as_svec(base_state(p(), d(), mom, AllSpin())) isa SVector
+            @test QEDbase._as_svec(base_state(p(), d(), mom, SpinUp())) isa SVector
+            @test QEDbase._as_svec(base_state(p(), d(), mom, SpinDown())) isa SVector
 
-            @test QEDbase._as_svec(base_state(P(), D(), mom, AllSpin())) ==
+            @test QEDbase._as_svec(base_state(p(), d(), mom, AllSpin())) ==
                 groundtruth_tuple
-            @test QEDbase._as_svec(base_state(P(), D(), mom, SpinUp()))[1] ==
+            @test QEDbase._as_svec(base_state(p(), d(), mom, SpinUp()))[1] ==
                 groundtruth_tuple[1]
-            @test QEDbase._as_svec(base_state(P(), D(), mom, SpinDown()))[1] ==
+            @test QEDbase._as_svec(base_state(p(), d(), mom, SpinDown()))[1] ==
                 groundtruth_tuple[2]
         end
     end
@@ -98,25 +118,36 @@ end
         @test is_particle(TestMajoranaBoson())
         @test is_anti_particle(TestMajoranaBoson())
     end
+end
 
-    @testset "photon" begin
-        @test !is_fermion(Photon())
-        @test is_boson(Photon())
-        @test is_particle(Photon())
-        @test is_anti_particle(Photon())
-        @test charge(Photon()) == 0.0
-        @test mass(Photon()) == 0.0
+@testset "photon" begin
+    @test !is_fermion(Photon())
+    @test is_boson(Photon())
+    @test is_particle(Photon())
+    @test is_anti_particle(Photon())
+    @test charge(Photon()) == 0.0
+    @test mass(Photon()) == 0.0
 
-        mom = SFourMomentum(sqrt(x^2 + y^2 + z^2 + mass(Photon())^2), x, y, z)
-        @testset "$D" for D in [Incoming, Outgoing]
+    @testset "$D" for D in [Incoming, Outgoing]
+        @testset "$om $cth $phi" for (om, cth, phi) in
+                                     Iterators.product(PHOTON_ENERGIES, COS_THETAS, PHIS)
+            #@testset "$x $y $z" for (x,y,z) in Iterators.product(X_arr,Y_arr,Z_arr)
+
+            mom = SFourMomentum(_cartesian_coordinates(om, om, cth, phi))
             both_photon_states = base_state(Photon(), D(), mom, AllPolarization())
 
             # property test the photon states
-            @test isapprox((both_photon_states[1] * mom), 0.0, atol=ATOL)
-            @test isapprox((both_photon_states[2] * mom), 0.0, atol=ATOL)
-            @test isapprox((both_photon_states[1] * both_photon_states[1]), -1.0)
-            @test isapprox((both_photon_states[2] * both_photon_states[2]), -1.0)
-            @test isapprox((both_photon_states[1] * both_photon_states[2]), 0.0, atol=ATOL)
+            @test isapprox((both_photon_states[1] * mom), 0.0, atol=ATOL, rtol=RTOL)
+            @test isapprox((both_photon_states[2] * mom), 0.0, atol=ATOL, rtol=RTOL)
+            @test isapprox(
+                (both_photon_states[1] * both_photon_states[1]), -1.0, atol=ATOL, rtol=RTOL
+            )
+            @test isapprox(
+                (both_photon_states[2] * both_photon_states[2]), -1.0, atol=ATOL, rtol=RTOL
+            )
+            @test isapprox(
+                (both_photon_states[1] * both_photon_states[2]), 0.0, atol=ATOL, rtol=RTOL
+            )
 
             # test the single polarization states
             @test base_state(Photon(), D(), mom, PolarizationX()) == both_photon_states[1]

test/runtests.jl

@@ -1,15 +1,32 @@
 using QEDbase
 using Test
+using SafeTestsets
 
-@testset "QEDbase.jl" begin
-    include("dirac_tensor.jl")
-    include("lorentz_vector.jl")
-    include("lorentz_interface.jl")
+begin
+    @time @safetestset "Dirac tensors" begin
+        include("dirac_tensor.jl")
+    end
 
-    include("gamma_matrices.jl")
+    @time @safetestset "Lorentz Vectors" begin
+        include("lorentz_vector.jl")
+    end
 
-    include("particle_spinors.jl")
-    include("four_momentum.jl")
+    @time @safetestset "Lorentz interface" begin
+        include("lorentz_interface.jl")
+    end
+    @time @safetestset "Gamma matrices" begin
+        include("gamma_matrices.jl")
+    end
 
-    include("particles.jl")
+    @time @safetestset "particle spinors" begin
+        include("particle_spinors.jl")
+    end
+
+    @time @safetestset "four momentum" begin
+        include("four_momentum.jl")
+    end
+
+    @time @safetestset "particles" begin
+        include("particles.jl")
+    end
 end

test/utils.jl

@@ -0,0 +1,11 @@
+"""
+
+Returns the cartesian coordinates (E,px,py,pz) for given spherical coordiantes
+(E, rho, cos_theta, phi), where rho denotes the length of the respective three-momentum, 
+theta is the polar- and phi the azimuthal angle. 
+"""
+function _cartesian_coordinates(E, rho, cth, phi)
+    sth = sqrt(1 - cth^2)
+    sphi, cphi = sincos(phi)
+    return (E, rho * sth * cphi, rho * sth * sphi, rho * cth)
+end