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" diff --git a/src/QEDprocesses.jl b/src/QEDprocesses.jl index 3bb38ef..60957e2 100644 --- a/src/QEDprocesses.jl +++ b/src/QEDprocesses.jl @@ -30,6 +30,8 @@ export propagator export AbstractCoordinateSystem, SphericalCoordinateSystem export AbstractFrameOfReference, CenterOfMomentumFrame, ElectronRestFrame export AbstractPhasespaceDefinition, PhasespaceDefinition +export ParticleStateful, PhaseSpacePoint +export spin, polarization, momentum, getindex using QEDbase 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 d35fe36..96a7f13 100644 --- a/src/phase_spaces.jl +++ b/src/phase_spaces.jl @@ -7,6 +7,30 @@ # 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 +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 @@ -32,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 @@ -40,10 +64,10 @@ 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))>", + "the number of incoming particles <$(number_incoming_particles(proc))> is inconsistent with input size <$(size(in_phase_space,1))>", ), ) end @@ -52,10 +76,10 @@ 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))>", + "the number of outgoing particles <$(number_outgoing_particles(proc))> is inconsistent with input size <$(size(out_phase_space,1))>", ), ) end @@ -67,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 @@ -79,7 +103,159 @@ 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 + +""" + ParticleStateful <: AbstractParticle + +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. +- `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`, `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 `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 + species::AbstractParticleType + mom::ElType + spin_or_pol::AbstractSpinOrPolarization + + function ParticleStateful( + 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 + + 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 + +@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 _polarization(::Species, particle::ParticleStateful) where {Species<:BosonLike} = + particle.spin_or_pol +@inline polarization(particle::ParticleStateful) = _polarization(particle.species, particle) + +""" + 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 an `Outgoing` direction, an error is thrown. +""" +struct PhaseSpacePoint{ + PROC<:AbstractProcessDefinition, + MODEL<:AbstractModelDefinition, + PSDEF<:AbstractPhasespaceDefinition, + PhaseSpaceElementType<:AbstractFourMomentum, + 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<:AbstractFourMomentum, + 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.species || throw( + InvalidInputError( + "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", + ), + ) + 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))", + ), + ) + 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 + +""" + Base.getindex(psp::PhaseSpacePoint, dir::Incoming, n::Int) + +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 + +""" + 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 + +""" + 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 momentum(psp::PhaseSpacePoint, dir::ParticleDirection, n::Int) + return psp[dir, n].mom +end 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 new file mode 100644 index 0000000..9bba7a7 --- /dev/null +++ b/test/phase_spaces.jl @@ -0,0 +1,133 @@ +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()] + SPECIES = [Electron(), Positron()] #=, Muon(), AntiMuon(), Tauon(), AntiTauon()=# + SPINANDPOLS = [AllSpin(), SpinUp(), SpinDown(), AllPol(), PolX(), PolY()] + + for (species, dir, spin_or_pol) in Iterators.product(SPECIES, DIRECTIONS, SPINANDPOLS) + momentum = rand(RNG, SFourMomentum) + + 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(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 polarization(particle_stateful) + else + @test polarization(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, species, momentum, spin_or_pol + ) + end + @test_throws MethodError ParticleStateful(dir, species, momentum, spin_or_pol) + end + end +end + +@testset "Phasespace Point" begin + 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) + + 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 + + psp = PhaseSpacePoint( + process, model, phasespace_def, in_particles_valid, out_particles_valid + ) + + @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 + @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 + @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 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( + process, model, phasespace_def, in_particles_valid, SVector(out_ph, out_el) + ) + end + + @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] + @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 + ) + + @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 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