From 48e42bf16509f240ba929c91d89cafe0bcf8c296 Mon Sep 17 00:00:00 2001
From: Anton Reinhard <anton.reinhard@proton.me>
Date: Fri, 3 May 2024 13:53:37 +0200
Subject: [PATCH 01/12] Add StaticArrays as dependency for the phase space
 points

---
 Project.toml | 1 +
 1 file changed, 1 insertion(+)

diff --git a/Project.toml b/Project.toml
index ebdec6e..e5cb4d2 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,6 +5,7 @@ version = "0.1.0"
 
 [deps]
 QEDbase = "10e22c08-3ccb-4172-bfcf-7d7aa3d04d93"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
 QEDbase = "0.1"

From 515798a8e034462355945db29f631aa42bb56fbb Mon Sep 17 00:00:00 2001
From: Anton Reinhard <anton.reinhard@proton.me>
Date: Fri, 3 May 2024 13:54:15 +0200
Subject: [PATCH 02/12] Add phase space point and stateful particle definitions
 and stateful particle tests

---
 src/QEDprocesses.jl  |  1 +
 src/phase_spaces.jl  | 72 ++++++++++++++++++++++++++++++++++++++++++++
 test/phase_spaces.jl | 33 ++++++++++++++++++++
 test/runtests.jl     |  3 ++
 4 files changed, 109 insertions(+)
 create mode 100644 test/phase_spaces.jl

diff --git a/src/QEDprocesses.jl b/src/QEDprocesses.jl
index 3bb38ef..15d8ad8 100644
--- a/src/QEDprocesses.jl
+++ b/src/QEDprocesses.jl
@@ -30,6 +30,7 @@ export propagator
 export AbstractCoordinateSystem, SphericalCoordinateSystem
 export AbstractFrameOfReference, CenterOfMomentumFrame, ElectronRestFrame
 export AbstractPhasespaceDefinition, PhasespaceDefinition
+export ParticleStateful, PhaseSpacePoint
 
 using QEDbase
 
diff --git a/src/phase_spaces.jl b/src/phase_spaces.jl
index d35fe36..7447bc9 100644
--- a/src/phase_spaces.jl
+++ b/src/phase_spaces.jl
@@ -7,6 +7,15 @@
 # TODO: ship to interfaces!
 #####################
 
+using StaticArrays
+
+import QEDbase.is_particle
+import QEDbase.is_anti_particle
+import QEDbase.is_fermion
+import QEDbase.is_boson
+import QEDbase.is_incoming
+import QEDbase.is_outgoing
+
 abstract type AbstractCoordinateSystem end
 struct SphericalCoordinateSystem <: AbstractCoordinateSystem end
 
@@ -83,3 +92,66 @@ end
         ),
     )
 end
+
+"""
+    ParticleStateful
+
+Representation of a particle with a state. It has four fields:
+- `mom::AbstractPhasespaceElement`: The momentum of the particle.
+- `type::AbstractParticle`: The type of the particle, [`Electron`](@ref)(), [`Positron`](@ref)() etc..
+- `spinorpol::AbstractSpinOrPolarization`: The spin or polarization of the particle, [`SpinUp`](@ref)(), [`PolX`](@ref)() etc.. Can only use spins with fermions and polarizations with bosons.
+- `dir::ParticleDirection`: The direction of the particle, [`Incoming`](@ref)() or [`Outgoing`](@ref)().
+
+"""
+struct ParticleStateful{ElType<:AbstractPhasespaceElement}
+    mom::ElType
+    type::QEDbase.AbstractParticle
+    spinorpol::AbstractSpinOrPolarization
+    dir::ParticleDirection
+
+    function ParticleStateful(
+        mom::ElType,
+        type::QEDbase.AbstractParticle,
+        spinorpol::AbstractSpinOrPolarization,
+        dir::ParticleDirection,
+    ) where {ElType<:AbstractPhasespaceElement}
+        if is_fermion(type) && !(spinorpol isa AbstractSpin)
+            throw(
+                InvalidInputError(
+                    "Tried to construct a stateful fermion with a $(spinorpol) instead of a spin",
+                ),
+            )
+        elseif is_boson(type) && !(spinorpol isa AbstractPolarization)
+            throw(
+                InvalidInputError(
+                    "Tried to construct a stateful boson with a $(spinorpol) instead of a polarization",
+                ),
+            )
+        end
+
+        return new{ElType}(mom, type, spinorpol, dir)
+    end
+end
+
+struct PhaseSpacePoint{
+    PROC<:AbstractProcessDefinition,
+    MODEL<:AbstractModelDefinition,
+    PSDEF<:AbstractPhasespaceDefinition,
+    PhaseSpaceElementType<:AbstractPhasespaceElement,
+    N_IN_PARTICLES,
+    N_OUT_PARTICLES,
+}
+    proc::PROC
+    model::MODEL
+    ps_def::PSDEF
+
+    in_particles::SVector{N_IN_PARTICLES,ParticleStateful{PhaseSpaceElementType}}
+    out_particles::SVector{N_OUT_PARTICLES,ParticleStateful{PhaseSpaceElementType}}
+end
+
+@inline is_incoming(particle::ParticleStateful) = is_incoming(particle.dir)
+@inline is_outgoing(particle::ParticleStateful) = is_outgoing(particle.dir)
+@inline is_fermion(particle::ParticleStateful) = is_fermion(particle.type)
+@inline is_boson(particle::ParticleStateful) = is_boson(particle.type)
+@inline is_particle(particle::ParticleStateful) = is_particle(particle.type)
+@inline is_anti_particle(particle::ParticleStateful) = is_anti_particle(particle.type)
diff --git a/test/phase_spaces.jl b/test/phase_spaces.jl
new file mode 100644
index 0000000..a8b4f6f
--- /dev/null
+++ b/test/phase_spaces.jl
@@ -0,0 +1,33 @@
+using Random
+using QEDbase
+using QEDprocesses
+
+# can be removed when QEDbase exports them
+import QEDbase.is_incoming, QEDbase.is_outgoing
+
+@testset "Phasespace Point" begin
+    DIRECTIONS = [Incoming(), Outgoing()]
+    PARTICLES = [Electron(), Positron()] #=, Muon(), AntiMuon(), Tauon(), AntiTauon()=#
+    SPINANDPOLS = [AllSpin(), SpinUp(), SpinDown(), AllPol(), PolX(), PolY()]
+
+    for (particle, dir, spinorpol) in Iterators.product(PARTICLES, DIRECTIONS, SPINANDPOLS)
+        momentum = rand(SFourMomentum)
+
+        if (is_fermion(particle) && (spinorpol isa AbstractSpin)) ||
+            (is_boson(particle) && (spinorpol isa AbstractPolarization))
+            particle_stateful = ParticleStateful(momentum, particle, spinorpol, dir)
+
+            @test particle_stateful.mom == momentum
+            @test is_fermion(particle_stateful) == is_fermion(particle)
+            @test is_boson(particle_stateful) == is_boson(particle)
+            @test is_particle(particle_stateful) == is_particle(particle)
+            @test is_anti_particle(particle_stateful) == is_anti_particle(particle)
+            @test is_incoming(particle_stateful) == is_incoming(dir)
+            @test is_outgoing(particle_stateful) == is_outgoing(dir)
+        else
+            @test_throws "Tried to construct a stateful $(is_fermion(particle) ? "fermion" : "boson") with a $(spinorpol) instead of a $(is_fermion(particle) ? "spin" : "polarization")" ParticleStateful(
+                momentum, particle, spinorpol, dir
+            )
+        end
+    end
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 470923a..9510013 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -21,4 +21,7 @@ begin
     @time @safetestset "cross section & probability" begin
         include("cross_sections.jl")
     end
+    @time @safetestset "phase spaces" begin
+        include("phase_spaces.jl")
+    end
 end

From 074e2b56095d02cfc88ce47158ea02915efdad29 Mon Sep 17 00:00:00 2001
From: Anton Reinhard <anton.reinhard@proton.me>
Date: Fri, 3 May 2024 16:31:30 +0200
Subject: [PATCH 03/12] Add PhaseSpacePoint definition and first test

---
 test/Project.toml    |  1 +
 test/phase_spaces.jl | 29 ++++++++++++++++++++++++++++-
 2 files changed, 29 insertions(+), 1 deletion(-)

diff --git a/test/Project.toml b/test/Project.toml
index c71dc1e..6ab0525 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -2,5 +2,6 @@
 QEDbase = "10e22c08-3ccb-4172-bfcf-7d7aa3d04d93"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/phase_spaces.jl b/test/phase_spaces.jl
