diff --git a/src/QEDprocesses.jl b/src/QEDprocesses.jl index 755a8e3..5aaf9cf 100644 --- a/src/QEDprocesses.jl +++ b/src/QEDprocesses.jl @@ -35,8 +35,7 @@ export AbstractCoordinateSystem, SphericalCoordinateSystem export AbstractFrameOfReference, CenterOfMomentumFrame, ElectronRestFrame export AbstractPhasespaceDefinition, PhasespaceDefinition export ParticleStateful, PhaseSpacePoint -export spin, polarization, particle_direction, particle_species, momentum, getindex -export generate_phase_space +export spin, polarization, particle_direction, particle_species, momentum, momenta, getindex # specific compute models export PerturbativeQED diff --git a/src/cross_section/diff_cross_section.jl b/src/cross_section/diff_cross_section.jl index 24f09bd..1430586 100644 --- a/src/cross_section/diff_cross_section.jl +++ b/src/cross_section/diff_cross_section.jl @@ -6,7 +6,7 @@ ######################## function _incident_flux(psp::PhaseSpacePoint) - return _incident_flux(psp.proc, psp.model, momentum.(psp.in_particles)) + return _incident_flux(psp.proc, psp.model, momenta(psp, Incoming())) end """ diff --git a/src/cross_section/internal.jl b/src/cross_section/internal.jl index 2dc2a1b..609602e 100644 --- a/src/cross_section/internal.jl +++ b/src/cross_section/internal.jl @@ -10,9 +10,7 @@ function _unsafe_differential_probability( in_phase_space::AbstractVector{T}, out_phase_space::AbstractVector{T}, ) where {T<:QEDbase.AbstractFourMomentum} - psp = generate_phase_space( - proc, model, phase_space_def, in_phase_space, out_phase_space - ) + psp = PhaseSpacePoint(proc, model, phase_space_def, in_phase_space, out_phase_space) return unsafe_differential_probability(psp) end @@ -45,9 +43,7 @@ function _differential_probability( in_phase_space::AbstractVector{T}, out_phase_space::AbstractVector{T}, ) where {T<:QEDbase.AbstractFourMomentum} - psp = generate_phase_space( - proc, model, phase_space_def, in_phase_space, out_phase_space - ) + psp = PhaseSpacePoint(proc, model, phase_space_def, in_phase_space, out_phase_space) return differential_probability(psp) end @@ -79,9 +75,7 @@ function _unsafe_differential_cross_section( in_phase_space::AbstractVector{T}, out_phase_space::AbstractVector{T}, ) where {T<:QEDbase.AbstractFourMomentum} - psp = generate_phase_space( - proc, model, phase_space_def, in_phase_space, out_phase_space - ) + psp = PhaseSpacePoint(proc, model, phase_space_def, in_phase_space, out_phase_space) return unsafe_differential_cross_section(psp) end @@ -114,9 +108,7 @@ function _differential_cross_section( in_phase_space::AbstractVector{T}, out_phase_space::AbstractVector{T}, ) where {T<:QEDbase.AbstractFourMomentum} - psp = generate_phase_space( - proc, model, phase_space_def, in_phase_space, out_phase_space - ) + psp = PhaseSpacePoint(proc, model, phase_space_def, in_phase_space, out_phase_space) return differential_cross_section(psp) end diff --git a/src/interfaces/process_interface.jl b/src/interfaces/process_interface.jl index c432c1f..b02d3ad 100644 --- a/src/interfaces/process_interface.jl +++ b/src/interfaces/process_interface.jl @@ -140,7 +140,6 @@ function _phase_space_factor end ####################### """ - number_incoming_particles(proc_def::AbstractProcessDefinition) Return the number of incoming particles of a given process. @@ -150,7 +149,6 @@ Return the number of incoming particles of a given process. end """ - number_outgoing_particles(proc_def::AbstractProcessDefinition) Return the number of outgoing particles of a given process. diff --git a/src/phase_spaces.jl b/src/phase_spaces.jl index 48e6c51..5010c4f 100644 --- a/src/phase_spaces.jl +++ b/src/phase_spaces.jl @@ -139,44 +139,34 @@ Representation of a particle with a state. It has four fields: - `dir::ParticleDirection`: The direction of the particle, `QEDbase.Incoming()` or `QEDbase.Outgoing()`. - `species::AbstractParticleType`: The species of the particle, `QEDbase.Electron()`, `QEDbase.Positron()` etc. - `mom::AbstractFourMomentum`: The momentum of the particle. -- `spin_or_pol::AbstractSpinOrPolarization`: The spin or polarization of the particle, `QEDbase.SpinUp()`, `QEDbase.PolX() etc. Can only use spins with fermions and polarizations with bosons. `AllSpin` or `AllPol` by default. Overloads for `QEDbase.is_fermion`, `QEDbase.is_boson`, `QEDbase.is_particle`, `QEDbase.is_anti_particle`, `QEDbase.is_incoming`, `QEDbase.is_outgoing`, `QEDbase.mass`, and `QEDbase.charge` are provided, delegating the call to the correct field and thus implementing the `QEDbase.AbstractParticle` interface. -The legality of the combination of `species` and `spin_or_pol` is checked on construction. If, for example, the construction of an `Electron()` with a polarization is attempted, an [`InvalidInputError`](@ref) is thrown. - ```jldoctest julia> using QEDbase; using QEDprocesses julia> ParticleStateful(Incoming(), Electron(), SFourMomentum(1, 0, 0, 0)) ParticleStateful: incoming electron - spin: all spins momentum: [1.0, 0.0, 0.0, 0.0] -julia> ParticleStateful(Outgoing(), Photon(), SFourMomentum(1, 0, 0, 0), PolX()) +julia> ParticleStateful(Outgoing(), Photon(), SFourMomentum(1, 0, 0, 0)) ParticleStateful: outgoing photon - polarization: x-polarized momentum: [1.0, 0.0, 0.0, 0.0] ``` """ -struct ParticleStateful{ElType<:AbstractFourMomentum} <: AbstractParticle - dir::ParticleDirection - species::AbstractParticleType - mom::ElType - spin_or_pol::AbstractSpinOrPolarization - - function ParticleStateful( - dir::ParticleDirection, species::Species, mom::ElType, spin::Spin=AllSpin() - ) where {Species<:FermionLike,ElType<:AbstractFourMomentum,Spin<:AbstractSpin} - # constructor for fermions with spin - return new{ElType}(dir, species, mom, spin) - end +struct ParticleStateful{ + DIR<:ParticleDirection,SPECIES<:AbstractParticleType,ELEMENT<:AbstractFourMomentum +} <: AbstractParticle + dir::DIR + species::SPECIES + mom::ELEMENT function ParticleStateful( - dir::ParticleDirection, species::Species, mom::ElType, pol::Pol=AllPol() - ) where {Species<:BosonLike,ElType<:AbstractFourMomentum,Pol<:AbstractPolarization} - # constructor for bosons with polarization - return new{ElType}(dir, species, mom, pol) + dir::DIR, species::SPECIES, mom::ELEMENT + ) where { + DIR<:ParticleDirection,SPECIES<:AbstractParticleType,ELEMENT<:AbstractFourMomentum + } + return new{DIR,SPECIES,ELEMENT}(dir, species, mom) end end @@ -195,27 +185,26 @@ particle_direction(part::ParticleStateful) = part.dir particle_species(part::ParticleStateful) = part.species momentum(part::ParticleStateful) = part.mom -@inline _spin(::Species, particle::ParticleStateful) where {Species<:FermionLike} = - particle.spin_or_pol -@inline spin(particle::ParticleStateful) = _spin(particle.species, particle) +# recursion termination: base case +@inline _assemble_tuple_type(::Tuple{}, ::ParticleDirection, ::Type) = () -@inline _polarization(::Species, particle::ParticleStateful) where {Species<:BosonLike} = - particle.spin_or_pol -@inline polarization(particle::ParticleStateful) = _polarization(particle.species, particle) +# 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.spin_or_pol)): $(particle.mom)" - ) + 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, - " $(particle.spin_or_pol isa AbstractSpin ? "spin" : "polarization"): $(particle.spin_or_pol)", - ) println(io, " momentum: $(particle.mom)") return nothing end @@ -234,94 +223,170 @@ 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 (all spins): [1.0, 0.0, 0.0, 0.0] - -> incoming photon (all polarizations): [1.0, 0.0, 0.0, 0.0] + -> incoming electron: [1.0, 0.0, 0.0, 0.0] + -> incoming photon: [1.0, 0.0, 0.0, 0.0] outgoing particles: - -> outgoing electron (all spins): [1.0, 0.0, 0.0, 0.0] - -> outgoing photon (all polarizations): [1.0, 0.0, 0.0, 0.0] + -> 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, - PhaseSpaceElementType<:AbstractFourMomentum, - N_IN_PARTICLES, - N_OUT_PARTICLES, + IN_PARTICLES<:Tuple{Vararg{ParticleStateful}}, + OUT_PARTICLES<:Tuple{Vararg{ParticleStateful}}, } proc::PROC model::MODEL ps_def::PSDEF - in_particles::SVector{N_IN_PARTICLES,ParticleStateful{PhaseSpaceElementType}} - out_particles::SVector{N_OUT_PARTICLES,ParticleStateful{PhaseSpaceElementType}} + 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_p::IN_P, out_p::OUT_P + proc::PROC, + model::MODEL, + ps_def::PSDEF, + in_ps::AbstractVector{ELEMENT}, + out_ps::AbstractVector{ELEMENT}, ) where { PROC<:AbstractProcessDefinition, MODEL<:AbstractModelDefinition, PSDEF<:AbstractPhasespaceDefinition, - PhaseSpaceElementType<:AbstractFourMomentum, - IN_P<:AbstractVector{ParticleStateful{PhaseSpaceElementType}}, - OUT_P<:AbstractVector{ParticleStateful{PhaseSpaceElementType}}, + ELEMENT<:AbstractFourMomentum, } - length(incoming_particles(proc)) == length(in_p) || throw( - InvalidInputError( - "the number of incoming particles given by the process ($(incoming_particles(proc))) mismatches the number of given stateful incoming particles ($(length(in_p)))", - ), + 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) ) - length(outgoing_particles(proc)) == length(out_p) || throw( - InvalidInputError( - "the number of outgoing particles given by the process ($(outgoing_particles(proc))) mismatches the number of given stateful outgoing particles ($(length(out_p)))", - ), + + 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) ) - for (proc_p, p) in zip(incoming_particles(proc), in_p) - proc_p == p.species || throw( - InvalidInputError( - "process given particle species ($(proc_p)) does not match stateful particle species ($(p.species))", - ), - ) - is_incoming(p) || throw( - InvalidInputError( - "stateful particle $(p) is given as an incoming particle but is outgoing", - ), - ) - end - for (proc_p, p) in zip(outgoing_particles(proc), out_p) - proc_p == p.species || throw( - InvalidInputError( - "process given particle species ($(proc_p)) does not match stateful particle species ($(p.species))", - ), - ) - is_outgoing(p) || throw( - InvalidInputError( - "stateful particle $(p) is given as an outgoing particle but is incoming", - ), - ) - end + # no need to check anything since it is constructed correctly above - return new{PROC,MODEL,PSDEF,PhaseSpaceElementType,length(in_p),length(out_p)}( - proc, model, ps_def, in_p, out_p + 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) @@ -349,45 +414,6 @@ function momentum(psp::PhaseSpacePoint, dir::ParticleDirection, n::Int) return psp[dir, n].mom end -""" - generate_phase_space( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - ps_def::AbstractPhasespaceDefinition, - in_ps::AbstractVector{ElType}, - out_ps::AbstractVector{ElType}, - ) where {ElType<:QEDbase.AbstractFourMomentum} - -Return the respective phase space point for given momenta of incoming and outgoing particles regarding a given process. -""" -function generate_phase_space( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - ps_def::AbstractPhasespaceDefinition, - in_ps::AbstractVector{ElType}, - out_ps::AbstractVector{ElType}, -) where {ElType<:QEDbase.AbstractFourMomentum} - in_particles = incoming_particles(proc) - in_n = number_incoming_particles(proc) - in_parts = SVector{in_n,ParticleStateful{SFourMomentum}}( - collect( - ParticleStateful(Incoming(), particle, mom) for - (particle, mom) in zip(in_particles, in_ps) - ), - ) - - out_particles = outgoing_particles(proc) - out_n = number_outgoing_particles(proc) - out_parts = SVector{out_n,ParticleStateful{SFourMomentum}}( - collect( - ParticleStateful(Outgoing(), particle, mom) for - (particle, mom) in zip(out_particles, out_ps) - ), - ) - - return PhaseSpacePoint(proc, model, ps_def, in_parts, out_parts) -end - function Base.show(io::IO, psp::PhaseSpacePoint) print(io, "PhaseSpacePoint of $(psp.proc)") return nothing diff --git a/src/processes/one_photon_compton/perturbative/cross_section.jl b/src/processes/one_photon_compton/perturbative/cross_section.jl index cb1173a..eaafa8c 100644 --- a/src/processes/one_photon_compton/perturbative/cross_section.jl +++ b/src/processes/one_photon_compton/perturbative/cross_section.jl @@ -10,8 +10,8 @@ function _incident_flux( end function _matrix_element(psp::PhaseSpacePoint{<:Compton,PerturbativeQED}) - in_ps = momentum.(psp.in_particles) - out_ps = momentum.(psp.out_particles) + in_ps = momenta(psp, Incoming()) + out_ps = momenta(psp, Outgoing()) return _pert_compton_matrix_element(psp.proc, in_ps, out_ps) end @@ -38,8 +38,8 @@ function _all_onshell( end function _is_in_phasespace(psp::PhaseSpacePoint{<:Compton,PerturbativeQED}) - in_ps = momentum.(psp.in_particles) - out_ps = momentum.(psp.out_particles) + in_ps = momenta(psp, Incoming()) + out_ps = momenta(psp, Outgoing()) if (!isapprox(sum(in_ps), sum(out_ps))) return false end @@ -47,8 +47,8 @@ function _is_in_phasespace(psp::PhaseSpacePoint{<:Compton,PerturbativeQED}) end @inline function _phase_space_factor(psp::PhaseSpacePoint{<:Compton,PerturbativeQED}) - in_ps = momentum.(psp.in_particles) - out_ps = momentum.(psp.out_particles) + in_ps = momenta(psp, Incoming()) + out_ps = momenta(psp, Outgoing()) return _pert_compton_ps_fac(psp.ps_def, in_ps[2], out_ps[2]) end diff --git a/test/cross_sections.jl b/test/cross_sections.jl index b051ca5..505514e 100644 --- a/test/cross_sections.jl +++ b/test/cross_sections.jl @@ -13,8 +13,8 @@ TESTPSDEF = TestImplementation.TestPhasespaceDef() @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, TestImplementation.PARTICLE_SET, N_INCOMING) - OUTGOING_PARTICLES = rand(RNG, TestImplementation.PARTICLE_SET, N_OUTGOING) + 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) @@ -35,7 +35,7 @@ TESTPSDEF = TestImplementation.TestPhasespaceDef() @testset "differential cross section" begin @testset "unsafe compute" begin - PS_POINT = generate_phase_space( + PS_POINT = PhaseSpacePoint( TESTPROC, TESTMODEL, TESTPSDEF, p_in_phys, p_out_phys ) @@ -49,7 +49,7 @@ TESTPSDEF = TestImplementation.TestPhasespaceDef() @testset "safe compute" begin for (P_IN, P_OUT) in p_combs - PS_POINT = generate_phase_space(TESTPROC, TESTMODEL, TESTPSDEF, P_IN, P_OUT) + PS_POINT = PhaseSpacePoint(TESTPROC, TESTMODEL, TESTPSDEF, P_IN, P_OUT) diffCS_on_psp = differential_cross_section(PS_POINT) groundtruth = TestImplementation._groundtruth_safe_diffCS( @@ -86,7 +86,7 @@ TESTPSDEF = TestImplementation.TestPhasespaceDef() @testset "differential probability" begin @testset "unsafe compute" begin - PS_POINT = generate_phase_space( + PS_POINT = PhaseSpacePoint( TESTPROC, TESTMODEL, TESTPSDEF, p_in_phys, p_out_phys ) prop_on_psp = unsafe_differential_probability(PS_POINT) @@ -98,7 +98,7 @@ TESTPSDEF = TestImplementation.TestPhasespaceDef() @testset "safe compute" begin for (P_IN, P_OUT) in p_combs - PS_POINT = generate_phase_space(TESTPROC, TESTMODEL, TESTPSDEF, P_IN, P_OUT) + PS_POINT = PhaseSpacePoint(TESTPROC, TESTMODEL, TESTPSDEF, P_IN, P_OUT) prop_on_psp = differential_probability(PS_POINT) groundtruth = TestImplementation._groundtruth_safe_probability( TESTPROC, P_IN, P_OUT diff --git a/test/interfaces/process_interface.jl b/test/interfaces/process_interface.jl index ccb86f4..6da09b1 100644 --- a/test/interfaces/process_interface.jl +++ b/test/interfaces/process_interface.jl @@ -11,15 +11,15 @@ 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 = rand(RNG, TestImplementation.PARTICLE_SET, N_INCOMING) - OUTGOING_PARTICLES = rand(RNG, TestImplementation.PARTICLE_SET, N_OUTGOING) + 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 = generate_phase_space(TESTPROC, TESTMODEL, TESTPSDEF, IN_PS, OUT_PS) + PSP = PhaseSpacePoint(TESTPROC, TESTMODEL, TESTPSDEF, IN_PS, OUT_PS) @testset "failed interface" begin TESTPROC_FAIL_ALL = TestImplementation.TestProcess_FAIL_ALL( @@ -40,7 +40,7 @@ include("../test_implementation/TestImplementation.jl") (TESTPROC, TESTPROC_FAIL_DIFFCS), (TESTMODEL, TESTMODEL_FAIL) ) if TestImplementation._any_fail(PROC, MODEL) - psp = generate_phase_space(PROC, MODEL, TESTPSDEF, IN_PS, OUT_PS) + psp = PhaseSpacePoint(PROC, MODEL, TESTPSDEF, IN_PS, OUT_PS) @test_throws MethodError QEDprocesses._incident_flux(psp) @test_throws MethodError QEDprocesses._averaging_norm(psp) @test_throws MethodError QEDprocesses._matrix_element(psp) @@ -48,7 +48,7 @@ include("../test_implementation/TestImplementation.jl") for PS_DEF in (TESTPSDEF, TESTPSDEF_FAIL) if TestImplementation._any_fail(PROC, MODEL, PS_DEF) - psp = generate_phase_space(PROC, MODEL, PS_DEF, IN_PS, OUT_PS) + psp = PhaseSpacePoint(PROC, MODEL, PS_DEF, IN_PS, OUT_PS) @test_throws MethodError QEDprocesses._phase_space_factor(psp) end end @@ -98,13 +98,13 @@ include("../test_implementation/TestImplementation.jl") IN_PS_unphysical[1] = SFourMomentum(zeros(4)) OUT_PS_unphysical = deepcopy(OUT_PS) OUT_PS_unphysical[end] = ones(SFourMomentum) - PSP_unphysical_in_ps = generate_phase_space( + PSP_unphysical_in_ps = PhaseSpacePoint( TESTPROC, TESTMODEL, TESTPSDEF, IN_PS_unphysical, OUT_PS ) - PSP_unphysical_out_ps = generate_phase_space( + PSP_unphysical_out_ps = PhaseSpacePoint( TESTPROC, TESTMODEL, TESTPSDEF, IN_PS, OUT_PS_unphysical ) - PSP_unphysical = generate_phase_space( + PSP_unphysical = PhaseSpacePoint( TESTPROC, TESTMODEL, TESTPSDEF, IN_PS_unphysical, OUT_PS_unphysical ) diff --git a/test/phase_spaces.jl b/test/phase_spaces.jl index 817d6d3..40407fd 100644 --- a/test/phase_spaces.jl +++ b/test/phase_spaces.jl @@ -21,56 +21,37 @@ end @testset "Stateful Particle" begin DIRECTIONS = [Incoming(), Outgoing()] SPECIES = [Electron(), Positron()] #=, Muon(), AntiMuon(), Tauon(), AntiTauon()=# - SPINANDPOLS = [AllSpin(), SpinUp(), SpinDown(), AllPol(), PolX(), PolY()] - for (species, dir, spin_or_pol) in Iterators.product(SPECIES, DIRECTIONS, SPINANDPOLS) + for (species, dir) in Iterators.product(SPECIES, DIRECTIONS) mom = rand(RNG, SFourMomentum) - if (is_fermion(species) && (spin_or_pol isa AbstractSpin)) || - (is_boson(species) && (spin_or_pol isa AbstractPolarization)) - particle_stateful = ParticleStateful(dir, species, mom, spin_or_pol) - - # particle interface - @test is_fermion(particle_stateful) == is_fermion(species) - @test is_boson(particle_stateful) == is_boson(species) - @test is_particle(particle_stateful) == is_particle(species) - @test is_anti_particle(particle_stateful) == is_anti_particle(species) - @test is_incoming(particle_stateful) == is_incoming(dir) - @test is_outgoing(particle_stateful) == is_outgoing(dir) - @test mass(particle_stateful) == mass(species) - @test charge(particle_stateful) == charge(species) - - # accessors - @test particle_stateful.dir == dir - @test particle_direction(particle_stateful) == particle_stateful.dir - @test particle_stateful.species == species - @test particle_species(particle_stateful) == particle_stateful.species - @test particle_stateful.mom == mom - @test momentum(particle_stateful) == mom - if (is_fermion(species)) - @test spin(particle_stateful) == spin_or_pol - @test_throws MethodError polarization(particle_stateful) - else - @test polarization(particle_stateful) == spin_or_pol - @test_throws MethodError spin(particle_stateful) - end - - # printing - print(BUF, particle_stateful) - @test String(take!(BUF)) == "$(dir) $(species) ($(spin_or_pol)): $(mom)" - - show(BUF, MIME"text/plain"(), particle_stateful) - @test String(take!(BUF)) == - "ParticleStateful: $(dir) $(species)\n $(spin_or_pol isa AbstractSpin ? "spin" : "polarization"): $(spin_or_pol)\n momentum: $(mom)\n" - else - if (VERSION >= v"1.8") - # julia versions before 1.8 did not have support for regex matching in @test_throws - @test_throws "MethodError: no method matching ParticleStateful" ParticleStateful( - dir, species, mom, spin_or_pol - ) - end - @test_throws MethodError ParticleStateful(dir, species, mom, spin_or_pol) - end + particle_stateful = ParticleStateful(dir, species, mom) + + # particle interface + @test is_fermion(particle_stateful) == is_fermion(species) + @test is_boson(particle_stateful) == is_boson(species) + @test is_particle(particle_stateful) == is_particle(species) + @test is_anti_particle(particle_stateful) == is_anti_particle(species) + @test is_incoming(particle_stateful) == is_incoming(dir) + @test is_outgoing(particle_stateful) == is_outgoing(dir) + @test mass(particle_stateful) == mass(species) + @test charge(particle_stateful) == charge(species) + + # accessors + @test particle_stateful.dir == dir + @test particle_direction(particle_stateful) == particle_stateful.dir + @test particle_stateful.species == species + @test particle_species(particle_stateful) == particle_stateful.species + @test particle_stateful.mom == mom + @test 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 @@ -85,29 +66,27 @@ end out_el = ParticleStateful(Outgoing(), Electron(), out_el_mom) out_ph = ParticleStateful(Outgoing(), Photon(), out_ph_mom) - in_particles_valid = SVector(in_el, in_ph) - in_particles_invalid = SVector(in_el, out_ph) + in_particles_valid = (in_el, in_ph) + in_particles_invalid = (in_el, out_ph) - out_particles_valid = SVector(out_el, out_ph) - out_particles_invalid = SVector(out_el, in_ph) + out_particles_valid = (out_el, out_ph) + out_particles_invalid = (out_el, in_ph) model = TESTMODEL - process = TestImplementation.TestProcess( - SVector{2,AbstractParticle}(Electron(), Photon()), - SVector{2,AbstractParticle}(Electron(), Photon()), - ) + 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 \(all spins\): (.*)\n -> incoming photon \(all polarizations\): (.*)\n outgoing particles:\n -> outgoing electron \(all spins\): (.*)\n -> outgoing photon \(all polarizations\): (.*)\n", + 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 @@ -126,20 +105,20 @@ 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"stateful particle (.*) is given as an incoming particle but is outgoing" PhaseSpacePoint( + @test_throws r"expected incoming photon but got outgoing photon" PhaseSpacePoint( process, model, phasespace_def, in_particles_invalid, out_particles_valid ) - @test_throws r"stateful particle (.*) is given as an outgoing particle but is incoming" PhaseSpacePoint( + @test_throws r"expected outgoing photon but got incoming photon" PhaseSpacePoint( process, model, phasespace_def, in_particles_valid, out_particles_invalid ) - @test_throws r"process given particle species \(electron\) does not match stateful particle species \(photon\)" PhaseSpacePoint( - process, model, phasespace_def, SVector(in_ph, in_el), out_particles_valid + @test_throws r"expected incoming electron but got incoming photon" PhaseSpacePoint( + process, model, phasespace_def, (in_ph, in_el), out_particles_valid ) - @test_throws r"process given particle species \(electron\) does not match stateful particle species \(photon\)" PhaseSpacePoint( - process, model, phasespace_def, in_particles_valid, SVector(out_ph, out_el) + @test_throws r"expected outgoing electron but got outgoing photon" PhaseSpacePoint( + process, model, phasespace_def, in_particles_valid, (out_ph, out_el) ) end @@ -162,17 +141,21 @@ end ) @test_throws InvalidInputError PhaseSpacePoint( - process, model, phasespace_def, SVector(in_ph, in_el), out_particles_valid + process, model, phasespace_def, (in_ph, in_el), out_particles_valid ) @test_throws InvalidInputError PhaseSpacePoint( - process, model, phasespace_def, in_particles_valid, SVector(out_ph, out_el) + process, model, phasespace_def, in_particles_valid, (out_ph, out_el) ) end @testset "Generation" begin - test_psp = generate_phase_space( - process, model, phasespace_def, [in_el_mom, in_ph_mom], [out_el_mom, out_ph_mom] + test_psp = PhaseSpacePoint( + process, + model, + phasespace_def, + SVector(in_el_mom, in_ph_mom), + SVector(out_el_mom, out_ph_mom), ) @test test_psp.proc == process diff --git a/test/processes/one_photon_compton/perturbative.jl b/test/processes/one_photon_compton/perturbative.jl index 1b61ac3..27f06ce 100644 --- a/test/processes/one_photon_compton/perturbative.jl +++ b/test/processes/one_photon_compton/perturbative.jl @@ -65,7 +65,7 @@ end omega, cos_theta, 1.0 ) - PSP = generate_phase_space(PROC, MODEL, PS_DEF, IN_PS, OUT_PS) + PSP = PhaseSpacePoint(PROC, MODEL, PS_DEF, IN_PS, OUT_PS) test_val = unsafe_differential_cross_section(PSP) @test isapprox(test_val, groundtruth, atol=ATOL, rtol=RTOL) @@ -87,7 +87,7 @@ end omega, cos_theta, phi, 1.0 ) - PSP = generate_phase_space(PROC, MODEL, PS_DEF, IN_PS, OUT_PS) + PSP = PhaseSpacePoint(PROC, MODEL, PS_DEF, IN_PS, OUT_PS) test_val = unsafe_differential_cross_section(PSP) @test isapprox(test_val, groundtruth, atol=ATOL, rtol=RTOL) @@ -109,7 +109,7 @@ end omega, cos_theta, phi, 1.0 ) - PSP = generate_phase_space(PROC, MODEL, PS_DEF, IN_PS, OUT_PS) + PSP = PhaseSpacePoint(PROC, MODEL, PS_DEF, IN_PS, OUT_PS) test_val = unsafe_differential_cross_section(PSP) @test isapprox(test_val, groundtruth, atol=ATOL, rtol=RTOL) diff --git a/test/test_implementation/test_process.jl b/test/test_implementation/test_process.jl index 8dfc209..003d14e 100644 --- a/test/test_implementation/test_process.jl +++ b/test/test_implementation/test_process.jl @@ -1,4 +1,3 @@ - # dummy particles struct TestParticleFermion <: FermionLike end struct TestParticleBoson <: BosonLike end @@ -10,7 +9,7 @@ const PARTICLE_SET = [TestParticleFermion(), TestParticleBoson()] TestProcess(rng,incoming_particles,outgoing_particles) """ -struct TestProcess{IP<:AbstractVector,OP<:AbstractVector} <: AbstractProcessDefinition +struct TestProcess{IP<:Tuple,OP<:Tuple} <: AbstractProcessDefinition incoming_particles::IP outgoing_particles::OP end @@ -24,6 +23,17 @@ end QEDprocesses.incoming_particles(proc::TestProcess) = proc.incoming_particles QEDprocesses.outgoing_particles(proc::TestProcess) = proc.outgoing_particles +struct TestProcess_FAIL{IP<:Tuple,OP<:Tuple} <: 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 QEDprocesses.in_phase_space_dimension(proc::TestProcess, ::TestModel) return number_incoming_particles(proc) * 4 end @@ -34,15 +44,14 @@ end """ Test process with no implemented interface. Should fail every usage except construction. """ -struct TestProcess_FAIL_ALL{IP<:AbstractVector,OP<:AbstractVector} <: - AbstractProcessDefinition +struct TestProcess_FAIL_ALL{IP<:Tuple,OP<:Tuple} <: AbstractProcessDefinition incoming_particles::IP outgoing_particles::OP end function TestProcess_FAIL_ALL(rng::AbstractRNG, N_in::Int, N_out::Int) - in_particles = rand(rng, PARTICLE_SET, N_in) - out_particles = rand(rng, PARTICLE_SET, N_out) + 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 @@ -50,15 +59,14 @@ 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<:AbstractVector,OP<:AbstractVector} <: - AbstractProcessDefinition +struct TestProcess_FAIL_DIFFCS{IP<:Tuple,OP<:Tuple} <: AbstractProcessDefinition incoming_particles::IP outgoing_particles::OP end function TestProcess_FAIL_DIFFCS(rng::AbstractRNG, N_in::Int, N_out::Int) - in_particles = rand(rng, PARTICLE_SET, N_in) - out_particles = rand(rng, PARTICLE_SET, N_out) + 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 @@ -82,22 +90,22 @@ function QEDprocesses._averaging_norm(proc::TestProcess) end function QEDprocesses._matrix_element(psp::PhaseSpacePoint{<:TestProcess,TestModel}) - in_ps = momentum.(psp.in_particles) - out_ps = momentum.(psp.out_particles) + in_ps = momenta(psp, Incoming()) + out_ps = momenta(psp, Outgoing()) return _groundtruth_matrix_element(in_ps, out_ps) end function QEDprocesses._is_in_phasespace(psp::PhaseSpacePoint{<:TestProcess,TestModel}) - in_ps = momentum.(psp.in_particles) - out_ps = momentum.(psp.out_particles) + in_ps = momenta(psp, Incoming()) + out_ps = momenta(psp, Outgoing()) return _groundtruth_is_in_phasespace(in_ps, out_ps) end function QEDprocesses._phase_space_factor( psp::PhaseSpacePoint{<:TestProcess,TestModel,TestPhasespaceDef} ) - in_ps = momentum.(psp.in_particles) - out_ps = momentum.(psp.out_particles) + in_ps = momenta(psp, Incoming()) + out_ps = momenta(psp, Outgoing()) return _groundtruth_phase_space_factor(in_ps, out_ps) end