diff --git a/docs/Project.toml b/docs/Project.toml index 2292a04..5ba0196 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -1,3 +1,5 @@ [deps] Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +QEDbase = "10e22c08-3ccb-4172-bfcf-7d7aa3d04d93" QEDcore = "35dc0263-cb5f-4c33-a114-1d7f54ab753e" +QEDprocesses = "46de9c38-1bb3-4547-a1ec-da24d767fdad" diff --git a/docs/make.jl b/docs/make.jl index 29867cd..979a5b0 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -1,7 +1,16 @@ +using Pkg + +# targeting the correct source code +# this asumes the make.jl script is located in QEDcore.jl/docs +project_path = Base.Filesystem.joinpath(Base.Filesystem.dirname(Base.source_path()), "..") +Pkg.develop(; path=project_path) + +Pkg.add(; url="https://github.com/QEDjl-project/QEDbase.jl/", rev="dev") + using QEDcore using Documenter -DocMeta.setdocmeta!(QEDcore, :DocTestSetup, :(using QEDcore); recursive=true) +# DocMeta.setdocmeta!(QEDcore, :DocTestSetup, :(using QEDcore); recursive=true) makedocs(; modules=[QEDcore], diff --git a/src/QEDcore.jl b/src/QEDcore.jl index 1503a5a..352b5c8 100644 --- a/src/QEDcore.jl +++ b/src/QEDcore.jl @@ -1,9 +1,5 @@ module QEDcore -import Base: * -import StaticArrays: similar_type -import QEDbase: base_state - # lorentz vectors export SLorentzVector, MLorentzVector @@ -32,6 +28,13 @@ export @valid_spinor_input # particle base states export base_state +# phase space +export SphericalCoordinateSystem +export CenterOfMomentumFrame, ElectronRestFrame +export PhasespaceDefinition +export ParticleStateful, PhaseSpacePoint, InPhaseSpacePoint, OutPhaseSpacePoint +export spin, polarization, particle_direction, particle_species, momentum, momenta, getindex + using QEDbase: QEDbase using DocStringExtensions using StaticArrays @@ -40,12 +43,19 @@ using SimpleTraits include("algebraic_objects/dirac_tensors/types.jl") include("algebraic_objects/dirac_tensors/multiplication.jl") +include("phase_spaces/types.jl") +include("phase_spaces/access.jl") +include("phase_spaces/create.jl") +include("phase_spaces/print.jl") +include("phase_spaces/utility.jl") + include("algebraic_objects/four_momentum.jl") include("algebraic_objects/lorentz_vector.jl") include("algebraic_objects/gamma_matrices.jl") include("particles/particle_types.jl") -include("particles/spinors.jl") +include("particles/propagators.jl") include("particles/states.jl") +include("particles/spinors.jl") end diff --git a/src/algebraic_objects/dirac_tensors/multiplication.jl b/src/algebraic_objects/dirac_tensors/multiplication.jl index 29d9e64..21c8e84 100644 --- a/src/algebraic_objects/dirac_tensors/multiplication.jl +++ b/src/algebraic_objects/dirac_tensors/multiplication.jl @@ -3,6 +3,9 @@ # Concrete implementation of multiplication for Dirac Tensors # ####### + +import Base: * + """ $(TYPEDSIGNATURES) diff --git a/src/algebraic_objects/four_momentum.jl b/src/algebraic_objects/four_momentum.jl index a3da97d..f39ddfe 100644 --- a/src/algebraic_objects/four_momentum.jl +++ b/src/algebraic_objects/four_momentum.jl @@ -5,6 +5,7 @@ ####### import QEDbase: getT, getX, getY, getZ, setT!, setX!, setY!, setZ! +import StaticArrays: similar_type """ $(TYPEDEF) diff --git a/src/algebraic_objects/lorentz_vector.jl b/src/algebraic_objects/lorentz_vector.jl index a0cfcfc..f83b422 100644 --- a/src/algebraic_objects/lorentz_vector.jl +++ b/src/algebraic_objects/lorentz_vector.jl @@ -5,6 +5,7 @@ ####### import QEDbase: getT, getX, getY, getZ, setT!, setX!, setY!, setZ! +import StaticArrays: similar_type ####### # diff --git a/src/particles/particle_types.jl b/src/particles/particle_types.jl index 2ef3c72..be51684 100644 --- a/src/particles/particle_types.jl +++ b/src/particles/particle_types.jl @@ -81,7 +81,7 @@ is_anti_particle(::MajoranaFermion) = true Concrete type for *electrons* as a particle species. Mostly used for dispatch. ```jldoctest -julia> using QEDbase +julia> using QEDcore julia> Electron() electron @@ -104,7 +104,7 @@ Base.show(io::IO, ::Electron) = print(io, "electron") Concrete type for *positrons* as a particle species. Mostly used for dispatch. ```jldoctest -julia> using QEDbase +julia> using QEDcore julia> Positron() positron @@ -182,7 +182,7 @@ is_anti_particle(::MajoranaBoson) = true Concrete type for the *photons* as a particle species. Mostly used for dispatch. ```jldoctest -julia> using QEDbase +julia> using QEDcore julia> Photon() photon diff --git a/src/particles/propagators.jl b/src/particles/propagators.jl new file mode 100644 index 0000000..1b69de5 --- /dev/null +++ b/src/particles/propagators.jl @@ -0,0 +1,29 @@ +import QEDbase: propagator + +function _scalar_propagator(K::QEDbase.AbstractFourMomentum, mass::Real) + return one(mass) / (K * K - mass^2) +end + +function _scalar_propagator(K::QEDbase.AbstractFourMomentum) + return one(getT(K)) / (K * K) +end + +function _fermion_propagator(P::QEDbase.AbstractFourMomentum, mass::Real) + return (slashed(P) + mass * one(DiracMatrix)) * _scalar_propagator(P, mass) +end + +function _fermion_propagator(P::QEDbase.AbstractFourMomentum) + return (slashed(P)) * _scalar_propagator(P) +end + +function QEDbase.propagator(particle_type::BosonLike, K::QEDbase.AbstractFourMomentum) + return _scalar_propagator(K, QEDbase.mass(particle_type)) +end + +function QEDbase.propagator(particle_type::Photon, K::QEDbase.AbstractFourMomentum) + return _scalar_propagator(K) +end + +function QEDbase.propagator(particle_type::FermionLike, P::QEDbase.AbstractFourMomentum) + return _fermion_propagator(P, QEDbase.mass(particle_type)) +end diff --git a/src/particles/states.jl b/src/particles/states.jl index a575b70..54affde 100644 --- a/src/particles/states.jl +++ b/src/particles/states.jl @@ -1,3 +1,5 @@ +import QEDbase: base_state + function _booster_fermion(mom::QEDbase.AbstractFourMomentum, mass::Real) return (slashed(mom) + mass * one(DiracMatrix)) / (sqrt(abs(QEDbase.getT(mom)) + mass)) end diff --git a/src/phase_spaces/access.jl b/src/phase_spaces/access.jl new file mode 100644 index 0000000..cc51c66 --- /dev/null +++ b/src/phase_spaces/access.jl @@ -0,0 +1,30 @@ +# accessor interface particle stateful +QEDbase.particle_direction(part::ParticleStateful) = part.dir +QEDbase.particle_species(part::ParticleStateful) = part.species +QEDbase.momentum(part::ParticleStateful) = part.mom + +# accessor interface phase space point +""" + 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, ::QEDbase.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, ::QEDbase.Outgoing, n::Int) + return psp.out_particles[n] +end + +QEDbase.process(psp::PhaseSpacePoint) = psp.proc +QEDbase.model(psp::PhaseSpacePoint) = psp.model +QEDbase.phase_space_definition(psp::PhaseSpacePoint) = psp.ps_def + +QEDbase.particles(psp::PhaseSpacePoint, ::QEDbase.Incoming) = psp.in_particles +QEDbase.particles(psp::PhaseSpacePoint, ::QEDbase.Outgoing) = psp.out_particles diff --git a/src/phase_spaces/create.jl b/src/phase_spaces/create.jl new file mode 100644 index 0000000..5f8f701 --- /dev/null +++ b/src/phase_spaces/create.jl @@ -0,0 +1,158 @@ +# PSP constructors from particle statefuls + +""" + InPhaseSpacePoint( + proc::QEDbase.AbstractProcessDefinition, + model::QEDbase.AbstractModelDefinition, + ps_def::QEDbase.AbstractPhasespaceDefinition, + in_ps::Tuple{ParticleStateful}, + ) + + Construct a [`PhaseSpacePoint`](@ref) with only input particles from [`ParticleStateful`](@ref)s. The result will be `<: InPhaseSpacePoint` but **not** `<: OutPhaseSpacePoint`. +""" +function InPhaseSpacePoint( + proc::PROC, model::MODEL, ps_def::PSDEF, in_ps::IN_PARTICLES +) where { + PROC<:QEDbase.AbstractProcessDefinition, + MODEL<:QEDbase.AbstractModelDefinition, + PSDEF<:QEDbase.AbstractPhasespaceDefinition, + IN_PARTICLES<:Tuple{Vararg{ParticleStateful}}, +} + return PhaseSpacePoint(proc, model, ps_def, in_ps, ()) +end + +""" + OutPhaseSpacePoint( + proc::QEDbase.AbstractProcessDefinition, + model::QEDbase.AbstractModelDefinition, + ps_def::QEDbase.AbstractPhasespaceDefinition, + out_ps::Tuple{ParticleStateful}, + ) + +Construct a [`PhaseSpacePoint`](@ref) with only output particles from [`ParticleStateful`](@ref)s. The result will be `<: OutPhaseSpacePoint` but **not** `<: InPhaseSpacePoint`. +""" +function OutPhaseSpacePoint( + proc::PROC, model::MODEL, ps_def::PSDEF, out_ps::OUT_PARTICLES +) where { + PROC<:QEDbase.AbstractProcessDefinition, + MODEL<:QEDbase.AbstractModelDefinition, + PSDEF<:QEDbase.AbstractPhasespaceDefinition, + OUT_PARTICLES<:Tuple{Vararg{ParticleStateful}}, +} + return PhaseSpacePoint(proc, model, ps_def, (), out_ps) +end + +# PSP constructors from momenta + +""" + PhaseSpacePoint( + proc::QEDbase.AbstractProcessDefinition, + model::QEDbase.AbstractModelDefinition, + ps_def::QEDbase.AbstractPhasespaceDefinition, + in_momenta::NTuple{N,QEDbase.AbstractFourMomentum}, + out_momenta::NTuple{M,QEDbase.AbstractFourMomentum}, + ) + +Construct the phase space point from given momenta of incoming and outgoing particles regarding a given process. +""" +function PhaseSpacePoint( + proc::QEDbase.AbstractProcessDefinition, + model::QEDbase.AbstractModelDefinition, + ps_def::QEDbase.AbstractPhasespaceDefinition, + in_momenta::NTuple{N,ELEMENT}, + out_momenta::NTuple{M,ELEMENT}, +) where {N,M,ELEMENT<:QEDbase.AbstractFourMomentum} + in_particles = _build_particle_statefuls(proc, in_momenta, QEDbase.Incoming()) + out_particles = _build_particle_statefuls(proc, out_momenta, QEDbase.Outgoing()) + + return PhaseSpacePoint(proc, model, ps_def, in_particles, out_particles) +end + +""" + InPhaseSpacePoint( + proc::QEDbase.AbstractProcessDefinition, + model::QEDbase.AbstractModelDefinition, + ps_def::QEDbase.AbstractPhasespaceDefinition, + in_momenta::NTuple{N,QEDbase.AbstractFourMomentum}, + ) + +Construct a [`PhaseSpacePoint`](@ref) with only input particles from given momenta. The result will be `<: InPhaseSpacePoint` but **not** `<: OutPhaseSpacePoint`. +""" +function InPhaseSpacePoint( + proc::QEDbase.AbstractProcessDefinition, + model::QEDbase.AbstractModelDefinition, + ps_def::QEDbase.AbstractPhasespaceDefinition, + in_momenta::NTuple{N,ELEMENT}, +) where {N,ELEMENT<:QEDbase.AbstractFourMomentum} + in_particles = _build_particle_statefuls(proc, in_momenta, QEDbase.Incoming()) + + return PhaseSpacePoint(proc, model, ps_def, in_particles, ()) +end + +""" + OutPhaseSpacePoint( + proc::QEDbase.AbstractProcessDefinition, + model::QEDbase.AbstractModelDefinition, + ps_def::QEDbase.AbstractPhasespaceDefinition, + out_momenta::NTuple{N,QEDbase.AbstractFourMomentum}, + ) + +Construct a [`PhaseSpacePoint`](@ref) with only output particles from given momenta. The result will be `<: OutPhaseSpacePoint` but **not** `<: InPhaseSpacePoint`. +""" +function OutPhaseSpacePoint( + proc::QEDbase.AbstractProcessDefinition, + model::QEDbase.AbstractModelDefinition, + ps_def::QEDbase.AbstractPhasespaceDefinition, + out_momenta::NTuple{N,ELEMENT}, +) where {N,ELEMENT<:QEDbase.AbstractFourMomentum} + out_particles = _build_particle_statefuls(proc, out_momenta, QEDbase.Outgoing()) + + return PhaseSpacePoint(proc, model, ps_def, (), out_particles) +end + +# PSP constructors from coordinates + +""" + PhaseSpacePoint( + proc::QEDbase.AbstractProcessDefinition, + model::QEDbase.AbstractModelDefinition, + ps_def::QEDbase.AbstractPhasespaceDefinition, + in_coords::NTuple{N,Real}, + out_coords::NTuple{M,Real}, + ) + +Construct a [`PhaseSpacePoint`](@ref) from given coordinates by using the `QEDbase._generate_momenta` interface. +""" +function PhaseSpacePoint( + proc::QEDbase.AbstractProcessDefinition, + model::QEDbase.AbstractModelDefinition, + ps_def::QEDbase.AbstractPhasespaceDefinition, + in_coords::NTuple{N,Real}, + out_coords::NTuple{M,Real}, +) where {N,M} + in_ps, out_ps = _generate_momenta(proc, model, ps_def, in_coords, out_coords) + return PhaseSpacePoint(proc, model, ps_def, in_ps, out_ps) +end + +""" + InPhaseSpacePoint( + proc::QEDbase.AbstractProcessDefinition, + model::QEDbase.AbstractModelDefinition, + ps_def::QEDbase.AbstractPhasespaceDefinition, + in_coords::NTuple{N,Real}, + ) + +Construct a [`PhaseSpacePoint`](@ref) from given coordinates by using the `QEDbase._generate_momenta` interface. The result will be `<: InPhaseSpacePoint` but **not** `<: OutPhaseSpacePoint`. + +!!! note + A similar function for [`OutPhaseSpacePoint`](@ref) does not exist from coordinates, only a full [`PhaseSpacePoint`](@ref). +""" +function InPhaseSpacePoint( + proc::QEDbase.AbstractProcessDefinition, + model::QEDbase.AbstractModelDefinition, + ps_def::QEDbase.AbstractPhasespaceDefinition, + in_coords::NTuple{N,Real}, +) where {N} + in_ps = _generate_incoming_momenta(proc, model, ps_def, in_coords) + return InPhaseSpacePoint(proc, model, ps_def, in_ps) +end diff --git a/src/phase_spaces/print.jl b/src/phase_spaces/print.jl new file mode 100644 index 0000000..e362218 --- /dev/null +++ b/src/phase_spaces/print.jl @@ -0,0 +1,58 @@ +function Base.show(io::IO, ::SphericalCoordinateSystem) + print(io, "spherical coordinates") + return nothing +end + +function Base.show(io::IO, ::CenterOfMomentumFrame) + print(io, "center-of-momentum frame") + return nothing +end + +function Base.show(io::IO, ::ElectronRestFrame) + print(io, "electron rest frame") + return nothing +end + +function Base.show(io::IO, m::MIME"text/plain", ps_def::PhasespaceDefinition) + println(io, "PhasespaceDefinition") + println(io, " coordinate system: $(ps_def.coord_sys)") + println(io, " frame: $(ps_def.frame)") + return nothing +end + +function Base.show(io::IO, ps_def::PhasespaceDefinition) + print(io, "$(ps_def.coord_sys) in $(ps_def.frame)") + return nothing +end + +function Base.show(io::IO, particle::ParticleStateful) + print(io, "$(particle.dir) $(particle.species): $(particle.mom)") + return nothing +end + +function Base.show(io::IO, m::MIME"text/plain", particle::ParticleStateful) + println(io, "ParticleStateful: $(particle.dir) $(particle.species)") + println(io, " momentum: $(particle.mom)") + return nothing +end + +function Base.show(io::IO, psp::PhaseSpacePoint) + print(io, "PhaseSpacePoint of $(psp.proc)") + return nothing +end + +function Base.show(io::IO, ::MIME"text/plain", psp::PhaseSpacePoint) + println(io, "PhaseSpacePoint:") + println(io, " process: $(psp.proc)") + println(io, " model: $(psp.model)") + println(io, " phasespace definition: $(psp.ps_def)") + println(io, " incoming particles:") + for p in psp.in_particles + println(io, " -> $(p)") + end + println(io, " outgoing particles:") + for p in psp.out_particles + println(io, " -> $(p)") + end + return nothing +end diff --git a/src/phase_spaces/types.jl b/src/phase_spaces/types.jl new file mode 100644 index 0000000..7f3fcb7 --- /dev/null +++ b/src/phase_spaces/types.jl @@ -0,0 +1,184 @@ +""" + SphericalCoordinateSystem <: AbstractCoordinateSystem + +TBW +""" +struct SphericalCoordinateSystem <: QEDbase.AbstractCoordinateSystem end + +""" + CenterOfMomentumFrame <: AbstractFrameOfReference + +TBW +""" +struct CenterOfMomentumFrame <: QEDbase.AbstractFrameOfReference end + +""" + ElectronRestFrame <: AbstractFrameOfReference + +TBW +""" +struct ElectronRestFrame <: QEDbase.AbstractFrameOfReference end + +""" + PhasespaceDefinition(coord_sys::QEDbase.AbstractCoordinateSystem, frame::QEDbase.AbstractFrameOfReference) + +Convenient type to dispatch on coordiante systems and frames of reference. Combines a `QEDbase.AbstractCoordinateSystem` with a `QEDbase.AbstractFrameOfReference`. +""" +struct PhasespaceDefinition{ + CS<:QEDbase.AbstractCoordinateSystem,F<:QEDbase.AbstractFrameOfReference +} <: QEDbase.AbstractPhasespaceDefinition + coord_sys::CS + frame::F +end + +""" + ParticleStateful <: QEDbase.AbstractParticle + +Representation of a particle with a state. It has four fields: +- `dir::QEDbase.ParticleDirection`: The direction of the particle, `QEDbase.Incoming()` or `QEDbase.Outgoing()`. +- `species::QEDbase.AbstractParticleType`: The species of the particle, `Electron()`, `Positron()` etc. +- `mom::QEDbase.AbstractFourMomentum`: The momentum of the particle. + +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. + +TODO: Turn this back into a `jldoctest` once refactoring is done. +```Julia +julia> using QEDcore; using QEDbase + +julia> ParticleStateful(Incoming(), Electron(), SFourMomentum(1, 0, 0, 0)) +ParticleStateful: incoming electron + momentum: [1.0, 0.0, 0.0, 0.0] + +julia> ParticleStateful(Outgoing(), Photon(), SFourMomentum(1, 0, 0, 0)) +ParticleStateful: outgoing photon + momentum: [1.0, 0.0, 0.0, 0.0] +``` +""" +struct ParticleStateful{ + DIR<:QEDbase.ParticleDirection, + SPECIES<:QEDbase.AbstractParticleType, + ELEMENT<:QEDbase.AbstractFourMomentum, +} <: QEDbase.AbstractParticleStateful{DIR,SPECIES,ELEMENT} + dir::DIR + species::SPECIES + mom::ELEMENT + + function ParticleStateful( + dir::DIR, species::SPECIES, mom::ELEMENT + ) where { + DIR<:QEDbase.ParticleDirection, + SPECIES<:QEDbase.AbstractParticleType, + ELEMENT<:QEDbase.AbstractFourMomentum, + } + return new{DIR,SPECIES,ELEMENT}(dir, species, mom) + end +end + +""" + PhaseSpacePoint + +Representation of a point in the phase space of a process. Contains the process (`QEDbase.AbstractProcessDefinition`), the model (`QEDbase.AbstractModelDefinition`), the phase space definition (`QEDbase.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 `QEDbase.Outgoing` direction, an error is thrown. + +TODO: Turn this back into a `jldoctest` once refactoring is done. +```Julia +julia> using QEDcore; using QEDbase; using QEDprocesses + +julia> PhaseSpacePoint( + Compton(), + PerturbativeQED(), + PhasespaceDefinition(SphericalCoordinateSystem(), ElectronRestFrame()), + ( + ParticleStateful(Incoming(), Electron(), SFourMomentum(1, 0, 0, 0)), + ParticleStateful(Incoming(), Photon(), SFourMomentum(1, 0, 0, 0)) + ), + ( + ParticleStateful(Outgoing(), Electron(), SFourMomentum(1, 0, 0, 0)), + ParticleStateful(Outgoing(), Photon(), SFourMomentum(1, 0, 0, 0)) + ) + ) +PhaseSpacePoint: + process: one-photon Compton scattering + model: perturbative QED + phasespace definition: spherical coordinates in electron rest frame + incoming particles: + -> incoming electron: [1.0, 0.0, 0.0, 0.0] + -> incoming photon: [1.0, 0.0, 0.0, 0.0] + outgoing particles: + -> outgoing electron: [1.0, 0.0, 0.0, 0.0] + -> outgoing photon: [1.0, 0.0, 0.0, 0.0] +``` + +!!! note + `PhaseSpacePoint`s can be constructed with only one of their in- or out-channel set. For this, see the special constructors [`InPhaseSpacePoint`](@ref) and [`OutPhaseSpacePoint`](@ref). + The [`InPhaseSpacePoint`](@ref) and [`OutPhaseSpacePoint`](@ref) type definitions can be used to dispatch on such `PhaseSpacePoint`s. Note that a full `PhaseSpacePoint` containing both its in- and out-channel matches both, .i.e. `psp isa InPhaseSpacePoint` and `psp isa OutPhaseSpacePoint` both evaluate to true if psp contains both channels. + A completely empty `PhaseSpacePoint` is not allowed. +""" +struct PhaseSpacePoint{ + PROC<:QEDbase.AbstractProcessDefinition, + MODEL<:QEDbase.AbstractModelDefinition, + PSDEF<:QEDbase.AbstractPhasespaceDefinition, + IN_PARTICLES<:Tuple{Vararg{ParticleStateful}}, + OUT_PARTICLES<:Tuple{Vararg{ParticleStateful}}, + ELEMENT<:QEDbase.AbstractFourMomentum, +} <: QEDbase.AbstractPhaseSpacePoint{PROC,MODEL,PSDEF,IN_PARTICLES,OUT_PARTICLES} + proc::PROC + model::MODEL + ps_def::PSDEF + + in_particles::IN_PARTICLES + out_particles::OUT_PARTICLES + + """ + PhaseSpacePoint( + proc::QEDbase.AbstractProcessDefinition, + model::QEDbase.AbstractModelDefinition, + ps_def::QEDbase.AbstractPhasespaceDefinition, + in_ps::Tuple{ParticleStateful}, + out_ps::Tuple{ParticleStateful}, + ) + + Construct a [`PhaseSpacePoint`](@ref) from a process, model, phasespace definition and a tuple of [`ParticleStateful`](@ref)s. + """ + function PhaseSpacePoint( + proc::PROC, model::MODEL, ps_def::PSDEF, in_p::IN_PARTICLES, out_p::OUT_PARTICLES + ) where { + PROC<:QEDbase.AbstractProcessDefinition, + MODEL<:QEDbase.AbstractModelDefinition, + PSDEF<:QEDbase.AbstractPhasespaceDefinition, + IN_PARTICLES<:Tuple{Vararg{ParticleStateful}}, + OUT_PARTICLES<:Tuple{Vararg{ParticleStateful}}, + } + # this entire check is compiled away every time, so there's no need to disable it for performance ever + ELEMENT = _check_psp( + QEDbase.incoming_particles(proc), QEDbase.outgoing_particles(proc), in_p, out_p + ) + + return new{PROC,MODEL,PSDEF,IN_PARTICLES,OUT_PARTICLES,ELEMENT}( + proc, model, ps_def, in_p, out_p + ) + end +end + +""" + InPhaseSpacePoint + +A partial type specialization on [`PhaseSpacePoint`](@ref) which can be used for dispatch in functions requiring only the in channel of the phase space to exist, for example implementations of `QEDbase._incident_flux`. No restrictions are imposed on the out-channel, which may or may not exist. + +See also: [`OutPhaseSpacePoint`](@ref) +""" +InPhaseSpacePoint{P,M,D,IN,OUT,E} = PhaseSpacePoint{ + P,M,D,IN,OUT,E +} where {IN<:Tuple{ParticleStateful,Vararg},OUT<:Tuple{Vararg}} + +""" + OutPhaseSpacePoint + +A partial type specialization on [`PhaseSpacePoint`](@ref) which can be used for dispatch in functions requiring only the out channel of the phase space to exist. No restrictions are imposed on the in-channel, which may or may not exist. + +See also: [`InPhaseSpacePoint`](@ref) +""" +OutPhaseSpacePoint{P,M,D,IN,OUT,E} = PhaseSpacePoint{ + P,M,D,IN,OUT,E +} where {IN<:Tuple{Vararg},OUT<:Tuple{ParticleStateful,Vararg}} diff --git a/src/phase_spaces/utility.jl b/src/phase_spaces/utility.jl new file mode 100644 index 0000000..02fb4d3 --- /dev/null +++ b/src/phase_spaces/utility.jl @@ -0,0 +1,150 @@ +# recursion termination: base case +@inline _assemble_tuple_type(::Tuple{}, ::QEDbase.ParticleDirection, ::Type) = () + +# function assembling the correct type information for the tuple of ParticleStatefuls in a phasespace point constructed from momenta +@inline function _assemble_tuple_type( + particle_types::Tuple{SPECIES_T,Vararg{QEDbase.AbstractParticleType}}, + dir::DIR_T, + ELTYPE::Type, +) where {SPECIES_T<:QEDbase.AbstractParticleType,DIR_T<:QEDbase.ParticleDirection} + return ( + ParticleStateful{DIR_T,SPECIES_T,ELTYPE}, + _assemble_tuple_type(particle_types[2:end], dir, ELTYPE)..., + ) +end + +# recursion termination: success +@inline _recursive_type_check(::Tuple{}, ::Tuple{}, ::QEDbase.ParticleDirection) = nothing + +# recursion termination: overload for unequal number of particles +@inline function _recursive_type_check( + ::Tuple{Vararg{ParticleStateful,N}}, + ::Tuple{Vararg{QEDbase.AbstractParticleType,M}}, + dir::QEDbase.ParticleDirection, +) where {N,M} + throw( + QEDbase.InvalidInputError( + "expected $(M) $(dir) particles for the process but got $(N)" + ), + ) + return nothing +end + +# recursion termination: overload for invalid types +@inline function _recursive_type_check( + ::Tuple{ParticleStateful{DIR_IN_T,SPECIES_IN_T},Vararg{ParticleStateful,N}}, + ::Tuple{SPECIES_T,Vararg{QEDbase.AbstractParticleType,N}}, + dir::DIR_T, +) where { + N, + DIR_IN_T<:QEDbase.ParticleDirection, + DIR_T<:QEDbase.ParticleDirection, + SPECIES_IN_T<:QEDbase.AbstractParticleType, + SPECIES_T<:QEDbase.AbstractParticleType, +} + throw( + QEDbase.InvalidInputError( + "expected $(dir) $(SPECIES_T()) but got $(DIR_IN_T()) $(SPECIES_IN_T())" + ), + ) + return nothing +end + +@inline function _recursive_type_check( + t::Tuple{ParticleStateful{DIR_T,SPECIES_T},Vararg{ParticleStateful,N}}, + p::Tuple{SPECIES_T,Vararg{QEDbase.AbstractParticleType,N}}, + dir::DIR_T, +) where {N,DIR_T<:QEDbase.ParticleDirection,SPECIES_T<:QEDbase.AbstractParticleType} + return _recursive_type_check(t[2:end], p[2:end], dir) +end + +@inline function _check_psp( + in_proc::P_IN_Ts, out_proc::P_OUT_Ts, in_p::IN_Ts, out_p::OUT_Ts +) where { + P_IN_Ts<:Tuple{Vararg{QEDbase.AbstractParticleType}}, + P_OUT_Ts<:Tuple{Vararg{QEDbase.AbstractParticleType}}, + IN_Ts<:Tuple{Vararg{ParticleStateful}}, + OUT_Ts<:Tuple{}, +} + # specific overload for InPhaseSpacePoint + _recursive_type_check(in_p, in_proc, QEDbase.Incoming()) + + return typeof(in_p[1].mom) +end + +@inline function _check_psp( + in_proc::P_IN_Ts, out_proc::P_OUT_Ts, in_p::IN_Ts, out_p::OUT_Ts +) where { + P_IN_Ts<:Tuple{Vararg{QEDbase.AbstractParticleType}}, + P_OUT_Ts<:Tuple{Vararg{QEDbase.AbstractParticleType}}, + IN_Ts<:Tuple{}, + OUT_Ts<:Tuple{Vararg{ParticleStateful}}, +} + # specific overload for OutPhaseSpacePoint + _recursive_type_check(out_p, out_proc, QEDbase.Outgoing()) + + return typeof(out_p[1].mom) +end + +@inline function _check_psp( + in_proc::P_IN_Ts, out_proc::P_OUT_Ts, in_p::IN_Ts, out_p::OUT_Ts +) where { + P_IN_Ts<:Tuple{Vararg{QEDbase.AbstractParticleType}}, + P_OUT_Ts<:Tuple{Vararg{QEDbase.AbstractParticleType}}, + IN_Ts<:Tuple{Vararg{ParticleStateful}}, + OUT_Ts<:Tuple{Vararg{ParticleStateful}}, +} + # in_proc/out_proc contain only species types + # in_p/out_p contain full ParticleStateful types + + _recursive_type_check(in_p, in_proc, QEDbase.Incoming()) + _recursive_type_check(out_p, out_proc, QEDbase.Outgoing()) + + return typeof(out_p[1].mom) +end + +""" + _momentum_type(psp::PhaseSpacePoint) + _momentum_type(type::Type{PhaseSpacePoint}) + +Returns the element type of the [`PhaseSpacePoint`](@ref) object or type, e.g. `SFourMomentum`. + +TODO: Turn this back into a `jldoctest` once refactoring is done. +```Julia +julia> using QEDcore; using QEDprocesses + +julia> psp = PhaseSpacePoint(Compton(), PerturbativeQED(), PhasespaceDefinition(SphericalCoordinateSystem(), ElectronRestFrame()), Tuple(rand(SFourMomentum) for _ in 1:2), Tuple(rand(SFourMomentum) for _ in 1:2)); + +julia> _momentum_type(psp) +SFourMomentum + +julia> _momentum_type(typeof(psp)) +SFourMomentum +``` +""" +@inline function _momentum_type( + ::Type{T} +) where {P,M,D,I,O,E,T<:PhaseSpacePoint{P,M,D,I,O,E}} + return E +end + +@inline _momentum_type(::T) where {T<:PhaseSpacePoint} = _momentum_type(T) + +# convenience function building a type stable tuple of ParticleStatefuls from the given process, momenta, and direction +function _build_particle_statefuls( + proc::QEDbase.AbstractProcessDefinition, + moms::NTuple{N,ELEMENT}, + dir::QEDbase.ParticleDirection, +) where {N,ELEMENT<:QEDbase.AbstractFourMomentum} + N == QEDbase.number_particles(proc, dir) || throw( + QEDbase.InvalidInputError( + "expected $(QEDbase.number_particles(proc, dir)) $(dir) particles for the process but got $(N)", + ), + ) + res = Tuple{_assemble_tuple_type(QEDbase.particles(proc, dir), dir, ELEMENT)...}( + ParticleStateful(dir, particle, mom) for + (particle, mom) in zip(QEDbase.particles(proc, dir), moms) + ) + + return res +end diff --git a/test/interfaces/process.jl b/test/interfaces/process.jl new file mode 100644 index 0000000..fae7118 --- /dev/null +++ b/test/interfaces/process.jl @@ -0,0 +1,130 @@ +using Random +using QEDbase: QEDbase +using QEDcore + +RNG = MersenneTwister(137137) +ATOL = 0.0 +RTOL = sqrt(eps()) + +include("../test_implementation/TestImplementation.jl") + +@testset "($N_INCOMING,$N_OUTGOING)" for (N_INCOMING, N_OUTGOING) in Iterators.product( + (1, rand(RNG, 2:8)), (1, rand(RNG, 2:8)) +) + INCOMING_PARTICLES = Tuple(rand(RNG, TestImplementation.PARTICLE_SET, N_INCOMING)) + OUTGOING_PARTICLES = Tuple(rand(RNG, TestImplementation.PARTICLE_SET, N_OUTGOING)) + + TESTPROC = TestImplementation.TestProcess(INCOMING_PARTICLES, OUTGOING_PARTICLES) + TESTMODEL = TestImplementation.TestModel() + TESTPSDEF = TestImplementation.TestPhasespaceDef() + IN_PS = TestImplementation._rand_momenta(RNG, N_INCOMING) + OUT_PS = TestImplementation._rand_momenta(RNG, N_OUTGOING) + PSP = PhaseSpacePoint(TESTPROC, TESTMODEL, TESTPSDEF, IN_PS, OUT_PS) + + @testset "failed interface" begin + TESTPROC_FAIL_ALL = TestImplementation.TestProcess_FAIL_ALL( + INCOMING_PARTICLES, OUTGOING_PARTICLES + ) + TESTPROC_FAIL_DIFFCS = TestImplementation.TestProcess_FAIL_DIFFCS( + INCOMING_PARTICLES, OUTGOING_PARTICLES + ) + TESTMODEL_FAIL = TestImplementation.TestModel_FAIL() + TESTPSDEF_FAIL = TestImplementation.TestPhasespaceDef_FAIL() + + @testset "failed process interface" begin + @test_throws MethodError QEDbase.incoming_particles(TESTPROC_FAIL_ALL) + @test_throws MethodError QEDbase.outgoing_particles(TESTPROC_FAIL_ALL) + end + + @testset "$PROC $MODEL" for (PROC, MODEL) in Iterators.product( + (TESTPROC, TESTPROC_FAIL_DIFFCS), (TESTMODEL, TESTMODEL_FAIL) + ) + if TestImplementation._any_fail(PROC, MODEL) + psp = PhaseSpacePoint(PROC, MODEL, TESTPSDEF, IN_PS, OUT_PS) + @test_throws MethodError QEDbase._incident_flux(psp) + @test_throws MethodError QEDbase._averaging_norm(psp) + @test_throws MethodError QEDbase._matrix_element(psp) + end + + for PS_DEF in (TESTPSDEF, TESTPSDEF_FAIL) + if TestImplementation._any_fail(PROC, MODEL, PS_DEF) + psp = PhaseSpacePoint(PROC, MODEL, PS_DEF, IN_PS, OUT_PS) + @test_throws MethodError QEDbase._phase_space_factor(psp) + end + end + end + end + + @testset "broadcast" begin + test_func(proc::QEDbase.AbstractProcessDefinition) = proc + @test test_func.(TESTPROC) == TESTPROC + + test_func(model::QEDbase.AbstractModelDefinition) = model + @test test_func.(TESTMODEL) == TESTMODEL + end + + @testset "incoming/outgoing particles" begin + @test QEDbase.incoming_particles(TESTPROC) == INCOMING_PARTICLES + @test QEDbase.outgoing_particles(TESTPROC) == OUTGOING_PARTICLES + @test QEDbase.number_incoming_particles(TESTPROC) == N_INCOMING + @test QEDbase.number_outgoing_particles(TESTPROC) == N_OUTGOING + end + + @testset "incident flux" begin + test_incident_flux = QEDbase._incident_flux( + InPhaseSpacePoint(TESTPROC, TESTMODEL, TESTPSDEF, IN_PS) + ) + groundtruth = TestImplementation._groundtruth_incident_flux(IN_PS) + @test isapprox(test_incident_flux, groundtruth, atol=ATOL, rtol=RTOL) + + test_incident_flux = QEDbase._incident_flux( + PhaseSpacePoint(TESTPROC, TESTMODEL, TESTPSDEF, IN_PS, OUT_PS) + ) + @test isapprox(test_incident_flux, groundtruth, atol=ATOL, rtol=RTOL) + + @test_throws MethodError QEDbase._incident_flux( + OutPhaseSpacePoint(TESTPROC, TESTMODEL, TESTPSDEF, OUT_PS) + ) + end + + @testset "averaging norm" begin + test_avg_norm = QEDbase._averaging_norm(TESTPROC) + groundtruth = TestImplementation._groundtruth_averaging_norm(TESTPROC) + @test isapprox(test_avg_norm, groundtruth, atol=ATOL, rtol=RTOL) + end + + @testset "matrix element" begin + test_matrix_element = QEDbase._matrix_element(PSP) + groundtruth = TestImplementation._groundtruth_matrix_element(IN_PS, OUT_PS) + @test length(test_matrix_element) == length(groundtruth) + for i in eachindex(test_matrix_element) + @test isapprox(test_matrix_element[i], groundtruth[i], atol=ATOL, rtol=RTOL) + end + end + + @testset "is in phasespace" begin + @test QEDbase._is_in_phasespace(PSP) + + IN_PS_unphysical = (zero(SFourMomentum), IN_PS[2:end]...) + OUT_PS_unphysical = (OUT_PS[1:(end - 1)]..., ones(SFourMomentum)) + PSP_unphysical_in_ps = PhaseSpacePoint( + TESTPROC, TESTMODEL, TESTPSDEF, IN_PS_unphysical, OUT_PS + ) + PSP_unphysical_out_ps = PhaseSpacePoint( + TESTPROC, TESTMODEL, TESTPSDEF, IN_PS, OUT_PS_unphysical + ) + PSP_unphysical = PhaseSpacePoint( + TESTPROC, TESTMODEL, TESTPSDEF, IN_PS_unphysical, OUT_PS_unphysical + ) + + @test !QEDbase._is_in_phasespace(PSP_unphysical_in_ps) + @test !QEDbase._is_in_phasespace(PSP_unphysical_out_ps) + @test !QEDbase._is_in_phasespace(PSP_unphysical) + end + + @testset "phase space factor" begin + test_phase_space_factor = QEDbase._phase_space_factor(PSP) + groundtruth = TestImplementation._groundtruth_phase_space_factor(IN_PS, OUT_PS) + @test isapprox(test_phase_space_factor, groundtruth, atol=ATOL, rtol=RTOL) + end +end diff --git a/test/interfaces/setup.jl b/test/interfaces/setup.jl new file mode 100644 index 0000000..067a8ed --- /dev/null +++ b/test/interfaces/setup.jl @@ -0,0 +1,194 @@ +using Random +using Suppressor +using QEDbase: QEDbase +using QEDcore + +RNG = MersenneTwister(137137) +ATOL = 0.0 +RTOL = sqrt(eps()) + +_groundtruth_compute(x) = x +_groundtruth_input_validation(x) = (x > 0) +struct TestException <: QEDbase.AbstractInvalidInputException end +function _groundtruth_valid_input_assert(x) + _groundtruth_input_validation(x) || throw(TestException()) + return nothing +end +_transform_to_invalid(x) = -abs(x) +_groundtruth_post_processing(x, y) = x + y + +# setups for which the interface is implemented +abstract type AbstractTestSetup <: QEDbase.AbstractComputationSetup end +QEDbase._compute(stp::AbstractTestSetup, x) = _groundtruth_compute(x) + +# setup with default implementations +struct TestSetupDefault <: AbstractTestSetup end + +# setup with custom _assert_valid_input +struct TestSetupCustomAssertValidInput <: AbstractTestSetup end +function QEDbase._assert_valid_input(stp::TestSetupCustomAssertValidInput, x) + return _groundtruth_valid_input_assert(x) +end + +# setup with custom post processing +struct TestSetupCustomPostProcessing <: AbstractTestSetup end +function QEDbase._post_processing(::TestSetupCustomPostProcessing, x, y) + return _groundtruth_post_processing(x, y) +end + +# setup with custom input validation and post processing +struct TestSetupCustom <: AbstractTestSetup end +function QEDbase._assert_valid_input(stp::TestSetupCustom, x) + return _groundtruth_valid_input_assert(x) +end +QEDbase._post_processing(::TestSetupCustom, x, y) = _groundtruth_post_processing(x, y) + +# setup which fail on computation with default implementations +struct TestSetupFAIL <: QEDbase.AbstractComputationSetup end + +# setup which fail on computation with custom input validation, where the +# invalid input will be caught before the computation. +struct TestSetupCustomValidationFAIL <: QEDbase.AbstractComputationSetup end +function QEDbase._assert_valid_input(stp::TestSetupCustomValidationFAIL, x) + return _groundtruth_valid_input_assert(x) +end + +# setup which fail on computation with custom post processing +struct TestSetupCustomPostProcessingFAIL <: QEDbase.AbstractComputationSetup end +function QEDbase._post_processing(::TestSetupCustomPostProcessingFAIL, x, y) + return _groundtruth_post_processing(x, y) +end +@testset "general computation setup interface" begin + @testset "interface fail" begin + rnd_input = rand(RNG) + + @test_throws MethodError QEDbase._compute(TestSetupFAIL(), rnd_input) + @test_throws MethodError QEDbase.compute(TestSetupFAIL(), rnd_input) + + @test_throws MethodError QEDbase._compute( + TestSetupCustomValidationFAIL(), rnd_input + ) + @test_throws MethodError QEDbase.compute(TestSetupCustomValidationFAIL(), rnd_input) + # invalid input should be caught without throwing a MethodError + @test_throws TestException QEDbase.compute( + TestSetupCustomValidationFAIL(), _transform_to_invalid(rnd_input) + ) + + @test_throws MethodError QEDbase._compute( + TestSetupCustomPostProcessingFAIL(), rnd_input + ) + @test_throws MethodError QEDbase.compute( + TestSetupCustomPostProcessingFAIL(), rnd_input + ) + end + + @testset "default interface" begin + stp = TestSetupDefault() + + rnd_input = rand(RNG) + rnd_output = rand(RNG) + @test QEDbase._post_processing(stp, rnd_input, rnd_output) == rnd_output + @test isapprox( + QEDbase._compute(stp, rnd_input), + _groundtruth_compute(rnd_input), + atol=ATOL, + rtol=RTOL, + ) + @test isapprox( + QEDbase.compute(stp, rnd_input), + _groundtruth_compute(rnd_input), + atol=ATOL, + rtol=RTOL, + ) + end + + @testset "custom input validation" begin + stp = TestSetupCustomAssertValidInput() + rnd_input = rand(RNG) + @test QEDbase._assert_valid_input(stp, rnd_input) == nothing + @test_throws TestException QEDbase._assert_valid_input( + stp, _transform_to_invalid(rnd_input) + ) + @test_throws TestException QEDbase.compute(stp, _transform_to_invalid(rnd_input)) + end + + @testset "custom post processing" begin + stp = TestSetupCustomPostProcessing() + rnd_input = rand(RNG) + rnd_output = rand(RNG) + @test isapprox( + QEDbase._post_processing(stp, rnd_input, rnd_output), + _groundtruth_post_processing(rnd_input, rnd_output), + ) + @test isapprox( + QEDbase.compute(stp, rnd_input), + _groundtruth_post_processing(rnd_input, _groundtruth_compute(rnd_input)), + ) + end + + @testset "custom input validation and post processing" begin + stp = TestSetupCustom() + rnd_input = rand(RNG) + rnd_output = rand(RNG) + + @test_throws TestException() QEDbase.compute(stp, _transform_to_invalid(rnd_input)) + @test isapprox( + QEDbase._post_processing(stp, rnd_input, rnd_output), + _groundtruth_post_processing(rnd_input, rnd_output), + ) + @test isapprox( + QEDbase.compute(stp, rnd_input), + _groundtruth_post_processing(rnd_input, _groundtruth_compute(rnd_input)), + ) + end +end +# process setup + +struct TestParticle1 <: QEDbase.AbstractParticle end +struct TestParticle2 <: QEDbase.AbstractParticle end +struct TestParticle3 <: QEDbase.AbstractParticle end +struct TestParticle4 <: QEDbase.AbstractParticle end + +PARTICLE_SET = [TestParticle1(), TestParticle2(), TestParticle3(), TestParticle4()] + +struct TestProcess <: QEDbase.AbstractProcessDefinition end +struct TestModel <: QEDbase.AbstractModelDefinition end + +struct TestProcessSetup <: QEDbase.AbstractProcessSetup end +QEDbase.scattering_process(::TestProcessSetup) = TestProcess() +QEDbase.physical_model(::TestProcessSetup) = TestModel() + +struct TestProcessSetupFAIL <: QEDbase.AbstractProcessSetup end + +@testset "process setup interface" begin + @testset "interface fail" begin + rnd_input = rand(RNG) + @test_throws MethodError QEDbase.scattering_process(TestProcessSetupFAIL()) + @test_throws MethodError QEDbase.physical_model(TestProcessSetupFAIL()) + @test_throws MethodError QEDbase._compute(TestProcessSetupFAIL(), rnd_input) + end + + @testset "hard interface" begin + stp = TestProcessSetup() + + @test QEDbase._is_computation_setup(stp) + @test QEDbase.scattering_process(stp) == TestProcess() + @test QEDbase.physical_model(stp) == TestModel() + end + + @testset "($N_INCOMING,$N_OUTGOING)" for (N_INCOMING, N_OUTGOING) in Iterators.product( + (1, rand(RNG, 2:8)), (1, rand(RNG, 2:8)) + ) + INCOMING_PARTICLES = rand(RNG, PARTICLE_SET, N_INCOMING) + OUTGOING_PARTICLES = rand(RNG, PARTICLE_SET, N_OUTGOING) + + @suppress QEDbase.incoming_particles(::TestProcess) = INCOMING_PARTICLES + @suppress QEDbase.outgoing_particles(::TestProcess) = OUTGOING_PARTICLES + + @testset "delegated functions" begin + stp = TestProcessSetup() + @test QEDbase.number_incoming_particles(stp) == N_INCOMING + @test QEDbase.number_outgoing_particles(stp) == N_OUTGOING + end + end +end diff --git a/test/phase_spaces.jl b/test/phase_spaces.jl new file mode 100644 index 0000000..7116a4b --- /dev/null +++ b/test/phase_spaces.jl @@ -0,0 +1,261 @@ +using Random +using StaticArrays +using QEDbase: QEDbase +using QEDcore + +# 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) +BUF = IOBuffer() + +@testset "broadcast" begin + test_func(ps_def) = ps_def + @test test_func.(TESTPSDEF) == TESTPSDEF +end + +@testset "Stateful Particle" begin + DIRECTIONS = [QEDbase.Incoming(), QEDbase.Outgoing()] + SPECIES = [Electron(), Positron()] #=, Muon(), AntiMuon(), Tauon(), AntiTauon()=# + + for (species, dir) in Iterators.product(SPECIES, DIRECTIONS) + mom = rand(RNG, SFourMomentum) + + particle_stateful = ParticleStateful(dir, species, mom) + + # particle interface + @test QEDbase.is_fermion(particle_stateful) == QEDbase.is_fermion(species) + @test QEDbase.is_boson(particle_stateful) == QEDbase.is_boson(species) + @test QEDbase.is_particle(particle_stateful) == QEDbase.is_particle(species) + @test QEDbase.is_anti_particle(particle_stateful) == + QEDbase.is_anti_particle(species) + @test QEDbase.is_incoming(particle_stateful) == QEDbase.is_incoming(dir) + @test QEDbase.is_outgoing(particle_stateful) == QEDbase.is_outgoing(dir) + @test QEDbase.mass(particle_stateful) == QEDbase.mass(species) + @test QEDbase.charge(particle_stateful) == QEDbase.charge(species) + + # accessors + @test particle_stateful.dir == dir + @test QEDbase.particle_direction(particle_stateful) == particle_stateful.dir + @test particle_stateful.species == species + @test QEDbase.particle_species(particle_stateful) == particle_stateful.species + @test particle_stateful.mom == mom + @test QEDbase.momentum(particle_stateful) == mom + + # printing + print(BUF, particle_stateful) + @test String(take!(BUF)) == "$(dir) $(species): $(mom)" + + show(BUF, MIME"text/plain"(), particle_stateful) + @test String(take!(BUF)) == + "ParticleStateful: $(dir) $(species)\n momentum: $(mom)\n" + end +end + +@testset "Phasespace Point" begin + in_el_mom = rand(RNG, SFourMomentum) + in_ph_mom = rand(RNG, SFourMomentum) + out_el_mom = rand(RNG, SFourMomentum) + out_ph_mom = rand(RNG, SFourMomentum) + + in_el = ParticleStateful(QEDbase.Incoming(), Electron(), in_el_mom) + in_ph = ParticleStateful(QEDbase.Incoming(), Photon(), in_ph_mom) + out_el = ParticleStateful(QEDbase.Outgoing(), Electron(), out_el_mom) + out_ph = ParticleStateful(QEDbase.Outgoing(), Photon(), out_ph_mom) + + in_particles_valid = (in_el, in_ph) + in_particles_invalid = (in_el, out_ph) + + out_particles_valid = (out_el, out_ph) + out_particles_invalid = (out_el, in_ph) + + model = TESTMODEL + process = TestImplementation.TestProcess((Electron(), Photon()), (Electron(), Photon())) + phasespace_def = TESTPSDEF + + psp = PhaseSpacePoint( + process, model, phasespace_def, in_particles_valid, out_particles_valid + ) + + take!(BUF) + print(BUF, psp) + @test String(take!(BUF)) == "PhaseSpacePoint of $(process)" + + show(BUF, MIME"text/plain"(), psp) + @test match( + r"PhaseSpacePoint:\n process: (.*)TestProcess(.*)\n model: (.*)TestModel(.*)\n phasespace definition: (.*)TestPhasespaceDef(.*)\n incoming particles:\n -> incoming electron: (.*)\n -> incoming photon: (.*)\n outgoing particles:\n -> outgoing electron: (.*)\n -> outgoing photon: (.*)\n", + String(take!(BUF)), + ) isa RegexMatch + + @testset "Accessor" begin + @test QEDbase.momentum(psp, QEDbase.Incoming(), 1) == in_el.mom + @test QEDbase.momentum(psp, QEDbase.Incoming(), 2) == in_ph.mom + @test QEDbase.momentum(psp, QEDbase.Outgoing(), 1) == out_el.mom + @test QEDbase.momentum(psp, QEDbase.Outgoing(), 2) == out_ph.mom + + @test psp[QEDbase.Incoming(), 1] == in_el + @test psp[QEDbase.Incoming(), 2] == in_ph + @test psp[QEDbase.Outgoing(), 1] == out_el + @test psp[QEDbase.Outgoing(), 2] == out_ph + end + + @testset "Error handling" begin + if (VERSION >= v"1.8") + # julia versions before 1.8 did not have support for regex matching in @test_throws + @test_throws r"expected incoming photon but got outgoing photon" PhaseSpacePoint( + process, model, phasespace_def, in_particles_invalid, out_particles_valid + ) + + @test_throws r"expected outgoing photon but got incoming photon" PhaseSpacePoint( + process, model, phasespace_def, in_particles_valid, out_particles_invalid + ) + + @test_throws r"expected incoming electron but got incoming photon" PhaseSpacePoint( + process, model, phasespace_def, (in_ph, in_el), out_particles_valid + ) + + @test_throws r"expected outgoing electron but got outgoing photon" PhaseSpacePoint( + process, model, phasespace_def, in_particles_valid, (out_ph, out_el) + ) + + @test_throws r"expected 2 outgoing particles for the process but got 1" PhaseSpacePoint( + process, model, phasespace_def, in_particles_valid, (out_el,) + ) + + @test_throws r"expected 2 incoming particles for the process but got 1" PhaseSpacePoint( + process, model, phasespace_def, (out_el,), out_particles_valid + ) + + @test_throws r"expected 2 outgoing particles for the process but got 3" PhaseSpacePoint( + process, model, phasespace_def, in_particles_valid, (out_el, out_el, out_ph) + ) + + @test_throws r"expected 2 incoming particles for the process but got 3" PhaseSpacePoint( + process, model, phasespace_def, (in_el, in_el, in_ph), out_particles_valid + ) + end + + @test_throws BoundsError QEDbase.momentum(psp, QEDbase.Incoming(), -1) + @test_throws BoundsError QEDbase.momentum(psp, QEDbase.Outgoing(), -1) + @test_throws BoundsError QEDbase.momentum(psp, QEDbase.Incoming(), 4) + @test_throws BoundsError QEDbase.momentum(psp, QEDbase.Outgoing(), 4) + + @test_throws BoundsError psp[QEDbase.Incoming(), -1] + @test_throws BoundsError psp[QEDbase.Outgoing(), -1] + @test_throws BoundsError psp[QEDbase.Incoming(), 4] + @test_throws BoundsError psp[QEDbase.Outgoing(), 4] + + @test_throws QEDbase.InvalidInputError PhaseSpacePoint( + process, model, phasespace_def, in_particles_invalid, out_particles_valid + ) + + @test_throws QEDbase.InvalidInputError PhaseSpacePoint( + process, model, phasespace_def, in_particles_valid, out_particles_invalid + ) + + @test_throws QEDbase.InvalidInputError PhaseSpacePoint( + process, model, phasespace_def, (in_ph, in_el), out_particles_valid + ) + + @test_throws QEDbase.InvalidInputError PhaseSpacePoint( + process, model, phasespace_def, in_particles_valid, (out_ph, out_el) + ) + end + + @testset "Generation from momenta" begin + test_psp = PhaseSpacePoint( + process, model, phasespace_def, (in_el_mom, in_ph_mom), (out_el_mom, out_ph_mom) + ) + + @test test_psp.proc == process + @test test_psp.model == model + @test test_psp.ps_def == phasespace_def + + @test test_psp[QEDbase.Incoming(), 1] == in_el + @test test_psp[QEDbase.Incoming(), 2] == in_ph + @test test_psp[QEDbase.Outgoing(), 1] == out_el + @test test_psp[QEDbase.Outgoing(), 2] == out_ph + end + + @testset "Error handling from momenta" for (i, o) in + Iterators.product([1, 3, 4, 5], [1, 3, 4, 5]) + @test_throws QEDbase.InvalidInputError PhaseSpacePoint( + process, + model, + phasespace_def, + TestImplementation._rand_momenta(RNG, i), + TestImplementation._rand_momenta(RNG, o), + ) + end + + @testset "Directional PhaseSpacePoint" begin + @test psp isa PhaseSpacePoint + @test psp isa InPhaseSpacePoint + @test psp isa OutPhaseSpacePoint + + in_psp = InPhaseSpacePoint(process, model, phasespace_def, in_particles_valid) + out_psp = OutPhaseSpacePoint(process, model, phasespace_def, out_particles_valid) + in_psp_from_moms = InPhaseSpacePoint( + process, model, phasespace_def, (in_el_mom, in_ph_mom) + ) + out_psp_from_moms = OutPhaseSpacePoint( + process, model, phasespace_def, (out_el_mom, out_ph_mom) + ) + + @test in_psp isa InPhaseSpacePoint + @test !(in_psp isa OutPhaseSpacePoint) + @test in_psp_from_moms isa InPhaseSpacePoint + @test !(in_psp_from_moms isa OutPhaseSpacePoint) + + @test out_psp isa OutPhaseSpacePoint + @test !(out_psp isa InPhaseSpacePoint) + @test out_psp_from_moms isa OutPhaseSpacePoint + @test !(out_psp_from_moms isa InPhaseSpacePoint) + + @test_throws QEDbase.InvalidInputError InPhaseSpacePoint( + process, model, phasespace_def, in_particles_invalid + ) + @test_throws QEDbase.InvalidInputError OutPhaseSpacePoint( + process, model, phasespace_def, out_particles_invalid + ) + + @testset "Error handling from momenta" for i in [1, 3, 4, 5] + @test_throws QEDbase.InvalidInputError InPhaseSpacePoint( + process, model, phasespace_def, TestImplementation._rand_momenta(RNG, i) + ) + @test_throws QEDbase.InvalidInputError OutPhaseSpacePoint( + process, model, phasespace_def, TestImplementation._rand_momenta(RNG, i) + ) + end + end +end + +@testset "Coordinate System" begin + @testset "Pretty printing" begin + print(BUF, SphericalCoordinateSystem()) + @test String(take!(BUF)) == "spherical coordinates" + end +end +@testset "Reference Frame" begin + @testset "Pretty printing" begin + print(BUF, ElectronRestFrame()) + @test String(take!(BUF)) == "electron rest frame" + print(BUF, CenterOfMomentumFrame()) + @test String(take!(BUF)) == "center-of-momentum frame" + end +end + +@testset "Phasespace Definition" for (coord_sys, frame) in Iterators.product( + [SphericalCoordinateSystem()], [ElectronRestFrame(), CenterOfMomentumFrame()] +) + ps_def = PhasespaceDefinition(coord_sys, frame) + + @testset "Pretty printing" begin + print(BUF, ps_def) + @test String(take!(BUF)) == "$coord_sys in $frame" + end +end diff --git a/test/runtests.jl b/test/runtests.jl index d204607..02a159f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,6 +3,10 @@ using Test using SafeTestsets begin + @time @safetestset "phase spaces" begin + include("phase_spaces.jl") + end + # algebraic objects @time @safetestset "four momentum" begin include("algebraic_objects/four_momentum.jl") @@ -25,11 +29,28 @@ begin include("particles/types.jl") end - @time @safetestset "particle types" begin + @time @safetestset "particle states" begin include("particles/states.jl") end @time @safetestset "particle spinors" begin include("particles/spinors.jl") end + + @time @safetestset "particle base states" begin + include("particles/states.jl") + end + + @time @safetestset "particle propagators" begin + include("particles/propagators.jl") + end + + # interfaces + @time @safetestset "process interface" begin + include("interfaces/process.jl") + end + + @time @safetestset "computation setup interface" begin + include("interfaces/setup.jl") + end end diff --git a/test/test_implementation/TestImplementation.jl b/test/test_implementation/TestImplementation.jl new file mode 100644 index 0000000..830e1ce --- /dev/null +++ b/test/test_implementation/TestImplementation.jl @@ -0,0 +1,39 @@ +""" +This module provides a full implementation of the model and process interface. Its purpose is only for testing and it does not reflect any +real-world physics. + +The module exports: + +``` +TestParticle1 # set of test particles without properties +TestParticle2 +TestParticle3 +TestParticle4 +TestModel # dummy compute model +TestModel_FAIL # failing compute model +TestProcess # dummy scattering process +TestProcess_FAIL # failing scattering process +TestPhasespaceDef # dummy phase space definition +TestPhasespaceDef_FAIL # failing phase space definition +``` +The respective groundtruth implementations for the interface functions are stored in `groundtruths.jl` +""" +module TestImplementation + +export TestParticle1, TestParticle2, TestParticle3, TestParticle4, PARTICLE_SET +export TestModel, TestModel_FAIL +export TestProcess, TestProcess_FAIL +export TestPhasespaceDef, TestPhasespaceDef_FAIL + +using Random +using QEDbase: QEDbase +using QEDcore +using StaticArrays + +include("groundtruths.jl") +include("test_model.jl") +include("test_process.jl") +include("random_momenta.jl") +include("utils.jl") + +end diff --git a/test/test_implementation/groundtruths.jl b/test/test_implementation/groundtruths.jl new file mode 100644 index 0000000..e6dc7e2 --- /dev/null +++ b/test/test_implementation/groundtruths.jl @@ -0,0 +1,267 @@ +import QEDbase.AbstractFourMomentum + +""" + _groundtruth_incident_flux(in_ps) + +Test implementation of the incident flux. Return the Minkowski square of the sum of the incoming momenta: + +```math +\\begin{align} +I = \\left(\\sum p_i\\right)^2, +\\end{align} +``` +where \$p_i\\in\\mathrm{ps_in}\$. +""" +function _groundtruth_incident_flux(in_ps) + s = sum(in_ps) + return s * s +end + +""" + _groundtruth_matrix_element(in_ps, out_ps) + +Test implementation for a matrix elements. Returns a list of three complex numbers without any physical meaning. +""" +function _groundtruth_matrix_element(in_ps, out_ps) + s_in = sum(in_ps) + s_out = sum(out_ps) + res = s_in * s_in + 1im * (s_out * s_out) + return (res, 2 * res, 3 * res) +end + +""" + _groundtruth_averaging_norm(proc) + +Test implementation of the averaging norm. Returns the inverse of the sum of all external particles of the passed process. +""" +function _groundtruth_averaging_norm(proc) + return 1.0 / ( + QEDbase.number_incoming_particles(proc) + QEDbase.number_outgoing_particles(proc) + ) +end + +""" + _groundtruth_is_in_phasespace(in_ps, out_ps) + +Test implementation of the phase space check. Return `false` if either the momentum of the first incoming particle is exactly `zero(SFourMomentum)`, or if the momentum of the last outgoing momentum is exactly `ones(SFourMomentum)`. Otherwise, return true. +""" +function _groundtruth_is_in_phasespace(in_ps, out_ps) + if in_ps[1] == SFourMomentum(zeros(4)) + return false + end + if out_ps[end] == ones(SFourMomentum) + return false + end + return true +end + +""" + _groundtruth_phase_space_factor(in_ps, out_ps) + +Test implementation of the phase space factor. Return the inverse of the product of the energies of all external particles. +""" +function _groundtruth_phase_space_factor(in_ps, out_ps) + en_in = QEDbase.getE.(in_ps) + en_out = QEDbase.getE.(out_ps) + return 1 / (prod(en_in) * prod(en_out)) +end + +function _groundtruth_generate_momenta(ps_coords) + moms = _furl_moms(ps_coords) + return moms +end +""" + _groundtruth_unsafe_probability(proc, in_ps, out_ps) + +Test implementation of the unsafe differential probability. Uses the test implementations of `_groundtruth_matrix_element`,`_groundtruth_averaging_norm` and `_groundtruth_phase_space_factor`. +""" +function _groundtruth_unsafe_probability(proc, in_ps, out_ps) + mat_el = _groundtruth_matrix_element(in_ps, out_ps) + mat_el_sq = abs2.(mat_el) + normalization = _groundtruth_averaging_norm(proc) + ps_fac = _groundtruth_phase_space_factor(in_ps, out_ps) + return sum(mat_el_sq) * ps_fac * normalization +end + +function _groundtruth_unsafe_probability( + proc, in_ps::AbstractVector, out_ps::AbstractMatrix +) + res = Vector{Float64}(undef, size(out_ps, 2)) + for i in 1:size(out_ps, 2) + res[i] = _groundtruth_unsafe_probability(proc, in_ps, view(out_ps, :, i)) + end + return res +end + +function _groundtruth_unsafe_probability( + proc, in_ps::AbstractMatrix, out_ps::AbstractVector +) + res = Vector{Float64}(undef, size(in_ps, 2)) + for i in 1:size(in_ps, 2) + res[i] = _groundtruth_unsafe_probability(proc, view(in_ps, :, i), out_ps) + end + return res +end + +function _groundtruth_unsafe_probability( + proc, in_ps::AbstractMatrix, out_ps::AbstractMatrix +) + res = Matrix{Float64}(undef, size(in_ps, 2), size(out_ps, 2)) + for i in 1:size(in_ps, 2) + for j in 1:size(out_ps, 2) + res[i, j] = _groundtruth_unsafe_probability( + proc, view(in_ps, :, i), view(out_ps, :, j) + ) + end + end + return res +end + +""" + _groundtruth_safe_probability(proc, in_ps, out_ps) + +Test implementation of the safe differential probability. Uses the test implementations of `_groundtruth_is_in_phasespace` and `_groundtruth_unsafe_probability`. +""" +function _groundtruth_safe_probability(proc, in_ps, out_ps) + if !_groundtruth_is_in_phasespace(in_ps, out_ps) + return zero(Float64) + end + return _groundtruth_unsafe_probability(proc, in_ps, out_ps) +end + +function _groundtruth_safe_probability(proc, in_ps::AbstractVector, out_ps::AbstractMatrix) + res = Vector{Float64}(undef, size(out_ps, 2)) + for i in 1:size(out_ps, 2) + res[i] = _groundtruth_safe_probability(proc, in_ps, view(out_ps, :, i)) + end + return res +end + +function _groundtruth_safe_probability(proc, in_ps::AbstractMatrix, out_ps::AbstractVector) + res = Vector{Float64}(undef, size(in_ps, 2)) + for i in 1:size(in_ps, 2) + res[i] = _groundtruth_safe_probability(proc, view(in_ps, :, i), out_ps) + end + return res +end + +function _groundtruth_safe_probability(proc, in_ps::AbstractMatrix, out_ps::AbstractMatrix) + res = Matrix{Float64}(undef, size(in_ps, 2), size(out_ps, 2)) + for i in 1:size(in_ps, 2) + for j in 1:size(out_ps, 2) + res[i, j] = _groundtruth_safe_probability( + proc, view(in_ps, :, i), view(out_ps, :, j) + ) + end + end + return res +end + +""" + _groundtruth_unsafe_diffCS(proc, in_ps, out_ps) + +Test implementation of the unsafe differential cross section. Uses the test implementations of `_groundtruth_incident_flux` and `_groundtruth_unsafe_probability`. +""" +function _groundtruth_unsafe_diffCS(proc, in_ps, out_ps) + init_flux = _groundtruth_incident_flux(in_ps) + return _groundtruth_unsafe_probability(proc, in_ps, out_ps) / (4 * init_flux) +end + +function _groundtruth_unsafe_diffCS(proc, in_ps::AbstractVector, out_ps::AbstractMatrix) + res = Vector{Float64}(undef, size(out_ps, 2)) + for i in 1:size(out_ps, 2) + res[i] = _groundtruth_unsafe_diffCS(proc, in_ps, view(out_ps, :, i)) + end + return res +end + +function _groundtruth_unsafe_diffCS(proc, in_ps::AbstractMatrix, out_ps::AbstractVector) + res = Vector{Float64}(undef, size(in_ps, 2)) + for i in 1:size(in_ps, 2) + res[i] = _groundtruth_unsafe_diffCS(proc, view(in_ps, :, i), out_ps) + end + return res +end + +function _groundtruth_unsafe_diffCS(proc, in_ps::AbstractMatrix, out_ps::AbstractMatrix) + res = Matrix{Float64}(undef, size(in_ps, 2), size(out_ps, 2)) + for i in 1:size(in_ps, 2) + for j in 1:size(out_ps, 2) + res[i, j] = _groundtruth_unsafe_diffCS( + proc, view(in_ps, :, i), view(out_ps, :, j) + ) + end + end + return res +end + +""" + _groundtruth_safe_diffCS(proc, in_ps, out_ps) + +Test implementation of the safe differential cross section. Uses the test implementations of `_groundtruth_is_in_phasespace` and `_groundtruth_unsafe_diffCS`. +""" +function _groundtruth_safe_diffCS(proc, in_ps, out_ps) + if !_groundtruth_is_in_phasespace(in_ps, out_ps) + return zero(Float64) + end + return _groundtruth_unsafe_diffCS(proc, in_ps, out_ps) +end + +function _groundtruth_safe_diffCS(proc, in_ps::AbstractVector, out_ps::AbstractMatrix) + res = Vector{Float64}(undef, size(out_ps, 2)) + for i in 1:size(out_ps, 2) + res[i] = _groundtruth_safe_diffCS(proc, in_ps, view(out_ps, :, i)) + end + return res +end + +function _groundtruth_safe_diffCS(proc, in_ps::AbstractMatrix, out_ps::AbstractVector) + res = Vector{Float64}(undef, size(in_ps, 2)) + for i in 1:size(in_ps, 2) + res[i] = _groundtruth_safe_diffCS(proc, view(in_ps, :, i), out_ps) + end + return res +end + +function _groundtruth_safe_diffCS(proc, in_ps::AbstractMatrix, out_ps::AbstractMatrix) + res = Matrix{Float64}(undef, size(in_ps, 2), size(out_ps, 2)) + for i in 1:size(in_ps, 2) + for j in 1:size(out_ps, 2) + res[i, j] = _groundtruth_safe_diffCS( + proc, view(in_ps, :, i), view(out_ps, :, j) + ) + end + end + return res +end + +""" + _groundtruth_total_probability(in_ps::AbstractVector) + +Test implementation of the total cross section. Return the Minkowski square of the sum the momenta of all incoming particles. +""" +function _groundtruth_total_probability( + in_ps::NTuple{N,T} +) where {N,T<:QEDbase.AbstractFourMomentum} + Ptot = sum(in_ps) + return Ptot * Ptot +end + +function _groundtruth_total_probability( + in_pss::Vector{NTuple{N,T}} +) where {N,T<:QEDbase.AbstractFourMomentum} + return _groundtruth_total_probability.(in_pss) +end + +function _groundtruth_total_cross_section( + in_ps::NTuple{N,T} +) where {N,T<:QEDbase.AbstractFourMomentum} + init_flux = _groundtruth_incident_flux(in_ps) + return _groundtruth_total_probability(in_ps) / (4 * init_flux) +end + +function _groundtruth_total_cross_section( + in_pss::Vector{NTuple{N,T}} +) where {N,T<:QEDbase.AbstractFourMomentum} + return _groundtruth_total_cross_section.(in_psps) +end diff --git a/test/test_implementation/random_momenta.jl b/test/test_implementation/random_momenta.jl new file mode 100644 index 0000000..b631551 --- /dev/null +++ b/test/test_implementation/random_momenta.jl @@ -0,0 +1,83 @@ + +""" +Return a tuple of random four momenta, i.e. a random phase space point. +""" +function _rand_momenta(rng::AbstractRNG, n) + return NTuple{n,SFourMomentum}(SFourMomentum(rand(rng, 4)) for _ in 1:n) +end + +""" +Return a vector of tuples of random four momenta, i.e. a collection of phase space points. +n1 is the size of the phase space point, n2 is the number of points. +""" +function _rand_momenta(rng::AbstractRNG, n1, n2) + moms = Vector{NTuple{n1,SFourMomentum}}(undef, n2) + for i in 1:n2 + moms[i] = _rand_momenta(rng, n1) + end + return moms +end + +""" +Return a random phase space point that is failing the incoming phase space constraint, +i.e. the first entry of the phase space is the null momentum. +""" +function _rand_in_momenta_failing(rng::AbstractRNG, n) + return (zero(SFourMomentum), _rand_momenta(rng, n - 1)...) +end + +""" +Return a random phase space point that is failing the outgoing phase space constraint, +i.e. the last entry of the phase space is the unit momentum. +""" +function _rand_out_momenta_failing(rng::AbstractRNG, n) + return (_rand_momenta(rng, n - 1)..., ones(SFourMomentum)) +end + +""" +Return a collection of incoming phase space points, where the first point is failing the phase space constraint, +i.e. the first entry of the vector is invalid but the others pass. +n1 is the size of the phase space point, n2 is the number of points. +""" +function _rand_in_momenta_failing_mix(rng::AbstractRNG, n1, n2) + moms = _rand_momenta(rng, n1, n2) + moms[1] = _rand_in_momenta_failing(rng, n1) + return moms +end + +""" +Return a collection of incoming phase space points, where all points are failing the phase space constraint, +i.e. their first entries are null momenta. +n1 is the size of the phase space point, n2 is the number of points. +""" +function _rand_in_momenta_failing_all(rng::AbstractRNG, n1, n2) + moms = Vector{NTuple{n1,SFourMomentum}}(undef, n2) + for i in 1:n2 + moms[i] = _rand_in_momenta_failing(rng, n1) + end + return moms +end + +""" +Return a vector of outgoing phase space points, where the first point is failing the phase space constraint, +i.e. the last entry of the vector is invalid but the others pass. +n1 is the size of the phase space point, n2 is the number of points. +""" +function _rand_out_momenta_failing_mix(rng::AbstractRNG, n1, n2) + moms = _rand_momenta(rng, n1, n2) + moms[end] = _rand_out_momenta_failing(rng, n1) + return moms +end + +""" +Return a vector of outgoing phase space points, where all points are failing the phase space constraint, +i.e. their last entries are unit momenta. +n1 is the size of the phase space point, n2 is the number of points. +""" +function _rand_out_momenta_failing_all(rng::AbstractRNG, n1, n2) + moms = Vector{NTuple{n1,SFourMomentum}}(undef, n2) + for i in 1:n2 + moms[i] = _rand_out_momenta_failing(rng, n1) + end + return moms +end diff --git a/test/test_implementation/test_model.jl b/test/test_implementation/test_model.jl new file mode 100644 index 0000000..6e5aaf2 --- /dev/null +++ b/test/test_implementation/test_model.jl @@ -0,0 +1,4 @@ +struct TestModel <: QEDbase.AbstractModelDefinition end +QEDbase.fundamental_interaction_type(::TestModel) = :test_interaction + +struct TestModel_FAIL <: QEDbase.AbstractModelDefinition end diff --git a/test/test_implementation/test_process.jl b/test/test_implementation/test_process.jl new file mode 100644 index 0000000..67c844a --- /dev/null +++ b/test/test_implementation/test_process.jl @@ -0,0 +1,131 @@ +# dummy particles +struct TestParticleFermion <: FermionLike end +struct TestParticleBoson <: BosonLike end + +const PARTICLE_SET = [TestParticleFermion(), TestParticleBoson()] + +""" + TestProcess(rng,incoming_particles,outgoing_particles) + +""" +struct TestProcess{IP<:Tuple,OP<:Tuple} <: QEDbase.AbstractProcessDefinition + incoming_particles::IP + outgoing_particles::OP +end + +function TestProcess(rng::AbstractRNG, N_in::Int, N_out::Int) + in_particles = rand(rng, PARTICLE_SET, N_in) + out_particles = rand(rng, PARTICLE_SET, N_out) + return TestProcess(in_particles, out_particles) +end + +QEDbase.incoming_particles(proc::TestProcess) = proc.incoming_particles +QEDbase.outgoing_particles(proc::TestProcess) = proc.outgoing_particles + +struct TestProcess_FAIL{IP<:Tuple,OP<:Tuple} <: QEDbase.AbstractProcessDefinition + incoming_particles::IP + outgoing_particles::OP +end + +function TestProcess_FAIL(rng::AbstractRNG, N_in::Int, N_out::Int) + in_particles = Tuple(rand(rng, PARTICLE_SET, N_in)) + out_particles = Tuple(rand(rng, PARTICLE_SET, N_out)) + return TestProcess_FAIL(in_particles, out_particles) +end + +function QEDbase.in_phase_space_dimension(proc::TestProcess, ::TestModel) + return QEDbase.number_incoming_particles(proc) * 4 +end +function QEDbase.out_phase_space_dimension(proc::TestProcess, ::TestModel) + return QEDbase.number_outgoing_particles(proc) * 4 +end + +""" +Test process with no implemented interface. Should fail every usage except construction. +""" +struct TestProcess_FAIL_ALL{IP<:Tuple,OP<:Tuple} <: QEDbase.AbstractProcessDefinition + incoming_particles::IP + outgoing_particles::OP +end + +function TestProcess_FAIL_ALL(rng::AbstractRNG, N_in::Int, N_out::Int) + in_particles = Tuple(rand(rng, PARTICLE_SET, N_in)) + out_particles = Tuple(rand(rng, PARTICLE_SET, N_out)) + return TestProcess_FAIL_ALL(in_particles, out_particles) +end + +""" +Test process with no implemented interface except the incoming and outgoing particles. +Should fail every usage except construction of itself and the respective phase space point for given four-momenta. +""" +struct TestProcess_FAIL_DIFFCS{IP<:Tuple,OP<:Tuple} <: QEDbase.AbstractProcessDefinition + incoming_particles::IP + outgoing_particles::OP +end + +function TestProcess_FAIL_DIFFCS(rng::AbstractRNG, N_in::Int, N_out::Int) + in_particles = Tuple(rand(rng, PARTICLE_SET, N_in)) + out_particles = Tuple(rand(rng, PARTICLE_SET, N_out)) + return TestProcess_FAIL_DIFFCS(in_particles, out_particles) +end + +QEDbase.incoming_particles(proc::TestProcess_FAIL_DIFFCS) = proc.incoming_particles +QEDbase.outgoing_particles(proc::TestProcess_FAIL_DIFFCS) = proc.outgoing_particles + +# dummy phase space definition + failing phase space definition +struct TestPhasespaceDef <: QEDbase.AbstractPhasespaceDefinition end +struct TestPhasespaceDef_FAIL <: QEDbase.AbstractPhasespaceDefinition end + +# dummy implementation of the process interface + +function QEDbase._incident_flux(in_psp::InPhaseSpacePoint{<:TestProcess,<:TestModel}) + return _groundtruth_incident_flux(QEDbase.momenta(in_psp, QEDbase.Incoming())) +end + +function QEDbase._averaging_norm(proc::TestProcess) + return _groundtruth_averaging_norm(proc) +end + +function QEDbase._matrix_element(psp::PhaseSpacePoint{<:TestProcess,TestModel}) + in_ps = QEDbase.momenta(psp, QEDbase.Incoming()) + out_ps = QEDbase.momenta(psp, QEDbase.Outgoing()) + return _groundtruth_matrix_element(in_ps, out_ps) +end + +function QEDbase._is_in_phasespace(psp::PhaseSpacePoint{<:TestProcess,TestModel}) + in_ps = QEDbase.momenta(psp, QEDbase.Incoming()) + out_ps = QEDbase.momenta(psp, QEDbase.Outgoing()) + return _groundtruth_is_in_phasespace(in_ps, out_ps) +end + +function QEDbase._phase_space_factor( + psp::PhaseSpacePoint{<:TestProcess,TestModel,TestPhasespaceDef} +) + in_ps = QEDbase.momenta(psp, QEDbase.Incoming()) + out_ps = QEDbase.momenta(psp, QEDbase.Outgoing()) + return _groundtruth_phase_space_factor(in_ps, out_ps) +end + +function QEDbase._generate_incoming_momenta( + proc::TestProcess, + model::TestModel, + phase_space_def::TestPhasespaceDef, + in_phase_space::NTuple{N,T}, +) where {N,T<:Real} + return _groundtruth_generate_momenta(in_phase_space) +end + +function QEDbase._generate_outgoing_momenta( + proc::TestProcess, + model::TestModel, + phase_space_def::TestPhasespaceDef, + out_phase_space::NTuple{N,T}, +) where {N,T<:Real} + return _groundtruth_generate_momenta(out_phase_space) +end + +function QEDbase._total_probability( + in_psp::InPhaseSpacePoint{<:TestProcess,<:TestModel,<:TestPhasespaceDef} +) + return _groundtruth_total_probability(QEDbase.momenta(in_psp, QEDbase.Incoming())) +end diff --git a/test/test_implementation/utils.jl b/test/test_implementation/utils.jl new file mode 100644 index 0000000..3a141d9 --- /dev/null +++ b/test/test_implementation/utils.jl @@ -0,0 +1,44 @@ + +# Check if any failed type is in the input +_any_fail(x...) = true +_any_fail(::TestProcess, ::TestModel) = false +_any_fail(::TestProcess, ::TestModel, ::TestPhasespaceDef) = false + +# unrolls all elements of a list of four-momenta into vector of coordinates +function _unroll_moms(ps_moms::AbstractVector{T}) where {T<:QEDbase.AbstractFourMomentum} + return collect(Iterators.flatten(ps_moms)) +end + +function _unroll_moms(ps_moms::AbstractMatrix{T}) where {T<:QEDbase.AbstractFourMomentum} + res = Matrix{eltype(T)}(undef, size(ps_moms, 1) * 4, size(ps_moms, 2)) + for i in 1:size(ps_moms, 2) + res[:, i] .= _unroll_moms(view(ps_moms, :, i)) + end + return res +end + +flat_components(moms::AbstractVecOrMat) = _unroll_moms(moms) +flat_components(moms::Tuple) = Tuple(_unroll_moms([moms...])) + +# collect components of four-momenta from a vector of coordinates +function __furl_moms(ps_coords::AbstractVector{T}) where {T<:Real} + return SFourMomentum.(eachcol(reshape(ps_coords, 4, :))) +end + +function _furl_moms(ps_coords::AbstractVector{T}) where {T<:Real} + @assert length(ps_coords) % 4 == 0 + return __furl_moms(ps_coords) +end + +function _furl_moms(ps_coords::AbstractMatrix{T}) where {T<:Real} + @assert size(ps_coords, 1) % 4 == 0 + res = Matrix{SFourMomentum}(undef, Int(size(ps_coords, 1)//4), size(ps_coords, 2)) + for i in 1:size(ps_coords, 2) + res[:, i] .= __furl_moms(view(ps_coords, :, i)) + end + return res +end + +function _furl_moms(moms::NTuple{N,Float64}) where {N} + return Tuple(_furl_moms(Vector{Float64}([moms...]))) +end