index a8b4f6f..47a7f2d 100644
--- a/test/phase_spaces.jl
+++ b/test/phase_spaces.jl
@@ -1,11 +1,16 @@
 using Random
+using StaticArrays
 using QEDbase
 using QEDprocesses
 
 # can be removed when QEDbase exports them
 import QEDbase.is_incoming, QEDbase.is_outgoing
 
-@testset "Phasespace Point" begin
+include("test_implementation/TestImplementation.jl")
+TESTMODEL = TestImplementation.TestModel()
+TESTPSDEF = TestImplementation.TestPhasespaceDef()
+
+@testset "Stateful Particle" begin
     DIRECTIONS = [Incoming(), Outgoing()]
     PARTICLES = [Electron(), Positron()] #=, Muon(), AntiMuon(), Tauon(), AntiTauon()=#
     SPINANDPOLS = [AllSpin(), SpinUp(), SpinDown(), AllPol(), PolX(), PolY()]
@@ -31,3 +36,25 @@ import QEDbase.is_incoming, QEDbase.is_outgoing
         end
     end
 end
+
+@testset "Phasespace Point" begin
+    in_el = ParticleStateful(rand(SFourMomentum), Electron(), AllSpin(), Incoming())
+    in_ph = ParticleStateful(rand(SFourMomentum), Photon(), AllPol(), Incoming())
+    out_el = ParticleStateful(rand(SFourMomentum), Electron(), AllSpin(), Outgoing())
+    out_ph = ParticleStateful(rand(SFourMomentum), Photon(), AllPol(), Outgoing())
+
+    in_particles_valid = SVector(in_el, in_ph)
+    in_particles_invalid = SVector(in_el, out_ph)
+
+    out_particles_valid = SVector(out_el, out_ph)
+    out_particles_invalid = SVector(out_el, in_ph)
+
+    model = TESTMODEL
+    process = TestImplementation.TestProcess(
+        SVector{2,AbstractParticle}(Electron(), Photon()),
+        SVector{2,AbstractParticle}(Electron(), Photon()),
+    )
+    phasespace_def = TESTPSDEF
+
+    PhaseSpacePoint(process, model, phasespace_def, in_particles_valid, out_particles_valid)
+end

From e8b30181d94a8baf95af059be2f44b828993fa66 Mon Sep 17 00:00:00 2001
From: Anton Reinhard <anton.reinhard@proton.me>
Date: Fri, 3 May 2024 17:20:51 +0200
Subject: [PATCH 04/12] Add phase space point validation and tests

---
 src/phase_spaces.jl  | 27 +++++++++++++++++++++++++++
 test/phase_spaces.jl | 16 ++++++++++++++++
 2 files changed, 43 insertions(+)

