From c768a574850d3135a9d5a5572e6881e3250b20a8 Mon Sep 17 00:00:00 2001 From: Anton Reinhard Date: Mon, 27 May 2024 22:21:42 +0200 Subject: [PATCH] Introduce InPhaseSpacePoint (#68) This PR adds a type `InPhaseSpacePoint`, which is a partially specialized version of a `PhaseSpacePoint`. It can be used to dispatch in places where only the incoming part of a phase space is required. Construction of PhaseSpacePoints now allows empty tuples for either incoming or outgoing phase spaces. --------- Co-authored-by: Uwe Hernandez Acosta Co-authored-by: Uwe Hernandez Acosta Co-authored-by: AntonReinhard --- src/QEDprocesses.jl | 11 +- src/cross_section/diff_cross_section.jl | 6 +- src/cross_section/internal.jl | 131 ------ src/cross_section/total_cross_section.jl | 83 +--- src/cross_section/total_probability.jl | 70 +-- src/interfaces/process_interface.jl | 84 ++-- src/momentum_generation.jl | 20 +- src/phase_spaces.jl | 436 ------------------ src/phase_spaces/access.jl | 59 +++ src/phase_spaces/create.jl | 158 +++++++ src/phase_spaces/print.jl | 58 +++ src/phase_spaces/types.jl | 177 +++++++ src/phase_spaces/utility.jl | 141 ++++++ .../perturbative/cross_section.jl | 46 +- .../perturbative/kinematics.jl | 22 +- .../perturbative/total_probability.jl | 13 +- test/cross_sections.jl | 46 +- test/interfaces/process_interface.jl | 19 +- test/phase_spaces.jl | 76 ++- .../one_photon_compton/perturbative.jl | 52 ++- test/test_implementation/groundtruths.jl | 29 +- test/test_implementation/random_momenta.jl | 87 ++-- test/test_implementation/test_process.jl | 21 +- test/test_implementation/utils.jl | 7 +- 24 files changed, 909 insertions(+), 943 deletions(-) delete mode 100644 src/cross_section/internal.jl delete mode 100644 src/phase_spaces.jl create mode 100644 src/phase_spaces/access.jl create mode 100644 src/phase_spaces/create.jl create mode 100644 src/phase_spaces/print.jl create mode 100644 src/phase_spaces/types.jl create mode 100644 src/phase_spaces/utility.jl diff --git a/src/QEDprocesses.jl b/src/QEDprocesses.jl index 5aaf9cf..608a96e 100644 --- a/src/QEDprocesses.jl +++ b/src/QEDprocesses.jl @@ -10,6 +10,7 @@ export AbstractModelDefinition, fundamental_interaction_type # Abstract process interface export AbstractProcessDefinition, incoming_particles, outgoing_particles export number_incoming_particles, number_outgoing_particles +export particles, number_particles # probabilities export differential_probability, unsafe_differential_probability @@ -34,7 +35,7 @@ export propagator export AbstractCoordinateSystem, SphericalCoordinateSystem export AbstractFrameOfReference, CenterOfMomentumFrame, ElectronRestFrame export AbstractPhasespaceDefinition, PhasespaceDefinition -export ParticleStateful, PhaseSpacePoint +export ParticleStateful, PhaseSpacePoint, InPhaseSpacePoint, OutPhaseSpacePoint export spin, polarization, particle_direction, particle_species, momentum, momenta, getindex # specific compute models @@ -54,7 +55,12 @@ include("interfaces/model_interface.jl") include("interfaces/process_interface.jl") include("interfaces/setup_interface.jl") -include("phase_spaces.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("momentum_generation.jl") include("propagators.jl") @@ -62,7 +68,6 @@ include("cross_section/diff_probability.jl") include("cross_section/diff_cross_section.jl") include("cross_section/total_probability.jl") include("cross_section/total_cross_section.jl") -include("cross_section/internal.jl") include("models/models.jl") include("processes/one_photon_compton/one_photon_compton.jl") diff --git a/src/cross_section/diff_cross_section.jl b/src/cross_section/diff_cross_section.jl index 1430586..0a32b1c 100644 --- a/src/cross_section/diff_cross_section.jl +++ b/src/cross_section/diff_cross_section.jl @@ -5,10 +5,6 @@ # cross sections based on the scattering process interface ######################## -function _incident_flux(psp::PhaseSpacePoint) - return _incident_flux(psp.proc, psp.model, momenta(psp, Incoming())) -end - """ unsafe_differential_cross_section(phase_space_point::PhaseSpacePoint) @@ -27,7 +23,7 @@ If the given phase spaces are physical, return differential cross section evalua """ function differential_cross_section(phase_space_point::PhaseSpacePoint) if !_is_in_phasespace(phase_space_point) - return zero(eltype(momentum(phase_space_point, Incoming(), 1))) + return zero(eltype(_momentum_type(phase_space_point))) end return unsafe_differential_cross_section(phase_space_point) diff --git a/src/cross_section/internal.jl b/src/cross_section/internal.jl deleted file mode 100644 index 609602e..0000000 --- a/src/cross_section/internal.jl +++ /dev/null @@ -1,131 +0,0 @@ - -# convenience function for internal use only -# differential probability without energy momentum conservation check -# single in phase space points/ single out phase space point -# based on four-momenta -function _unsafe_differential_probability( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - phase_space_def::AbstractPhasespaceDefinition, - in_phase_space::AbstractVector{T}, - out_phase_space::AbstractVector{T}, -) where {T<:QEDbase.AbstractFourMomentum} - psp = PhaseSpacePoint(proc, model, phase_space_def, in_phase_space, out_phase_space) - return unsafe_differential_probability(psp) -end - -# convenience function for internal use only -# differential probability without energy momentum conservation check -# based on coordinates -function _unsafe_differential_probability( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - phase_space_def::AbstractPhasespaceDefinition, - in_phase_space::AbstractVector{T}, - out_phase_space::AbstractVector{T}, -) where {T<:Real} - in_momenta, out_momenta = _generate_momenta( - proc, model, phase_space_def, in_phase_space, out_phase_space - ) - return _unsafe_differential_probability( - proc, model, phase_space_def, in_momenta, out_momenta - ) -end - -# convenience function for internal use only -# differential probability with energy momentum conservation check -# one in phase space point/ one out phase space point -# based on four-momenta -function _differential_probability( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - phase_space_def::AbstractPhasespaceDefinition, - in_phase_space::AbstractVector{T}, - out_phase_space::AbstractVector{T}, -) where {T<:QEDbase.AbstractFourMomentum} - psp = PhaseSpacePoint(proc, model, phase_space_def, in_phase_space, out_phase_space) - return differential_probability(psp) -end - -# convenience function for internal use only -# differential probability with energy momentum conservation check -# one in phase space point/ one out phase space point -# based on coordinates -function _differential_probability( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - phase_space_def::AbstractPhasespaceDefinition, - in_phase_space::AbstractVector{T}, - out_phase_space::AbstractVector{T}, -) where {T<:Real} - in_momenta, out_momenta = _generate_momenta( - proc, model, phase_space_def, in_phase_space, out_phase_space - ) - return _differential_probability(proc, model, phase_space_def, in_momenta, out_momenta) -end - -# convenience function for internal use only -# differential cross sections without energy momentum conservation check -# single in phase space point/ single out phase space point -# based on four-momenta -function _unsafe_differential_cross_section( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - phase_space_def::AbstractPhasespaceDefinition, - in_phase_space::AbstractVector{T}, - out_phase_space::AbstractVector{T}, -) where {T<:QEDbase.AbstractFourMomentum} - psp = PhaseSpacePoint(proc, model, phase_space_def, in_phase_space, out_phase_space) - return unsafe_differential_cross_section(psp) -end - -# convenience function for internal use only -# differential cross sections without energy momentum conservation check -# single in phase space point/ single out phase space point -# based on coordinates -function _unsafe_differential_cross_section( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - phase_space_def::AbstractPhasespaceDefinition, - in_phase_space::AbstractVector{T}, - out_phase_space::AbstractVector{T}, -) where {T<:Real} - in_momenta, out_momenta = _generate_momenta( - proc, model, phase_space_def, in_phase_space, out_phase_space - ) - return _unsafe_differential_cross_section( - proc, model, phase_space_def, in_momenta, out_momenta - ) -end - -# convenience function for internal use only -# differential cross sections with energy momentum conservation check -# single in phase space point/ single out phase space point -function _differential_cross_section( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - phase_space_def::AbstractPhasespaceDefinition, - in_phase_space::AbstractVector{T}, - out_phase_space::AbstractVector{T}, -) where {T<:QEDbase.AbstractFourMomentum} - psp = PhaseSpacePoint(proc, model, phase_space_def, in_phase_space, out_phase_space) - return differential_cross_section(psp) -end - -# convenience function for internal use only -# differential cross sections with energy momentum conservation check -# single in phase space point/ several out phase space points -function _differential_cross_section( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - phase_space_def::AbstractPhasespaceDefinition, - in_phase_space::AbstractVector{T}, - out_phase_space::AbstractVector{T}, -)::Float64 where {T<:Real} - in_momenta, out_momenta = _generate_momenta( - proc, model, phase_space_def, in_phase_space, out_phase_space - ) - return _differential_cross_section( - proc, model, phase_space_def, in_momenta, out_momenta - ) -end diff --git a/src/cross_section/total_cross_section.jl b/src/cross_section/total_cross_section.jl index a7574dd..f4a5ffa 100644 --- a/src/cross_section/total_cross_section.jl +++ b/src/cross_section/total_cross_section.jl @@ -3,85 +3,12 @@ # Total cross sections ############ -# total cross section on single phase space point -# based on four-momenta -function _total_cross_section( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - phase_space_def::AbstractPhasespaceDefinition, - in_phase_space::AbstractVector{T}, -) where {T<:QEDbase.AbstractFourMomentum} - I = 1 / (4 * _incident_flux(proc, model, in_phase_space)) - - return I * _total_probability(proc, model, phase_space_def, in_phase_space) -end - -# total cross section on single phase space point -# based on coordinates -function _total_cross_section( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - phase_space_def::AbstractPhasespaceDefinition, - in_phase_space::AbstractVector{T}, -) where {T<:Real} - in_momenta = _generate_incoming_momenta(proc, model, phase_space_def, in_phase_space) - return _total_cross_section(proc, model, phase_space_def, in_momenta) -end - -# total cross section on several phase space points -function _total_cross_section( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - in_phase_space_def::AbstractPhasespaceDefinition, - in_phase_space::AbstractMatrix{T}, -) where {T<:AbstractPhasespaceElement} - res = Vector{eltype(T)}(undef, size(in_phase_space, 2)) - for i in 1:size(in_phase_space, 2) - res[i] = _total_cross_section( - proc, model, in_phase_space_def, view(in_phase_space, :, i) - ) - end - return res -end - """ - total_cross_section( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - in_phase_space_def::AbstractPhasespaceDefinition, - in_phase_space::AbstractVecOrMat{T}, - ) where {T<:QEDbase.AbstractFourMomentum} + total_cross_section(in_psp::InPhaseSpacePoint) -Return the total cross section for a given combination of scattering process and compute model, evaluated at the particle momenta. +Return the total cross section for a given [`InPhaseSpacePoint`](@ref). """ -function total_cross_section( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - in_phase_space_def::AbstractPhasespaceDefinition, - in_phase_space::AbstractVecOrMat{T}, -) where {T<:QEDbase.AbstractFourMomentum} - _check_in_phase_space_dimension(proc, model, in_phase_space) - - return _total_cross_section(proc, model, in_phase_space_def, in_phase_space) -end - -""" - total_cross_section( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - in_phase_space_def::AbstractPhasespaceDefinition, - in_phase_space::AbstractVecOrMat{T}, - ) where {T<:Real} - -Return the total cross section for a given combination of scattering process and compute model, evaluated at the coordinates. -""" -function total_cross_section( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - in_phase_space_def::AbstractPhasespaceDefinition, - in_phase_space::AbstractVecOrMat{T}, -) where {T<:Real} - _check_in_phase_space_dimension(proc, model, in_phase_space) - - return _total_cross_section(proc, model, in_phase_space_def, in_phase_space) +function total_cross_section(in_psp::InPhaseSpacePoint) + I = 1 / (4 * _incident_flux(in_psp)) + return I * _total_probability(in_psp) end diff --git a/src/cross_section/total_probability.jl b/src/cross_section/total_probability.jl index 3296ba0..446900f 100644 --- a/src/cross_section/total_probability.jl +++ b/src/cross_section/total_probability.jl @@ -1,74 +1,12 @@ - ########### # Total probability ########### -# total probability on a phase space point -# based on coordinates -function _total_probability( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - phase_space_def::AbstractPhasespaceDefinition, - in_phase_space::AbstractVector{T}, -) where {T<:Real} - in_momenta = _generate_incoming_momenta(proc, model, phase_space_def, in_phase_space) - return _total_probability(proc, model, phase_space_def, in_momenta) -end - -# total probability on several phase space points -function _total_probability( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - phase_space_def::AbstractPhasespaceDefinition, - in_phase_space::AbstractMatrix{T}, -) where {T<:AbstractPhasespaceElement} - res = Vector{eltype(T)}(undef, size(in_phase_space, 2)) - for i in 1:size(in_phase_space, 2) - res[i] = _total_probability( - proc, model, phase_space_def, view(in_phase_space, :, i) - ) - end - return res -end - -""" - total_probability( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - phase_space_def::AbstractPhasespaceDefinition, - in_phase_space::AbstractMatrix{T}, - ) where {T<:QEDbase.AbstractFourMomentum} - -Return the total probability of a given model and process combination, evaluated at the particle momenta. -""" -function total_probability( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - phase_space_def::AbstractPhasespaceDefinition, - in_phase_space::AbstractVecOrMat{T}, -) where {T<:QEDbase.AbstractFourMomentum} - _check_in_phase_space_dimension(proc, model, in_phase_space) - - return _total_probability(proc, model, phase_space_def, in_phase_space) -end - """ - total_probability( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - phase_space_def::AbstractPhasespaceDefinition, - in_phase_space::AbstractMatrix{T}, - ) where {T<:Real} + total_probability(in_psp::InPhaseSpacePoint) -Return the total probability of a given model and process combination, evaluated at the coordinates. +Return the total probability of a given [`InPhaseSpacePoint`](@ref). """ -function total_probability( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - phase_space_def::AbstractPhasespaceDefinition, - in_phase_space::AbstractVecOrMat{T}, -) where {T<:Real} - _check_in_phase_space_dimension(proc, model, in_phase_space) - - return _total_probability(proc, model, phase_space_def, in_phase_space) +function total_probability(in_psp::InPhaseSpacePoint) + return _total_probability(in_psp) end diff --git a/src/interfaces/process_interface.jl b/src/interfaces/process_interface.jl index b02d3ad..e54513a 100644 --- a/src/interfaces/process_interface.jl +++ b/src/interfaces/process_interface.jl @@ -22,11 +22,7 @@ interface functions need to be implemented for every combination of `CustomProce `CustomModel<:AbstractModelDefinition`, and `CustomPhasespaceDefinition<:AbstractPhasespaceDefinition`. ```Julia - _incident_flux( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - in_phase_space::AbstractVector{T} - ) where {T<:QEDbase.AbstractFourMomentum} + _incident_flux(psp::InPhaseSpacePoint{CustomProcess,CustomModel}) _matrix_element(psp::PhaseSpacePoint{CustomProcess,CustomModel}) @@ -71,13 +67,12 @@ This function needs to be given to implement the scattering process interface. function outgoing_particles end """ - _incident_flux( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - in_phase_space::AbstractVector{T} - ) where {T<:QEDbase.AbstractFourMomentum} + _incident_flux(in_psp::InPhaseSpacePoint{PROC,MODEL}) where { + PROC <: AbstractProcessDefinition, + MODEL <: AbstractModelDefinition, + } -Interface function, which returns the incident flux of the given scattering process for a given incoming phase-space. +Interface function which returns the incident flux of the given scattering process for a given [`InPhaseSpacePoint`](@ref). """ function _incident_flux end @@ -88,14 +83,12 @@ function _incident_flux end MODEL <: AbstractModelDefinition, } -Interface function which returns a tuple of scattering matrix elements for each spin and polarisation combination of `proc`. +Interface function which returns a tuple of scattering matrix elements for each spin and polarization combination of `proc`. """ function _matrix_element end """ - _averaging_norm( - proc::AbstractProcessDefinition - ) + _averaging_norm(proc::AbstractProcessDefinition) Interface function, which returns a normalization for the averaging of the squared matrix elements over spins and polarizations. """ @@ -157,11 +150,31 @@ Return the number of outgoing particles of a given process. return length(outgoing_particles(proc_def)) end +""" + particles(proc_def::AbstractProcessDefinition, ::ParticleDirection) + +Convenience function dispatching to [`incoming_particles`](@ref) or [`outgoing_particles`](@ref) depending on the given direction. +""" +@inline particles(proc_def::AbstractProcessDefinition, ::Incoming) = + incoming_particles(proc_def) +@inline particles(proc_def::AbstractProcessDefinition, ::Outgoing) = + outgoing_particles(proc_def) + +""" + number_particles(proc_def::AbstractProcessDefinition, ::ParticleDirection) + +Convenience function dispatching to [`number_incoming_particles`](@ref) or [`number_outgoing_particles`](@ref) depending on the given direction. +""" +@inline number_particles(proc_def::AbstractProcessDefinition, ::Incoming) = + number_incoming_particles(proc_def) +@inline number_particles(proc_def::AbstractProcessDefinition, ::Outgoing) = + number_outgoing_particles(proc_def) + """ in_phase_space_dimension( proc::AbstractProcessDefinition, model::AbstractModelDefinition, - ) + ) TBW """ function in_phase_space_dimension end @@ -170,48 +183,23 @@ function in_phase_space_dimension end out_phase_space_dimension( proc::AbstractProcessDefinition, model::AbstractModelDefinition, - ) + ) TBW """ function out_phase_space_dimension end """ - _total_probability( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - ps_def::AbstractPhasespaceDefinition, - in_phase_space::AbstractVector{T} - ) where {T<: QEDbase.AbstractFourMomentum} + _total_probability(in_psp::InPhaseSpacePoint{PROC,MODEL}) where { + PROC <: AbstractProcessDefinition, + MODEL <: AbstractModelDefinition, + } Interface function for the combination of a scattering process and a physical model. Return the total of a -given process and model for a passed initial phase space definition and point. The elements of `in_phase_space`, -which represent the initial phase space, are the momenta of the respective particles. -The implementation of this function for a concrete process and model must not check if the length of -the passed initial phase spaces match the number of incoming particles. - -!!! note "probability interface" - - Given an implementation of this method, the following generic implementation without input check is provided: - - ```julia - - _total_probability(proc_def,model_def,in_phasespace::AbstractMatrix{T}) - - ``` - - where `T<:QEDbase.AbstractFourMomentum`, i.e. `_total_probability` is also implemented for a vector of initial phase space points. - Furthermore, a safe version of `_total_probability` is also implemented: [`total_probability`](@ref). +given process and model for a passed [`InPhaseSpacePoint`](@ref). !!! note "total cross section" - Given an implementaion of this method and [`_incident_flux`](@ref), the respective functions for the total cross section are available, - i.e. `_total_cross_section` (unsafe and not exported), and [`total_cross_section`](@ref), respectively. - -!!! note - - Each instance of this function does not check the validity of the input. - This function is not exported and should be used with caution. To add a method in order to implement the cross section interface, - it is recommented to directly use `QEDprocesses._total_cross_section` instead of globally `using QEDprocesses: _total_cross_section`. + Given an implementation of this method and [`_incident_flux`](@ref), the respective function for the total cross section [`total_cross_section`](@ref) is also available. """ function _total_probability end diff --git a/src/momentum_generation.jl b/src/momentum_generation.jl index 355bd4c..7e7d278 100644 --- a/src/momentum_generation.jl +++ b/src/momentum_generation.jl @@ -10,8 +10,8 @@ proc::AbstractProcessDefinition, model::AbstractModelDefinition, phase_space_def::AbstractPhasespaceDefinition, - in_phase_space::AbstractVector{T}, - ) where {T<:Real} + in_phase_space::NTuple{N,T}, + ) where {N,T<:Real} Interface function to generate the four-momenta of the incoming particles from coordinates for a given phase-space definition. """ @@ -22,8 +22,8 @@ function _generate_incoming_momenta end proc::AbstractProcessDefinition, model::AbstractModelDefinition, phase_space_def::AbstractPhasespaceDefinition, - out_phase_space::AbstractVector{T}, - ) where {T<:Real} + out_phase_space::NTuple{N,T}, + ) where {N,T<:Real} Interface function to generate the four-momenta of the outgoing particles from coordinates for a given phase-space definition. """ @@ -34,9 +34,9 @@ function _generate_outgoing_momenta end proc::AbstractProcessDefinition, model::AbstractModelDefinition, phase_space_def::AbstractPhasespaceDefinition, - in_phase_space::AbstractVector{T}, - out_phase_space::AbstractVector{T}, -) where {T<:Real} + in_phase_space::NTuple{N,T}, + out_phase_space::NTuple{M,T}, +) where {N,M,T<:Real} Return four-momenta for incoming and outgoing particles for given coordinate based phase-space points. """ @@ -44,9 +44,9 @@ function _generate_momenta( proc::AbstractProcessDefinition, model::AbstractModelDefinition, phase_space_def::AbstractPhasespaceDefinition, - in_phase_space::AbstractVector{T}, - out_phase_space::AbstractVector{T}, -) where {T<:Real} + in_phase_space::NTuple{N,T}, + out_phase_space::NTuple{M,T}, +) where {N,M,T<:Real} in_momenta = _generate_incoming_momenta(proc, model, phase_space_def, in_phase_space) out_momenta = _generate_outgoing_momenta(proc, model, phase_space_def, out_phase_space) diff --git a/src/phase_spaces.jl b/src/phase_spaces.jl deleted file mode 100644 index 5010c4f..0000000 --- a/src/phase_spaces.jl +++ /dev/null @@ -1,436 +0,0 @@ -##################### -# phase spaces -# -# This file contains a collection of types and functions to handle phase spaces -# for scattering processes. -# -# 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 - -abstract type AbstractFrameOfReference end -struct CenterOfMomentumFrame <: AbstractFrameOfReference end -struct ElectronRestFrame <: AbstractFrameOfReference end - -abstract type AbstractPhasespaceDefinition end - -Broadcast.broadcastable(ps_def::AbstractPhasespaceDefinition) = Ref(ps_def) - -""" - PhasespaceDefinition(coord_sys::AbstractCoordinateSystem, frame::AbstractFrameOfReference) - -Convenient type to dispatch on coordiante systems and frames of reference. -""" -struct PhasespaceDefinition{CS<:AbstractCoordinateSystem,F<:AbstractFrameOfReference} <: - AbstractPhasespaceDefinition - coord_sys::CS - frame::F -end - -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 - -# abstract type for generic phase spaces -# -# Currently, elements can be either four-momenta, or real numbers, -# i.e. coordinates. -AbstractPhasespaceElement = Union{AbstractFourMomentum,Real} - -# utility functions - -@inline function _check_in_phase_space_dimension( - proc::AbstractProcessDefinition, - ::AbstractModelDefinition, - in_phase_space::AbstractVecOrMat{T}, -) 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))>", - ), - ) -end - -@inline function _check_out_phase_space_dimension( - proc::AbstractProcessDefinition, - ::AbstractModelDefinition, - out_phase_space::AbstractVecOrMat{T}, -) 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))>", - ), - ) -end - -@inline function _check_in_phase_space_dimension( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - in_phase_space::AbstractVecOrMat{T}, -) 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))>", - ), - ) -end - -@inline function _check_out_phase_space_dimension( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - out_phase_space::AbstractVecOrMat{T}, -) 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))>", - ), - ) -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. - -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. - -```jldoctest -julia> using QEDbase; using QEDprocesses - -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<:ParticleDirection,SPECIES<:AbstractParticleType,ELEMENT<:AbstractFourMomentum -} <: AbstractParticle - dir::DIR - species::SPECIES - mom::ELEMENT - - function ParticleStateful( - dir::DIR, species::SPECIES, mom::ELEMENT - ) where { - DIR<:ParticleDirection,SPECIES<:AbstractParticleType,ELEMENT<:AbstractFourMomentum - } - return new{DIR,SPECIES,ELEMENT}(dir, species, mom) - end -end - -# particle interface -@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) - -# accessors -particle_direction(part::ParticleStateful) = part.dir -particle_species(part::ParticleStateful) = part.species -momentum(part::ParticleStateful) = part.mom - -# recursion termination: base case -@inline _assemble_tuple_type(::Tuple{}, ::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{AbstractParticleType}}, dir::DIR_T, ELTYPE::Type -) where {SPECIES_T<:AbstractParticleType,DIR_T<:ParticleDirection} - return ( - ParticleStateful{DIR_T,SPECIES_T,ELTYPE}, - _assemble_tuple_type(particle_types[2:end], dir, ELTYPE)..., - ) -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 - -""" - 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. - -```jldoctest -julia> using QEDprocesses; using QEDbase - -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] -``` -""" -struct PhaseSpacePoint{ - PROC<:AbstractProcessDefinition, - MODEL<:AbstractModelDefinition, - PSDEF<:AbstractPhasespaceDefinition, - IN_PARTICLES<:Tuple{Vararg{ParticleStateful}}, - OUT_PARTICLES<:Tuple{Vararg{ParticleStateful}}, -} - proc::PROC - model::MODEL - ps_def::PSDEF - - in_particles::IN_PARTICLES - out_particles::OUT_PARTICLES - - function PhaseSpacePoint( - proc::PROC, model::MODEL, ps_def::PSDEF, in_p::IN_PARTICLES, out_p::OUT_PARTICLES - ) where { - PROC<:AbstractProcessDefinition, - MODEL<:AbstractModelDefinition, - PSDEF<:AbstractPhasespaceDefinition, - IN_PARTICLES<:Tuple{Vararg{ParticleStateful}}, - OUT_PARTICLES<:Tuple{Vararg{ParticleStateful}}, - } - _check_psp(incoming_particles(proc), outgoing_particles(proc), in_p, out_p) - - return new{PROC,MODEL,PSDEF,IN_PARTICLES,OUT_PARTICLES}( - proc, model, ps_def, in_p, out_p - ) - end - - """ - PhaseSpacePoint( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - ps_def::AbstractPhasespaceDefinition, - in_ps::AbstractVector{ELEMENT}, - out_ps::AbstractVector{ELEMENT}, - ) where {ELEMENT<:QEDbase.AbstractFourMomentum} - - Construct the phase space point from given momenta of incoming and outgoing particles regarding a given process. - """ - function PhaseSpacePoint( - proc::PROC, - model::MODEL, - ps_def::PSDEF, - in_ps::AbstractVector{ELEMENT}, - out_ps::AbstractVector{ELEMENT}, - ) where { - PROC<:AbstractProcessDefinition, - MODEL<:AbstractModelDefinition, - PSDEF<:AbstractPhasespaceDefinition, - ELEMENT<:AbstractFourMomentum, - } - in_particles = Tuple{ - _assemble_tuple_type(incoming_particles(proc), Incoming(), ELEMENT)... - }( - ParticleStateful(Incoming(), particle, mom) for - (particle, mom) in zip(incoming_particles(proc), in_ps) - ) - - out_particles = Tuple{ - _assemble_tuple_type(outgoing_particles(proc), Outgoing(), ELEMENT)... - }( - ParticleStateful(Outgoing(), particle, mom) for - (particle, mom) in zip(outgoing_particles(proc), out_ps) - ) - - # no need to check anything since it is constructed correctly above - - return new{PROC,MODEL,PSDEF,typeof(in_particles),typeof(out_particles)}( - proc, model, ps_def, in_particles, out_particles - ) - end -end - -function momenta(psp::PhaseSpacePoint, ::Incoming) - return QEDbase._as_svec(momentum.(psp.in_particles)) -end - -function momenta(psp::PhaseSpacePoint, ::Outgoing) - return QEDbase._as_svec(momentum.(psp.out_particles)) -end - -# recursion termination: success -@inline _recursive_type_check(t::Tuple{}, p::Tuple{}, dir::ParticleDirection) = nothing - -# recursion termination: overload for unequal number of particles -@inline function _recursive_type_check( - ::Tuple{Vararg{ParticleStateful,N}}, - ::Tuple{Vararg{AbstractParticleType,M}}, - dir::ParticleDirection, -) where {N,M} - return throw( - InvalidInputError( - "the number of $(dir) particles in the process $(M) does not match number of particles in the input $(N)", - ), - ) -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{AbstractParticleType,M}}, - dir::DIR_T, -) where { - N, - M, - DIR_IN_T<:ParticleDirection, - DIR_T<:ParticleDirection, - SPECIES_IN_T<:AbstractParticleType, - SPECIES_T<:AbstractParticleType, -} - return throw( - InvalidInputError( - "expected $(dir) $(SPECIES_T()) but got $(DIR_IN_T()) $(SPECIES_IN_T())" - ), - ) -end - -@inline function _recursive_type_check( - t::Tuple{ParticleStateful{DIR_T,SPECIES_T},Vararg{ParticleStateful,N}}, - p::Tuple{SPECIES_T,Vararg{AbstractParticleType,N}}, - dir::DIR_T, -) where {N,DIR_T<:ParticleDirection,SPECIES_T<: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{AbstractParticleType}}, - P_OUT_Ts<:Tuple{Vararg{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, Incoming()) - _recursive_type_check(out_p, out_proc, Outgoing()) - return nothing -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 - -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/access.jl b/src/phase_spaces/access.jl new file mode 100644 index 0000000..6436be0 --- /dev/null +++ b/src/phase_spaces/access.jl @@ -0,0 +1,59 @@ +import QEDbase: + is_particle, + is_anti_particle, + is_fermion, + is_boson, + is_incoming, + is_outgoing, + mass, + charge + +# particle interface +@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) + +# accessors +particle_direction(part::ParticleStateful) = part.dir +particle_species(part::ParticleStateful) = part.species +momentum(part::ParticleStateful) = part.mom + +""" + momenta(psp::PhaseSpacePoint, ::ParticleDirection) + +Return a `Tuple` of all the particles' momenta for the given `ParticleDirection`. +""" +momenta(psp::PhaseSpacePoint, ::Incoming) = momentum.(psp.in_particles) +momenta(psp::PhaseSpacePoint, ::Outgoing) = momentum.(psp.out_particles) + +""" + 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/src/phase_spaces/create.jl b/src/phase_spaces/create.jl new file mode 100644 index 0000000..c8fa26f --- /dev/null +++ b/src/phase_spaces/create.jl @@ -0,0 +1,158 @@ +# PSP constructors from particle statefuls + +""" + InPhaseSpacePoint( + proc::AbstractProcessDefinition, + model::AbstractModelDefinition, + ps_def::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<:AbstractProcessDefinition, + MODEL<:AbstractModelDefinition, + PSDEF<:AbstractPhasespaceDefinition, + IN_PARTICLES<:Tuple{Vararg{ParticleStateful}}, +} + return PhaseSpacePoint(proc, model, ps_def, in_ps, ()) +end + +""" + OutPhaseSpacePoint( + proc::AbstractProcessDefinition, + model::AbstractModelDefinition, + ps_def::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<:AbstractProcessDefinition, + MODEL<:AbstractModelDefinition, + PSDEF<:AbstractPhasespaceDefinition, + OUT_PARTICLES<:Tuple{Vararg{ParticleStateful}}, +} + return PhaseSpacePoint(proc, model, ps_def, (), out_ps) +end + +# PSP constructors from momenta + +""" + PhaseSpacePoint( + proc::AbstractProcessDefinition, + model::AbstractModelDefinition, + ps_def::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::AbstractProcessDefinition, + model::AbstractModelDefinition, + ps_def::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, Incoming()) + out_particles = _build_particle_statefuls(proc, out_momenta, Outgoing()) + + return PhaseSpacePoint(proc, model, ps_def, in_particles, out_particles) +end + +""" + InPhaseSpacePoint( + proc::AbstractProcessDefinition, + model::AbstractModelDefinition, + ps_def::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::AbstractProcessDefinition, + model::AbstractModelDefinition, + ps_def::AbstractPhasespaceDefinition, + in_momenta::NTuple{N,ELEMENT}, +) where {N,ELEMENT<:QEDbase.AbstractFourMomentum} + in_particles = _build_particle_statefuls(proc, in_momenta, Incoming()) + + return PhaseSpacePoint(proc, model, ps_def, in_particles, ()) +end + +""" + OutPhaseSpacePoint( + proc::AbstractProcessDefinition, + model::AbstractModelDefinition, + ps_def::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::AbstractProcessDefinition, + model::AbstractModelDefinition, + ps_def::AbstractPhasespaceDefinition, + out_momenta::NTuple{N,ELEMENT}, +) where {N,ELEMENT<:QEDbase.AbstractFourMomentum} + out_particles = _build_particle_statefuls(proc, out_momenta, Outgoing()) + + return PhaseSpacePoint(proc, model, ps_def, (), out_particles) +end + +# PSP constructors from coordinates + +""" + PhaseSpacePoint( + proc::AbstractProcessDefinition, + model::AbstractModelDefinition, + ps_def::AbstractPhasespaceDefinition, + in_coords::NTuple{N,Real}, + out_coords::NTuple{M,Real}, + ) + +Construct a [`PhaseSpacePoint`](@ref) from given coordinates by using the [`_generate_momenta`](@ref) interface. +""" +function PhaseSpacePoint( + proc::AbstractProcessDefinition, + model::AbstractModelDefinition, + ps_def::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::AbstractProcessDefinition, + model::AbstractModelDefinition, + ps_def::AbstractPhasespaceDefinition, + in_coords::NTuple{N,Real}, + ) + +Construct a [`PhaseSpacePoint`](@ref) from given coordinates by using the [`_generate_momenta`](@ref) 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::AbstractProcessDefinition, + model::AbstractModelDefinition, + ps_def::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..0645975 --- /dev/null +++ b/src/phase_spaces/types.jl @@ -0,0 +1,177 @@ +abstract type AbstractCoordinateSystem end +struct SphericalCoordinateSystem <: AbstractCoordinateSystem end + +abstract type AbstractFrameOfReference end +struct CenterOfMomentumFrame <: AbstractFrameOfReference end +struct ElectronRestFrame <: AbstractFrameOfReference end + +abstract type AbstractPhasespaceDefinition end + +""" + PhasespaceDefinition(coord_sys::AbstractCoordinateSystem, frame::AbstractFrameOfReference) + +Convenient type to dispatch on coordiante systems and frames of reference. +""" +struct PhasespaceDefinition{CS<:AbstractCoordinateSystem,F<:AbstractFrameOfReference} <: + AbstractPhasespaceDefinition + coord_sys::CS + frame::F +end + +Broadcast.broadcastable(ps_def::AbstractPhasespaceDefinition) = Ref(ps_def) + +# abstract type for generic phase spaces +# +# Currently, elements can be either four-momenta, or real numbers, +# i.e. coordinates. +AbstractPhasespaceElement = Union{QEDbase.AbstractFourMomentum,Real} + +""" + 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::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. + +```jldoctest +julia> using QEDbase; using QEDprocesses + +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<:ParticleDirection, + SPECIES<:AbstractParticleType, + ELEMENT<:QEDbase.AbstractFourMomentum, +} <: AbstractParticle + dir::DIR + species::SPECIES + mom::ELEMENT + + function ParticleStateful( + dir::DIR, species::SPECIES, mom::ELEMENT + ) where { + DIR<:ParticleDirection, + SPECIES<: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 ([`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. + +```jldoctest +julia> using QEDprocesses; using QEDbase + +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<:AbstractProcessDefinition, + MODEL<:AbstractModelDefinition, + PSDEF<:AbstractPhasespaceDefinition, + IN_PARTICLES<:Tuple{Vararg{ParticleStateful}}, + OUT_PARTICLES<:Tuple{Vararg{ParticleStateful}}, + ELEMENT<:QEDbase.AbstractFourMomentum, +} + proc::PROC + model::MODEL + ps_def::PSDEF + + in_particles::IN_PARTICLES + out_particles::OUT_PARTICLES + + """ + PhaseSpacePoint( + proc::AbstractProcessDefinition, + model::AbstractModelDefinition, + ps_def::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<:AbstractProcessDefinition, + MODEL<:AbstractModelDefinition, + PSDEF<: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( + incoming_particles(proc), 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 [`_incident_flux`](@ref). 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..cba0710 --- /dev/null +++ b/src/phase_spaces/utility.jl @@ -0,0 +1,141 @@ +# recursion termination: base case +@inline _assemble_tuple_type(::Tuple{}, ::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{AbstractParticleType}}, dir::DIR_T, ELTYPE::Type +) where {SPECIES_T<:AbstractParticleType,DIR_T<: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{}, ::ParticleDirection) = nothing + +# recursion termination: overload for unequal number of particles +@inline function _recursive_type_check( + ::Tuple{Vararg{ParticleStateful,N}}, + ::Tuple{Vararg{AbstractParticleType,M}}, + dir::ParticleDirection, +) where {N,M} + throw(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{AbstractParticleType,N}}, + dir::DIR_T, +) where { + N, + DIR_IN_T<:ParticleDirection, + DIR_T<:ParticleDirection, + SPECIES_IN_T<:AbstractParticleType, + SPECIES_T<:AbstractParticleType, +} + throw( + 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{AbstractParticleType,N}}, + dir::DIR_T, +) where {N,DIR_T<:ParticleDirection,SPECIES_T<: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{AbstractParticleType}}, + P_OUT_Ts<:Tuple{Vararg{AbstractParticleType}}, + IN_Ts<:Tuple{Vararg{ParticleStateful}}, + OUT_Ts<:Tuple{}, +} + # specific overload for InPhaseSpacePoint + _recursive_type_check(in_p, in_proc, 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{AbstractParticleType}}, + P_OUT_Ts<:Tuple{Vararg{AbstractParticleType}}, + IN_Ts<:Tuple{}, + OUT_Ts<:Tuple{Vararg{ParticleStateful}}, +} + # specific overload for OutPhaseSpacePoint + _recursive_type_check(out_p, out_proc, 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{AbstractParticleType}}, + P_OUT_Ts<:Tuple{Vararg{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, Incoming()) + _recursive_type_check(out_p, out_proc, 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`. + +```jldoctest +julia> using QEDprocesses; using QEDbase; + +julia> psp = PhaseSpacePoint(Compton(), PerturbativeQED(), PhasespaceDefinition(SphericalCoordinateSystem(), ElectronRestFrame()), Tuple(rand(SFourMomentum) for _ in 1:2), Tuple(rand(SFourMomentum) for _ in 1:2)); + +julia> QEDprocesses._momentum_type(psp) +SFourMomentum + +julia> QEDprocesses._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::AbstractProcessDefinition, moms::NTuple{N,ELEMENT}, dir::ParticleDirection +) where {N,ELEMENT<:QEDbase.AbstractFourMomentum} + N == number_particles(proc, dir) || throw( + InvalidInputError( + "expected $(number_particles(proc, dir)) $(dir) particles for the process but got $(N)", + ), + ) + res = Tuple{_assemble_tuple_type(particles(proc, dir), dir, ELEMENT)...}( + ParticleStateful(dir, particle, mom) for + (particle, mom) in zip(particles(proc, dir), moms) + ) + + return res +end diff --git a/src/processes/one_photon_compton/perturbative/cross_section.jl b/src/processes/one_photon_compton/perturbative/cross_section.jl index eaafa8c..a451da9 100644 --- a/src/processes/one_photon_compton/perturbative/cross_section.jl +++ b/src/processes/one_photon_compton/perturbative/cross_section.jl @@ -3,10 +3,8 @@ # Implementation of the cross section interface ##### -function _incident_flux( - proc::Compton, model::PerturbativeQED, in_ps::AbstractVector{T} -) where {T<:QEDbase.AbstractFourMomentum} - return @inbounds in_ps[1] * in_ps[2] +function _incident_flux(in_psp::InPhaseSpacePoint{<:Compton,<:PerturbativeQED}) + return momentum(in_psp, Incoming(), 1) * momentum(in_psp, Incoming(), 2) end function _matrix_element(psp::PhaseSpacePoint{<:Compton,PerturbativeQED}) @@ -16,6 +14,7 @@ function _matrix_element(psp::PhaseSpacePoint{<:Compton,PerturbativeQED}) end """ + _averaging_norm(proc::Compton) !!! note "Convention" @@ -26,24 +25,31 @@ function _averaging_norm(proc::Compton) return inv(prod(normalizations)) end -function _all_onshell( - proc::Compton, in_ps::AbstractVector{T}, out_ps::AbstractVector{T} -) where {T<:QEDbase.AbstractFourMomentum} - sq_in_moms = getMass2.(in_ps) - sq_out_moms = getMass2.(out_ps) - sq_in_masses = mass.(incoming_particles(proc)) .^ 2 - sq_out_masses = mass.(outgoing_particles(proc)) .^ 2 - return isapprox(sq_in_moms, SVector(sq_in_masses)) && - isapprox(sq_out_moms, SVector(sq_out_masses)) +@inline function _all_onshell(psp::PhaseSpacePoint{<:Compton}) + return @inbounds isapprox( + getMass2(momentum(psp, Incoming(), 1)), mass(incoming_particles(psp.proc)[1])^2 + ) && + isapprox( + getMass2(momentum(psp, Incoming(), 2)), mass(incoming_particles(psp.proc)[2])^2 + ) && + isapprox( + getMass2(momentum(psp, Outgoing(), 1)), mass(outgoing_particles(psp.proc)[1])^2 + ) && + isapprox( + getMass2(momentum(psp, Outgoing(), 2)), mass(outgoing_particles(psp.proc)[2])^2 + ) end -function _is_in_phasespace(psp::PhaseSpacePoint{<:Compton,PerturbativeQED}) - in_ps = momenta(psp, Incoming()) - out_ps = momenta(psp, Outgoing()) - if (!isapprox(sum(in_ps), sum(out_ps))) +@inline function _is_in_phasespace(psp::PhaseSpacePoint{<:Compton,<:PerturbativeQED}) + @inbounds if ( + !isapprox( + momentum(psp, Incoming(), 1) + momentum(psp, Incoming(), 2), + momentum(psp, Outgoing(), 1) + momentum(psp, Outgoing(), 2), + ) + ) return false end - return _all_onshell(psp.proc, in_ps, out_ps) + return _all_onshell(psp) end @inline function _phase_space_factor(psp::PhaseSpacePoint{<:Compton,PerturbativeQED}) @@ -57,8 +63,8 @@ end ####### @inline function _pert_compton_matrix_element( - proc::Compton, in_ps::AbstractVector{T}, out_ps::AbstractVector{T} -) where {T<:QEDbase.AbstractFourMomentum} + proc::Compton, in_ps::NTuple{N,T}, out_ps::NTuple{M,T} +) where {N,M,T<:QEDbase.AbstractFourMomentum} in_electron_mom = in_ps[1] in_photon_mom = in_ps[2] out_electron_mom = out_ps[1] diff --git a/src/processes/one_photon_compton/perturbative/kinematics.jl b/src/processes/one_photon_compton/perturbative/kinematics.jl index 75fa0fb..a299726 100644 --- a/src/processes/one_photon_compton/perturbative/kinematics.jl +++ b/src/processes/one_photon_compton/perturbative/kinematics.jl @@ -6,9 +6,9 @@ function generate_momenta( proc::Compton, model::PerturbativeQED, in_ps_def::PhasespaceDefinition{SphericalCoordinateSystem,ElectronRestFrame}, - in_ps::AbstractVector{T}, - out_ps::AbstractVector{T}, -) where {T<:Real} + in_ps::NTuple{N,T}, + out_ps::NTuple{M,T}, +) where {N,M,T<:Real} return _generate_momenta(proc, model, in_ps_def, in_ps, out_ps) end @@ -16,29 +16,29 @@ function _generate_incoming_momenta( proc::Compton, model::PerturbativeQED, in_ps_def::PhasespaceDefinition{SphericalCoordinateSystem,ElectronRestFrame}, - in_ps::AbstractVector{T}, -) where {T<:Real} + in_ps::NTuple{N,T}, +) where {N,T<:Real} om = in_ps[1] P = SFourMomentum(one(om), zero(om), zero(om), zero(om)) K = SFourMomentum(om, zero(om), zero(om), om) - return [P, K] + return P, K end function _generate_momenta( proc::Compton, model::PerturbativeQED, in_ps_def::PhasespaceDefinition{SphericalCoordinateSystem,ElectronRestFrame}, - in_ps::AbstractVector{T}, - out_ps::AbstractVector{T}, -) where {T<:Real} + in_ps::NTuple{N,T}, + out_ps::NTuple{M,T}, +) where {N,M,T<:Real} omega = in_ps[1] cth = out_ps[1] phi = out_ps[2] P, K, Pp, Kp = _generate_momenta_elab_sph(omega, cth, phi) # TODO: do this coord and frame dependent - in_moms = SVector(P, K) - out_moms = SVector(Pp, Kp) + in_moms = (P, K) + out_moms = (Pp, Kp) return in_moms, out_moms end diff --git a/src/processes/one_photon_compton/perturbative/total_probability.jl b/src/processes/one_photon_compton/perturbative/total_probability.jl index 4ec00ed..6699b17 100644 --- a/src/processes/one_photon_compton/perturbative/total_probability.jl +++ b/src/processes/one_photon_compton/perturbative/total_probability.jl @@ -1,15 +1,10 @@ -function _total_probability( - proc::Compton, - model::PerturbativeQED, - in_phase_space_def::AbstractPhasespaceDefinition, - in_phase_space::AbstractVector{T}, -) where {T<:QEDbase.AbstractFourMomentum} - omega = getE(in_phase_space[2]) +function _total_probability(in_psp::InPhaseSpacePoint{<:Compton,<:PerturbativeQED}) + omega = getE(momentum(in_psp[Incoming(), 2])) function func(x) - return _unsafe_differential_probability( - proc, model, in_phase_space_def, [omega], [x, 0.0] + return unsafe_differential_probability( + PhaseSpacePoint(in_psp.proc, in_psp.model, in_psp.ps_def, (omega,), (x, 0.0)) ) end diff --git a/test/cross_sections.jl b/test/cross_sections.jl index 505514e..bee260c 100644 --- a/test/cross_sections.jl +++ b/test/cross_sections.jl @@ -63,24 +63,19 @@ TESTPSDEF = TestImplementation.TestPhasespaceDef() @testset "total cross section" begin @testset "compute" begin - groundtruth = TestImplementation._groundtruth_total_cross_section(p_in_phys) - totCS_on_moms = total_cross_section(TESTPROC, TESTMODEL, TESTPSDEF, p_in_phys) - COORDS_IN = TestImplementation.flat_components(p_in_phys) - totCS_on_coords = total_cross_section(TESTPROC, TESTMODEL, TESTPSDEF, COORDS_IN) - - @test isapprox(totCS_on_moms, groundtruth, atol=ATOL, rtol=RTOL) - @test isapprox(totCS_on_coords, groundtruth, atol=ATOL, rtol=RTOL) - end - @testset "invalid input" begin - @test_throws DimensionMismatch total_cross_section( - TESTPROC, TESTMODEL, TESTPSDEF, p_in_phys_invalid - ) - COORDS_IN = TestImplementation.flat_components(p_in_phys_invalid) - @test_throws DimensionMismatch total_cross_section( + IN_PS_POINT = InPhaseSpacePoint(TESTPROC, TESTMODEL, TESTPSDEF, p_in_phys) + IN_PS_POINT_COORDS = InPhaseSpacePoint( TESTPROC, TESTMODEL, TESTPSDEF, COORDS_IN ) + + groundtruth = TestImplementation._groundtruth_total_cross_section(p_in_phys) + totCS_on_moms = total_cross_section(IN_PS_POINT) + totCS_on_coords = total_cross_section(IN_PS_POINT_COORDS) + + @test isapprox(totCS_on_moms, groundtruth, atol=ATOL, rtol=RTOL) + @test isapprox(totCS_on_coords, groundtruth, atol=ATOL, rtol=RTOL) end end @@ -110,24 +105,19 @@ TESTPSDEF = TestImplementation.TestPhasespaceDef() @testset "total probability" begin @testset "compute" begin - groundtruth = TestImplementation._groundtruth_total_probability(p_in_phys) - totCS_on_moms = total_probability(TESTPROC, TESTMODEL, TESTPSDEF, p_in_phys) - COORDS_IN = TestImplementation.flat_components(p_in_phys) - totCS_on_coords = total_probability(TESTPROC, TESTMODEL, TESTPSDEF, COORDS_IN) - - @test isapprox(totCS_on_moms, groundtruth, atol=ATOL, rtol=RTOL) - @test isapprox(totCS_on_coords, groundtruth, atol=ATOL, rtol=RTOL) - end - @testset "invalid input" begin - @test_throws DimensionMismatch total_probability( - TESTPROC, TESTMODEL, TESTPSDEF, p_in_phys_invalid - ) - COORDS_IN = TestImplementation.flat_components(p_in_phys_invalid) - @test_throws DimensionMismatch total_probability( + IN_PS_POINT = InPhaseSpacePoint(TESTPROC, TESTMODEL, TESTPSDEF, p_in_phys) + IN_PS_POINT_COORDS = InPhaseSpacePoint( TESTPROC, TESTMODEL, TESTPSDEF, COORDS_IN ) + + groundtruth = TestImplementation._groundtruth_total_probability(p_in_phys) + totCS_on_moms = total_probability(IN_PS_POINT) + totCS_on_coords = total_probability(IN_PS_POINT_COORDS) + + @test isapprox(totCS_on_moms, groundtruth, atol=ATOL, rtol=RTOL) + @test isapprox(totCS_on_coords, groundtruth, atol=ATOL, rtol=RTOL) end end end diff --git a/test/interfaces/process_interface.jl b/test/interfaces/process_interface.jl index 6da09b1..28b29db 100644 --- a/test/interfaces/process_interface.jl +++ b/test/interfaces/process_interface.jl @@ -71,9 +71,20 @@ include("../test_implementation/TestImplementation.jl") end @testset "incident flux" begin - test_incident_flux = QEDprocesses._incident_flux(TESTPROC, TESTMODEL, IN_PS) + test_incident_flux = QEDprocesses._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 = QEDprocesses._incident_flux( + PhaseSpacePoint(TESTPROC, TESTMODEL, TESTPSDEF, IN_PS, OUT_PS) + ) + @test isapprox(test_incident_flux, groundtruth, atol=ATOL, rtol=RTOL) + + @test_throws MethodError QEDprocesses._incident_flux( + OutPhaseSpacePoint(TESTPROC, TESTMODEL, TESTPSDEF, OUT_PS) + ) end @testset "averaging norm" begin @@ -94,10 +105,8 @@ include("../test_implementation/TestImplementation.jl") @testset "is in phasespace" begin @test QEDprocesses._is_in_phasespace(PSP) - IN_PS_unphysical = deepcopy(IN_PS) - IN_PS_unphysical[1] = SFourMomentum(zeros(4)) - OUT_PS_unphysical = deepcopy(OUT_PS) - OUT_PS_unphysical[end] = ones(SFourMomentum) + 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 ) diff --git a/test/phase_spaces.jl b/test/phase_spaces.jl index 40407fd..0bb290d 100644 --- a/test/phase_spaces.jl +++ b/test/phase_spaces.jl @@ -120,6 +120,22 @@ end @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 momentum(psp, Incoming(), -1) @@ -149,13 +165,9 @@ end ) end - @testset "Generation" begin + @testset "Generation from momenta" begin test_psp = PhaseSpacePoint( - process, - model, - phasespace_def, - SVector(in_el_mom, in_ph_mom), - SVector(out_el_mom, out_ph_mom), + process, model, phasespace_def, (in_el_mom, in_ph_mom), (out_el_mom, out_ph_mom) ) @test test_psp.proc == process @@ -167,6 +179,58 @@ end @test test_psp[Outgoing(), 1] == out_el @test test_psp[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 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 InvalidInputError InPhaseSpacePoint( + process, model, phasespace_def, in_particles_invalid + ) + @test_throws InvalidInputError OutPhaseSpacePoint( + process, model, phasespace_def, out_particles_invalid + ) + + @testset "Error handling from momenta" for i in [1, 3, 4, 5] + @test_throws InvalidInputError InPhaseSpacePoint( + process, model, phasespace_def, TestImplementation._rand_momenta(RNG, i) + ) + @test_throws InvalidInputError OutPhaseSpacePoint( + process, model, phasespace_def, TestImplementation._rand_momenta(RNG, i) + ) + end + end end @testset "Coordinate System" begin diff --git a/test/processes/one_photon_compton/perturbative.jl b/test/processes/one_photon_compton/perturbative.jl index 27f06ce..248e2b1 100644 --- a/test/processes/one_photon_compton/perturbative.jl +++ b/test/processes/one_photon_compton/perturbative.jl @@ -32,8 +32,8 @@ end @testset "momentum generation" begin @testset "$om, $cth, $phi" for (om, cth, phi) in Iterators.product(OMEGAS, COS_THETAS, PHIS) - IN_COORDS = [om] - OUT_COORDS = [cth, phi] + IN_COORDS = (om,) + OUT_COORDS = (cth, phi) IN_PS, OUT_PS = QEDprocesses._generate_momenta( PROC, MODEL, PS_DEF, IN_COORDS, OUT_COORDS ) @@ -41,8 +41,10 @@ end out_mom_square = getMass2.(OUT_PS) in_masses = mass.(incoming_particles(PROC)) .^ 2 out_masses = mass.(outgoing_particles(PROC)) .^ 2 - @test isapprox(in_mom_square, SVector(in_masses)) - @test isapprox(out_mom_square, SVector(out_masses)) + + # we need a larger ATOL than eps() here because the error is accumulated over several additions + @test all(isapprox.(in_mom_square, in_masses, atol=4 * ATOL, rtol=RTOL)) + @test all(isapprox.(out_mom_square, out_masses, atol=4 * ATOL, rtol=RTOL)) end end end @@ -55,8 +57,8 @@ end @testset "$cos_theta $phi" for (cos_theta, phi) in Iterators.product(COS_THETAS, PHIS) - IN_COORDS = [omega] - OUT_COORDS = [cos_theta, phi] + IN_COORDS = (omega,) + OUT_COORDS = (cos_theta, phi) IN_PS, OUT_PS = QEDprocesses._generate_momenta( PROC, MODEL, PS_DEF, IN_COORDS, OUT_COORDS ) @@ -77,8 +79,8 @@ end @testset "$cos_theta $phi" for (cos_theta, phi) in Iterators.product(COS_THETAS, PHIS) - IN_COORDS = [omega] - OUT_COORDS = [cos_theta, phi] + IN_COORDS = (omega,) + OUT_COORDS = (cos_theta, phi) IN_PS, OUT_PS = QEDprocesses._generate_momenta( PROC, MODEL, PS_DEF, IN_COORDS, OUT_COORDS ) @@ -99,8 +101,8 @@ end @testset "$cos_theta $phi" for (cos_theta, phi) in Iterators.product(COS_THETAS, PHIS) - IN_COORDS = [omega] - OUT_COORDS = [cos_theta, phi] + IN_COORDS = (omega,) + OUT_COORDS = (cos_theta, phi) IN_PS, OUT_PS = QEDprocesses._generate_momenta( PROC, MODEL, PS_DEF, IN_COORDS, OUT_COORDS ) @@ -122,8 +124,10 @@ end # Klein-Nishina: total cross section function klein_nishina_total_cross_section(in_ps) function func(x) - return QEDprocesses._unsafe_differential_cross_section( - Compton(), PerturbativeQED(), PS_DEF, in_ps, [x, 0.0] + return unsafe_differential_cross_section( + PhaseSpacePoint( + Compton(), PerturbativeQED(), PS_DEF, in_ps, (x, 0.0) + ), ) end res, err = quadgk(func, -1, 1) @@ -132,10 +136,30 @@ end return 2 * pi * res end - IN_COORDS = [omega] + IN_COORDS = (omega,) groundtruth = klein_nishina_total_cross_section(IN_COORDS) - test_val = @inferred total_cross_section(PROC, MODEL, PS_DEF, IN_COORDS) + test_val = @inferred total_cross_section( + InPhaseSpacePoint(PROC, MODEL, PS_DEF, IN_COORDS) + ) @test isapprox(test_val, groundtruth, atol=ATOL, rtol=RTOL) + + @testset "$cos_theta $phi" for (cos_theta, phi) in + Iterators.product(COS_THETAS, PHIS) + OUT_COORDS = (cos_theta, phi) + + test_val = @inferred total_cross_section( + PhaseSpacePoint(PROC, MODEL, PS_DEF, IN_COORDS, OUT_COORDS) + ) + @test isapprox(test_val, groundtruth, atol=ATOL, rtol=RTOL) + + out_moms = momenta( + PhaseSpacePoint(PROC, MODEL, PS_DEF, IN_COORDS, OUT_COORDS), + Outgoing(), + ) + @test_throws MethodError total_cross_section( + OutPhaseSpacePoint(PROC, MODEL, PS_DEF, out_moms) + ) + end end end end diff --git a/test/test_implementation/groundtruths.jl b/test/test_implementation/groundtruths.jl index 78731ba..f05a010 100644 --- a/test/test_implementation/groundtruths.jl +++ b/test/test_implementation/groundtruths.jl @@ -1,3 +1,4 @@ +import QEDbase.AbstractFourMomentum """ _groundtruth_incident_flux(in_ps) @@ -237,28 +238,28 @@ end 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::AbstractVector) +function _groundtruth_total_probability( + in_ps::NTuple{N,T} +) where {N,T<:AbstractFourMomentum} Ptot = sum(in_ps) return Ptot * Ptot end -function _groundtruth_total_probability(in_ps::AbstractMatrix) - res = Vector{Float64}(undef, size(in_ps, 2)) - for i in 1:size(in_ps, 2) - res[i] = _groundtruth_total_probability(view(in_ps, :, i)) - end - return res +function _groundtruth_total_probability( + in_pss::Vector{NTuple{N,T}} +) where {N,T<:AbstractFourMomentum} + return _groundtruth_total_probability.(in_pss) end -function _groundtruth_total_cross_section(in_ps) +function _groundtruth_total_cross_section( + in_ps::NTuple{N,T} +) where {N,T<:AbstractFourMomentum} init_flux = _groundtruth_incident_flux(in_ps) return _groundtruth_total_probability(in_ps) / (4 * init_flux) end -function _groundtruth_total_cross_section(in_ps::AbstractMatrix) - res = Vector{Float64}(undef, size(in_ps, 2)) - for i in 1:size(in_ps, 2) - res[i] = _groundtruth_total_cross_section(view(in_ps, :, i)) - end - return res +function _groundtruth_total_cross_section( + in_pss::Vector{NTuple{N,T}} +) where {N,T<: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 index 1972c16..b631551 100644 --- a/test/test_implementation/random_momenta.jl +++ b/test/test_implementation/random_momenta.jl @@ -1,88 +1,83 @@ """ -Return a vector of random four momenta, i.e. a random phase space point +Return a tuple of random four momenta, i.e. a random phase space point. """ -function _rand_momenta(rng::AbstractRNG, N) - moms = Vector{SFourMomentum}(undef, N) - for i in 1:N - moms[i] = SFourMomentum(rand(rng, 4)) - end - return moms +function _rand_momenta(rng::AbstractRNG, n) + return NTuple{n,SFourMomentum}(SFourMomentum(rand(rng, 4)) for _ in 1:n) end """ -Return a matrix of random four momenta, i.e. a collection of phase space points +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 = Matrix{SFourMomentum}(undef, N1, N2) - for i in 1:N1 - for j in 1:N2 - moms[i, j] = SFourMomentum(rand(rng, 4)) - end +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, which is failing the incoming phase space constraint, -i.e. the first entry of the vector is the null momentum. +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) - moms = _rand_momenta(rng, N) - moms[1] = zero(SFourMomentum) - return moms +function _rand_in_momenta_failing(rng::AbstractRNG, n) + return (zero(SFourMomentum), _rand_momenta(rng, n - 1)...) end """ -Return a random phase space point, which is failing the outgoing phase space constraint, -i.e. the last entry of the vector is the unit momentum. +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) - moms = _rand_momenta(rng, N) - moms[end] = ones(SFourMomentum) - return moms +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 matrix is the null momentum, but the others pass. +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, 1] = zero(SFourMomentum) +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, +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 = _rand_momenta(rng, N1, N2) - for n in 1:N2 - moms[1, n] = zero(SFourMomentum) +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 collection of outgoing phase space points, where the first point is failing the phase space constraint, -i.e. the last entry of the matrix is the unit momentum, but the others pass. +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, 1] = ones(SFourMomentum) +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 collection of outgoing phase space points, where all points are failing the phase space constraint, +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 = _rand_momenta(rng, N1, N2) - for n in 1:N2 - moms[end, n] = ones(SFourMomentum) +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_process.jl b/test/test_implementation/test_process.jl index 003d14e..9c7fd5d 100644 --- a/test/test_implementation/test_process.jl +++ b/test/test_implementation/test_process.jl @@ -5,7 +5,6 @@ struct TestParticleBoson <: BosonLike end const PARTICLE_SET = [TestParticleFermion(), TestParticleBoson()] """ - TestProcess(rng,incoming_particles,outgoing_particles) """ @@ -79,10 +78,8 @@ struct TestPhasespaceDef_FAIL <: AbstractPhasespaceDefinition end # dummy implementation of the process interface -function QEDprocesses._incident_flux( - ::TestProcess, ::TestModel, in_ps::AbstractVector{T} -) where {T<:QEDbase.AbstractFourMomentum} - return _groundtruth_incident_flux(in_ps) +function QEDprocesses._incident_flux(in_psp::InPhaseSpacePoint{<:TestProcess,<:TestModel}) + return _groundtruth_incident_flux(momenta(in_psp, Incoming())) end function QEDprocesses._averaging_norm(proc::TestProcess) @@ -113,8 +110,8 @@ function QEDprocesses._generate_incoming_momenta( proc::TestProcess, model::TestModel, phase_space_def::TestPhasespaceDef, - in_phase_space::AbstractVector{T}, -) where {T<:Real} + in_phase_space::NTuple{N,T}, +) where {N,T<:Real} return _groundtruth_generate_momenta(in_phase_space) end @@ -122,13 +119,13 @@ function QEDprocesses._generate_outgoing_momenta( proc::TestProcess, model::TestModel, phase_space_def::TestPhasespaceDef, - out_phase_space::AbstractVector{T}, -) where {T<:Real} + out_phase_space::NTuple{N,T}, +) where {N,T<:Real} return _groundtruth_generate_momenta(out_phase_space) end function QEDprocesses._total_probability( - proc::TestProcess, model::TestModel, ps_def::TestPhasespaceDef, in_ps::AbstractVector{T} -) where {T<:QEDbase.AbstractFourMomentum} - return _groundtruth_total_probability(in_ps) + in_psp::InPhaseSpacePoint{<:TestProcess,<:TestModel,<:TestPhasespaceDef} +) + return _groundtruth_total_probability(momenta(in_psp, Incoming())) end diff --git a/test/test_implementation/utils.jl b/test/test_implementation/utils.jl index 194ca9e..3a141d9 100644 --- a/test/test_implementation/utils.jl +++ b/test/test_implementation/utils.jl @@ -17,7 +17,8 @@ function _unroll_moms(ps_moms::AbstractMatrix{T}) where {T<:QEDbase.AbstractFour return res end -flat_components(moms) = _unroll_moms(moms) +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} @@ -37,3 +38,7 @@ function _furl_moms(ps_coords::AbstractMatrix{T}) where {T<:Real} end return res end + +function _furl_moms(moms::NTuple{N,Float64}) where {N} + return Tuple(_furl_moms(Vector{Float64}([moms...]))) +end