Add PhaseSpacePoint definition

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"

src/QEDprocesses.jl

@@ -30,6 +30,7 @@ export propagator
 export AbstractCoordinateSystem, SphericalCoordinateSystem
 export AbstractFrameOfReference, CenterOfMomentumFrame, ElectronRestFrame
 export AbstractPhasespaceDefinition, PhasespaceDefinition
+export ParticleStateful, PhaseSpacePoint
 
 using QEDbase
 

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,127 @@ 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, `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
+    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
+
+"""
+    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,
+    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}}
+
+    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}},
+    }
+        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)
+            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)
+            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)}(
+            proc, model, ps_def, in_p, out_p
+        )
+    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.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)

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"

test/phase_spaces.jl

@@ -0,0 +1,103 @@
+using Random
+using StaticArrays
+using QEDbase
+using QEDprocesses
+
+# can be removed when QEDbase exports them
+import QEDbase.is_incoming, QEDbase.is_outgoing
+
+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(RNG, 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
+            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
+    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_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)
+
+    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 InvalidInputError PhaseSpacePoint(
+        process, model, phasespace_def, in_particles_valid, out_particles_invalid
+    )
+
+    @test_throws InvalidInputError PhaseSpacePoint(
+        process, model, phasespace_def, SVector(in_ph, in_el), out_particles_valid
+    )
+
+    @test_throws InvalidInputError PhaseSpacePoint(
+        process, model, phasespace_def, in_particles_valid, SVector(out_ph, out_el)
+    )
+end

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