diff --git a/src/phase_spaces.jl b/src/phase_spaces.jl
index 7447bc9..28046a2 100644
--- a/src/phase_spaces.jl
+++ b/src/phase_spaces.jl
@@ -147,6 +147,33 @@ struct PhaseSpacePoint{
 
     in_particles::SVector{N_IN_PARTICLES,ParticleStateful{PhaseSpaceElementType}}
     out_particles::SVector{N_OUT_PARTICLES,ParticleStateful{PhaseSpaceElementType}}
+
+    function PhaseSpacePoint(
+        proc::PROC, model::MODEL, ps_def::PSDEF, in_p::IN_P, out_p::OUT_P
+    ) where {
+        PROC<:AbstractProcessDefinition,
+        MODEL<:AbstractModelDefinition,
+        PSDEF<:AbstractPhasespaceDefinition,
+        PhaseSpaceElementType<:AbstractPhasespaceElement,
+        IN_P<:AbstractVector{ParticleStateful{PhaseSpaceElementType}},
+        OUT_P<:AbstractVector{ParticleStateful{PhaseSpaceElementType}},
+    }
+        @assert length(incoming_particles(proc)) == length(in_p)
+        @assert length(outgoing_particles(proc)) == length(out_p)
+
+        for (proc_p, p) in zip(incoming_particles(proc), in_p)
+            @assert proc_p == p.type "Process given particle type ($(proc_p)) does not match stateful particle type ($(p.type))"
+            @assert is_incoming(p) "Stateful particle $(p) is given as an incoming particle but is outgoing"
+        end
+        for (proc_p, p) in zip(outgoing_particles(proc), out_p)
+            @assert proc_p == p.type "Process given particle type ($(proc_p)) does not match stateful particle type ($(p.type))"
+            @assert is_outgoing(p) "Stateful particle $(p) is given as an outgoing particle but is incoming"
+        end
+
+        return new{PROC,MODEL,PSDEF,PhaseSpaceElementType,length(in_p),length(out_p)}(
+            proc, model, ps_def, in_p, out_p
+        )
+    end
 end
 
 @inline is_incoming(particle::ParticleStateful) = is_incoming(particle.dir)
diff --git a/test/phase_spaces.jl b/test/phase_spaces.jl
index 47a7f2d..eaf9f43 100644
--- a/test/phase_spaces.jl
+++ b/test/phase_spaces.jl
@@ -57,4 +57,20 @@ end
     phasespace_def = TESTPSDEF
 
     PhaseSpacePoint(process, model, phasespace_def, in_particles_valid, out_particles_valid)
+
+    @test_throws r"Stateful particle (.*) is given as an incoming particle but is outgoing" PhaseSpacePoint(
+        process, model, phasespace_def, in_particles_invalid, out_particles_valid
+    )
+
+    @test_throws r"Stateful particle (.*) is given as an outgoing particle but is incoming" PhaseSpacePoint(
+        process, model, phasespace_def, in_particles_valid, out_particles_invalid
+    )
+
+    @test_throws r"Process given particle type \((.*)Electron\(\)\) does not match stateful particle type \((.*)Photon\(\)\)" PhaseSpacePoint(
+        process, model, phasespace_def, SVector(in_ph, in_el), out_particles_valid
+    )
+
+    @test_throws r"Process given particle type \((.*)Electron\(\)\) does not match stateful particle type \((.*)Photon\(\)\)" PhaseSpacePoint(
+        process, model, phasespace_def, in_particles_valid, SVector(out_ph, out_el)
+    )
 end

From c84f1aed89d340783284c9d7716ae80174b34fde Mon Sep 17 00:00:00 2001
From: Anton Reinhard <anton.reinhard@proton.me>
Date: Fri, 3 May 2024 18:31:59 +0200
Subject: [PATCH 05/12] Add some more documentation

---
 src/phase_spaces.jl | 52 +++++++++++++++++++++++++++++++++++++--------
 1 file changed, 43 insertions(+), 9 deletions(-)

diff --git a/src/phase_spaces.jl b/src/phase_spaces.jl
index 28046a2..fd31563 100644
--- a/src/phase_spaces.jl
+++ b/src/phase_spaces.jl
@@ -98,10 +98,13 @@ end
 
 Representation of a particle with a state. It has four fields:
 - `mom::AbstractPhasespaceElement`: The momentum of the particle.
-- `type::AbstractParticle`: The type of the particle, [`Electron`](@ref)(), [`Positron`](@ref)() etc..
-- `spinorpol::AbstractSpinOrPolarization`: The spin or polarization of the particle, [`SpinUp`](@ref)(), [`PolX`](@ref)() etc.. Can only use spins with fermions and polarizations with bosons.
-- `dir::ParticleDirection`: The direction of the particle, [`Incoming`](@ref)() or [`Outgoing`](@ref)().
+- `type::AbstractParticle`: The type of the particle, `QEDbase.Electron()`, `QEDbase.Positron()` etc..
+- `spinorpol::AbstractSpinOrPolarization`: The spin or polarization of the particle, `QEDbase.SpinUp()`, `QEDbase.PolX() etc.. Can only use spins with fermions and polarizations with bosons.
+- `dir::ParticleDirection`: The direction of the particle, `QEDbase.Incoming()` or `QEDbase.Outgoing()`.
 
+Overloads for `QEDbase.is_fermion`, `QEDbase.is_boson`, `QEDbase.is_particle`, `QEDbase.is_anti_particle`, `QEDbase.is_incoming`, and `QEDbase.is_outgoing` are provided, delegating the call to the correct field.
+
+The legality of the combination of `type` and `spinorpol` is checked on construction. If, for example, the construction of an `Electron()` with a polarization is attempted, an [`InvalidInputError`](@ref) is thrown.
 """
 struct ParticleStateful{ElType<:AbstractPhasespaceElement}
     mom::ElType
@@ -133,6 +136,13 @@ struct ParticleStateful{ElType<:AbstractPhasespaceElement}
     end
 end
 
+"""
+    PhaseSpacePoint
+
+Representation of a point in the phase space of a process. Contains the process ([`AbstractProcessDefinition`](@ref)), the model ([`AbstractModelDefinition`](@ref)), the phase space definition ([`AbstractPhasespaceDefinition`]), and stateful incoming and outgoing particles ([`ParticleStateful`](@ref)).
+
+The legality of the combination of the given process and the incoming and outgoing particles is checked on construction. If the numbers of particles mismatch, the types of particles mismatch (note that order is important), or incoming particles have a `Outgoing` direction, an error is thrown.
+"""
 struct PhaseSpacePoint{
     PROC<:AbstractProcessDefinition,
     MODEL<:AbstractModelDefinition,
@@ -158,16 +168,40 @@ struct PhaseSpacePoint{
         IN_P<:AbstractVector{ParticleStateful{PhaseSpaceElementType}},
         OUT_P<:AbstractVector{ParticleStateful{PhaseSpaceElementType}},
     }
-        @assert length(incoming_particles(proc)) == length(in_p)
-        @assert length(outgoing_particles(proc)) == length(out_p)
+        length(incoming_particles(proc)) == length(in_p) || throw(
+            InvalidInputError(
+                "The number of incoming particles given by the process ($(incoming_particles(proc))) mismatches the number of given stateful incoming particles ($(length(in_p)))",
+            ),
+        )
+        length(outgoing_particles(proc)) == length(out_p) || throw(
+            InvalidInputError(
+                "The number of outgoing particles given by the process ($(outgoing_particles(proc))) mismatches the number of given stateful outgoing particles ($(length(out_p)))",
+            ),
+        )
 
         for (proc_p, p) in zip(incoming_particles(proc), in_p)
-            @assert proc_p == p.type "Process given particle type ($(proc_p)) does not match stateful particle type ($(p.type))"
-            @assert is_incoming(p) "Stateful particle $(p) is given as an incoming particle but is outgoing"
+            proc_p == p.type || throw(
+                InvalidInputError(
+                    "Process given particle type ($(proc_p)) does not match stateful particle type ($(p.type))",
+                ),
+            )
+            is_incoming(p) || throw(
+                InvalidInputError(
+                    "Stateful particle $(p) is given as an incoming particle but is outgoing",
+                ),
+            )
         end
         for (proc_p, p) in zip(outgoing_particles(proc), out_p)
-            @assert proc_p == p.type "Process given particle type ($(proc_p)) does not match stateful particle type ($(p.type))"
-            @assert is_outgoing(p) "Stateful particle $(p) is given as an outgoing particle but is incoming"
+            proc_p == p.type || throw(
+                InvalidInputError(
+                    "Process given particle type ($(proc_p)) does not match stateful particle type ($(p.type))",
+                ),
+            )
+            is_outgoing(p) || throw(
+                InvalidInputError(
+                    "Stateful particle $(p) is given as an outgoing particle but is incoming",
+                ),
+            )
         end
 
         return new{PROC,MODEL,PSDEF,PhaseSpaceElementType,length(in_p),length(out_p)}(

From e7f6f90d9eb894fd121bc2eb6c6942e3aa5dded5 Mon Sep 17 00:00:00 2001
From: Anton Reinhard <anton.reinhard@proton.me>
Date: Fri, 3 May 2024 21:41:51 +0200
Subject: [PATCH 06/12] seed test randomness and fix tests for julia versions
 <= 1.7

---
 test/phase_spaces.jl | 47 ++++++++++++++++++++++++++++++++++----------
 1 file changed, 37 insertions(+), 10 deletions(-)

diff --git a/test/phase_spaces.jl b/test/phase_spaces.jl
index eaf9f43..7ebc3d8 100644
--- a/test/phase_spaces.jl
+++ b/test/phase_spaces.jl
@@ -10,13 +10,15 @@ include("test_implementation/TestImplementation.jl")
 TESTMODEL = TestImplementation.TestModel()
 TESTPSDEF = TestImplementation.TestPhasespaceDef()
 
+RNG = Random.MersenneTwister(727)
+
 @testset "Stateful Particle" begin
     DIRECTIONS = [Incoming(), Outgoing()]
     PARTICLES = [Electron(), Positron()] #=, Muon(), AntiMuon(), Tauon(), AntiTauon()=#
     SPINANDPOLS = [AllSpin(), SpinUp(), SpinDown(), AllPol(), PolX(), PolY()]
 
     for (particle, dir, spinorpol) in Iterators.product(PARTICLES, DIRECTIONS, SPINANDPOLS)
-        momentum = rand(SFourMomentum)
+        momentum = rand(RNG, SFourMomentum)
 
         if (is_fermion(particle) && (spinorpol isa AbstractSpin)) ||
             (is_boson(particle) && (spinorpol isa AbstractPolarization))
@@ -30,7 +32,13 @@ TESTPSDEF = TestImplementation.TestPhasespaceDef()
             @test is_incoming(particle_stateful) == is_incoming(dir)
             @test is_outgoing(particle_stateful) == is_outgoing(dir)
         else
-            @test_throws "Tried to construct a stateful $(is_fermion(particle) ? "fermion" : "boson") with a $(spinorpol) instead of a $(is_fermion(particle) ? "spin" : "polarization")" ParticleStateful(
+            if (VERSION >= v"1.8")
+                # julia versions before 1.8 did not have support for regex matching in @test_throws
+                @test_throws "Tried to construct a stateful $(is_fermion(particle) ? "fermion" : "boson") with a $(spinorpol) instead of a $(is_fermion(particle) ? "spin" : "polarization")" ParticleStateful(
+                    momentum, particle, spinorpol, dir
+                )
+            end
+            @test_throws InvalidInputError ParticleStateful(
                 momentum, particle, spinorpol, dir
             )
         end
@@ -38,10 +46,10 @@ TESTPSDEF = TestImplementation.TestPhasespaceDef()
 end
 
 @testset "Phasespace Point" begin
-    in_el = ParticleStateful(rand(SFourMomentum), Electron(), AllSpin(), Incoming())
-    in_ph = ParticleStateful(rand(SFourMomentum), Photon(), AllPol(), Incoming())
-    out_el = ParticleStateful(rand(SFourMomentum), Electron(), AllSpin(), Outgoing())
-    out_ph = ParticleStateful(rand(SFourMomentum), Photon(), AllPol(), Outgoing())
+    in_el = ParticleStateful(rand(RNG, SFourMomentum), Electron(), AllSpin(), Incoming())
+    in_ph = ParticleStateful(rand(RNG, SFourMomentum), Photon(), AllPol(), Incoming())
+    out_el = ParticleStateful(rand(RNG, SFourMomentum), Electron(), AllSpin(), Outgoing())
+    out_ph = ParticleStateful(rand(RNG, SFourMomentum), Photon(), AllPol(), Outgoing())
 
     in_particles_valid = SVector(in_el, in_ph)
     in_particles_invalid = SVector(in_el, out_ph)
@@ -58,19 +66,38 @@ end
 
     PhaseSpacePoint(process, model, phasespace_def, in_particles_valid, out_particles_valid)
 
-    @test_throws r"Stateful particle (.*) is given as an incoming particle but is outgoing" PhaseSpacePoint(
+    if (VERSION >= v"1.8")
+        # julia versions before 1.8 did not have support for regex matching in @test_throws
+        @test_throws r"Stateful particle (.*) is given as an incoming particle but is outgoing" PhaseSpacePoint(
+            process, model, phasespace_def, in_particles_invalid, out_particles_valid
+        )
+
+        @test_throws r"Stateful particle (.*) is given as an outgoing particle but is incoming" PhaseSpacePoint(
+            process, model, phasespace_def, in_particles_valid, out_particles_invalid
+        )
+
+        @test_throws r"Process given particle type \((.*)Electron\(\)\) does not match stateful particle type \((.*)Photon\(\)\)" PhaseSpacePoint(
+            process, model, phasespace_def, SVector(in_ph, in_el), out_particles_valid
+        )
+
+        @test_throws r"Process given particle type \((.*)Electron\(\)\) does not match stateful particle type \((.*)Photon\(\)\)" PhaseSpacePoint(
+            process, model, phasespace_def, in_particles_valid, SVector(out_ph, out_el)
+        )
+    end
+
+    @test_throws InvalidInputError PhaseSpacePoint(
         process, model, phasespace_def, in_particles_invalid, out_particles_valid
     )
 
-    @test_throws r"Stateful particle (.*) is given as an outgoing particle but is incoming" PhaseSpacePoint(
+    @test_throws InvalidInputError PhaseSpacePoint(
         process, model, phasespace_def, in_particles_valid, out_particles_invalid
     )
 
-    @test_throws r"Process given particle type \((.*)Electron\(\)\) does not match stateful particle type \((.*)Photon\(\)\)" PhaseSpacePoint(
+    @test_throws InvalidInputError PhaseSpacePoint(
         process, model, phasespace_def, SVector(in_ph, in_el), out_particles_valid
     )
 
-    @test_throws r"Process given particle type \((.*)Electron\(\)\) does not match stateful particle type \((.*)Photon\(\)\)" PhaseSpacePoint(
+    @test_throws InvalidInputError PhaseSpacePoint(
         process, model, phasespace_def, in_particles_valid, SVector(out_ph, out_el)
     )
 end

From 65951def6222b7a0a0809c4ac2f3c2a7e5cd7a85 Mon Sep 17 00:00:00 2001
From: Anton Reinhard <anton.reinhard@proton.me>
Date: Mon, 6 May 2024 18:45:23 +0200
Subject: [PATCH 07/12] Review changes

---
 src/phase_spaces.jl  | 90 ++++++++++++++++++++++++--------------------
 test/phase_spaces.jl | 29 +++++++-------
 2 files changed, 63 insertions(+), 56 deletions(-)

diff --git a/src/phase_spaces.jl b/src/phase_spaces.jl
index fd31563..5131f92 100644
--- a/src/phase_spaces.jl
+++ b/src/phase_spaces.jl
@@ -15,6 +15,21 @@ import QEDbase.is_fermion
 import QEDbase.is_boson
 import QEDbase.is_incoming
 import QEDbase.is_outgoing
+import QEDbase.mass
+import QEDbase.charge
+
+import QEDbase:
+    is_particle,
+    is_anti_particle,
+    is_fermion,
+    is_boson,
+    is_incoming,
+    is_outgoing,
+    mass,
+    charge,
+    AbstractFourMomentum,
+    AbstractParticleType,
+    AbstractParticle
 
 abstract type AbstractCoordinateSystem end
 struct SphericalCoordinateSystem <: AbstractCoordinateSystem end
@@ -41,7 +56,7 @@ end
 #
 # Currently, elements can be either four-momenta, or real numbers,
 # i.e. coordinates.
-AbstractPhasespaceElement = Union{QEDbase.AbstractFourMomentum,Real}
+AbstractPhasespaceElement = Union{AbstractFourMomentum,Real}
 
 # utility functions 
 
@@ -49,7 +64,7 @@ AbstractPhasespaceElement = Union{QEDbase.AbstractFourMomentum,Real}
     proc::AbstractProcessDefinition,
     ::AbstractModelDefinition,
     in_phase_space::AbstractVecOrMat{T},
-) where {T<:QEDbase.AbstractFourMomentum}
+) where {T<:AbstractFourMomentum}
     return size(in_phase_space, 1) == number_incoming_particles(proc) || throw(
         DimensionMismatch(
             "The number of incoming particles <$(number_incoming_particles(proc))> is inconsistent with input size <$(size(in_phase_space,1))>",
@@ -61,7 +76,7 @@ end
     proc::AbstractProcessDefinition,
     ::AbstractModelDefinition,
     out_phase_space::AbstractVecOrMat{T},
-) where {T<:QEDbase.AbstractFourMomentum}
+) where {T<:AbstractFourMomentum}
     return size(out_phase_space, 1) == number_outgoing_particles(proc) || throw(
         DimensionMismatch(
             "The number of outgoing particles <$(number_outgoing_particles(proc))> is inconsistent with input size <$(size(out_phase_space,1))>",
@@ -94,45 +109,36 @@ end
 end
 
 """
-    ParticleStateful
+    ParticleStateful <: AbstractParticle
 
 Representation of a particle with a state. It has four fields:
-- `mom::AbstractPhasespaceElement`: The momentum of the particle.
-- `type::AbstractParticle`: The type of the particle, `QEDbase.Electron()`, `QEDbase.Positron()` etc..
-- `spinorpol::AbstractSpinOrPolarization`: The spin or polarization of the particle, `QEDbase.SpinUp()`, `QEDbase.PolX() etc.. Can only use spins with fermions and polarizations with bosons.
 - `dir::ParticleDirection`: The direction of the particle, `QEDbase.Incoming()` or `QEDbase.Outgoing()`.
+- `species::AbstractParticleType`: The species of the particle, `QEDbase.Electron()`, `QEDbase.Positron()` etc..
+- `mom::AbstractFourMomentum`: The momentum of the particle.
+- `spin_or_pol::AbstractSpinOrPolarization`: The spin or polarization of the particle, `QEDbase.SpinUp()`, `QEDbase.PolX() etc.. Can only use spins with fermions and polarizations with bosons.
 
-Overloads for `QEDbase.is_fermion`, `QEDbase.is_boson`, `QEDbase.is_particle`, `QEDbase.is_anti_particle`, `QEDbase.is_incoming`, and `QEDbase.is_outgoing` are provided, delegating the call to the correct field.
+Overloads for `QEDbase.is_fermion`, `QEDbase.is_boson`, `QEDbase.is_particle`, `QEDbase.is_anti_particle`, `QEDbase.is_incoming`, `QEDbase.is_outgoing`, `QEDbase.mass`, and `QEDbase.charge` are provided, delegating the call to the correct field and thus implementing the `QEDbase.AbstractParticle` interface.
 
 The legality of the combination of `type` and `spinorpol` is checked on construction. If, for example, the construction of an `Electron()` with a polarization is attempted, an [`InvalidInputError`](@ref) is thrown.
 """
-struct ParticleStateful{ElType<:AbstractPhasespaceElement}
-    mom::ElType
-    type::QEDbase.AbstractParticle
-    spinorpol::AbstractSpinOrPolarization
+struct ParticleStateful{ElType<:AbstractFourMomentum} <: AbstractParticle
     dir::ParticleDirection
+    species::AbstractParticleType
+    mom::ElType
+    spin_or_pol::AbstractSpinOrPolarization
 
     function ParticleStateful(
-        mom::ElType,
-        type::QEDbase.AbstractParticle,
-        spinorpol::AbstractSpinOrPolarization,
-        dir::ParticleDirection,
-    ) where {ElType<:AbstractPhasespaceElement}
-        if is_fermion(type) && !(spinorpol isa AbstractSpin)
-            throw(
-                InvalidInputError(
-                    "Tried to construct a stateful fermion with a $(spinorpol) instead of a spin",
-                ),
-            )
-        elseif is_boson(type) && !(spinorpol isa AbstractPolarization)
-            throw(
-                InvalidInputError(
-                    "Tried to construct a stateful boson with a $(spinorpol) instead of a polarization",
-                ),
-            )
-        end
+        dir::ParticleDirection, species::Species, mom::ElType, spin::Spin=AllSpin()
+    ) where {Species<:FermionLike,ElType<:AbstractFourMomentum,Spin<:AbstractSpin}
+        # constructor for fermions with spin
+        return new{ElType}(dir, species, mom, spin)
+    end
 
-        return new{ElType}(mom, type, spinorpol, dir)
+    function ParticleStateful(
+        dir::ParticleDirection, species::Species, mom::ElType, pol::Pol=AllPol()
+    ) where {Species<:BosonLike,ElType<:AbstractFourMomentum,Pol<:AbstractPolarization}
+        # constructor for bosons with polarization
+        return new{ElType}(dir, species, mom, pol)
     end
 end
 
@@ -147,7 +153,7 @@ struct PhaseSpacePoint{
     PROC<:AbstractProcessDefinition,
     MODEL<:AbstractModelDefinition,
     PSDEF<:AbstractPhasespaceDefinition,
-    PhaseSpaceElementType<:AbstractPhasespaceElement,
+    PhaseSpaceElementType<:AbstractFourMomentum,
     N_IN_PARTICLES,
     N_OUT_PARTICLES,
 }
@@ -164,7 +170,7 @@ struct PhaseSpacePoint{
         PROC<:AbstractProcessDefinition,
         MODEL<:AbstractModelDefinition,
         PSDEF<:AbstractPhasespaceDefinition,
-        PhaseSpaceElementType<:AbstractPhasespaceElement,
+        PhaseSpaceElementType<:AbstractFourMomentum,
         IN_P<:AbstractVector{ParticleStateful{PhaseSpaceElementType}},
         OUT_P<:AbstractVector{ParticleStateful{PhaseSpaceElementType}},
     }
@@ -180,9 +186,9 @@ struct PhaseSpacePoint{
         )
 
         for (proc_p, p) in zip(incoming_particles(proc), in_p)
-            proc_p == p.type || throw(
+            proc_p == p.species || throw(
                 InvalidInputError(
-                    "Process given particle type ($(proc_p)) does not match stateful particle type ($(p.type))",
+                    "Process given particle species ($(proc_p)) does not match stateful particle species ($(p.species))",
                 ),
             )
             is_incoming(p) || throw(
@@ -192,9 +198,9 @@ struct PhaseSpacePoint{
             )
         end
         for (proc_p, p) in zip(outgoing_particles(proc), out_p)
-            proc_p == p.type || throw(
+            proc_p == p.species || throw(
                 InvalidInputError(
-                    "Process given particle type ($(proc_p)) does not match stateful particle type ($(p.type))",
+                    "Process given particle species ($(proc_p)) does not match stateful particle species ($(p.species))",
                 ),
             )
             is_outgoing(p) || throw(
@@ -212,7 +218,9 @@ end
 
 @inline is_incoming(particle::ParticleStateful) = is_incoming(particle.dir)
 @inline is_outgoing(particle::ParticleStateful) = is_outgoing(particle.dir)
-@inline is_fermion(particle::ParticleStateful) = is_fermion(particle.type)
-@inline is_boson(particle::ParticleStateful) = is_boson(particle.type)
-@inline is_particle(particle::ParticleStateful) = is_particle(particle.type)
-@inline is_anti_particle(particle::ParticleStateful) = is_anti_particle(particle.type)
+@inline is_fermion(particle::ParticleStateful) = is_fermion(particle.species)
+@inline is_boson(particle::ParticleStateful) = is_boson(particle.species)
+@inline is_particle(particle::ParticleStateful) = is_particle(particle.species)
+@inline is_anti_particle(particle::ParticleStateful) = is_anti_particle(particle.species)
+@inline mass(particle::ParticleStateful) = mass(particle.species)
+@inline charge(particle::ParticleStateful) = charge(particle.species)
diff --git a/test/phase_spaces.jl b/test/phase_spaces.jl
index 7ebc3d8..f3b281f 100644
--- a/test/phase_spaces.jl
+++ b/test/phase_spaces.jl
@@ -17,12 +17,13 @@ RNG = Random.MersenneTwister(727)
     PARTICLES = [Electron(), Positron()] #=, Muon(), AntiMuon(), Tauon(), AntiTauon()=#
     SPINANDPOLS = [AllSpin(), SpinUp(), SpinDown(), AllPol(), PolX(), PolY()]
 
-    for (particle, dir, spinorpol) in Iterators.product(PARTICLES, DIRECTIONS, SPINANDPOLS)
+    for (particle, dir, spin_or_pol) in
+        Iterators.product(PARTICLES, DIRECTIONS, SPINANDPOLS)
         momentum = rand(RNG, SFourMomentum)
 
-        if (is_fermion(particle) && (spinorpol isa AbstractSpin)) ||
-            (is_boson(particle) && (spinorpol isa AbstractPolarization))
-            particle_stateful = ParticleStateful(momentum, particle, spinorpol, dir)
+        if (is_fermion(particle) && (spin_or_pol isa AbstractSpin)) ||
+            (is_boson(particle) && (spin_or_pol isa AbstractPolarization))
+            particle_stateful = ParticleStateful(dir, particle, momentum, spin_or_pol)
 
             @test particle_stateful.mom == momentum
             @test is_fermion(particle_stateful) == is_fermion(particle)
@@ -34,22 +35,20 @@ RNG = Random.MersenneTwister(727)
         else
             if (VERSION >= v"1.8")
                 # julia versions before 1.8 did not have support for regex matching in @test_throws
-                @test_throws "Tried to construct a stateful $(is_fermion(particle) ? "fermion" : "boson") with a $(spinorpol) instead of a $(is_fermion(particle) ? "spin" : "polarization")" ParticleStateful(
-                    momentum, particle, spinorpol, dir
+                @test_throws "MethodError: no method matching ParticleStateful" ParticleStateful(
+                    dir, particle, momentum, spin_or_pol
                 )
             end
-            @test_throws InvalidInputError ParticleStateful(
-                momentum, particle, spinorpol, dir
-            )
+            @test_throws MethodError ParticleStateful(dir, particle, momentum, spin_or_pol)
         end
     end
 end
 
 @testset "Phasespace Point" begin
-    in_el = ParticleStateful(rand(RNG, SFourMomentum), Electron(), AllSpin(), Incoming())
-    in_ph = ParticleStateful(rand(RNG, SFourMomentum), Photon(), AllPol(), Incoming())
-    out_el = ParticleStateful(rand(RNG, SFourMomentum), Electron(), AllSpin(), Outgoing())
-    out_ph = ParticleStateful(rand(RNG, SFourMomentum), Photon(), AllPol(), Outgoing())
+    in_el = ParticleStateful(Incoming(), Electron(), rand(RNG, SFourMomentum))
+    in_ph = ParticleStateful(Incoming(), Photon(), rand(RNG, SFourMomentum))
+    out_el = ParticleStateful(Outgoing(), Electron(), rand(RNG, SFourMomentum))
+    out_ph = ParticleStateful(Outgoing(), Photon(), rand(RNG, SFourMomentum))
 
     in_particles_valid = SVector(in_el, in_ph)
     in_particles_invalid = SVector(in_el, out_ph)
@@ -76,11 +75,11 @@ end
             process, model, phasespace_def, in_particles_valid, out_particles_invalid
         )
 
-        @test_throws r"Process given particle type \((.*)Electron\(\)\) does not match stateful particle type \((.*)Photon\(\)\)" PhaseSpacePoint(
+        @test_throws r"Process given particle species \((.*)Electron\(\)\) does not match stateful particle species \((.*)Photon\(\)\)" PhaseSpacePoint(
             process, model, phasespace_def, SVector(in_ph, in_el), out_particles_valid
         )
 
-        @test_throws r"Process given particle type \((.*)Electron\(\)\) does not match stateful particle type \((.*)Photon\(\)\)" PhaseSpacePoint(
+        @test_throws r"Process given particle species \((.*)Electron\(\)\) does not match stateful particle species \((.*)Photon\(\)\)" PhaseSpacePoint(
             process, model, phasespace_def, in_particles_valid, SVector(out_ph, out_el)
         )
     end

From b9b43696b479625256b08c8bb74904992ae0c60d Mon Sep 17 00:00:00 2001
From: Anton Reinhard <anton.reinhard@proton.me>
Date: Mon, 6 May 2024 18:56:37 +0200
Subject: [PATCH 08/12] Add spin and pol convenience functions for
 ParticleStateful

---
 src/QEDprocesses.jl  |  1 +
 src/phase_spaces.jl  |  8 ++++++++
 test/phase_spaces.jl | 33 +++++++++++++++++++++------------
 3 files changed, 30 insertions(+), 12 deletions(-)

diff --git a/src/QEDprocesses.jl b/src/QEDprocesses.jl
index 15d8ad8..36681ea 100644
--- a/src/QEDprocesses.jl
+++ b/src/QEDprocesses.jl
@@ -31,6 +31,7 @@ export AbstractCoordinateSystem, SphericalCoordinateSystem
 export AbstractFrameOfReference, CenterOfMomentumFrame, ElectronRestFrame
 export AbstractPhasespaceDefinition, PhasespaceDefinition
 export ParticleStateful, PhaseSpacePoint
+export spin, pol
 
 using QEDbase
 
diff --git a/src/phase_spaces.jl b/src/phase_spaces.jl
index 5131f92..816b60c 100644
--- a/src/phase_spaces.jl
+++ b/src/phase_spaces.jl
@@ -224,3 +224,11 @@ end
 @inline is_anti_particle(particle::ParticleStateful) = is_anti_particle(particle.species)
 @inline mass(particle::ParticleStateful) = mass(particle.species)
 @inline charge(particle::ParticleStateful) = charge(particle.species)
+
+@inline _spin(::Species, particle::ParticleStateful) where {Species<:FermionLike} =
+    particle.spin_or_pol
+@inline spin(particle::ParticleStateful) = _spin(particle.species, particle)
+
+@inline _pol(::Species, particle::ParticleStateful) where {Species<:BosonLike} =
+    particle.spin_or_pol
+@inline pol(particle::ParticleStateful) = _pol(particle.species, particle)
diff --git a/test/phase_spaces.jl b/test/phase_spaces.jl
index f3b281f..c11b1a9 100644
--- a/test/phase_spaces.jl
+++ b/test/phase_spaces.jl
@@ -14,32 +14,41 @@ RNG = Random.MersenneTwister(727)
 
 @testset "Stateful Particle" begin
     DIRECTIONS = [Incoming(), Outgoing()]
-    PARTICLES = [Electron(), Positron()] #=, Muon(), AntiMuon(), Tauon(), AntiTauon()=#
+    SPECIES = [Electron(), Positron()] #=, Muon(), AntiMuon(), Tauon(), AntiTauon()=#
     SPINANDPOLS = [AllSpin(), SpinUp(), SpinDown(), AllPol(), PolX(), PolY()]
 
-    for (particle, dir, spin_or_pol) in
-        Iterators.product(PARTICLES, DIRECTIONS, SPINANDPOLS)
+    for (species, dir, spin_or_pol) in Iterators.product(SPECIES, DIRECTIONS, SPINANDPOLS)
         momentum = rand(RNG, SFourMomentum)
 
-        if (is_fermion(particle) && (spin_or_pol isa AbstractSpin)) ||
-            (is_boson(particle) && (spin_or_pol isa AbstractPolarization))
-            particle_stateful = ParticleStateful(dir, particle, momentum, spin_or_pol)
+        if (is_fermion(species) && (spin_or_pol isa AbstractSpin)) ||
+            (is_boson(species) && (spin_or_pol isa AbstractPolarization))
+            particle_stateful = ParticleStateful(dir, species, momentum, spin_or_pol)
 
             @test particle_stateful.mom == momentum
-            @test is_fermion(particle_stateful) == is_fermion(particle)
-            @test is_boson(particle_stateful) == is_boson(particle)
-            @test is_particle(particle_stateful) == is_particle(particle)
-            @test is_anti_particle(particle_stateful) == is_anti_particle(particle)
+            @test is_fermion(particle_stateful) == is_fermion(species)
+            @test is_boson(particle_stateful) == is_boson(species)
+            @test is_particle(particle_stateful) == is_particle(species)
+            @test is_anti_particle(particle_stateful) == is_anti_particle(species)
             @test is_incoming(particle_stateful) == is_incoming(dir)
             @test is_outgoing(particle_stateful) == is_outgoing(dir)
+            @test mass(particle_stateful) == mass(species)
+            @test charge(particle_stateful) == charge(species)
+
+            if (is_fermion(species))
+                @test spin(particle_stateful) == spin_or_pol
+                @test_throws MethodError pol(particle_stateful)
+            else
+                @test pol(particle_stateful) == spin_or_pol
+                @test_throws MethodError spin(particle_stateful)
+            end
         else
             if (VERSION >= v"1.8")
                 # julia versions before 1.8 did not have support for regex matching in @test_throws
                 @test_throws "MethodError: no method matching ParticleStateful" ParticleStateful(
-                    dir, particle, momentum, spin_or_pol
+                    dir, species, momentum, spin_or_pol
                 )
             end
-            @test_throws MethodError ParticleStateful(dir, particle, momentum, spin_or_pol)
+            @test_throws MethodError ParticleStateful(dir, species, momentum, spin_or_pol)
         end
     end
 end

From 5a5f85596a613a670d5f256ba4cc3dca3431ff7a Mon Sep 17 00:00:00 2001
From: Anton Reinhard <anton.reinhard@proton.me>
Date: Mon, 6 May 2024 19:31:34 +0200
Subject: [PATCH 09/12] Add convenient accessor functions for phase space point

---
 src/QEDprocesses.jl  |  2 +-
 src/phase_spaces.jl  | 55 +++++++++++++++++++++++++++++++++-----------
 test/phase_spaces.jl | 24 ++++++++++++++++++-
 3 files changed, 65 insertions(+), 16 deletions(-)

diff --git a/src/QEDprocesses.jl b/src/QEDprocesses.jl
index 36681ea..a811b74 100644
--- a/src/QEDprocesses.jl
+++ b/src/QEDprocesses.jl
@@ -31,7 +31,7 @@ export AbstractCoordinateSystem, SphericalCoordinateSystem
 export AbstractFrameOfReference, CenterOfMomentumFrame, ElectronRestFrame
 export AbstractPhasespaceDefinition, PhasespaceDefinition
 export ParticleStateful, PhaseSpacePoint
-export spin, pol
+export spin, pol, nth_momentum, getindex
 
 using QEDbase
 
diff --git a/src/phase_spaces.jl b/src/phase_spaces.jl
index 816b60c..44bf7f3 100644
--- a/src/phase_spaces.jl
+++ b/src/phase_spaces.jl
@@ -142,6 +142,23 @@ struct ParticleStateful{ElType<:AbstractFourMomentum} <: AbstractParticle
     end
 end
 
+@inline is_incoming(particle::ParticleStateful) = is_incoming(particle.dir)
+@inline is_outgoing(particle::ParticleStateful) = is_outgoing(particle.dir)
+@inline is_fermion(particle::ParticleStateful) = is_fermion(particle.species)
+@inline is_boson(particle::ParticleStateful) = is_boson(particle.species)
+@inline is_particle(particle::ParticleStateful) = is_particle(particle.species)
+@inline is_anti_particle(particle::ParticleStateful) = is_anti_particle(particle.species)
+@inline mass(particle::ParticleStateful) = mass(particle.species)
+@inline charge(particle::ParticleStateful) = charge(particle.species)
+
+@inline _spin(::Species, particle::ParticleStateful) where {Species<:FermionLike} =
+    particle.spin_or_pol
+@inline spin(particle::ParticleStateful) = _spin(particle.species, particle)
+
+@inline _pol(::Species, particle::ParticleStateful) where {Species<:BosonLike} =
+    particle.spin_or_pol
+@inline pol(particle::ParticleStateful) = _pol(particle.species, particle)
+
 """
     PhaseSpacePoint
 
@@ -216,19 +233,29 @@ struct PhaseSpacePoint{
     end
 end
 
-@inline is_incoming(particle::ParticleStateful) = is_incoming(particle.dir)
-@inline is_outgoing(particle::ParticleStateful) = is_outgoing(particle.dir)
-@inline is_fermion(particle::ParticleStateful) = is_fermion(particle.species)
-@inline is_boson(particle::ParticleStateful) = is_boson(particle.species)
-@inline is_particle(particle::ParticleStateful) = is_particle(particle.species)
-@inline is_anti_particle(particle::ParticleStateful) = is_anti_particle(particle.species)
-@inline mass(particle::ParticleStateful) = mass(particle.species)
-@inline charge(particle::ParticleStateful) = charge(particle.species)
+"""
+    Base.getindex(psp::PhaseSpacePoint, dir::Incoming, n::Int)
 
-@inline _spin(::Species, particle::ParticleStateful) where {Species<:FermionLike} =
-    particle.spin_or_pol
-@inline spin(particle::ParticleStateful) = _spin(particle.species, particle)
+Overload for the array indexing operator `[]`. Returns the nth incoming particle in this phase space point.
+"""
+function Base.getindex(psp::PhaseSpacePoint, ::Incoming, n::Int)
+    return psp.in_particles[n]
+end
 
-@inline _pol(::Species, particle::ParticleStateful) where {Species<:BosonLike} =
-    particle.spin_or_pol
-@inline pol(particle::ParticleStateful) = _pol(particle.species, particle)
+"""
+    Base.getindex(psp::PhaseSpacePoint, dir::Outgoing, n::Int)
+
+Overload for the array indexing operator `[]`. Returns the nth outgoing particle in this phase space point.
+"""
+function Base.getindex(psp::PhaseSpacePoint, ::Outgoing, n::Int)
+    return psp.out_particles[n]
+end
+
+"""
+    nth_momentum(psp::PhaseSpacePoint, dir::ParticleDirection, n::Int)
+
+Returns the momentum of the `n`th particle in the given [`PhaseSpacePoint`](@ref) which has direction `dir`. If `n` is outside the valid range for this phase space point, an `BoundsError` is thrown.
+"""
+function nth_momentum(psp::PhaseSpacePoint, dir::ParticleDirection, n::Int)
+    return psp[dir, n].mom
+end
diff --git a/test/phase_spaces.jl b/test/phase_spaces.jl
index c11b1a9..0d5aafc 100644
--- a/test/phase_spaces.jl
+++ b/test/phase_spaces.jl
@@ -72,7 +72,19 @@ end
     )
     phasespace_def = TESTPSDEF
 
-    PhaseSpacePoint(process, model, phasespace_def, in_particles_valid, out_particles_valid)
+    psp = PhaseSpacePoint(
+        process, model, phasespace_def, in_particles_valid, out_particles_valid
+    )
+
+    @test nth_momentum(psp, Incoming(), 1) == in_el.mom
+    @test nth_momentum(psp, Incoming(), 2) == in_ph.mom
+    @test nth_momentum(psp, Outgoing(), 1) == out_el.mom
+    @test nth_momentum(psp, Outgoing(), 2) == out_ph.mom
+
+    @test psp[Incoming(), 1] == in_el
+    @test psp[Incoming(), 2] == in_ph
+    @test psp[Outgoing(), 1] == out_el
+    @test psp[Outgoing(), 2] == out_ph
 
     if (VERSION >= v"1.8")
         # julia versions before 1.8 did not have support for regex matching in @test_throws
@@ -93,6 +105,16 @@ end
         )
     end
 
+    @test_throws BoundsError nth_momentum(psp, Incoming(), -1)
+    @test_throws BoundsError nth_momentum(psp, Outgoing(), -1)
+    @test_throws BoundsError nth_momentum(psp, Incoming(), 4)
+    @test_throws BoundsError nth_momentum(psp, Outgoing(), 4)
+
+    @test_throws BoundsError psp[Incoming(), -1]
+    @test_throws BoundsError psp[Outgoing(), -1]
+    @test_throws BoundsError psp[Incoming(), 4]
+    @test_throws BoundsError psp[Outgoing(), 4]
+
     @test_throws InvalidInputError PhaseSpacePoint(
         process, model, phasespace_def, in_particles_invalid, out_particles_valid
     )

From c8e3c51cf1d4f64c08cb16a764453060e95b1886 Mon Sep 17 00:00:00 2001
From: Anton Reinhard <anton.reinhard@proton.me>
Date: Mon, 6 May 2024 21:53:09 +0200
Subject: [PATCH 10/12] Apply suggestions from code review

Co-authored-by: Uwe Hernandez Acosta <szabo137@users.noreply.github.com>
---
 src/phase_spaces.jl | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/phase_spaces.jl b/src/phase_spaces.jl
index 44bf7f3..0acb8d5 100644
--- a/src/phase_spaces.jl
+++ b/src/phase_spaces.jl
@@ -113,13 +113,13 @@ end
 
 Representation of a particle with a state. It has four fields:
 - `dir::ParticleDirection`: The direction of the particle, `QEDbase.Incoming()` or `QEDbase.Outgoing()`.
-- `species::AbstractParticleType`: The species of the particle, `QEDbase.Electron()`, `QEDbase.Positron()` etc..
+- `species::AbstractParticleType`: The species of the particle, `QEDbase.Electron()`, `QEDbase.Positron()` etc.
 - `mom::AbstractFourMomentum`: The momentum of the particle.
-- `spin_or_pol::AbstractSpinOrPolarization`: The spin or polarization of the particle, `QEDbase.SpinUp()`, `QEDbase.PolX() etc.. Can only use spins with fermions and polarizations with bosons.
+- `spin_or_pol::AbstractSpinOrPolarization`: The spin or polarization of the particle, `QEDbase.SpinUp()`, `QEDbase.PolX() etc. Can only use spins with fermions and polarizations with bosons.
 
 Overloads for `QEDbase.is_fermion`, `QEDbase.is_boson`, `QEDbase.is_particle`, `QEDbase.is_anti_particle`, `QEDbase.is_incoming`, `QEDbase.is_outgoing`, `QEDbase.mass`, and `QEDbase.charge` are provided, delegating the call to the correct field and thus implementing the `QEDbase.AbstractParticle` interface.
 
-The legality of the combination of `type` and `spinorpol` is checked on construction. If, for example, the construction of an `Electron()` with a polarization is attempted, an [`InvalidInputError`](@ref) is thrown.
+The legality of the combination of `species` and `spin_or_pol` is checked on construction. If, for example, the construction of an `Electron()` with a polarization is attempted, an [`InvalidInputError`](@ref) is thrown.
 """
 struct ParticleStateful{ElType<:AbstractFourMomentum} <: AbstractParticle
     dir::ParticleDirection
@@ -164,7 +164,7 @@ end
 
 Representation of a point in the phase space of a process. Contains the process ([`AbstractProcessDefinition`](@ref)), the model ([`AbstractModelDefinition`](@ref)), the phase space definition ([`AbstractPhasespaceDefinition`]), and stateful incoming and outgoing particles ([`ParticleStateful`](@ref)).
 
-The legality of the combination of the given process and the incoming and outgoing particles is checked on construction. If the numbers of particles mismatch, the types of particles mismatch (note that order is important), or incoming particles have a `Outgoing` direction, an error is thrown.
+The legality of the combination of the given process and the incoming and outgoing particles is checked on construction. If the numbers of particles mismatch, the types of particles mismatch (note that order is important), or incoming particles have an `Outgoing` direction, an error is thrown.
 """
 struct PhaseSpacePoint{
     PROC<:AbstractProcessDefinition,
@@ -254,7 +254,7 @@ end
 """
     nth_momentum(psp::PhaseSpacePoint, dir::ParticleDirection, n::Int)
 
-Returns the momentum of the `n`th particle in the given [`PhaseSpacePoint`](@ref) which has direction `dir`. If `n` is outside the valid range for this phase space point, an `BoundsError` is thrown.
+Returns the momentum of the `n`th particle in the given [`PhaseSpacePoint`](@ref) which has direction `dir`. If `n` is outside the valid range for this phase space point, a `BoundsError` is thrown.
 """
 function nth_momentum(psp::PhaseSpacePoint, dir::ParticleDirection, n::Int)
     return psp[dir, n].mom

From c6fd5ce01e407e7638c0560ea9cf8d9c61a2250a Mon Sep 17 00:00:00 2001
From: Anton Reinhard <anton.reinhard@proton.me>
Date: Mon, 6 May 2024 21:59:03 +0200
Subject: [PATCH 11/12] Renaming according to review

---
 src/QEDprocesses.jl  |  2 +-
 src/phase_spaces.jl  |  8 ++++----
 test/phase_spaces.jl | 20 ++++++++++----------
 3 files changed, 15 insertions(+), 15 deletions(-)

diff --git a/src/QEDprocesses.jl b/src/QEDprocesses.jl
index a811b74..60957e2 100644
--- a/src/QEDprocesses.jl
+++ b/src/QEDprocesses.jl
@@ -31,7 +31,7 @@ export AbstractCoordinateSystem, SphericalCoordinateSystem
 export AbstractFrameOfReference, CenterOfMomentumFrame, ElectronRestFrame
 export AbstractPhasespaceDefinition, PhasespaceDefinition
 export ParticleStateful, PhaseSpacePoint
-export spin, pol, nth_momentum, getindex
+export spin, polarization, momentum, getindex
 
 using QEDbase
 
diff --git a/src/phase_spaces.jl b/src/phase_spaces.jl
index 0acb8d5..7fef7a1 100644
--- a/src/phase_spaces.jl
+++ b/src/phase_spaces.jl
@@ -155,9 +155,9 @@ end
     particle.spin_or_pol
 @inline spin(particle::ParticleStateful) = _spin(particle.species, particle)
 
-@inline _pol(::Species, particle::ParticleStateful) where {Species<:BosonLike} =
+@inline _polarization(::Species, particle::ParticleStateful) where {Species<:BosonLike} =
     particle.spin_or_pol
-@inline pol(particle::ParticleStateful) = _pol(particle.species, particle)
+@inline polarization(particle::ParticleStateful) = _polarization(particle.species, particle)
 
 """
     PhaseSpacePoint
@@ -252,10 +252,10 @@ function Base.getindex(psp::PhaseSpacePoint, ::Outgoing, n::Int)
 end
 
 """
-    nth_momentum(psp::PhaseSpacePoint, dir::ParticleDirection, n::Int)
+    momentum(psp::PhaseSpacePoint, dir::ParticleDirection, n::Int)
 
 Returns the momentum of the `n`th particle in the given [`PhaseSpacePoint`](@ref) which has direction `dir`. If `n` is outside the valid range for this phase space point, a `BoundsError` is thrown.
 """
-function nth_momentum(psp::PhaseSpacePoint, dir::ParticleDirection, n::Int)
+function momentum(psp::PhaseSpacePoint, dir::ParticleDirection, n::Int)
     return psp[dir, n].mom
 end
diff --git a/test/phase_spaces.jl b/test/phase_spaces.jl
index 0d5aafc..743f443 100644
--- a/test/phase_spaces.jl
+++ b/test/phase_spaces.jl
@@ -36,9 +36,9 @@ RNG = Random.MersenneTwister(727)
 
             if (is_fermion(species))
                 @test spin(particle_stateful) == spin_or_pol
-                @test_throws MethodError pol(particle_stateful)
+                @test_throws MethodError polarization(particle_stateful)
             else
-                @test pol(particle_stateful) == spin_or_pol
+                @test polarization(particle_stateful) == spin_or_pol
                 @test_throws MethodError spin(particle_stateful)
             end
         else
@@ -76,10 +76,10 @@ end
         process, model, phasespace_def, in_particles_valid, out_particles_valid
     )
 
-    @test nth_momentum(psp, Incoming(), 1) == in_el.mom
-    @test nth_momentum(psp, Incoming(), 2) == in_ph.mom
-    @test nth_momentum(psp, Outgoing(), 1) == out_el.mom
-    @test nth_momentum(psp, Outgoing(), 2) == out_ph.mom
+    @test momentum(psp, Incoming(), 1) == in_el.mom
+    @test momentum(psp, Incoming(), 2) == in_ph.mom
+    @test momentum(psp, Outgoing(), 1) == out_el.mom
+    @test momentum(psp, Outgoing(), 2) == out_ph.mom
 
     @test psp[Incoming(), 1] == in_el
     @test psp[Incoming(), 2] == in_ph
@@ -105,10 +105,10 @@ end
         )
     end
 
-    @test_throws BoundsError nth_momentum(psp, Incoming(), -1)
-    @test_throws BoundsError nth_momentum(psp, Outgoing(), -1)
-    @test_throws BoundsError nth_momentum(psp, Incoming(), 4)
-    @test_throws BoundsError nth_momentum(psp, Outgoing(), 4)
+    @test_throws BoundsError momentum(psp, Incoming(), -1)
+    @test_throws BoundsError momentum(psp, Outgoing(), -1)
+    @test_throws BoundsError momentum(psp, Incoming(), 4)
+    @test_throws BoundsError momentum(psp, Outgoing(), 4)
 
     @test_throws BoundsError psp[Incoming(), -1]
     @test_throws BoundsError psp[Outgoing(), -1]

From c62aa75ff3b568afd8c3dea2d63baf6b13406f5d Mon Sep 17 00:00:00 2001
From: Anton Reinhard <anton.reinhard@proton.me>
Date: Mon, 6 May 2024 22:04:18 +0200
Subject: [PATCH 12/12] Error message dcapitalization

---
 src/interfaces/setup_interface.jl |  2 +-
 src/phase_spaces.jl               | 20 ++++++++++----------
 test/phase_spaces.jl              |  8 ++++----
 3 files changed, 15 insertions(+), 15 deletions(-)

diff --git a/src/interfaces/setup_interface.jl b/src/interfaces/setup_interface.jl
index fbae838..ac406f4 100644
--- a/src/interfaces/setup_interface.jl
+++ b/src/interfaces/setup_interface.jl
@@ -82,7 +82,7 @@ struct InvalidInputError <: AbstractInvalidInputException
     msg::String
 end
 function Base.showerror(io::IO, err::InvalidInputError)
-    return println(io, "InvalidInputError: $(err.msg).")
+    return println(io, "InvalidInputError: $(err.msg)")
 end
 
 """
diff --git a/src/phase_spaces.jl b/src/phase_spaces.jl
index 7fef7a1..96a7f13 100644
--- a/src/phase_spaces.jl
+++ b/src/phase_spaces.jl
@@ -67,7 +67,7 @@ AbstractPhasespaceElement = Union{AbstractFourMomentum,Real}
 ) where {T<:AbstractFourMomentum}
     return size(in_phase_space, 1) == number_incoming_particles(proc) || throw(
         DimensionMismatch(
-            "The number of incoming particles <$(number_incoming_particles(proc))> is inconsistent with input size <$(size(in_phase_space,1))>",
+            "the number of incoming particles <$(number_incoming_particles(proc))> is inconsistent with input size <$(size(in_phase_space,1))>",
         ),
     )
 end
@@ -79,7 +79,7 @@ end
 ) where {T<:AbstractFourMomentum}
     return size(out_phase_space, 1) == number_outgoing_particles(proc) || throw(
         DimensionMismatch(
-            "The number of outgoing particles <$(number_outgoing_particles(proc))> is inconsistent with input size <$(size(out_phase_space,1))>",
+            "the number of outgoing particles <$(number_outgoing_particles(proc))> is inconsistent with input size <$(size(out_phase_space,1))>",
         ),
     )
 end
@@ -91,7 +91,7 @@ end
 ) where {T<:Real}
     return size(in_phase_space, 1) == in_phase_space_dimension(proc, model) || throw(
         DimensionMismatch(
-            "The dimension of the in-phase-space <$(in_phase_space_dimension(proc,model))> is inconsistent with input size <$(size(in_phase_space,1))>",
+            "the dimension of the in-phase-space <$(in_phase_space_dimension(proc,model))> is inconsistent with input size <$(size(in_phase_space,1))>",
         ),
     )
 end
@@ -103,7 +103,7 @@ end
 ) where {T<:Real}
     return size(out_phase_space, 1) == out_phase_space_dimension(proc, model) || throw(
         DimensionMismatch(
-            "The dimension of the out-phase-space <$(out_phase_space_dimension(proc,model))> is inconsistent with input size <$(size(out_phase_space,1))>",
+            "the dimension of the out-phase-space <$(out_phase_space_dimension(proc,model))> is inconsistent with input size <$(size(out_phase_space,1))>",
         ),
     )
 end
@@ -193,36 +193,36 @@ struct PhaseSpacePoint{
     }
         length(incoming_particles(proc)) == length(in_p) || throw(
             InvalidInputError(
-                "The number of incoming particles given by the process ($(incoming_particles(proc))) mismatches the number of given stateful incoming particles ($(length(in_p)))",
+                "the number of incoming particles given by the process ($(incoming_particles(proc))) mismatches the number of given stateful incoming particles ($(length(in_p)))",
             ),
         )
         length(outgoing_particles(proc)) == length(out_p) || throw(
             InvalidInputError(
-                "The number of outgoing particles given by the process ($(outgoing_particles(proc))) mismatches the number of given stateful outgoing particles ($(length(out_p)))",
+                "the number of outgoing particles given by the process ($(outgoing_particles(proc))) mismatches the number of given stateful outgoing particles ($(length(out_p)))",
             ),
         )
 
         for (proc_p, p) in zip(incoming_particles(proc), in_p)
             proc_p == p.species || throw(
                 InvalidInputError(
-                    "Process given particle species ($(proc_p)) does not match stateful particle species ($(p.species))",
+                    "process given particle species ($(proc_p)) does not match stateful particle species ($(p.species))",
                 ),
             )
             is_incoming(p) || throw(
                 InvalidInputError(
-                    "Stateful particle $(p) is given as an incoming particle but is outgoing",
+                    "stateful particle $(p) is given as an incoming particle but is outgoing",
                 ),
             )
         end
         for (proc_p, p) in zip(outgoing_particles(proc), out_p)
             proc_p == p.species || throw(
                 InvalidInputError(
-                    "Process given particle species ($(proc_p)) does not match stateful particle species ($(p.species))",
+                    "process given particle species ($(proc_p)) does not match stateful particle species ($(p.species))",
                 ),
             )
             is_outgoing(p) || throw(
                 InvalidInputError(
-                    "Stateful particle $(p) is given as an outgoing particle but is incoming",
+                    "stateful particle $(p) is given as an outgoing particle but is incoming",
                 ),
             )
         end
diff --git a/test/phase_spaces.jl b/test/phase_spaces.jl
index 743f443..9bba7a7 100644
--- a/test/phase_spaces.jl
+++ b/test/phase_spaces.jl
@@ -88,19 +88,19 @@ end
 
     if (VERSION >= v"1.8")
         # julia versions before 1.8 did not have support for regex matching in @test_throws
-        @test_throws r"Stateful particle (.*) is given as an incoming particle but is outgoing" PhaseSpacePoint(
+        @test_throws r"stateful particle (.*) is given as an incoming particle but is outgoing" PhaseSpacePoint(
             process, model, phasespace_def, in_particles_invalid, out_particles_valid
         )
 
-        @test_throws r"Stateful particle (.*) is given as an outgoing particle but is incoming" PhaseSpacePoint(
+        @test_throws r"stateful particle (.*) is given as an outgoing particle but is incoming" PhaseSpacePoint(
             process, model, phasespace_def, in_particles_valid, out_particles_invalid
         )
 
-        @test_throws r"Process given particle species \((.*)Electron\(\)\) does not match stateful particle species \((.*)Photon\(\)\)" PhaseSpacePoint(
+        @test_throws r"process given particle species \((.*)Electron\(\)\) does not match stateful particle species \((.*)Photon\(\)\)" PhaseSpacePoint(
             process, model, phasespace_def, SVector(in_ph, in_el), out_particles_valid
         )
 
-        @test_throws r"Process given particle species \((.*)Electron\(\)\) does not match stateful particle species \((.*)Photon\(\)\)" PhaseSpacePoint(
+        @test_throws r"process given particle species \((.*)Electron\(\)\) does not match stateful particle species \((.*)Photon\(\)\)" PhaseSpacePoint(
             process, model, phasespace_def, in_particles_valid, SVector(out_ph, out_el)
         )
     end