From 008c6ee9bd673cd21caabeeaaddb7e3de4643be6 Mon Sep 17 00:00:00 2001 From: Anton Reinhard Date: Thu, 6 Jun 2024 11:32:52 +0200 Subject: [PATCH 01/21] Add interfaces from QEDprocesses.jl --- src/QEDbase.jl | 24 +++ src/interfaces/model_interface.jl | 26 +++ .../momentum_generation_interface.jl | 54 +++++ src/interfaces/process_interface.jl | 202 ++++++++++++++++++ src/interfaces/setup_interface.jl | 188 ++++++++++++++++ src/particles/particle_types.jl | 27 +++ test/utils.jl | 1 - 7 files changed, 521 insertions(+), 1 deletion(-) create mode 100644 src/interfaces/model_interface.jl create mode 100644 src/interfaces/momentum_generation_interface.jl create mode 100644 src/interfaces/process_interface.jl create mode 100644 src/interfaces/setup_interface.jl diff --git a/src/QEDbase.jl b/src/QEDbase.jl index 8439461..ccedd89 100644 --- a/src/QEDbase.jl +++ b/src/QEDbase.jl @@ -56,6 +56,7 @@ export BosonLike, Boson, AntiBoson, MajoranaBoson export Electron, Positron, Photon export ParticleDirection, Incoming, Outgoing export is_incoming, is_outgoing +export propagator # polarizations and spins export AbstractSpinOrPolarization, AbstractPolarization, AbstractSpin @@ -69,6 +70,26 @@ using StaticArrays using LinearAlgebra using DocStringExtensions +# probabilities +export differential_probability, unsafe_differential_probability +export total_probability + +# differential cross section +export differential_cross_section, unsafe_differential_cross_section +export total_cross_section + +# Abstract model interface +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 + +# Abstract setup interface +export AbstractComputationSetup, InvalidInputError, compute +export AbstractProcessSetup, scattering_process, physical_model + include("dirac_tensors.jl") include("lorentz_interface.jl") include("lorentz_vector.jl") @@ -77,6 +98,9 @@ include("gamma_matrices.jl") include("four_momentum.jl") # maybe go to a kinematics module!! include("interfaces/particle_interface.jl") +include("interfaces/model_interface.jl") +include("interfaces/setup_interface.jl") + include("particles/particle_types.jl") include("particles/particle_direction.jl") include("particles/particle_spin_pol.jl") diff --git a/src/interfaces/model_interface.jl b/src/interfaces/model_interface.jl new file mode 100644 index 0000000..7b98325 --- /dev/null +++ b/src/interfaces/model_interface.jl @@ -0,0 +1,26 @@ +############### +# The model interface +# +# In this file, we define the interface of working with compute models in +# general. +############### +# root type for models +""" +Abstract base type for all compute model definitions in the context of scattering processes. Every subtype of `AbstractModelDefinition` is associated with a fundamental interaction. +Therefore, one needs to implement the following soft interface function + +```Julia +fundamental_interaction_type(::AbstractModelDefinition) +``` +""" +abstract type AbstractModelDefinition end + +# broadcast every model as a scalar +Broadcast.broadcastable(model::AbstractModelDefinition) = Ref(model) + +""" + fundamental_interaction_type(models_def::AbstractModelDefinition) + +Return the fundamental interaction associated with the passed model definition. +""" +function fundamental_interaction_type end diff --git a/src/interfaces/momentum_generation_interface.jl b/src/interfaces/momentum_generation_interface.jl new file mode 100644 index 0000000..7e7d278 --- /dev/null +++ b/src/interfaces/momentum_generation_interface.jl @@ -0,0 +1,54 @@ +##### +# Generation of four-momenta from coordinates +# +# This file contains the interface and functionality to compute momenta from +# given coordinates. +##### + +""" + _generate_incoming_momenta + proc::AbstractProcessDefinition, + model::AbstractModelDefinition, + phase_space_def::AbstractPhasespaceDefinition, + 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. +""" +function _generate_incoming_momenta end + +""" + _generate_outgoing_momenta + proc::AbstractProcessDefinition, + model::AbstractModelDefinition, + phase_space_def::AbstractPhasespaceDefinition, + 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. +""" +function _generate_outgoing_momenta end + +""" + _generate_momenta( + proc::AbstractProcessDefinition, + model::AbstractModelDefinition, + phase_space_def::AbstractPhasespaceDefinition, + 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. +""" +function _generate_momenta( + proc::AbstractProcessDefinition, + model::AbstractModelDefinition, + phase_space_def::AbstractPhasespaceDefinition, + 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) + + return in_momenta, out_momenta +end diff --git a/src/interfaces/process_interface.jl b/src/interfaces/process_interface.jl new file mode 100644 index 0000000..44291ae --- /dev/null +++ b/src/interfaces/process_interface.jl @@ -0,0 +1,202 @@ +############### +# The process interface +# +# In this file, we define the interface for working with scattering processes in +# general. +############### + +""" +Abstract base type for definitions of scattering processes. It is the root type for the +process interface, which assumes that every subtype of `AbstractProcessDefinition` +implements at least + +```Julia +incoming_particles(proc_def::AbstractProcessDefinition) +outgoing_particles(proc_def::AbstractProcessDefinition) +``` + +which return a tuple of the incoming and outgoing particles, respectively. + +Furthermore, to calculate scattering probabilities and differential cross sections, the following +interface functions need to be implemented for every combination of `CustomProcess<:AbstractProcessDefinition`, +`CustomModel<:AbstractModelDefinition`, and `CustomPhasespaceDefinition<:AbstractPhasespaceDefinition`. + +```Julia + _incident_flux(psp::InPhaseSpacePoint{CustomProcess,CustomModel}) + + _matrix_element(psp::PhaseSpacePoint{CustomProcess,CustomModel}) + + _averaging_norm(proc::CustomProcess) + + _is_in_phasespace(psp::PhaseSpacePoint{CustomProcess,CustomModel}) + + _phase_space_factor(psp::PhaseSpacePoint{CustomProcess,CustomModel,CustomPhasespaceDefinition}) +``` + +Optional is the implementation of + +```Julia + + _total_probability(psp::PhaseSpacePoint{CustomProcess,CustomModel,CustomPhasespaceDefinition}) + +``` +to enable the calculation of total probabilities and cross sections. + +""" +abstract type AbstractProcessDefinition end + +# broadcast every model as a scalar +Broadcast.broadcastable(proc::AbstractProcessDefinition) = Ref(proc) + +""" + incoming_particles(proc_def::AbstractProcessDefinition) + +Interface function for scattering processes. Return a tuple of the incoming particles for the given process definition. +This function needs to be given to implement the scattering process interface. +""" +function incoming_particles end + +""" + outgoing_particles(proc_def::AbstractProcessDefinition) + +Interface function for scattering processes. Return the tuple of outgoing particles for the given process definition. +This function needs to be given to implement the scattering process interface. +""" +function outgoing_particles end + +""" + _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 [`InPhaseSpacePoint`](@ref). + +""" +function _incident_flux end + +""" + _matrix_element(PhaseSpacePoint{PROC,MODEL}) where { + PROC <: AbstractProcessDefinition, + MODEL <: AbstractModelDefinition, + } + +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) + +Interface function, which returns a normalization for the averaging of the squared matrix elements over spins and polarizations. +""" +function _averaging_norm end + +""" + _is_in_phasespace(PhaseSpacePoint{PROC,MODEL}) where { + PROC <: AbstractProcessDefinition, + MODEL <: AbstractModelDefinition, + } + +Interface function which returns `true` if the combination of the given incoming and outgoing phase space +is physical, i.e. all momenta are on-shell and some sort of energy-momentum conservation holds. +""" +function _is_in_phasespace end + +""" + _phase_space_factor(PhaseSpacePoint{PROC,MODEL,PSDEF}) where { + PROC <: AbstractProcessDefinition, + MODEL <: AbstractModelDefinition + PSDEF <: AbstractPhasespaceDefinition, + } + +Interface function, which returns the pre-differential factor of the invariant phase space intergral measure. + +!!! note "Convention" + + It is assumed, that this function returns the value of + + ```math + \\mathrm{d}\\Pi_n:= \\prod_{i=1}^N \\frac{\\mathrm{d}^3p_i}{(2\\pi)^3 2 p_i^0} H(P_t, p_1, \\dots, p_N), + ``` +where ``H(\\dots)`` is a characteristic function (or distribution) which constrains the phase space, e.g. ``\\delta^{(4)}(P_t - \\sum_i p_i)``. +""" +function _phase_space_factor end + +""" + in_phase_space_dimension( + proc::AbstractProcessDefinition, + model::AbstractModelDefinition, + ) +TBW +""" +function in_phase_space_dimension end + +""" + out_phase_space_dimension( + proc::AbstractProcessDefinition, + model::AbstractModelDefinition, + ) +TBW +""" +function out_phase_space_dimension end + +""" + _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 [`InPhaseSpacePoint`](@ref). + +!!! note "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 + +####################### +# +# utility functions +# +####################### + +""" + number_incoming_particles(proc_def::AbstractProcessDefinition) + +Return the number of incoming particles of a given process. +""" +@inline function number_incoming_particles(proc_def::AbstractProcessDefinition) + return length(incoming_particles(proc_def)) +end + +""" + number_outgoing_particles(proc_def::AbstractProcessDefinition) + +Return the number of outgoing particles of a given process. +""" +@inline function number_outgoing_particles(proc_def::AbstractProcessDefinition) + 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) diff --git a/src/interfaces/setup_interface.jl b/src/interfaces/setup_interface.jl new file mode 100644 index 0000000..e495549 --- /dev/null +++ b/src/interfaces/setup_interface.jl @@ -0,0 +1,188 @@ +############### +# The process setup +# +# In this file, we define the interface for general computation and process setups. +############### + +""" +Abstract base type for computation setups. A *setup* means +a collection of setup data needed to evaluate a dedicated quantity of given +running data. Therefore, each setup is associated with a single quantity, which one may compute using the setup data and the running data. +Despite that, the decomposition into setup and running data is +arbitrary, and this can be used for cases where a subset of the variables a +quantity depends on is kept constant. + +!!! note "Computation setup interface" + + The computation performed using a computation setup is separated into three steps: + + 1. input validation + 2. actual computation + 3. post processing + + where every step has its own interface function (see [`compute`](@ref) for details). + + ## Input validation + + Every subtype of `AbstractComputationSetup` should implement the interface function + + ```Julia + _assert_valid_input(stp::AbstractComputationSetup, input) + ``` + + which should throw and an exception subtyped from [`AbstractInvalidInputException`](@ref) if the `input` is not valid for the computation of the associated quantity (see [`_assert_valid_input`](@ref) for more details). + The default implementation does nothing, i.e. every input is valid by default. Provide a custom implementation if a different behavior is required. + + ## Actual computation + + Every subtype of `AbstractComputationSetup` must at least implement the required interface function + + ```Julia + _compute(stp::AbstractComputationSetup, input) + ``` + + which computes the value of the associated quantity for a given `input` (see [`_compute`](@ref) for more details). + + + ## Post processing + + Every subtype of `AbstractComputationSetup` should implement the interface function + + ```Julia + _post_processing(stp::AbstractComputationSetup, input, result) + ``` + + which performs task after the actual computation, e.g. conversions or normalizations (see [`_post_processing`](@ref) for more details). + +""" +abstract type AbstractComputationSetup end + +# convenience function to check if an object is a computation setup +_is_computation_setup(::AbstractComputationSetup) = true + +""" +Abstract base type for exceptions indicating invalid input. See [`InvalidInputError`](@ref) for a simple concrete implementation. +Concrete implementations should at least implement + +```Julia + +Base.showerror(io::IO, err::CustomInvalidError) where {CustomInvalidError<:AbstractInvalidInputException} + +``` +""" +abstract type AbstractInvalidInputException <: Exception end + +""" + InvalidInputError(msg::String) + +Exception which is thrown if a given input is invalid, e.g. passed to [`_assert_valid_input`](@ref). +""" +struct InvalidInputError <: AbstractInvalidInputException + msg::String +end +function Base.showerror(io::IO, err::InvalidInputError) + return println(io, "InvalidInputError: $(err.msg)") +end + +""" + _assert_valid_input(stp::AbstractComputationSetup, input::Any) + +Interface function, which asserts that the given `input` is valid, and throws an [`InvalidInputError`](@ref) if not. + +!!! note "default implementation" + + By default, every input is assumed to be valid. Therefore, this function does nothing. + To customize this behavior, add your own implementation of + + ```Julia + _assert_valid_input(stp::YourCustomSetup,input) + ``` + which should throw an exception, which is a subtype of [`AbstractInvalidInputException`](@ref). One may also use the concrete implementation [`InvalidInputError`](@ref) if the input is invalid instead of writing a custom exception type. + +""" +@inline function _assert_valid_input(stp::AbstractComputationSetup, input) + return nothing +end + +""" + function _post_processing(stp::AbstractComputationSetup, input::Any, result::Any) + +Interface function, which is called in [`compute`](@ref) after [`_compute`](@ref) has been called. This function is dedicated to +finalize the result of a computation. + +!!! note "default implementation" + + Since in the case of no post processing the result of [`_compute`](@ref) is unchanged, this function returns `result` by default. + +""" +@inline function _post_processing(stp::AbstractComputationSetup, input, result) + return result +end + +""" + _compute(stp::AbstractComputationSetup, input::Any) + +Interface function that returns the value of the associated quantity evaluated on `input`, which can be anything the associated quantity is defined to be feasible for. + +!!! note "unsafe implementation" + + This function must be implemented for any subtype of [`AbstractComputationSetup`](@ref). It should not do any input validation or post processing (see [`_assert_valid_input`](@ref) and [`_post_processing`](@ref)), as those two are performed while calling + the safe version of this function [`compute`](@ref). + +""" +function _compute end + +""" + compute(stp::AbstractComputationSetup, input::Any) + +Return the value of the quantity associated with `stp` for a given `input`. +In addition to the actual call of the associated unsafe version [`_compute`](@ref), +input validation ([`_assert_valid_input`]) and post processing +(using [`_post_processing`](@ref)) are wrapped around the calculation (see [`AbstractComputationSetup`](@ref) for details). +""" +function compute(stp::AbstractComputationSetup, input) + _assert_valid_input(stp, input) + raw_result = _compute(stp, input) + return _post_processing(stp, input, raw_result) +end + +""" +Abstract base type for setups related to combining scattering processes and physical models. +Every subtype of `AbstractProcessSetup` must implement at least the following +interface functions: + +```Julia +scattering_process(::AbstractProcessSetup) +physical_model(::AbstractProcessSetup) +``` + +Derived from these interface functions, the following delegations are provided: + +```Julia +number_incoming_particles(::AbstractProcessSetup) +number_outgoing_particles(::AbstractProcessSetup) +``` + +""" +abstract type AbstractProcessSetup <: AbstractComputationSetup end + +""" + scattering_process(stp::AbstractProcessSetup) + +Interface function that returns the scattering process associated with `stp`, +i.e. an object which is a subtype of [`AbstractProcessDefinition`](@ref). +""" +function scattering_process end + +""" + physical_model(stp::AbstractProcessSetup) + +Interface function that returns the physical model associated with `stp`, i.e. +an object which is a subtype of [`AbstractModelDefinition`](@ref). +""" +function physical_model end + +@inline number_incoming_particles(stp::AbstractProcessSetup) = + number_incoming_particles(scattering_process(stp)) +@inline number_outgoing_particles(stp::AbstractProcessSetup) = + number_outgoing_particles(scattering_process(stp)) diff --git a/src/particles/particle_types.jl b/src/particles/particle_types.jl index 16282e1..01573cd 100644 --- a/src/particles/particle_types.jl +++ b/src/particles/particle_types.jl @@ -207,3 +207,30 @@ struct Photon <: MajoranaBoson end mass(::Photon) = 0.0 charge(::Photon) = 0.0 Base.show(io::IO, ::Photon) = print(io, "photon") + +""" + propagator(particle::AbstractParticleType, mom::QEDbase.AbstractFourMomentum, [mass::Real]) + +Return the propagator of a particle for a given four-momentum. If `mass` is passed, the respective propagator for massive particles is used, if not, it is assumed the particle passed in is massless. + +!!! note "Convention" + + There are two types of implementations for propagators given in `QEDProcesses`: + For a `BosonLike` particle with four-momentum ``k`` and mass ``m``, the propagator is given as + + ```math + D(k) = \\frac{1}{k^2 - m^2}. + ``` + + For a `FermionLike` particle with four-momentum ``p`` and mass ``m``, the propagator is given as + + ```math + S(p) = \\frac{\\gamma^\\mu p_\\mu + mass}{p^2 - m^2}. + ``` + +!!! warning + + This function does not throw when the given particle is off-shell. If an off-shell particle is passed, the function `propagator` returns `Inf`. + +""" +function propagator end diff --git a/test/utils.jl b/test/utils.jl index b4d4f8e..522a9bd 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -1,5 +1,4 @@ """ - Returns the cartesian coordinates (E,px,py,pz) for given spherical coordiantes (E, rho, cos_theta, phi), where rho denotes the length of the respective three-momentum, theta is the polar- and phi the azimuthal angle. From c5fba5121c61122cb34ee039eb8376b5d6b42c6d Mon Sep 17 00:00:00 2001 From: Anton Reinhard Date: Thu, 6 Jun 2024 12:39:53 +0200 Subject: [PATCH 02/21] Export missing interfaces --- src/QEDbase.jl | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/src/QEDbase.jl b/src/QEDbase.jl index ccedd89..8e9ec7a 100644 --- a/src/QEDbase.jl +++ b/src/QEDbase.jl @@ -90,15 +90,10 @@ export particles, number_particles export AbstractComputationSetup, InvalidInputError, compute export AbstractProcessSetup, scattering_process, physical_model -include("dirac_tensors.jl") -include("lorentz_interface.jl") -include("lorentz_vector.jl") -include("gamma_matrices.jl") - -include("four_momentum.jl") # maybe go to a kinematics module!! - -include("interfaces/particle_interface.jl") include("interfaces/model_interface.jl") +include("interfaces/momentum_generation_interface.jl") +include("interfaces/particle_interface.jl") +include("interfaces/process_interface.jl") include("interfaces/setup_interface.jl") include("particles/particle_types.jl") @@ -107,4 +102,11 @@ include("particles/particle_spin_pol.jl") include("particles/particle_spinors.jl") include("particles/particle_states.jl") +include("dirac_tensors.jl") +include("four_momentum.jl") # maybe go to a kinematics module!! +include("four_polarisation.jl") +include("gamma_matrices.jl") +include("lorentz_interface.jl") +include("lorentz_vector.jl") + end #QEDbase From 98103944b688d85e33e884b5a70859ce4dd5e2f9 Mon Sep 17 00:00:00 2001 From: Anton Reinhard Date: Thu, 6 Jun 2024 12:42:36 +0200 Subject: [PATCH 03/21] Move momentum generation interface --- src/QEDbase.jl | 1 - .../momentum_generation_interface.jl | 54 ------------------ src/interfaces/process_interface.jl | 55 +++++++++++++++++++ 3 files changed, 55 insertions(+), 55 deletions(-) delete mode 100644 src/interfaces/momentum_generation_interface.jl diff --git a/src/QEDbase.jl b/src/QEDbase.jl index 8e9ec7a..8e68ade 100644 --- a/src/QEDbase.jl +++ b/src/QEDbase.jl @@ -91,7 +91,6 @@ export AbstractComputationSetup, InvalidInputError, compute export AbstractProcessSetup, scattering_process, physical_model include("interfaces/model_interface.jl") -include("interfaces/momentum_generation_interface.jl") include("interfaces/particle_interface.jl") include("interfaces/process_interface.jl") include("interfaces/setup_interface.jl") diff --git a/src/interfaces/momentum_generation_interface.jl b/src/interfaces/momentum_generation_interface.jl deleted file mode 100644 index 7e7d278..0000000 --- a/src/interfaces/momentum_generation_interface.jl +++ /dev/null @@ -1,54 +0,0 @@ -##### -# Generation of four-momenta from coordinates -# -# This file contains the interface and functionality to compute momenta from -# given coordinates. -##### - -""" - _generate_incoming_momenta - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - phase_space_def::AbstractPhasespaceDefinition, - 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. -""" -function _generate_incoming_momenta end - -""" - _generate_outgoing_momenta - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - phase_space_def::AbstractPhasespaceDefinition, - 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. -""" -function _generate_outgoing_momenta end - -""" - _generate_momenta( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - phase_space_def::AbstractPhasespaceDefinition, - 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. -""" -function _generate_momenta( - proc::AbstractProcessDefinition, - model::AbstractModelDefinition, - phase_space_def::AbstractPhasespaceDefinition, - 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) - - return in_momenta, out_momenta -end diff --git a/src/interfaces/process_interface.jl b/src/interfaces/process_interface.jl index 44291ae..469a1cd 100644 --- a/src/interfaces/process_interface.jl +++ b/src/interfaces/process_interface.jl @@ -200,3 +200,58 @@ Convenience function dispatching to [`number_incoming_particles`](@ref) or [`num number_incoming_particles(proc_def) @inline number_particles(proc_def::AbstractProcessDefinition, ::Outgoing) = number_outgoing_particles(proc_def) + +##### +# Generation of four-momenta from coordinates +# +# This file contains the interface and functionality to compute momenta from +# given coordinates. +##### + +""" + _generate_incoming_momenta + proc::AbstractProcessDefinition, + model::AbstractModelDefinition, + phase_space_def::AbstractPhasespaceDefinition, + 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. +""" +function _generate_incoming_momenta end + +""" + _generate_outgoing_momenta + proc::AbstractProcessDefinition, + model::AbstractModelDefinition, + phase_space_def::AbstractPhasespaceDefinition, + 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. +""" +function _generate_outgoing_momenta end + +""" + _generate_momenta( + proc::AbstractProcessDefinition, + model::AbstractModelDefinition, + phase_space_def::AbstractPhasespaceDefinition, + 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. +""" +function _generate_momenta( + proc::AbstractProcessDefinition, + model::AbstractModelDefinition, + phase_space_def::AbstractPhasespaceDefinition, + 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) + + return in_momenta, out_momenta +end From 1cac751e12601e17e5440ceebfaedb07db5b869a Mon Sep 17 00:00:00 2001 From: Anton Reinhard Date: Thu, 6 Jun 2024 13:05:34 +0200 Subject: [PATCH 04/21] Rename and move includes around --- src/QEDbase.jl | 31 ++++++++++--------- src/four_polarisation.jl | 20 ------------ .../lorentz.jl} | 0 .../{model_interface.jl => model.jl} | 0 .../{particle_interface.jl => particle.jl} | 0 src/interfaces/phase_space.jl | 22 +++++++++++++ .../{process_interface.jl => process.jl} | 0 .../{setup_interface.jl => setup.jl} | 0 .../{particle_direction.jl => direction.jl} | 0 .../{particle_spin_pol.jl => spin_pol.jl} | 0 .../{particle_spinors.jl => spinors.jl} | 0 .../{particle_states.jl => states.jl} | 0 src/particles/{particle_types.jl => types.jl} | 0 13 files changed, 39 insertions(+), 34 deletions(-) delete mode 100644 src/four_polarisation.jl rename src/{lorentz_interface.jl => interfaces/lorentz.jl} (100%) rename src/interfaces/{model_interface.jl => model.jl} (100%) rename src/interfaces/{particle_interface.jl => particle.jl} (100%) create mode 100644 src/interfaces/phase_space.jl rename src/interfaces/{process_interface.jl => process.jl} (100%) rename src/interfaces/{setup_interface.jl => setup.jl} (100%) rename src/particles/{particle_direction.jl => direction.jl} (100%) rename src/particles/{particle_spin_pol.jl => spin_pol.jl} (100%) rename src/particles/{particle_spinors.jl => spinors.jl} (100%) rename src/particles/{particle_states.jl => states.jl} (100%) rename src/particles/{particle_types.jl => types.jl} (100%) diff --git a/src/QEDbase.jl b/src/QEDbase.jl index 8e68ade..7738939 100644 --- a/src/QEDbase.jl +++ b/src/QEDbase.jl @@ -90,22 +90,25 @@ export particles, number_particles export AbstractComputationSetup, InvalidInputError, compute export AbstractProcessSetup, scattering_process, physical_model -include("interfaces/model_interface.jl") -include("interfaces/particle_interface.jl") -include("interfaces/process_interface.jl") -include("interfaces/setup_interface.jl") - -include("particles/particle_types.jl") -include("particles/particle_direction.jl") -include("particles/particle_spin_pol.jl") -include("particles/particle_spinors.jl") -include("particles/particle_states.jl") +include("interfaces/phase_space.jl") +include("interfaces/lorentz.jl") include("dirac_tensors.jl") -include("four_momentum.jl") # maybe go to a kinematics module!! -include("four_polarisation.jl") -include("gamma_matrices.jl") -include("lorentz_interface.jl") include("lorentz_vector.jl") +include("gamma_matrices.jl") +include("four_momentum.jl") # maybe go to a kinematics module!! + +include("interfaces/particle.jl") +include("particles/types.jl") +include("particles/direction.jl") +include("particles/spin_pol.jl") +include("particles/spinors.jl") +include("particles/states.jl") + +include("interfaces/model.jl") + +include("interfaces/process.jl") + +include("interfaces/setup.jl") end #QEDbase diff --git a/src/four_polarisation.jl b/src/four_polarisation.jl deleted file mode 100644 index 7b6d265..0000000 --- a/src/four_polarisation.jl +++ /dev/null @@ -1,20 +0,0 @@ -""" -FourPolarisation type - -TODO: -- type promotion into complex fields - -""" - -struct FourPolarisation <: AbstractLorentzVector{ComplexF64} - t::ComplexF64 - x::ComplexF64 - y::ComplexF64 - z::ComplexF64 -end - -FourPolarisation(t, x, y, z) = FourPolarisation(promote(t, x, y, z)...) - -function FourPolarisation(t::T, x::T, y::T, z::T) where {T<:Number} - return FourPolarisation(complex(t), x, y, z) -end diff --git a/src/lorentz_interface.jl b/src/interfaces/lorentz.jl similarity index 100% rename from src/lorentz_interface.jl rename to src/interfaces/lorentz.jl diff --git a/src/interfaces/model_interface.jl b/src/interfaces/model.jl similarity index 100% rename from src/interfaces/model_interface.jl rename to src/interfaces/model.jl diff --git a/src/interfaces/particle_interface.jl b/src/interfaces/particle.jl similarity index 100% rename from src/interfaces/particle_interface.jl rename to src/interfaces/particle.jl diff --git a/src/interfaces/phase_space.jl b/src/interfaces/phase_space.jl new file mode 100644 index 0000000..7f6d76e --- /dev/null +++ b/src/interfaces/phase_space.jl @@ -0,0 +1,22 @@ +""" + AbstractCoordinateSystem + +TBW +""" +abstract type AbstractCoordinateSystem end + +""" + AbstractFrameOfReference + +TBW +""" +abstract type AbstractFrameOfReference end + +""" + AbstractPhasespaceDefinition + +TBW +""" +abstract type AbstractPhasespaceDefinition end + +Broadcast.broadcastable(ps_def::AbstractPhasespaceDefinition) = Ref(ps_def) diff --git a/src/interfaces/process_interface.jl b/src/interfaces/process.jl similarity index 100% rename from src/interfaces/process_interface.jl rename to src/interfaces/process.jl diff --git a/src/interfaces/setup_interface.jl b/src/interfaces/setup.jl similarity index 100% rename from src/interfaces/setup_interface.jl rename to src/interfaces/setup.jl diff --git a/src/particles/particle_direction.jl b/src/particles/direction.jl similarity index 100% rename from src/particles/particle_direction.jl rename to src/particles/direction.jl diff --git a/src/particles/particle_spin_pol.jl b/src/particles/spin_pol.jl similarity index 100% rename from src/particles/particle_spin_pol.jl rename to src/particles/spin_pol.jl diff --git a/src/particles/particle_spinors.jl b/src/particles/spinors.jl similarity index 100% rename from src/particles/particle_spinors.jl rename to src/particles/spinors.jl diff --git a/src/particles/particle_states.jl b/src/particles/states.jl similarity index 100% rename from src/particles/particle_states.jl rename to src/particles/states.jl diff --git a/src/particles/particle_types.jl b/src/particles/types.jl similarity index 100% rename from src/particles/particle_types.jl rename to src/particles/types.jl From 398b6f187cbdb8db79ef11c71c2ecc80079580db Mon Sep 17 00:00:00 2001 From: Anton Reinhard Date: Thu, 6 Jun 2024 14:32:58 +0200 Subject: [PATCH 05/21] WIP --- Project.toml | 3 +- src/QEDbase.jl | 53 ++-- src/dirac_tensors.jl | 160 ----------- src/four_momentum.jl | 205 -------------- src/gamma_matrices.jl | 72 ----- src/interfaces/dirac_tensors.jl | 18 ++ src/interfaces/four_momentum.jl | 9 + src/interfaces/gamma_matrices.jl | 9 + src/interfaces/lorentz.jl | 57 +--- .../particle_functions.jl} | 139 +-------- .../types.jl => interfaces/particle_types.jl} | 7 + src/lorentz_vector.jl | 154 ---------- src/particles/spinors.jl | 134 --------- src/utils.jl | 12 + test/four_momentum.jl | 212 -------------- test/gamma_matrices.jl | 112 -------- test/interfaces/model.jl | 19 ++ test/interfaces/process.jl | 129 +++++++++ test/interfaces/setup.jl | 188 +++++++++++++ test/runtests.jl | 26 +- .../test_implementation/TestImplementation.jl | 38 +++ test/test_implementation/groundtruths.jl | 265 ++++++++++++++++++ test/test_implementation/random_momenta.jl | 83 ++++++ test/test_implementation/test_model.jl | 5 + test/test_implementation/test_process.jl | 131 +++++++++ test/test_implementation/utils.jl | 44 +++ 26 files changed, 1012 insertions(+), 1272 deletions(-) delete mode 100644 src/dirac_tensors.jl delete mode 100644 src/four_momentum.jl delete mode 100644 src/gamma_matrices.jl create mode 100644 src/interfaces/dirac_tensors.jl create mode 100644 src/interfaces/four_momentum.jl create mode 100644 src/interfaces/gamma_matrices.jl rename src/{particles/states.jl => interfaces/particle_functions.jl} (58%) rename src/{particles/types.jl => interfaces/particle_types.jl} (98%) delete mode 100644 src/lorentz_vector.jl delete mode 100644 src/particles/spinors.jl create mode 100644 src/utils.jl delete mode 100644 test/four_momentum.jl delete mode 100644 test/gamma_matrices.jl create mode 100644 test/interfaces/model.jl create mode 100644 test/interfaces/process.jl create mode 100644 test/interfaces/setup.jl create mode 100644 test/test_implementation/TestImplementation.jl create mode 100644 test/test_implementation/groundtruths.jl create mode 100644 test/test_implementation/random_momenta.jl create mode 100644 test/test_implementation/test_model.jl create mode 100644 test/test_implementation/test_process.jl create mode 100644 test/test_implementation/utils.jl diff --git a/Project.toml b/Project.toml index 75dc58d..948200c 100644 --- a/Project.toml +++ b/Project.toml @@ -26,6 +26,7 @@ julia = "1.6" [extras] SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb" [targets] -test = ["SafeTestsets","Test"] +test = ["SafeTestsets", "Test", "Suppressor"] diff --git a/src/QEDbase.jl b/src/QEDbase.jl index 7738939..c12efc0 100644 --- a/src/QEDbase.jl +++ b/src/QEDbase.jl @@ -1,9 +1,5 @@ module QEDbase -using SimpleTraits -using ArgCheck -using ConstructionBase - import Base: * import StaticArrays: similar_type @@ -28,20 +24,14 @@ export setTransverseMomentum!, setPerp!, setPt! export setTransverseMass!, setMt! export setRapidity! -export AbstractLorentzVector, SLorentzVector, MLorentzVector, dot -export AbstractFourMomentum, SFourMomentum, MFourMomentum +export AbstractLorentzVector, dot +export AbstractFourMomentum export isonshell, assert_onshell -export BiSpinor, AdjointBiSpinor, DiracMatrix, mul export AbstractDiracVector, AbstractDiracMatrix +export mul -export gamma, GAMMA, AbstractGammaRepresentation, DiracGammaRepresentation, slashed - -export BASE_PARTICLE_SPINOR, BASE_ANTIPARTICLE_SPINOR -export IncomingFermionSpinor, - OutgoingFermionSpinor, IncomingAntiFermionSpinor, OutgoingAntiFermionSpinor -export SpinorU, SpinorUbar, SpinorV, SpinorVbar -export @valid_spinor_input +export AbstractGammaRepresentation # particle interface export AbstractParticle @@ -66,10 +56,6 @@ export AbstractDefiniteSpin, AbstractIndefiniteSpin export SpinUp, SpinDown, AllSpin export multiplicity -using StaticArrays -using LinearAlgebra -using DocStringExtensions - # probabilities export differential_probability, unsafe_differential_probability export total_probability @@ -90,25 +76,34 @@ export particles, number_particles export AbstractComputationSetup, InvalidInputError, compute export AbstractProcessSetup, scattering_process, physical_model -include("interfaces/phase_space.jl") +# Abstract phase space interface +export AbstractCoordinateSystem, AbstractFrameOfReference, AbstractPhasespaceDefinition +using StaticArrays +using LinearAlgebra +using DocStringExtensions + +using SimpleTraits +using ArgCheck +using ConstructionBase + +include("utils.jl") + +include("interfaces/dirac_tensors.jl") +include("interfaces/gamma_matrices.jl") include("interfaces/lorentz.jl") -include("dirac_tensors.jl") -include("lorentz_vector.jl") -include("gamma_matrices.jl") -include("four_momentum.jl") # maybe go to a kinematics module!! +include("interfaces/four_momentum.jl") +include("interfaces/model.jl") include("interfaces/particle.jl") -include("particles/types.jl") +include("interfaces/particle_types.jl") +include("interfaces/particle_functions.jl") + include("particles/direction.jl") include("particles/spin_pol.jl") -include("particles/spinors.jl") -include("particles/states.jl") - -include("interfaces/model.jl") +include("interfaces/phase_space.jl") include("interfaces/process.jl") - include("interfaces/setup.jl") end #QEDbase diff --git a/src/dirac_tensors.jl b/src/dirac_tensors.jl deleted file mode 100644 index 774631f..0000000 --- a/src/dirac_tensors.jl +++ /dev/null @@ -1,160 +0,0 @@ -""" -Dirac Tensors -""" -####### -# -# Abstract Dirac Tensor types -# -####### -""" -$(TYPEDEF) - -Abstract type for Dirac vectors, e.g. four dimensional vectors from a spinor space. -""" -abstract type AbstractDiracVector{T} <: FieldVector{4,T} end - -""" -$(TYPEDEF) - -Abstract type for Dirac matrices, i.e. matrix representations for linear mappings from a spinor space into another. -""" -abstract type AbstractDiracMatrix{T} <: FieldMatrix{4,4,T} end - -####### -# -# Concrete Dirac Tensor types -# -####### -""" -$(TYPEDEF) - -Concrete type to model a Dirac four-spinor with complex-valued components. These are the elements of an actual spinor space. -""" -struct BiSpinor <: AbstractDiracVector{ComplexF64} - el1::ComplexF64 - el2::ComplexF64 - el3::ComplexF64 - el4::ComplexF64 -end - -""" -$(TYPEDEF) - -Concrete type to model an adjoint Dirac four-spinor with complex-valued components. These are the elements of the dual spinor space. -""" -struct AdjointBiSpinor <: AbstractDiracVector{ComplexF64} - el1::ComplexF64 - el2::ComplexF64 - el3::ComplexF64 - el4::ComplexF64 -end - -#interface -AdjointBiSpinor(spn::BiSpinor) = AdjointBiSpinor(conj(SVector(spn))) -BiSpinor(spn::AdjointBiSpinor) = BiSpinor(conj(SVector(spn))) - -""" -$(TYPEDEF) - -Concrete type to model Dirac matrices, i.e. matrix representations of linear mappings between two spinor spaces. -""" -struct DiracMatrix <: AbstractDiracMatrix{ComplexF64} - el11::ComplexF64 - el12::ComplexF64 - el13::ComplexF64 - el14::ComplexF64 - el21::ComplexF64 - el22::ComplexF64 - el23::ComplexF64 - el24::ComplexF64 - el31::ComplexF64 - el32::ComplexF64 - el33::ComplexF64 - el34::ComplexF64 - el41::ComplexF64 - el42::ComplexF64 - el43::ComplexF64 - el44::ComplexF64 -end - -####### -# -# Concrete implementation of multiplication for Dirac Tensors -# -####### -""" -$(TYPEDSIGNATURES) - -Tensor product of an adjoint with a standard bi-spinor resulting in a scalar. - -!!! note "Multiplication operator" - This also overloads the `*` operator for this types. - -""" -function mul(aBS::AdjointBiSpinor, BS::BiSpinor)::ComplexF64 - return aBS.el1 * BS.el1 + aBS.el2 * BS.el2 + aBS.el3 * BS.el3 + aBS.el4 * BS.el4 -end -@inline *(aBS::AdjointBiSpinor, BS::BiSpinor) = mul(aBS::AdjointBiSpinor, BS::BiSpinor) - -""" -$(TYPEDSIGNATURES) - -Tensor product of a standard with an adjoint bi-spinor resulting in a Dirac matrix. - -!!! note "Multiplication operator" - This also overloads the `*` operator for this types. - -""" -function mul(BS::BiSpinor, aBS::AdjointBiSpinor)::DiracMatrix - return DiracMatrix(BS * transpose(aBS)) -end -@inline *(BS::BiSpinor, aBS::AdjointBiSpinor) = mul(BS::BiSpinor, aBS::AdjointBiSpinor) - -""" -$(TYPEDSIGNATURES) - -Tensor product of an Dirac matrix with a standard bi-spinor resulting in another standard bi-spinor. - -!!! note "Multiplication operator" - This also overloads the `*` operator for this types. - -""" -function mul(DM::DiracMatrix, BS::BiSpinor)::BiSpinor - return DM * BS -end - -""" -$(TYPEDSIGNATURES) - -Tensor product of an adjoint bi-spinor with a Dirac matrix resulting in another adjoint bi-spinor. - -!!! note "Multiplication operator" - This also overloads the `*` operator for this types. - -""" -function mul(aBS::AdjointBiSpinor, DM::DiracMatrix)::AdjointBiSpinor - return transpose(aBS) * DM -end -@inline *(aBS::AdjointBiSpinor, DM::DiracMatrix) = mul(aBS, DM) - -""" -$(TYPEDSIGNATURES) - -Tensor product two Dirac matrices resulting in another Dirac matrix. - -!!! note "Multiplication operator" - This also overloads the `*` operator for this types. - -""" -function mul(DM1::DiracMatrix, DM2::DiracMatrix)::DiracMatrix - return DM1 * DM2 -end - -""" -$(TYPEDSIGNATURES) - -Tensor product of Dirac matrix sandwiched between an adjoint and a standard bi-spinor resulting in a scalar. -""" -function mul(aBS::AdjointBiSpinor, DM::DiracMatrix, BS::BiSpinor)::ComplexF64 - return transpose(aBS) * DM * BS -end diff --git a/src/four_momentum.jl b/src/four_momentum.jl deleted file mode 100644 index 60bfb68..0000000 --- a/src/four_momentum.jl +++ /dev/null @@ -1,205 +0,0 @@ -""" - AbstractFourMomentum - -Abstract base type for four-momentas, representing one energy and three spacial components. - -Also see: [`SFourMomentum`](@ref) -""" - -abstract type AbstractFourMomentum <: AbstractLorentzVector{Float64} end - -####### -# -# Concrete SFourMomentum type -# -####### -""" -$(TYPEDEF) - -Builds a static LorentzVectorLike with real components used to statically model the four-momentum of a particle or field. - -# Fields -$(TYPEDFIELDS) -""" -struct SFourMomentum <: AbstractFourMomentum - "energy component" - E::Float64 - - "`x` component" - px::Float64 - - "`y` component" - py::Float64 - - "`z` component" - pz::Float64 -end - -""" -$(SIGNATURES) - -The interface transforms each number-like input to float64: - -$(TYPEDSIGNATURES) -""" -function SFourMomentum(t::T, x::T, y::T, z::T) where {T<:Union{Integer,Rational,Irrational}} - return SFourMomentum(float(t), x, y, z) -end - -function similar_type(::Type{A}, ::Type{T}, ::Size{S}) where {A<:SFourMomentum,T<:Real,S} - return SFourMomentum -end -function similar_type(::Type{A}, ::Type{T}, ::Size{S}) where {A<:SFourMomentum,T,S} - return SLorentzVector{T} -end - -@inline getT(p::SFourMomentum) = p.E -@inline getX(p::SFourMomentum) = p.px -@inline getY(p::SFourMomentum) = p.py -@inline getZ(p::SFourMomentum) = p.pz - -register_LorentzVectorLike(SFourMomentum) - -####### -# -# Concrete MFourMomentum type -# -####### -""" -$(TYPEDEF) - -Builds a mutable LorentzVector with real components used to statically model the four-momentum of a particle or field. - -# Fields -$(TYPEDFIELDS) -""" -mutable struct MFourMomentum <: AbstractFourMomentum - "energy component" - E::Float64 - - "`x` component" - px::Float64 - - "`y` component" - py::Float64 - - "`z` component" - pz::Float64 -end - -""" -$(SIGNATURES) - -The interface transforms each number-like input to float64: - -$(TYPEDSIGNATURES) -""" -function MFourMomentum(t::T, x::T, y::T, z::T) where {T<:Union{Integer,Rational,Irrational}} - return MFourMomentum(float(t), x, y, z) -end - -function similar_type(::Type{A}, ::Type{T}, ::Size{S}) where {A<:MFourMomentum,T<:Real,S} - return MFourMomentum -end -function similar_type(::Type{A}, ::Type{T}, ::Size{S}) where {A<:MFourMomentum,T,S} - return MLorentzVector{T} -end - -@inline getT(p::MFourMomentum) = p.E -@inline getX(p::MFourMomentum) = p.px -@inline getY(p::MFourMomentum) = p.py -@inline getZ(p::MFourMomentum) = p.pz - -function QEDbase.setT!(lv::MFourMomentum, value::Float64) - return lv.E = value -end - -function QEDbase.setX!(lv::MFourMomentum, value::Float64) - return lv.px = value -end - -function QEDbase.setY!(lv::MFourMomentum, value::Float64) - return lv.py = value -end - -function QEDbase.setZ!(lv::MFourMomentum, value::Float64) - return lv.pz = value -end - -register_LorentzVectorLike(MFourMomentum) - -function Base.getproperty(P::TM, sym::Symbol) where {TM<:AbstractFourMomentum} - if sym == :t - return P.E - elseif sym == :x - return P.px - elseif sym == :y - return P.py - elseif sym == :z - return P.pz - else - return getfield(P, sym) - end -end - -####### -# -# Utility functions on FourMomenta -# -####### - -function isonshell(mom::QEDbase.AbstractLorentzVector{T}) where {T<:Real} - mag2 = getMag2(mom) - E = getE(mom) - return isapprox(E^2, mag2; rtol=eps(T)) -end - -""" -$(SIGNATURES) - -On-shell check of a given four-momentum `mom` w.r.t. a given mass `mass`. - -!!! note "Precision" - For `AbstactFourMomentum`, the element type is fixed to `Float64`, limiting the precision of comparisons between elements. - The current implementation has been tested within the boundaries for energy scales `E` with `1e-9 <= E <= 1e5`. - In those bounds, the mass error, which is correctly detected as off-shell, is `1e-4` times the mean value of the components, but has at most the value `0.01`, e.g. at the high energy end. - The energy scales correspond to `0.5 meV` for the lower bound and `50 GeV` for the upper bound. - - -!!! todo "FourMomenta with real entries" - * if `AbstractFourMomentum` is updated to elementtypes `T<:Real`, the `AbstractLorentzVector` should be updated with the `AbstractFourMomentum`. -""" -function isonshell(mom::QEDbase.AbstractLorentzVector{T}, mass::Real) where {T<:Real} - if iszero(mass) - return isonshell(mom) - end - mag2 = getMag2(mom) - E = getE(mom) - return isapprox(E^2, (mass)^2 + mag2; atol=2 * eps(T), rtol=eps(T)) -end - -struct OnshellError{M,T} <: Exception - mom::M - mass::T -end - -function Base.showerror(io::IO, e::OnshellError) - return print( - io, - "OnshellError: The momentum $(e.mom) is not onshell w.r.t. the mass $(e.mass).\n mom*mom = $(e.mom*e.mom)", - ) -end - -""" -$(SIGNATURES) - -Assertion if a FourMomentum `mom` is on-shell w.r.t a given mass `mass`. - -!!! note "See also" - The precision of this functions is explained in [`isonshell`](@ref). - -""" -function assert_onshell(mom::QEDbase.AbstractLorentzVector, mass::Real) - isonshell(mom, mass) || throw(OnshellError(mom, mass)) - return nothing -end diff --git a/src/gamma_matrices.jl b/src/gamma_matrices.jl deleted file mode 100644 index 9e936f4..0000000 --- a/src/gamma_matrices.jl +++ /dev/null @@ -1,72 +0,0 @@ -""" -submodule to model Dirac's gamma matrices. -""" - -abstract type AbstractGammaRepresentation end - -#### -# generic definition of the gamma matrices -#### - -function gamma(::Type{T})::SLorentzVector where {T<:AbstractGammaRepresentation} - return SLorentzVector(_gamma0(T), _gamma1(T), _gamma2(T), _gamma3(T)) -end - -#### -# concrete implementatio of gamma matrices in Diracs representation -# -# Note: lower-index version of the gamma matrices are used -# e.g. see https://en.wikipedia.org/wiki/Gamma_matrices -# Note: caused by the column major construction of matrices in Julia, -# the definition below looks *transposed*. -# -#### - -struct DiracGammaRepresentation <: AbstractGammaRepresentation end - -#! format: off -function _gamma0(::Type{DiracGammaRepresentation})::DiracMatrix - return DiracMatrix(1, 0, 0, 0, - 0, 1, 0, 0, - 0, 0, -1, 0, - 0, 0, 0, -1) -end - -function _gamma1(::Type{DiracGammaRepresentation})::DiracMatrix - return DiracMatrix(0, 0, 0, 1, - 0, 0, 1, 0, - 0, -1, 0, 0, - -1, 0, 0, 0) -end - -function _gamma2(::Type{DiracGammaRepresentation})::DiracMatrix - return DiracMatrix( 0,0,0,1im, - 0,0,-1im,0, - 0,-1im,0,0, - 1im,0,0,0) -end - -function _gamma3(::Type{DiracGammaRepresentation})::DiracMatrix - return DiracMatrix(0, 0, 1, 0, - 0, 0, 0, -1, - -1, 0, 0, 0, - 0, 1, 0, 0) -end -#! format: on - -# default gamma matrix is in Dirac's representation -gamma() = gamma(DiracGammaRepresentation) - -const GAMMA = gamma() - -# feynman slash notation - -function slashed( - ::Type{TG}, LV::TV -) where {TG<:AbstractGammaRepresentation,TV<:AbstractLorentzVector} - return gamma(TG) * LV -end - -function slashed(LV::T) where {T<:AbstractLorentzVector} - return GAMMA * LV -end diff --git a/src/interfaces/dirac_tensors.jl b/src/interfaces/dirac_tensors.jl new file mode 100644 index 0000000..88d84b3 --- /dev/null +++ b/src/interfaces/dirac_tensors.jl @@ -0,0 +1,18 @@ +####### +# +# Abstract Dirac Tensor types +# +####### +""" +$(TYPEDEF) + +Abstract type for Dirac vectors, e.g. four dimensional vectors from a spinor space. +""" +abstract type AbstractDiracVector{T} <: FieldVector{4,T} end + +""" +$(TYPEDEF) + +Abstract type for Dirac matrices, i.e. matrix representations for linear mappings from a spinor space into another. +""" +abstract type AbstractDiracMatrix{T} <: FieldMatrix{4,4,T} end diff --git a/src/interfaces/four_momentum.jl b/src/interfaces/four_momentum.jl new file mode 100644 index 0000000..ed05534 --- /dev/null +++ b/src/interfaces/four_momentum.jl @@ -0,0 +1,9 @@ +""" + AbstractFourMomentum + +Abstract base type for four-momentas, representing one energy and three spacial components. + +Also see: [`SFourMomentum`](@ref) +""" + +abstract type AbstractFourMomentum <: AbstractLorentzVector{Float64} end diff --git a/src/interfaces/gamma_matrices.jl b/src/interfaces/gamma_matrices.jl new file mode 100644 index 0000000..1505710 --- /dev/null +++ b/src/interfaces/gamma_matrices.jl @@ -0,0 +1,9 @@ +#### +# generic definition of the gamma matrices +#### + +abstract type AbstractGammaRepresentation end + +function gamma(::Type{T})::SLorentzVector where {T<:AbstractGammaRepresentation} + return SLorentzVector(_gamma0(T), _gamma1(T), _gamma2(T), _gamma3(T)) +end diff --git a/src/interfaces/lorentz.jl b/src/interfaces/lorentz.jl index 8123e00..7456a93 100644 --- a/src/interfaces/lorentz.jl +++ b/src/interfaces/lorentz.jl @@ -1,3 +1,15 @@ +####### +# +# Abstract types +# +####### +""" +$(TYPEDEF) + +Abstract type to model generic Lorentz vectors, i.e. elements of a minkowski-like space, where the component space is arbitray. +""" +abstract type AbstractLorentzVector{T} <: FieldVector{4,T} end + """ ## Definition of LorentzVector interface. @@ -72,7 +84,6 @@ end @traitdef IsMutableLorentzVectorLike{T} """ - getT(lv) Return the 0-component of a given `LorentzVectorLike`. @@ -85,7 +96,6 @@ Return the 0-component of a given `LorentzVectorLike`. function getT end """ - getX(lv) Return the 1-component of a given `LorentzVectorLike`. @@ -98,7 +108,6 @@ Return the 1-component of a given `LorentzVectorLike`. function getX end """ - getY(lv) Return the 2-component of a given `LorentzVectorLike`. @@ -111,7 +120,6 @@ Return the 2-component of a given `LorentzVectorLike`. function getY end """ - getZ(lv) Return the 3-component of a given `LorentzVectorLike`. @@ -124,7 +132,6 @@ Return the 3-component of a given `LorentzVectorLike`. function getZ end """ - setT!(lv,value) Sets the 0-component of a given `LorentzVectorLike` to the given `value`. @@ -132,7 +139,6 @@ Sets the 0-component of a given `LorentzVectorLike` to the given `value`. function setT! end """ - setX!(lv,value) Sets the 1-component of a given `LorentzVectorLike` to the given `value`. @@ -140,7 +146,6 @@ Sets the 1-component of a given `LorentzVectorLike` to the given `value`. function setX! end """ - setY!(lv,value) Sets the 2-component of a given `LorentzVectorLike` to the given `value`. @@ -148,7 +153,6 @@ Sets the 2-component of a given `LorentzVectorLike` to the given `value`. function setY! end """ - setZ!(lv,value) Sets the 3-component of a given `LorentzVectorLike` to the given `value`. @@ -183,7 +187,6 @@ function register_LorentzVectorLike(T::Type) end """ - minkowski_dot(v1,v2) Return the Minkowski dot product of two `LorentzVectorLike`. @@ -216,7 +219,6 @@ const mdot = minkowski_dot # getter ######## """ - getMagnitude2(lv) Return the square of the magnitude of a given `LorentzVectorLike`, i.e. the sum of the squared spatial components. @@ -241,7 +243,6 @@ Functiom alias for [`getMagnitude2`](@ref). const getMag2 = getMagnitude2 """ - getMagnitude(lv) Return the magnitude of a given `LorentzVectorLike`, i.e. the euklidian norm spatial components. @@ -266,7 +267,6 @@ Functiom alias for [`getMagnitude`](@ref). const getMag = getMagnitude """ - getInvariantMass2(lv) Return the squared invariant mass of a given `LorentzVectorLike`, i.e. the minkowski dot with itself. @@ -285,7 +285,6 @@ end const getMass2 = getInvariantMass2 """ - getInvariantMass(lv) Return the the invariant mass of a given `LorentzVectorLike`, i.e. the square root of the minkowski dot with itself. @@ -319,7 +318,6 @@ const getMass = getInvariantMass ########################## """ - getE(lv) Return the energy component of a given `LorentzVectorLike`, i.e. its 0-component. @@ -335,7 +333,6 @@ Return the energy component of a given `LorentzVectorLike`, i.e. its 0-component const getEnergy = getE """ - getPx(lv) Return the ``p_x`` component of a given `LorentzVectorLike`, i.e. its 1-component. @@ -348,7 +345,6 @@ Return the ``p_x`` component of a given `LorentzVectorLike`, i.e. its 1-componen @inline @traitfn getPx(lv::T) where {T; IsLorentzVectorLike{T}} = getX(lv) """ - getPy(lv) Return the ``p_y`` component of a given `LorentzVectorLike`, i.e. its 2-component. @@ -361,7 +357,6 @@ Return the ``p_y`` component of a given `LorentzVectorLike`, i.e. its 2-componen @inline @traitfn getPy(lv::T) where {T; IsLorentzVectorLike{T}} = getY(lv) """ - getPz(lv) Return the ``p_z`` component of a given `LorentzVectorLike`, i.e. its 3-component. @@ -374,7 +369,6 @@ Return the ``p_z`` component of a given `LorentzVectorLike`, i.e. its 3-componen @inline @traitfn getPz(lv::T) where {T; IsLorentzVectorLike{T}} = getZ(lv) """ - getBeta(lv) Return magnitude of the beta vector for a given `LorentzVectorLike`, i.e. the magnitude of the `LorentzVectorLike` divided by its 0-component. @@ -399,7 +393,6 @@ Return magnitude of the beta vector for a given `LorentzVectorLike`, i.e. the ma end """ - getGamma(lv) Return the relativistic gamma factor for a given `LorentzVectorLike`, i.e. the inverse square root of one minus the beta vector squared. @@ -417,7 +410,6 @@ end # transverse coordinates ######################## """ - getTransverseMomentum2(lv) Return the squared transverse momentum for a given `LorentzVectorLike`, i.e. the sum of its squared 1- and 2-component. @@ -441,7 +433,6 @@ const getPt2 = getTransverseMomentum2 const getPerp2 = getTransverseMomentum2 """ - getTransverseMomentum(lv) Return the transverse momentum for a given `LorentzVectorLike`, i.e. the magnitude of its transverse components. @@ -466,7 +457,6 @@ const getPt = getTransverseMomentum const getPerp = getTransverseMomentum """ - getTransverseMass2(lv) Return the squared transverse mass for a given `LorentzVectorLike`, i.e. the difference of its squared 0- and 3-component. @@ -480,7 +470,6 @@ Return the squared transverse mass for a given `LorentzVectorLike`, i.e. the dif The transverse components are defined w.r.t. to the 3-axis. - """ @inline @traitfn function getTransverseMass2(lv::T) where {T; IsLorentzVectorLike{T}} return getT(lv)^2 - getZ(lv)^2 @@ -489,7 +478,6 @@ end const getMt2 = getTransverseMass2 """ - getTransverseMass(lv) Return the transverse momentum for a given `LorentzVectorLike`, i.e. the square root of its squared transverse mass. @@ -521,7 +509,6 @@ end const getMt = getTransverseMass """ - getRapidity(lv) Return the [rapidity](https://en.wikipedia.org/wiki/Rapidity) for a given `LorentzVectorLike`. @@ -550,7 +537,6 @@ const getRho2 = getMagnitude2 const getRho = getMagnitude """ - getTheta(lv) Return the theta angle for a given `LorentzVectorLike`, i.e. the polar angle of its spatial components in [spherical coordinates](https://en.wikipedia.org/wiki/Spherical_coordinate_system). @@ -575,7 +561,6 @@ Return the theta angle for a given `LorentzVectorLike`, i.e. the polar angle of end """ - getCosTheta(lv) Return the cosine of the theta angle for a given `LorentzVectorLike`. @@ -591,7 +576,6 @@ Return the cosine of the theta angle for a given `LorentzVectorLike`. end """ - getPhi(lv) Return the phi angle for a given `LorentzVectorLike`, i.e. the azimuthal angle of its spatial components in [spherical coordinates](https://en.wikipedia.org/wiki/Spherical_coordinate_system). @@ -611,7 +595,6 @@ Return the phi angle for a given `LorentzVectorLike`, i.e. the azimuthal angle o end """ - getCosPhi(lv) Return the cosine of the phi angle for a given `LorentzVectorLike`. @@ -627,7 +610,6 @@ Return the cosine of the phi angle for a given `LorentzVectorLike`. end """ - getSinPhi(lv) Return the sine of the phi angle for a given `LorentzVectorLike`. @@ -646,7 +628,6 @@ end # light cone coordinates ######################## """ - getPlus(lv) Return the plus component for a given `LorentzVectorLike` in [light-cone coordinates](https://en.wikipedia.org/wiki/Light-cone_coordinates). @@ -669,7 +650,6 @@ Return the plus component for a given `LorentzVectorLike` in [light-cone coordin end """ - getMinus(lv) Return the minus component for a given `LorentzVectorLike` in [light-cone coordinates](https://en.wikipedia.org/wiki/Light-cone_coordinates). @@ -698,7 +678,6 @@ end #### """ - setE!(lv,value) Sets the energy component of a given `LorentzVectorLike` to a given `value`. @@ -715,7 +694,6 @@ end const setEnergy! = setE! """ - setPx!(lv,value) Sets the 1-component of a given `LorentzVectorLike` to a given `value`. @@ -732,7 +710,6 @@ Sets the 1-component of a given `LorentzVectorLike` to a given `value`. end """ - setPy!(lv,value) Sets the 2-component of a given `LorentzVectorLike` to a given `value`. @@ -749,7 +726,6 @@ Sets the 2-component of a given `LorentzVectorLike` to a given `value`. end """ - setPz!(lv,value) Sets the 3-component of a given `LorentzVectorLike` to a given `value`. @@ -768,7 +744,6 @@ end # setter spherical coordinates """ - setTheta!(lv,value) Sets the theta angle of a `LorentzVectorLike` to a given `value`. @@ -791,7 +766,6 @@ Sets the theta angle of a `LorentzVectorLike` to a given `value`. end """ - setCosTheta!(lv,value) Sets the cosine of the theta angle of a `LorentzVectorLike` to a given `value`. @@ -816,7 +790,6 @@ Sets the cosine of the theta angle of a `LorentzVectorLike` to a given `value`. end """ - setPhi!(lv,value) Sets the phi angle of a `LorentzVectorLike` to a given `value`. @@ -840,7 +813,6 @@ Sets the phi angle of a `LorentzVectorLike` to a given `value`. end """ - setRho!(lv,value) Sets the magnitude of a `LorentzVectorLike` to a given `value`. @@ -864,7 +836,6 @@ end # setter light cone coordinates """ - setPlus!(lv,value) Sets the plus component of a `LorentzVectorLike` to a given `value`. @@ -882,7 +853,6 @@ Sets the plus component of a `LorentzVectorLike` to a given `value`. end """ - setMinus!(lv,value) Sets the minus component of a `LorentzVectorLike` to a given `value`. @@ -901,7 +871,6 @@ end # transverse coordinates """ - setTransverseMomentum!(lv,value) Sets the transverse momentum of a `LorentzVectorLike` to a given `value`. @@ -928,7 +897,6 @@ const setPerp! = setTransverseMomentum! const setPt! = setTransverseMomentum! """ - setTransverseMass!(lv,value) Sets the transverse mass of a `LorentzVectorLike` to a given `value`. @@ -953,7 +921,6 @@ end const setMt! = setTransverseMass! """ - setRapidity!(lv,value) Sets the rapidity of a `LorentzVectorLike` to a given `value`. diff --git a/src/particles/states.jl b/src/interfaces/particle_functions.jl similarity index 58% rename from src/particles/states.jl rename to src/interfaces/particle_functions.jl index 4729939..9d89943 100644 --- a/src/particles/states.jl +++ b/src/interfaces/particle_functions.jl @@ -1,10 +1,3 @@ -function _booster_fermion(mom::QEDbase.AbstractFourMomentum, mass::Real) - return (slashed(mom) + mass * one(DiracMatrix)) / (sqrt(abs(getT(mom)) + mass)) -end - -function _booster_antifermion(mom::QEDbase.AbstractFourMomentum, mass::Real) - return (mass * one(DiracMatrix) - slashed(mom)) / (sqrt(abs(getT(mom)) + mass)) -end """ ```julia @@ -68,7 +61,7 @@ electron_state = base_state(Electron(), Incoming(), mom, SpinUp()) ``` ```jldoctest -julia> using QEDbase +julia> using QEDbase; using QEDcore julia> mass = 1.0; px,py,pz = (0.1, 0.2, 0.3); E = sqrt(px^2 + py^2 + pz^2 + mass^2); mom = SFourMomentum(E, px, py, pz) 4-element SFourMomentum with indices SOneTo(4): @@ -150,133 +143,3 @@ julia> electron_states = base_state(Electron(), Incoming(), mom, AllSpin()) In the current implementation there are **no checks** built-in, which verify the passed arguments whether they describe on-shell particles, i.e. `p*p≈mass^2`. Using `base_state` with off-shell particles will cause unpredictable behavior. """ function base_state end - -function base_state( - particle::Fermion, - ::Incoming, - mom::QEDbase.AbstractFourMomentum, - spin::AbstractDefiniteSpin, -) - booster = _booster_fermion(mom, mass(particle)) - return BiSpinor(booster[:, _spin_index(spin)]) -end - -function base_state( - particle::Fermion, ::Incoming, mom::QEDbase.AbstractFourMomentum, spin::AllSpin -) - booster = _booster_fermion(mom, mass(particle)) - return SVector(BiSpinor(booster[:, 1]), BiSpinor(booster[:, 2])) -end - -function base_state( - particle::AntiFermion, - ::Incoming, - mom::QEDbase.AbstractFourMomentum, - spin::AbstractDefiniteSpin, -) - booster = _booster_antifermion(mom, mass(particle)) - return AdjointBiSpinor(BiSpinor(booster[:, _spin_index(spin) + 2])) * GAMMA[1] -end - -function base_state( - particle::AntiFermion, ::Incoming, mom::QEDbase.AbstractFourMomentum, spin::AllSpin -) - booster = _booster_antifermion(mom, mass(particle)) - return SVector( - AdjointBiSpinor(BiSpinor(booster[:, 3])) * GAMMA[1], - AdjointBiSpinor(BiSpinor(booster[:, 4])) * GAMMA[1], - ) -end - -function base_state( - particle::Fermion, - ::Outgoing, - mom::QEDbase.AbstractFourMomentum, - spin::AbstractDefiniteSpin, -) - booster = _booster_fermion(mom, mass(particle)) - return AdjointBiSpinor(BiSpinor(booster[:, _spin_index(spin)])) * GAMMA[1] -end - -function base_state( - particle::Fermion, ::Outgoing, mom::QEDbase.AbstractFourMomentum, spin::AllSpin -) - booster = _booster_fermion(mom, mass(particle)) - return SVector( - AdjointBiSpinor(BiSpinor(booster[:, 1])) * GAMMA[1], - AdjointBiSpinor(BiSpinor(booster[:, 2])) * GAMMA[1], - ) -end - -function base_state( - particle::AntiFermion, - ::Outgoing, - mom::QEDbase.AbstractFourMomentum, - spin::AbstractDefiniteSpin, -) - booster = _booster_antifermion(mom, mass(particle)) - return BiSpinor(booster[:, _spin_index(spin) + 2]) -end - -function base_state( - particle::AntiFermion, ::Outgoing, mom::QEDbase.AbstractFourMomentum, spin::AllSpin -) - booster = _booster_antifermion(mom, mass(particle)) - return SVector(BiSpinor(booster[:, 3]), BiSpinor(booster[:, 4])) -end - -function _photon_state(pol::AllPolarization, mom::QEDbase.AbstractFourMomentum) - cth = getCosTheta(mom) - sth = sqrt(1 - cth^2) - cos_phi = getCosPhi(mom) - sin_phi = getSinPhi(mom) - return SVector( - SLorentzVector{Float64}(0.0, cth * cos_phi, cth * sin_phi, -sth), - SLorentzVector{Float64}(0.0, -sin_phi, cos_phi, 0.0), - ) -end - -function _photon_state(pol::PolarizationX, mom::QEDbase.AbstractFourMomentum) - cth = getCosTheta(mom) - sth = sqrt(1 - cth^2) - cos_phi = getCosPhi(mom) - sin_phi = getSinPhi(mom) - return SLorentzVector{Float64}(0.0, cth * cos_phi, cth * sin_phi, -sth) -end - -function _photon_state(pol::PolarizationY, mom::QEDbase.AbstractFourMomentum) - cos_phi = getCosPhi(mom) - sin_phi = getSinPhi(mom) - return SLorentzVector{Float64}(0.0, -sin_phi, cos_phi, 0.0) -end - -@inline function base_state( - particle::Photon, - ::ParticleDirection, - mom::QEDbase.AbstractFourMomentum, - pol::AllPolarization, -) - return _photon_state(pol, mom) -end - -@inline function base_state( - particle::Photon, - ::ParticleDirection, - mom::QEDbase.AbstractFourMomentum, - pol::AbstractPolarization, -) - return _photon_state(pol, mom) -end - -""" - _as_svec(x) - -Accepts a single object, an `SVector` of objects or a tuple of objects, and returns them in a single "layer" of SVector. - -Intended for usage with [`base_state`](@ref). -""" -function _as_svec end - -@inline _as_svec(x) = SVector((x,)) -@inline _as_svec(x::SVector{N,T}) where {N,T} = x -@inline _as_svec(x::NTuple) = SVector(x) diff --git a/src/particles/types.jl b/src/interfaces/particle_types.jl similarity index 98% rename from src/particles/types.jl rename to src/interfaces/particle_types.jl index 01573cd..121327d 100644 --- a/src/particles/types.jl +++ b/src/interfaces/particle_types.jl @@ -5,6 +5,13 @@ # implement the abstact particle interface accordingly. ############### +""" + AbstractParticleSpinor + +TBW +""" +abstract type AbstractParticleSpinor end + """ AbstractParticleType <: AbstractParticle diff --git a/src/lorentz_vector.jl b/src/lorentz_vector.jl deleted file mode 100644 index 5e26e14..0000000 --- a/src/lorentz_vector.jl +++ /dev/null @@ -1,154 +0,0 @@ -""" -LorentzVectors -""" - -####### -# -# Abstract types -# -####### -""" -$(TYPEDEF) - -Abstract type to model generic Lorentz vectors, i.e. elements of a minkowski-like space, where the component space is arbitray. -""" -abstract type AbstractLorentzVector{T} <: FieldVector{4,T} end - -#interface with dirac tensors - -""" -$(TYPEDSIGNATURES) - -Product of generic Lorentz vector with a Dirac tensor from the left. Basically, the multiplication is piped to the components from the Lorentz vector. - -!!! note "Multiplication operator" - This also overloads the `*` operator for this types. -""" -function mul( - DM::T, L::TL -) where {T<:Union{AbstractDiracMatrix,AbstractDiracVector},TL<:AbstractLorentzVector} - return constructorof(TL)(DM * L[1], DM * L[2], DM * L[3], DM * L[4]) -end -@inline function *( - DM::T, L::TL -) where {T<:Union{AbstractDiracMatrix,AbstractDiracVector},TL<:AbstractLorentzVector} - return mul(DM, L) -end - -""" -$(TYPEDSIGNATURES) - -Product of generic Lorentz vector with a Dirac tensor from the right. Basically, the multiplication is piped to the components from the Lorentz vector. - -!!! note "Multiplication operator" - This also overloads the `*` operator for this types. - -""" -function mul( - L::TL, DM::T -) where {TL<:AbstractLorentzVector,T<:Union{AbstractDiracMatrix,AbstractDiracVector}} - return constructorof(TL)(L[1] * DM, L[2] * DM, L[3] * DM, L[4] * DM) -end -@inline function *( - L::TL, DM::T -) where {TL<:AbstractLorentzVector,T<:Union{AbstractDiracMatrix,AbstractDiracVector}} - return mul(L, DM) -end - -####### -# -# Concrete LorentzVector types -# -####### -""" -$(TYPEDEF) - -Concrete implementation of a generic static Lorentz vector. Each manipulation of an concrete implementation which is not self-contained (i.e. produces the same Lorentz vector type) will result in this type. - -# Fields -$(TYPEDFIELDS) -""" -struct SLorentzVector{T} <: AbstractLorentzVector{T} - "`t` component" - t::T - - "`x` component" - x::T - - "`y` component" - y::T - - "`z` component" - z::T -end -SLorentzVector(t, x, y, z) = SLorentzVector(promote(t, x, y, z)...) - -function similar_type(::Type{A}, ::Type{T}, ::Size{S}) where {A<:SLorentzVector,T,S} - return SLorentzVector{T} -end - -@inline getT(lv::SLorentzVector) = lv.t -@inline getX(lv::SLorentzVector) = lv.x -@inline getY(lv::SLorentzVector) = lv.y -@inline getZ(lv::SLorentzVector) = lv.z - -register_LorentzVectorLike(SLorentzVector) - -""" -$(TYPEDEF) - -Concrete implementation of a generic mutable Lorentz vector. Each manipulation of an concrete implementation which is not self-contained (i.e. produces the same Lorentz vector type) will result in this type. - -# Fields -$(TYPEDFIELDS) -""" -mutable struct MLorentzVector{T} <: AbstractLorentzVector{T} - "`t` component" - t::T - - "`x` component" - x::T - - "`y` component" - y::T - - "`z` component" - z::T -end -MLorentzVector(t, x, y, z) = MLorentzVector(promote(t, x, y, z)...) - -function similar_type(::Type{A}, ::Type{T}, ::Size{S}) where {A<:MLorentzVector,T,S} - return MLorentzVector{T} -end - -@inline getT(lv::MLorentzVector) = lv.t -@inline getX(lv::MLorentzVector) = lv.x -@inline getY(lv::MLorentzVector) = lv.y -@inline getZ(lv::MLorentzVector) = lv.z - -function QEDbase.setT!(lv::MLorentzVector, value::T) where {T} - return lv.t = value -end - -function QEDbase.setX!(lv::MLorentzVector, value::T) where {T} - return lv.x = value -end - -function QEDbase.setY!(lv::MLorentzVector, value::T) where {T} - return lv.y = value -end - -function QEDbase.setZ!(lv::MLorentzVector, value::T) where {T} - return lv.z = value -end - -register_LorentzVectorLike(MLorentzVector) - -function dot(p1::T1, p2::T2) where {T1<:AbstractLorentzVector,T2<:AbstractLorentzVector} - return mdot(p1, p2) -end -@inline function *( - p1::T1, p2::T2 -) where {T1<:AbstractLorentzVector,T2<:AbstractLorentzVector} - return dot(p1, p2) -end diff --git a/src/particles/spinors.jl b/src/particles/spinors.jl deleted file mode 100644 index 56fb1c3..0000000 --- a/src/particles/spinors.jl +++ /dev/null @@ -1,134 +0,0 @@ - -import Base.getindex - -const SPINOR_VALIDITY_CHECK = Ref(true) - -macro valid_spinor_input(ex) - return quote - SPINOR_VALIDITY_CHECK.x = false - local val = $(esc(ex)) - SPINOR_VALIDITY_CHECK.x = true - val - end -end - -const BASE_PARTICLE_SPINOR_UP = BiSpinor(1.0, 0.0, 0.0, 0.0) -const BASE_PARTICLE_SPINOR_DOWN = BiSpinor(0.0, 1.0, 0.0, 0.0) -const BASE_PARTICLE_SPINOR = [BASE_PARTICLE_SPINOR_UP, BASE_PARTICLE_SPINOR_DOWN] -const BASE_ANTIPARTICLE_SPINOR_UP = BiSpinor(0.0, 0.0, 1.0, 0.0) -const BASE_ANTIPARTICLE_SPINOR_DOWN = BiSpinor(0.0, 0.0, 0.0, 1.0) -const BASE_ANTIPARTICLE_SPINOR = [ - BASE_ANTIPARTICLE_SPINOR_UP, BASE_ANTIPARTICLE_SPINOR_DOWN -] - -mutable struct SpinorConstructionError <: Exception - var::String -end - -function Base.showerror(io::IO, e::SpinorConstructionError) - return print(io, "SpinorConstructionError: ", e.var) -end - -@inline function _check_spinor_input( - mom::T, mass::Float64 -) where {T<:AbstractLorentzVector{TE}} where {TE<:Real} - if SPINOR_VALIDITY_CHECK[] && !isonshell(mom, mass) - throw( - SpinorConstructionError( - "P^2 = $(getMass2(mom)) needs to be equal to mass^2=$(mass^2)" - ), - ) - end -end - -abstract type AbstractParticleSpinor end - -# -# fermion spinors -# - -function _build_particle_booster( - mom::T, mass::Float64 -) where {T<:AbstractLorentzVector{TE}} where {TE<:Real} - _check_spinor_input(mom, mass) - return (slashed(mom) + mass * one(DiracMatrix)) / (sqrt(abs(mom.t) + mass)) -end - -struct IncomingFermionSpinor <: AbstractParticleSpinor - booster::DiracMatrix -end - -function IncomingFermionSpinor( - mom::T, mass::Float64 -) where {T<:AbstractLorentzVector{TE}} where {TE<:Real} - return IncomingFermionSpinor(_build_particle_booster(mom, mass)) -end - -function (SP::IncomingFermionSpinor)(spin::Int64) - return SP.booster * BASE_PARTICLE_SPINOR[spin] -end - -const SpinorU = IncomingFermionSpinor - -struct OutgoingFermionSpinor <: AbstractParticleSpinor - booster::DiracMatrix -end - -function OutgoingFermionSpinor( - mom::T, mass::Float64 -) where {T<:AbstractLorentzVector{TE}} where {TE<:Real} - return OutgoingFermionSpinor(_build_particle_booster(mom, mass)) -end - -function (SP::OutgoingFermionSpinor)(spin::Int64) - return AdjointBiSpinor(SP.booster * BASE_PARTICLE_SPINOR[spin]) * GAMMA[1] -end - -const SpinorUbar = OutgoingFermionSpinor - -# -# Anti fermion spinors -# - -function _build_antiparticle_booster( - mom::T, mass::Float64 -) where {T<:AbstractLorentzVector{TE}} where {TE<:Real} - _check_spinor_input(mom, mass) - return (mass * one(DiracMatrix) - slashed(mom)) / (sqrt(abs(mom.t) + mass)) -end - -struct OutgoingAntiFermionSpinor <: AbstractParticleSpinor - booster::DiracMatrix -end - -function OutgoingAntiFermionSpinor( - mom::T, mass::Float64 -) where {T<:AbstractLorentzVector{TE}} where {TE<:Real} - return OutgoingAntiFermionSpinor(_build_antiparticle_booster(mom, mass)) -end - -function (SP::OutgoingAntiFermionSpinor)(spin::Int64) - return SP.booster * BASE_ANTIPARTICLE_SPINOR[spin] -end - -const SpinorV = OutgoingAntiFermionSpinor - -struct IncomingAntiFermionSpinor <: AbstractParticleSpinor - booster::DiracMatrix -end - -function IncomingAntiFermionSpinor( - mom::T, mass::Float64 -) where {T<:AbstractLorentzVector{TE}} where {TE<:Real} - return IncomingAntiFermionSpinor(_build_antiparticle_booster(mom, mass)) -end - -function (SP::IncomingAntiFermionSpinor)(spin::Int64) - return AdjointBiSpinor(SP.booster * BASE_ANTIPARTICLE_SPINOR[spin]) * GAMMA[1] -end - -const SpinorVbar = IncomingAntiFermionSpinor - -function getindex(SP::T, idx) where {T<:AbstractParticleSpinor} - return idx in (1, 2) ? SP(idx) : throw(BoundsError()) -end diff --git a/src/utils.jl b/src/utils.jl new file mode 100644 index 0000000..763abe5 --- /dev/null +++ b/src/utils.jl @@ -0,0 +1,12 @@ +""" + _as_svec(x) + +Accepts a single object, an `SVector` of objects or a tuple of objects, and returns them in a single "layer" of SVector. + +Useful with [`base_state`](@ref). +""" +function _as_svec end + +@inline _as_svec(x) = SVector((x,)) +@inline _as_svec(x::SVector{N,T}) where {N,T} = x +@inline _as_svec(x::NTuple) = SVector(x) diff --git a/test/four_momentum.jl b/test/four_momentum.jl deleted file mode 100644 index 2ee008b..0000000 --- a/test/four_momentum.jl +++ /dev/null @@ -1,212 +0,0 @@ - -using QEDbase -using Random - -const ATOL = 1e-15 - -@testset "FourMomentum getter" for MomentumType in [SFourMomentum, MFourMomentum] - rng = MersenneTwister(12345) - x, y, z = rand(rng, 3) - mass = rand(rng) - E = sqrt(x^2 + y^2 + z^2 + mass^2) - mom_onshell = MomentumType(E, x, y, z) - mom_zero = MomentumType(0.0, 0.0, 0.0, 0.0) - mom_offshell = MomentumType(0.0, 0.0, 0.0, mass) - - @testset "magnitude consistence" for mom in [mom_onshell, mom_offshell, mom_zero] - @test getMagnitude2(mom) == getMag2(mom) - @test getMagnitude(mom) == getMag(mom) - @test isapprox(getMagnitude(mom), sqrt(getMagnitude2(mom))) - end - - @testset "magnitude values" begin - @test isapprox(getMagnitude2(mom_onshell), x^2 + y^2 + z^2) - @test isapprox(getMagnitude(mom_onshell), sqrt(x^2 + y^2 + z^2)) - end - - @testset "mass consistence" for mom_on in [mom_onshell, mom_zero] - @test getInvariantMass2(mom_on) == getMass2(mom_on) - @test getInvariantMass(mom_on) == getMass(mom_on) - @test isapprox(getInvariantMass(mom_on), sqrt(getInvariantMass2(mom_on))) - end - - @testset "mass value" begin - @test isapprox(getInvariantMass2(mom_onshell), E^2 - (x^2 + y^2 + z^2)) - @test isapprox(getInvariantMass(mom_onshell), sqrt(E^2 - (x^2 + y^2 + z^2))) - - @test isapprox(getInvariantMass(mom_onshell), mass) - @test isapprox(getInvariantMass(mom_offshell), -mass) - @test isapprox(getInvariantMass(mom_zero), 0.0) - end - - @testset "momentum components" begin - @test getE(mom_onshell) == E - @test getEnergy(mom_onshell) == getE(mom_onshell) - @test getPx(mom_onshell) == x - @test getPy(mom_onshell) == y - @test getPz(mom_onshell) == z - - @test isapprox(getBeta(mom_onshell), sqrt(x^2 + y^2 + z^2) / E) - @test isapprox(getGamma(mom_onshell), 1 / sqrt(1.0 - getBeta(mom_onshell)^2)) - - @test getE(mom_zero) == 0.0 - @test getEnergy(mom_zero) == 0.0 - @test getPx(mom_zero) == 0.0 - @test getPy(mom_zero) == 0.0 - @test getPz(mom_zero) == 0.0 - - @test isapprox(getBeta(mom_zero), 0.0) - @test isapprox(getGamma(mom_zero), 1.0) - end - - @testset "transverse coordinates" for mom_on in [mom_onshell, mom_zero] - @test getTransverseMomentum2(mom_on) == getPt2(mom_on) - @test getTransverseMomentum2(mom_on) == getPerp2(mom_on) - @test getTransverseMomentum(mom_on) == getPt(mom_on) - @test getTransverseMomentum(mom_on) == getPerp(mom_on) - - @test isapprox(getPt(mom_on), sqrt(getPt2(mom_on))) - - @test getTransverseMass2(mom_on) == getMt2(mom_on) - @test getTransverseMass(mom_on) == getMt(mom_on) - end - - @testset "transverse coordiantes value" begin - @test isapprox(getTransverseMomentum2(mom_onshell), x^2 + y^2) - @test isapprox(getTransverseMomentum(mom_onshell), sqrt(x^2 + y^2)) - @test isapprox(getTransverseMass2(mom_onshell), E^2 - z^2) - @test isapprox(getTransverseMass(mom_onshell), sqrt(E^2 - z^2)) - @test isapprox(getMt(mom_offshell), -mass) - @test isapprox(getRapidity(mom_onshell), 0.5 * log((E + z) / (E - z))) - - @test isapprox(getTransverseMomentum2(mom_zero), 0.0) - @test isapprox(getTransverseMomentum(mom_zero), 0.0) - @test isapprox(getTransverseMass2(mom_zero), 0.0) - @test isapprox(getTransverseMass(mom_zero), 0.0) - @test isapprox(getMt(mom_zero), 0.0) - end - - @testset "spherical coordiantes consistence" for mom_on in [mom_onshell, mom_zero] - @test getRho2(mom_on) == getMagnitude2(mom_on) - @test getRho(mom_on) == getMagnitude(mom_on) - - @test isapprox(getCosTheta(mom_on), cos(getTheta(mom_on))) - @test isapprox(getCosPhi(mom_on), cos(getPhi(mom_on))) - @test isapprox(getSinPhi(mom_on), sin(getPhi(mom_on))) - end - - @testset "spherical coordiantes values" begin - @test isapprox(getTheta(mom_onshell), atan(getPt(mom_onshell), z)) - @test isapprox(getTheta(mom_zero), 0.0) - - @test isapprox(getPhi(mom_onshell), atan(y, x)) - @test isapprox(getPhi(mom_zero), 0.0) - end - - @testset "light-cone coordiantes" begin - @test isapprox(getPlus(mom_onshell), 0.5 * (E + z)) - @test isapprox(getMinus(mom_onshell), 0.5 * (E - z)) - - @test isapprox(getPlus(mom_zero), 0.0) - @test isapprox(getMinus(mom_zero), 0.0) - end -end # FourMomentum getter - -function test_get_set(rng, setter, getter; value=rand(rng)) - x, y, z = rand(rng, 3) - mass = rand(rng) - E = sqrt(x^2 + y^2 + z^2 + mass^2) - mom = MFourMomentum(E, x, y, z) - setter(mom, value) - return isapprox(getter(mom), value) -end - -@testset "FourMomentum setter" begin - rng = MersenneTwister(123456) - - @testset "Momentum components" begin - @test test_get_set(rng, setE!, getE) - @test test_get_set(rng, setEnergy!, getE) - @test test_get_set(rng, setPx!, getPx) - @test test_get_set(rng, setPy!, getPy) - @test test_get_set(rng, setPz!, getPz) - end - - @testset "spherical coordiantes" begin - @test test_get_set(rng, setTheta!, getTheta) - @test test_get_set(rng, setTheta!, getTheta, value=0.0) - @test test_get_set(rng, setCosTheta!, getCosTheta) - @test test_get_set(rng, setCosTheta!, getCosTheta, value=1.0) - @test test_get_set(rng, setPhi!, getPhi) - @test test_get_set(rng, setPhi!, getPhi, value=0.0) - @test test_get_set(rng, setRho!, getRho) - @test test_get_set(rng, setRho!, getRho, value=0.0) - end - - @testset "light-cone coordiantes" begin - @test test_get_set(rng, setPlus!, getPlus) - @test test_get_set(rng, setPlus!, getPlus, value=0.0) - @test test_get_set(rng, setMinus!, getMinus) - @test test_get_set(rng, setMinus!, getMinus, value=0.0) - end - - @testset "transverse coordinates" begin - @test test_get_set(rng, setTransverseMomentum!, getTransverseMomentum) - @test test_get_set(rng, setTransverseMomentum!, getTransverseMomentum, value=0.0) - @test test_get_set(rng, setPerp!, getTransverseMomentum) - @test test_get_set(rng, setPt!, getTransverseMomentum) - @test test_get_set(rng, setTransverseMass!, getTransverseMass) - @test test_get_set(rng, setTransverseMass!, getTransverseMass, value=0.0) - @test test_get_set(rng, setMt!, getTransverseMass) - @test test_get_set(rng, setRapidity!, getRapidity) - @test test_get_set(rng, setRapidity!, getRapidity, value=0.0) - end -end # FourMomentum setter - -const SCALE = 10.0 .^ [-9, 0, 5] -const M_MASSIVE = 1.0 -const M_MASSLESS = 0.0 - -const M_ABSERR = 0.01 -const M_RELERR = 0.0001 - -@testset "isonshell" begin - rng = MersenneTwister(42) - x_base, y_base, z_base = rand(rng, 3) - - @testset "correct onshell" begin - @testset "($x_scale, $y_scale, $z_scale)" for (x_scale, y_scale, z_scale) in - Iterators.product(SCALE, SCALE, SCALE) - x, y, z = x_base * x_scale, y_base * y_scale, z_base * z_scale - E_massless = sqrt(x^2 + y^2 + z^2 + M_MASSLESS^2) - E_massive = sqrt(x^2 + y^2 + z^2 + M_MASSIVE^2) - mom_massless = SFourMomentum(E_massless, x, y, z) - mom_massive = SFourMomentum(E_massive, x, y, z) - @test isonshell(mom_massless, M_MASSLESS) - @test isonshell(mom_massive, M_MASSIVE) - - @test assert_onshell(mom_massless, M_MASSLESS) == nothing - @test assert_onshell(mom_massive, M_MASSIVE) == nothing - end - end - - @testset "correct not onshell" begin - @testset "$x_scale, $y_scale, $z_scale" for (x_scale, y_scale, z_scale) in - Iterators.product(SCALE, SCALE, SCALE) - x, y, z = x_base * x_scale, y_base * y_scale, z_base * z_scale - m_err = min(M_ABSERR, M_RELERR * sum([x, y, z]) / 3.0) # mass error is M_RELERR of the mean of the components - # but has at most the value M_ABSERR - - E_massless = sqrt(x^2 + y^2 + z^2 + (M_MASSLESS + m_err)^2) - E_massive = sqrt(x^2 + y^2 + z^2 + (M_MASSIVE + m_err)^2) - mom_massless = SFourMomentum(E_massless, x, y, z) - mom_massive = SFourMomentum(E_massive, x, y, z) - - @test !isonshell(mom_massless, M_MASSLESS) - @test !isonshell(mom_massive, M_MASSIVE) - - @test_throws QEDbase.OnshellError assert_onshell(mom_massless, M_MASSLESS) - @test_throws QEDbase.OnshellError assert_onshell(mom_massive, M_MASSIVE) - end - end -end diff --git a/test/gamma_matrices.jl b/test/gamma_matrices.jl deleted file mode 100644 index e6f91a7..0000000 --- a/test/gamma_matrices.jl +++ /dev/null @@ -1,112 +0,0 @@ -using QEDbase -using LinearAlgebra -using Random -using SparseArrays - -function _gamma_anticommutator(i, j) - return GAMMA[i] * GAMMA[j] + GAMMA[j] * GAMMA[i] -end - -function _gamma_sandwich(nu::Int64) - return GAMMA * (GAMMA[nu] * GAMMA) -end - -METRIC = Diagonal([1, -1, -1, -1]) -EYE = one(DiracMatrix) -GROUNDTRUTH_GAMMA0_DIRAC = spzeros(4, 4) -GROUNDTRUTH_GAMMA0_DIRAC[1, 1] = 1 -GROUNDTRUTH_GAMMA0_DIRAC[2, 2] = 1 -GROUNDTRUTH_GAMMA0_DIRAC[3, 3] = -1 -GROUNDTRUTH_GAMMA0_DIRAC[4, 4] = -1 - -GROUNDTRUTH_GAMMA1_DIRAC = spzeros(4, 4) -GROUNDTRUTH_GAMMA1_DIRAC[4, 1] = -1 -GROUNDTRUTH_GAMMA1_DIRAC[3, 2] = -1 -GROUNDTRUTH_GAMMA1_DIRAC[2, 3] = 1 -GROUNDTRUTH_GAMMA1_DIRAC[1, 4] = 1 - -GROUNDTRUTH_GAMMA2_DIRAC = spzeros(ComplexF64, 4, 4) -GROUNDTRUTH_GAMMA2_DIRAC[4, 1] = -1im -GROUNDTRUTH_GAMMA2_DIRAC[3, 2] = 1im -GROUNDTRUTH_GAMMA2_DIRAC[2, 3] = 1im -GROUNDTRUTH_GAMMA2_DIRAC[1, 4] = -1im - -GROUNDTRUTH_GAMMA3_DIRAC = spzeros(4, 4) -GROUNDTRUTH_GAMMA3_DIRAC[3, 1] = -1 -GROUNDTRUTH_GAMMA3_DIRAC[4, 2] = 1 -GROUNDTRUTH_GAMMA3_DIRAC[1, 3] = 1 -GROUNDTRUTH_GAMMA3_DIRAC[2, 4] = -1 - -@testset "gamma matrices" begin - rng = MersenneTwister(42) - - @testset "commutator" begin - for (i, j) in Iterators.product(1:4, 1:4) - @test isapprox(_gamma_anticommutator(i, j), 2 * METRIC[i, j] * EYE) - end - end #commutator - - @testset "Minkowski product" begin - @test GAMMA * GAMMA == 4 * one(DiracMatrix) - end # Minkowski product - - @testset "gamma sandwich" begin - @test SLorentzVector([_gamma_sandwich(nu) for nu in 1:4]) == -2 * GAMMA - end # gamma sandwich - - @testset "adjoints" begin - @test adjoint(GAMMA[1]) == GAMMA[1] - for i in 2:4 - @test adjoint(GAMMA[i]) == -GAMMA[i] - end - - @test SLorentzVector([adjoint(GAMMA[nu]) for nu in 1:4]) == GAMMA[1] * (GAMMA * GAMMA[1]) - end # adjoints - - @testset "interface SLorentzVector" begin - a = SLorentzVector(rand(rng, 4)) - b = SLorentzVector(rand(rng, 4)) - - a_slash = GAMMA * a - b_slash = GAMMA * b - - @test isapprox(tr(a_slash * b_slash), 4 * a * b) - @test isapprox(a_slash * a_slash, (a * a) * one(DiracMatrix)) - @test isapprox(GAMMA * (a_slash * GAMMA), -2 * a_slash) - end # interface SLorentzVector - - @testset "Feynman slash" begin - a = SLorentzVector(rand(rng, 4) + 1im * rand(rng, 4)) - - @test isapprox(slashed(a), GAMMA * a) - @test isapprox(slashed(a), slashed(DiracGammaRepresentation, a)) - end - - @testset "Dirac representation" begin - # check the components of the gamma matrices against the - # Dirac representations, e.g. from https://en.wikipedia.org/wiki/Gamma_matrices - # note: we use the convention of lower indices for the gamma matrix definition. - # This motivates the minus sign in front of the spatial components - comps = 1:4 - @testset "gamma_0" begin - @testset "($col,$row)" for (row, col) in Iterators.product(comps, comps) - @test isapprox(GAMMA[1][row, col], GROUNDTRUTH_GAMMA0_DIRAC[row, col]) - end - end - @testset "gamma_1" begin - @testset "($col,$row)" for (row, col) in Iterators.product(comps, comps) - @test isapprox(GAMMA[2][row, col], -GROUNDTRUTH_GAMMA1_DIRAC[row, col]) - end - end - @testset "gamma_2" begin - @testset "($col,$row)" for (row, col) in Iterators.product(comps, comps) - @test isapprox(GAMMA[3][row, col], -GROUNDTRUTH_GAMMA2_DIRAC[row, col]) - end - end - @testset "gamma_3" begin - @testset "($col,$row)" for (row, col) in Iterators.product(comps, comps) - @test isapprox(GAMMA[4][row, col], -GROUNDTRUTH_GAMMA3_DIRAC[row, col]) - end - end - end -end diff --git a/test/interfaces/model.jl b/test/interfaces/model.jl new file mode 100644 index 0000000..cdb8887 --- /dev/null +++ b/test/interfaces/model.jl @@ -0,0 +1,19 @@ +using QEDbase + +include("../test_implementation/TestImplementation.jl") + +@testset "hard interface" begin + TESTMODEL = TestImplementation.TestModel() + @test fundamental_interaction_type(TESTMODEL) == :test_interaction +end + +@testset "interface fail" begin + TESTMODEL_FAIL = TestImplementation.TestModel_FAIL() + @test_throws MethodError fundamental_interaction_type(TESTMODEL_FAIL) +end + +@testset "broadcast" begin + test_func(model) = model + TESTMODEL = TestImplementation.TestModel() + @test test_func.(TESTMODEL) == TESTMODEL +end diff --git a/test/interfaces/process.jl b/test/interfaces/process.jl new file mode 100644 index 0000000..3ecc139 --- /dev/null +++ b/test/interfaces/process.jl @@ -0,0 +1,129 @@ +using Random +using QEDbase + +RNG = MersenneTwister(137137) +ATOL = 0.0 +RTOL = sqrt(eps()) + +include("../test_implementation/TestImplementation.jl") + +@testset "($N_INCOMING,$N_OUTGOING)" for (N_INCOMING, N_OUTGOING) in Iterators.product( + (1, rand(RNG, 2:8)), (1, rand(RNG, 2:8)) +) + INCOMING_PARTICLES = Tuple(rand(RNG, TestImplementation.PARTICLE_SET, N_INCOMING)) + OUTGOING_PARTICLES = Tuple(rand(RNG, TestImplementation.PARTICLE_SET, N_OUTGOING)) + + TESTPROC = TestImplementation.TestProcess(INCOMING_PARTICLES, OUTGOING_PARTICLES) + TESTMODEL = TestImplementation.TestModel() + TESTPSDEF = TestImplementation.TestPhasespaceDef() + IN_PS = TestImplementation._rand_momenta(RNG, N_INCOMING) + OUT_PS = TestImplementation._rand_momenta(RNG, N_OUTGOING) + PSP = PhaseSpacePoint(TESTPROC, TESTMODEL, TESTPSDEF, IN_PS, OUT_PS) + + @testset "failed interface" begin + TESTPROC_FAIL_ALL = TestImplementation.TestProcess_FAIL_ALL( + INCOMING_PARTICLES, OUTGOING_PARTICLES + ) + TESTPROC_FAIL_DIFFCS = TestImplementation.TestProcess_FAIL_DIFFCS( + INCOMING_PARTICLES, OUTGOING_PARTICLES + ) + TESTMODEL_FAIL = TestImplementation.TestModel_FAIL() + TESTPSDEF_FAIL = TestImplementation.TestPhasespaceDef_FAIL() + + @testset "failed process interface" begin + @test_throws MethodError incoming_particles(TESTPROC_FAIL_ALL) + @test_throws MethodError outgoing_particles(TESTPROC_FAIL_ALL) + end + + @testset "$PROC $MODEL" for (PROC, MODEL) in Iterators.product( + (TESTPROC, TESTPROC_FAIL_DIFFCS), (TESTMODEL, TESTMODEL_FAIL) + ) + if TestImplementation._any_fail(PROC, MODEL) + psp = PhaseSpacePoint(PROC, MODEL, TESTPSDEF, IN_PS, OUT_PS) + @test_throws MethodError QEDbase._incident_flux(psp) + @test_throws MethodError QEDbase._averaging_norm(psp) + @test_throws MethodError QEDbase._matrix_element(psp) + end + + for PS_DEF in (TESTPSDEF, TESTPSDEF_FAIL) + if TestImplementation._any_fail(PROC, MODEL, PS_DEF) + psp = PhaseSpacePoint(PROC, MODEL, PS_DEF, IN_PS, OUT_PS) + @test_throws MethodError QEDbase._phase_space_factor(psp) + end + end + end + end + + @testset "broadcast" begin + test_func(proc::AbstractProcessDefinition) = proc + @test test_func.(TESTPROC) == TESTPROC + + test_func(model::AbstractModelDefinition) = model + @test test_func.(TESTMODEL) == TESTMODEL + end + + @testset "incoming/outgoing particles" begin + @test incoming_particles(TESTPROC) == INCOMING_PARTICLES + @test outgoing_particles(TESTPROC) == OUTGOING_PARTICLES + @test number_incoming_particles(TESTPROC) == N_INCOMING + @test number_outgoing_particles(TESTPROC) == N_OUTGOING + end + + @testset "incident flux" begin + test_incident_flux = QEDbase._incident_flux( + InPhaseSpacePoint(TESTPROC, TESTMODEL, TESTPSDEF, IN_PS) + ) + groundtruth = TestImplementation._groundtruth_incident_flux(IN_PS) + @test isapprox(test_incident_flux, groundtruth, atol=ATOL, rtol=RTOL) + + test_incident_flux = QEDbase._incident_flux( + PhaseSpacePoint(TESTPROC, TESTMODEL, TESTPSDEF, IN_PS, OUT_PS) + ) + @test isapprox(test_incident_flux, groundtruth, atol=ATOL, rtol=RTOL) + + @test_throws MethodError QEDbase._incident_flux( + OutPhaseSpacePoint(TESTPROC, TESTMODEL, TESTPSDEF, OUT_PS) + ) + end + + @testset "averaging norm" begin + test_avg_norm = QEDbase._averaging_norm(TESTPROC) + groundtruth = TestImplementation._groundtruth_averaging_norm(TESTPROC) + @test isapprox(test_avg_norm, groundtruth, atol=ATOL, rtol=RTOL) + end + + @testset "matrix element" begin + test_matrix_element = QEDbase._matrix_element(PSP) + groundtruth = TestImplementation._groundtruth_matrix_element(IN_PS, OUT_PS) + @test length(test_matrix_element) == length(groundtruth) + for i in eachindex(test_matrix_element) + @test isapprox(test_matrix_element[i], groundtruth[i], atol=ATOL, rtol=RTOL) + end + end + + @testset "is in phasespace" begin + @test QEDbase._is_in_phasespace(PSP) + + IN_PS_unphysical = (zero(SFourMomentum), IN_PS[2:end]...) + OUT_PS_unphysical = (OUT_PS[1:(end - 1)]..., ones(SFourMomentum)) + PSP_unphysical_in_ps = PhaseSpacePoint( + TESTPROC, TESTMODEL, TESTPSDEF, IN_PS_unphysical, OUT_PS + ) + PSP_unphysical_out_ps = PhaseSpacePoint( + TESTPROC, TESTMODEL, TESTPSDEF, IN_PS, OUT_PS_unphysical + ) + PSP_unphysical = PhaseSpacePoint( + TESTPROC, TESTMODEL, TESTPSDEF, IN_PS_unphysical, OUT_PS_unphysical + ) + + @test !QEDbase._is_in_phasespace(PSP_unphysical_in_ps) + @test !QEDbase._is_in_phasespace(PSP_unphysical_out_ps) + @test !QEDbase._is_in_phasespace(PSP_unphysical) + end + + @testset "phase space factor" begin + test_phase_space_factor = QEDbase._phase_space_factor(PSP) + groundtruth = TestImplementation._groundtruth_phase_space_factor(IN_PS, OUT_PS) + @test isapprox(test_phase_space_factor, groundtruth, atol=ATOL, rtol=RTOL) + end +end diff --git a/test/interfaces/setup.jl b/test/interfaces/setup.jl new file mode 100644 index 0000000..44d5157 --- /dev/null +++ b/test/interfaces/setup.jl @@ -0,0 +1,188 @@ +using Random +using Suppressor +using QEDbase + +RNG = MersenneTwister(137137) +ATOL = 0.0 +RTOL = sqrt(eps()) + +_groundtruth_compute(x) = x +_groundtruth_input_validation(x) = (x > 0) +struct TestException <: QEDbase.AbstractInvalidInputException end +function _groundtruth_valid_input_assert(x) + _groundtruth_input_validation(x) || throw(TestException()) + return nothing +end +_transform_to_invalid(x) = -abs(x) +_groundtruth_post_processing(x, y) = x + y + +# setups for which the interface is implemented +abstract type AbstractTestSetup <: AbstractComputationSetup end +QEDbase._compute(stp::AbstractTestSetup, x) = _groundtruth_compute(x) + +# setup with default implementations +struct TestSetupDefault <: AbstractTestSetup end + +# setup with custom _assert_valid_input +struct TestSetupCustomAssertValidInput <: AbstractTestSetup end +function QEDbase._assert_valid_input(stp::TestSetupCustomAssertValidInput, x) + return _groundtruth_valid_input_assert(x) +end + +# setup with custom post processing +struct TestSetupCustomPostProcessing <: AbstractTestSetup end +function QEDbase._post_processing(::TestSetupCustomPostProcessing, x, y) + return _groundtruth_post_processing(x, y) +end + +# setup with custom input validation and post processing +struct TestSetupCustom <: AbstractTestSetup end +function QEDbase._assert_valid_input(stp::TestSetupCustom, x) + return _groundtruth_valid_input_assert(x) +end +QEDbase._post_processing(::TestSetupCustom, x, y) = _groundtruth_post_processing(x, y) + +# setup which fail on computation with default implementations +struct TestSetupFAIL <: AbstractComputationSetup end + +# setup which fail on computation with custom input validation, where the +# invalid input will be caught before the computation. +struct TestSetupCustomValidationFAIL <: AbstractComputationSetup end +function QEDbase._assert_valid_input(stp::TestSetupCustomValidationFAIL, x) + return _groundtruth_valid_input_assert(x) +end + +# setup which fail on computation with custom post processing +struct TestSetupCustomPostProcessingFAIL <: AbstractComputationSetup end +function QEDbase._post_processing(::TestSetupCustomPostProcessingFAIL, x, y) + return _groundtruth_post_processing(x, y) +end +@testset "general computation setup interface" begin + @testset "interface fail" begin + rnd_input = rand(RNG) + + @test_throws MethodError QEDbase._compute(TestSetupFAIL(), rnd_input) + @test_throws MethodError compute(TestSetupFAIL(), rnd_input) + + @test_throws MethodError QEDbase._compute( + TestSetupCustomValidationFAIL(), rnd_input + ) + @test_throws MethodError compute(TestSetupCustomValidationFAIL(), rnd_input) + # invalid input should be caught without throwing a MethodError + @test_throws TestException compute( + TestSetupCustomValidationFAIL(), _transform_to_invalid(rnd_input) + ) + + @test_throws MethodError QEDbase._compute( + TestSetupCustomPostProcessingFAIL(), rnd_input + ) + @test_throws MethodError compute(TestSetupCustomPostProcessingFAIL(), rnd_input) + end + + @testset "default interface" begin + stp = TestSetupDefault() + + rnd_input = rand(RNG) + rnd_output = rand(RNG) + @test QEDbase._post_processing(stp, rnd_input, rnd_output) == rnd_output + @test isapprox( + QEDbase._compute(stp, rnd_input), + _groundtruth_compute(rnd_input), + atol=ATOL, + rtol=RTOL, + ) + @test isapprox( + compute(stp, rnd_input), _groundtruth_compute(rnd_input), atol=ATOL, rtol=RTOL + ) + end + + @testset "custom input validation" begin + stp = TestSetupCustomAssertValidInput() + rnd_input = rand(RNG) + @test QEDbase._assert_valid_input(stp, rnd_input) == nothing + @test_throws TestException QEDbase._assert_valid_input( + stp, _transform_to_invalid(rnd_input) + ) + @test_throws TestException compute(stp, _transform_to_invalid(rnd_input)) + end + + @testset "custom post processing" begin + stp = TestSetupCustomPostProcessing() + rnd_input = rand(RNG) + rnd_output = rand(RNG) + @test isapprox( + QEDbase._post_processing(stp, rnd_input, rnd_output), + _groundtruth_post_processing(rnd_input, rnd_output), + ) + @test isapprox( + compute(stp, rnd_input), + _groundtruth_post_processing(rnd_input, _groundtruth_compute(rnd_input)), + ) + end + + @testset "custom input validation and post processing" begin + stp = TestSetupCustom() + rnd_input = rand(RNG) + rnd_output = rand(RNG) + + @test_throws TestException() compute(stp, _transform_to_invalid(rnd_input)) + @test isapprox( + QEDbase._post_processing(stp, rnd_input, rnd_output), + _groundtruth_post_processing(rnd_input, rnd_output), + ) + @test isapprox( + compute(stp, rnd_input), + _groundtruth_post_processing(rnd_input, _groundtruth_compute(rnd_input)), + ) + end +end +# process setup + +struct TestParticle1 <: AbstractParticle end +struct TestParticle2 <: AbstractParticle end +struct TestParticle3 <: AbstractParticle end +struct TestParticle4 <: AbstractParticle end + +PARTICLE_SET = [TestParticle1(), TestParticle2(), TestParticle3(), TestParticle4()] + +struct TestProcess <: AbstractProcessDefinition end +struct TestModel <: AbstractModelDefinition end + +struct TestProcessSetup <: AbstractProcessSetup end +QEDbase.scattering_process(::TestProcessSetup) = TestProcess() +QEDbase.physical_model(::TestProcessSetup) = TestModel() + +struct TestProcessSetupFAIL <: AbstractProcessSetup end + +@testset "process setup interface" begin + @testset "interface fail" begin + rnd_input = rand(RNG) + @test_throws MethodError scattering_process(TestProcessSetupFAIL()) + @test_throws MethodError physical_model(TestProcessSetupFAIL()) + @test_throws MethodError QEDbase._compute(TestProcessSetupFAIL(), rnd_input) + end + + @testset "hard interface" begin + stp = TestProcessSetup() + + @test QEDbase._is_computation_setup(stp) + @test scattering_process(stp) == TestProcess() + @test physical_model(stp) == TestModel() + end + + @testset "($N_INCOMING,$N_OUTGOING)" for (N_INCOMING, N_OUTGOING) in Iterators.product( + (1, rand(RNG, 2:8)), (1, rand(RNG, 2:8)) + ) + INCOMING_PARTICLES = rand(RNG, PARTICLE_SET, N_INCOMING) + OUTGOING_PARTICLES = rand(RNG, PARTICLE_SET, N_OUTGOING) + + @suppress QEDbase.incoming_particles(::TestProcess) = INCOMING_PARTICLES + @suppress QEDbase.outgoing_particles(::TestProcess) = OUTGOING_PARTICLES + + @testset "delegated functions" begin + stp = TestProcessSetup() + @test number_incoming_particles(stp) == N_INCOMING + @test number_outgoing_particles(stp) == N_OUTGOING + end + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 3860a86..67b8aa1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,27 +3,33 @@ using Test using SafeTestsets begin - @time @safetestset "Dirac tensors" begin - include("dirac_tensor.jl") + # # Interfaces + @time @safetestset "model interface" begin + include("interfaces/model.jl") end - @time @safetestset "Lorentz Vectors" begin - include("lorentz_vector.jl") + @time @safetestset "process interface" begin + include("interfaces/process.jl") + end + + @time @safetestset "computation setup interface" begin + include("interfaces/setup.jl") end @time @safetestset "Lorentz interface" begin include("lorentz_interface.jl") end - @time @safetestset "Gamma matrices" begin - include("gamma_matrices.jl") + + @time @safetestset "Dirac tensors" begin + include("dirac_tensor.jl") end - @time @safetestset "particle spinors" begin - include("particle_spinors.jl") + @time @safetestset "Lorentz Vectors" begin + include("lorentz_vector.jl") end - @time @safetestset "four momentum" begin - include("four_momentum.jl") + @time @safetestset "particle spinors" begin + include("particle_spinors.jl") end @time @safetestset "particles" begin diff --git a/test/test_implementation/TestImplementation.jl b/test/test_implementation/TestImplementation.jl new file mode 100644 index 0000000..4cd4aac --- /dev/null +++ b/test/test_implementation/TestImplementation.jl @@ -0,0 +1,38 @@ +""" +This module provides a full implementation of the model and process interface. Its purpose is only for testing and it does not reflect any +real-world physics. + +The module exports: + +``` +TestParticle1 # set of test particles without properties +TestParticle2 +TestParticle3 +TestParticle4 +TestModel # dummy compute model +TestModel_FAIL # failing compute model +TestProcess # dummy scattering process +TestProcess_FAIL # failing scattering process +TestPhasespaceDef # dummy phase space definition +TestPhasespaceDef_FAIL # failing phase space definition +``` +The respective groundtruth implementations for the interface functions are stored in `groundtruths.jl` +""" +module TestImplementation + +export TestParticle1, TestParticle2, TestParticle3, TestParticle4, PARTICLE_SET +export TestModel, TestModel_FAIL +export TestProcess, TestProcess_FAIL +export TestPhasespaceDef, TestPhasespaceDef_FAIL + +using Random +using QEDbase +using StaticArrays + +include("groundtruths.jl") +include("test_model.jl") +include("test_process.jl") +include("random_momenta.jl") +include("utils.jl") + +end diff --git a/test/test_implementation/groundtruths.jl b/test/test_implementation/groundtruths.jl new file mode 100644 index 0000000..f05a010 --- /dev/null +++ b/test/test_implementation/groundtruths.jl @@ -0,0 +1,265 @@ +import QEDbase.AbstractFourMomentum + +""" + _groundtruth_incident_flux(in_ps) + +Test implementation of the incident flux. Return the Minkowski square of the sum of the incoming momenta: + +```math +\\begin{align} +I = \\left(\\sum p_i\\right)^2, +\\end{align} +``` +where \$p_i\\in\\mathrm{ps_in}\$. +""" +function _groundtruth_incident_flux(in_ps) + s = sum(in_ps) + return s * s +end + +""" + _groundtruth_matrix_element(in_ps, out_ps) + +Test implementation for a matrix elements. Returns a list of three complex numbers without any physical meaning. +""" +function _groundtruth_matrix_element(in_ps, out_ps) + s_in = sum(in_ps) + s_out = sum(out_ps) + res = s_in * s_in + 1im * (s_out * s_out) + return (res, 2 * res, 3 * res) +end + +""" + _groundtruth_averaging_norm(proc) + +Test implementation of the averaging norm. Returns the inverse of the sum of all external particles of the passed process. +""" +function _groundtruth_averaging_norm(proc) + return 1.0 / (number_incoming_particles(proc) + number_outgoing_particles(proc)) +end + +""" + _groundtruth_is_in_phasespace(in_ps, out_ps) + +Test implementation of the phase space check. Return `false` if either the momentum of the first incoming particle is exactly `zero(SFourMomentum)`, or if the momentum of the last outgoing momentum is exactly `ones(SFourMomentum)`. Otherwise, return true. +""" +function _groundtruth_is_in_phasespace(in_ps, out_ps) + if in_ps[1] == SFourMomentum(zeros(4)) + return false + end + if out_ps[end] == ones(SFourMomentum) + return false + end + return true +end + +""" + _groundtruth_phase_space_factor(in_ps, out_ps) + +Test implementation of the phase space factor. Return the inverse of the product of the energies of all external particles. +""" +function _groundtruth_phase_space_factor(in_ps, out_ps) + en_in = getE.(in_ps) + en_out = getE.(out_ps) + return 1 / (prod(en_in) * prod(en_out)) +end + +function _groundtruth_generate_momenta(ps_coords) + moms = _furl_moms(ps_coords) + return moms +end +""" + _groundtruth_unsafe_probability(proc, in_ps, out_ps) + +Test implementation of the unsafe differential probability. Uses the test implementations of `_groundtruth_matrix_element`,`_groundtruth_averaging_norm` and `_groundtruth_phase_space_factor`. +""" +function _groundtruth_unsafe_probability(proc, in_ps, out_ps) + mat_el = _groundtruth_matrix_element(in_ps, out_ps) + mat_el_sq = abs2.(mat_el) + normalization = _groundtruth_averaging_norm(proc) + ps_fac = _groundtruth_phase_space_factor(in_ps, out_ps) + return sum(mat_el_sq) * ps_fac * normalization +end + +function _groundtruth_unsafe_probability( + proc, in_ps::AbstractVector, out_ps::AbstractMatrix +) + res = Vector{Float64}(undef, size(out_ps, 2)) + for i in 1:size(out_ps, 2) + res[i] = _groundtruth_unsafe_probability(proc, in_ps, view(out_ps, :, i)) + end + return res +end + +function _groundtruth_unsafe_probability( + proc, in_ps::AbstractMatrix, out_ps::AbstractVector +) + res = Vector{Float64}(undef, size(in_ps, 2)) + for i in 1:size(in_ps, 2) + res[i] = _groundtruth_unsafe_probability(proc, view(in_ps, :, i), out_ps) + end + return res +end + +function _groundtruth_unsafe_probability( + proc, in_ps::AbstractMatrix, out_ps::AbstractMatrix +) + res = Matrix{Float64}(undef, size(in_ps, 2), size(out_ps, 2)) + for i in 1:size(in_ps, 2) + for j in 1:size(out_ps, 2) + res[i, j] = _groundtruth_unsafe_probability( + proc, view(in_ps, :, i), view(out_ps, :, j) + ) + end + end + return res +end + +""" + _groundtruth_safe_probability(proc, in_ps, out_ps) + +Test implementation of the safe differential probability. Uses the test implementations of `_groundtruth_is_in_phasespace` and `_groundtruth_unsafe_probability`. +""" +function _groundtruth_safe_probability(proc, in_ps, out_ps) + if !_groundtruth_is_in_phasespace(in_ps, out_ps) + return zero(Float64) + end + return _groundtruth_unsafe_probability(proc, in_ps, out_ps) +end + +function _groundtruth_safe_probability(proc, in_ps::AbstractVector, out_ps::AbstractMatrix) + res = Vector{Float64}(undef, size(out_ps, 2)) + for i in 1:size(out_ps, 2) + res[i] = _groundtruth_safe_probability(proc, in_ps, view(out_ps, :, i)) + end + return res +end + +function _groundtruth_safe_probability(proc, in_ps::AbstractMatrix, out_ps::AbstractVector) + res = Vector{Float64}(undef, size(in_ps, 2)) + for i in 1:size(in_ps, 2) + res[i] = _groundtruth_safe_probability(proc, view(in_ps, :, i), out_ps) + end + return res +end + +function _groundtruth_safe_probability(proc, in_ps::AbstractMatrix, out_ps::AbstractMatrix) + res = Matrix{Float64}(undef, size(in_ps, 2), size(out_ps, 2)) + for i in 1:size(in_ps, 2) + for j in 1:size(out_ps, 2) + res[i, j] = _groundtruth_safe_probability( + proc, view(in_ps, :, i), view(out_ps, :, j) + ) + end + end + return res +end + +""" + _groundtruth_unsafe_diffCS(proc, in_ps, out_ps) + +Test implementation of the unsafe differential cross section. Uses the test implementations of `_groundtruth_incident_flux` and `_groundtruth_unsafe_probability`. +""" +function _groundtruth_unsafe_diffCS(proc, in_ps, out_ps) + init_flux = _groundtruth_incident_flux(in_ps) + return _groundtruth_unsafe_probability(proc, in_ps, out_ps) / (4 * init_flux) +end + +function _groundtruth_unsafe_diffCS(proc, in_ps::AbstractVector, out_ps::AbstractMatrix) + res = Vector{Float64}(undef, size(out_ps, 2)) + for i in 1:size(out_ps, 2) + res[i] = _groundtruth_unsafe_diffCS(proc, in_ps, view(out_ps, :, i)) + end + return res +end + +function _groundtruth_unsafe_diffCS(proc, in_ps::AbstractMatrix, out_ps::AbstractVector) + res = Vector{Float64}(undef, size(in_ps, 2)) + for i in 1:size(in_ps, 2) + res[i] = _groundtruth_unsafe_diffCS(proc, view(in_ps, :, i), out_ps) + end + return res +end + +function _groundtruth_unsafe_diffCS(proc, in_ps::AbstractMatrix, out_ps::AbstractMatrix) + res = Matrix{Float64}(undef, size(in_ps, 2), size(out_ps, 2)) + for i in 1:size(in_ps, 2) + for j in 1:size(out_ps, 2) + res[i, j] = _groundtruth_unsafe_diffCS( + proc, view(in_ps, :, i), view(out_ps, :, j) + ) + end + end + return res +end + +""" + _groundtruth_safe_diffCS(proc, in_ps, out_ps) + +Test implementation of the safe differential cross section. Uses the test implementations of `_groundtruth_is_in_phasespace` and `_groundtruth_unsafe_diffCS`. +""" +function _groundtruth_safe_diffCS(proc, in_ps, out_ps) + if !_groundtruth_is_in_phasespace(in_ps, out_ps) + return zero(Float64) + end + return _groundtruth_unsafe_diffCS(proc, in_ps, out_ps) +end + +function _groundtruth_safe_diffCS(proc, in_ps::AbstractVector, out_ps::AbstractMatrix) + res = Vector{Float64}(undef, size(out_ps, 2)) + for i in 1:size(out_ps, 2) + res[i] = _groundtruth_safe_diffCS(proc, in_ps, view(out_ps, :, i)) + end + return res +end + +function _groundtruth_safe_diffCS(proc, in_ps::AbstractMatrix, out_ps::AbstractVector) + res = Vector{Float64}(undef, size(in_ps, 2)) + for i in 1:size(in_ps, 2) + res[i] = _groundtruth_safe_diffCS(proc, view(in_ps, :, i), out_ps) + end + return res +end + +function _groundtruth_safe_diffCS(proc, in_ps::AbstractMatrix, out_ps::AbstractMatrix) + res = Matrix{Float64}(undef, size(in_ps, 2), size(out_ps, 2)) + for i in 1:size(in_ps, 2) + for j in 1:size(out_ps, 2) + res[i, j] = _groundtruth_safe_diffCS( + proc, view(in_ps, :, i), view(out_ps, :, j) + ) + end + end + return res +end + +""" + _groundtruth_total_probability(in_ps::AbstractVector) + +Test implementation of the total cross section. Return the Minkowski square of the sum the momenta of all incoming particles. +""" +function _groundtruth_total_probability( + in_ps::NTuple{N,T} +) where {N,T<:AbstractFourMomentum} + Ptot = sum(in_ps) + return Ptot * Ptot +end + +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::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_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 new file mode 100644 index 0000000..b631551 --- /dev/null +++ b/test/test_implementation/random_momenta.jl @@ -0,0 +1,83 @@ + +""" +Return a tuple of random four momenta, i.e. a random phase space point. +""" +function _rand_momenta(rng::AbstractRNG, n) + return NTuple{n,SFourMomentum}(SFourMomentum(rand(rng, 4)) for _ in 1:n) +end + +""" +Return a vector of tuples of random four momenta, i.e. a collection of phase space points. +n1 is the size of the phase space point, n2 is the number of points. +""" +function _rand_momenta(rng::AbstractRNG, n1, n2) + moms = Vector{NTuple{n1,SFourMomentum}}(undef, n2) + for i in 1:n2 + moms[i] = _rand_momenta(rng, n1) + end + return moms +end + +""" +Return a random phase space point that is failing the incoming phase space constraint, +i.e. the first entry of the phase space is the null momentum. +""" +function _rand_in_momenta_failing(rng::AbstractRNG, n) + return (zero(SFourMomentum), _rand_momenta(rng, n - 1)...) +end + +""" +Return a random phase space point that is failing the outgoing phase space constraint, +i.e. the last entry of the phase space is the unit momentum. +""" +function _rand_out_momenta_failing(rng::AbstractRNG, n) + return (_rand_momenta(rng, n - 1)..., ones(SFourMomentum)) +end + +""" +Return a collection of incoming phase space points, where the first point is failing the phase space constraint, +i.e. the first entry of the vector is invalid but the others pass. +n1 is the size of the phase space point, n2 is the number of points. +""" +function _rand_in_momenta_failing_mix(rng::AbstractRNG, n1, n2) + moms = _rand_momenta(rng, n1, n2) + moms[1] = _rand_in_momenta_failing(rng, n1) + return moms +end + +""" +Return a collection of incoming phase space points, where all points are failing the phase space constraint, +i.e. their first entries are null momenta. +n1 is the size of the phase space point, n2 is the number of points. +""" +function _rand_in_momenta_failing_all(rng::AbstractRNG, n1, n2) + moms = Vector{NTuple{n1,SFourMomentum}}(undef, n2) + for i in 1:n2 + moms[i] = _rand_in_momenta_failing(rng, n1) + end + return moms +end + +""" +Return a vector of outgoing phase space points, where the first point is failing the phase space constraint, +i.e. the last entry of the vector is invalid but the others pass. +n1 is the size of the phase space point, n2 is the number of points. +""" +function _rand_out_momenta_failing_mix(rng::AbstractRNG, n1, n2) + moms = _rand_momenta(rng, n1, n2) + moms[end] = _rand_out_momenta_failing(rng, n1) + return moms +end + +""" +Return a vector of outgoing phase space points, where all points are failing the phase space constraint, +i.e. their last entries are unit momenta. +n1 is the size of the phase space point, n2 is the number of points. +""" +function _rand_out_momenta_failing_all(rng::AbstractRNG, n1, n2) + moms = Vector{NTuple{n1,SFourMomentum}}(undef, n2) + for i in 1:n2 + moms[i] = _rand_out_momenta_failing(rng, n1) + end + return moms +end diff --git a/test/test_implementation/test_model.jl b/test/test_implementation/test_model.jl new file mode 100644 index 0000000..12b961f --- /dev/null +++ b/test/test_implementation/test_model.jl @@ -0,0 +1,5 @@ + +struct TestModel <: AbstractModelDefinition end +QEDbase.fundamental_interaction_type(::TestModel) = :test_interaction + +struct TestModel_FAIL <: AbstractModelDefinition end diff --git a/test/test_implementation/test_process.jl b/test/test_implementation/test_process.jl new file mode 100644 index 0000000..5a81146 --- /dev/null +++ b/test/test_implementation/test_process.jl @@ -0,0 +1,131 @@ +# dummy particles +struct TestParticleFermion <: FermionLike end +struct TestParticleBoson <: BosonLike end + +const PARTICLE_SET = [TestParticleFermion(), TestParticleBoson()] + +""" + TestProcess(rng,incoming_particles,outgoing_particles) + +""" +struct TestProcess{IP<:Tuple,OP<:Tuple} <: AbstractProcessDefinition + incoming_particles::IP + outgoing_particles::OP +end + +function TestProcess(rng::AbstractRNG, N_in::Int, N_out::Int) + in_particles = rand(rng, PARTICLE_SET, N_in) + out_particles = rand(rng, PARTICLE_SET, N_out) + return TestProcess(in_particles, out_particles) +end + +QEDbase.incoming_particles(proc::TestProcess) = proc.incoming_particles +QEDbase.outgoing_particles(proc::TestProcess) = proc.outgoing_particles + +struct TestProcess_FAIL{IP<:Tuple,OP<:Tuple} <: AbstractProcessDefinition + incoming_particles::IP + outgoing_particles::OP +end + +function TestProcess_FAIL(rng::AbstractRNG, N_in::Int, N_out::Int) + in_particles = Tuple(rand(rng, PARTICLE_SET, N_in)) + out_particles = Tuple(rand(rng, PARTICLE_SET, N_out)) + return TestProcess_FAIL(in_particles, out_particles) +end + +function QEDbase.in_phase_space_dimension(proc::TestProcess, ::TestModel) + return number_incoming_particles(proc) * 4 +end +function QEDbase.out_phase_space_dimension(proc::TestProcess, ::TestModel) + return number_outgoing_particles(proc) * 4 +end + +""" +Test process with no implemented interface. Should fail every usage except construction. +""" +struct TestProcess_FAIL_ALL{IP<:Tuple,OP<:Tuple} <: AbstractProcessDefinition + incoming_particles::IP + outgoing_particles::OP +end + +function TestProcess_FAIL_ALL(rng::AbstractRNG, N_in::Int, N_out::Int) + in_particles = Tuple(rand(rng, PARTICLE_SET, N_in)) + out_particles = Tuple(rand(rng, PARTICLE_SET, N_out)) + return TestProcess_FAIL_ALL(in_particles, out_particles) +end + +""" +Test process with no implemented interface except the incoming and outgoing particles. +Should fail every usage except construction of itself and the respective phase space point for given four-momenta. +""" +struct TestProcess_FAIL_DIFFCS{IP<:Tuple,OP<:Tuple} <: AbstractProcessDefinition + incoming_particles::IP + outgoing_particles::OP +end + +function TestProcess_FAIL_DIFFCS(rng::AbstractRNG, N_in::Int, N_out::Int) + in_particles = Tuple(rand(rng, PARTICLE_SET, N_in)) + out_particles = Tuple(rand(rng, PARTICLE_SET, N_out)) + return TestProcess_FAIL_DIFFCS(in_particles, out_particles) +end + +QEDbase.incoming_particles(proc::TestProcess_FAIL_DIFFCS) = proc.incoming_particles +QEDbase.outgoing_particles(proc::TestProcess_FAIL_DIFFCS) = proc.outgoing_particles + +# dummy phase space definition + failing phase space definition +struct TestPhasespaceDef <: AbstractPhasespaceDefinition end +struct TestPhasespaceDef_FAIL <: AbstractPhasespaceDefinition end + +# dummy implementation of the process interface + +function QEDbase._incident_flux(in_psp::InPhaseSpacePoint{<:TestProcess,<:TestModel}) + return _groundtruth_incident_flux(momenta(in_psp, Incoming())) +end + +function QEDbase._averaging_norm(proc::TestProcess) + return _groundtruth_averaging_norm(proc) +end + +function QEDbase._matrix_element(psp::PhaseSpacePoint{<:TestProcess,TestModel}) + in_ps = momenta(psp, Incoming()) + out_ps = momenta(psp, Outgoing()) + return _groundtruth_matrix_element(in_ps, out_ps) +end + +function QEDbase._is_in_phasespace(psp::PhaseSpacePoint{<:TestProcess,TestModel}) + in_ps = momenta(psp, Incoming()) + out_ps = momenta(psp, Outgoing()) + return _groundtruth_is_in_phasespace(in_ps, out_ps) +end + +function QEDbase._phase_space_factor( + psp::PhaseSpacePoint{<:TestProcess,TestModel,TestPhasespaceDef} +) + in_ps = momenta(psp, Incoming()) + out_ps = momenta(psp, Outgoing()) + return _groundtruth_phase_space_factor(in_ps, out_ps) +end + +function QEDbase._generate_incoming_momenta( + proc::TestProcess, + model::TestModel, + phase_space_def::TestPhasespaceDef, + in_phase_space::NTuple{N,T}, +) where {N,T<:Real} + return _groundtruth_generate_momenta(in_phase_space) +end + +function QEDbase._generate_outgoing_momenta( + proc::TestProcess, + model::TestModel, + phase_space_def::TestPhasespaceDef, + out_phase_space::NTuple{N,T}, +) where {N,T<:Real} + return _groundtruth_generate_momenta(out_phase_space) +end + +function QEDbase._total_probability( + in_psp::InPhaseSpacePoint{<:TestProcess,<:TestModel,<:TestPhasespaceDef} +) + return _groundtruth_total_probability(momenta(in_psp, Incoming())) +end diff --git a/test/test_implementation/utils.jl b/test/test_implementation/utils.jl new file mode 100644 index 0000000..3a141d9 --- /dev/null +++ b/test/test_implementation/utils.jl @@ -0,0 +1,44 @@ + +# Check if any failed type is in the input +_any_fail(x...) = true +_any_fail(::TestProcess, ::TestModel) = false +_any_fail(::TestProcess, ::TestModel, ::TestPhasespaceDef) = false + +# unrolls all elements of a list of four-momenta into vector of coordinates +function _unroll_moms(ps_moms::AbstractVector{T}) where {T<:QEDbase.AbstractFourMomentum} + return collect(Iterators.flatten(ps_moms)) +end + +function _unroll_moms(ps_moms::AbstractMatrix{T}) where {T<:QEDbase.AbstractFourMomentum} + res = Matrix{eltype(T)}(undef, size(ps_moms, 1) * 4, size(ps_moms, 2)) + for i in 1:size(ps_moms, 2) + res[:, i] .= _unroll_moms(view(ps_moms, :, i)) + end + return res +end + +flat_components(moms::AbstractVecOrMat) = _unroll_moms(moms) +flat_components(moms::Tuple) = Tuple(_unroll_moms([moms...])) + +# collect components of four-momenta from a vector of coordinates +function __furl_moms(ps_coords::AbstractVector{T}) where {T<:Real} + return SFourMomentum.(eachcol(reshape(ps_coords, 4, :))) +end + +function _furl_moms(ps_coords::AbstractVector{T}) where {T<:Real} + @assert length(ps_coords) % 4 == 0 + return __furl_moms(ps_coords) +end + +function _furl_moms(ps_coords::AbstractMatrix{T}) where {T<:Real} + @assert size(ps_coords, 1) % 4 == 0 + res = Matrix{SFourMomentum}(undef, Int(size(ps_coords, 1)//4), size(ps_coords, 2)) + for i in 1:size(ps_coords, 2) + res[:, i] .= __furl_moms(view(ps_coords, :, i)) + end + return res +end + +function _furl_moms(moms::NTuple{N,Float64}) where {N} + return Tuple(_furl_moms(Vector{Float64}([moms...]))) +end From 52b80696332e7f07d9b8c59d841d225d6d23ae25 Mon Sep 17 00:00:00 2001 From: Anton Reinhard Date: Thu, 6 Jun 2024 16:36:49 +0200 Subject: [PATCH 06/21] Remove gamma() function --- src/QEDbase.jl | 1 - src/interfaces/gamma_matrices.jl | 4 ---- src/interfaces/lorentz.jl | 2 +- test/interfaces/model.jl | 11 +++++++---- test/test_implementation/test_model.jl | 1 - 5 files changed, 8 insertions(+), 11 deletions(-) diff --git a/src/QEDbase.jl b/src/QEDbase.jl index c12efc0..f81a694 100644 --- a/src/QEDbase.jl +++ b/src/QEDbase.jl @@ -1,6 +1,5 @@ module QEDbase -import Base: * import StaticArrays: similar_type export minkowski_dot, mdot, register_LorentzVectorLike diff --git a/src/interfaces/gamma_matrices.jl b/src/interfaces/gamma_matrices.jl index 1505710..2bdc5e5 100644 --- a/src/interfaces/gamma_matrices.jl +++ b/src/interfaces/gamma_matrices.jl @@ -3,7 +3,3 @@ #### abstract type AbstractGammaRepresentation end - -function gamma(::Type{T})::SLorentzVector where {T<:AbstractGammaRepresentation} - return SLorentzVector(_gamma0(T), _gamma1(T), _gamma2(T), _gamma3(T)) -end diff --git a/src/interfaces/lorentz.jl b/src/interfaces/lorentz.jl index 7456a93..6495621 100644 --- a/src/interfaces/lorentz.jl +++ b/src/interfaces/lorentz.jl @@ -38,7 +38,7 @@ With this done, we can aleady register the `CustomType` as a `LorentzVectorLike` ```julia register_LorentzVectorLike(CustomType) ``` -If something goes wrong, this function call will raise an `RegistryError` indicating, what is missing. With this done, `CustomType` is ready to be used as a LorentzVectorLike +If something goes wrong, this function call will raise an `RegistryError` indicating what is missing. With this done, `CustomType` is ready to be used as a `LorentzVectorLike` ```julia L = CustomType(2.0,1.0,0.0,-1.0) diff --git a/test/interfaces/model.jl b/test/interfaces/model.jl index cdb8887..42ec606 100644 --- a/test/interfaces/model.jl +++ b/test/interfaces/model.jl @@ -1,19 +1,22 @@ using QEDbase -include("../test_implementation/TestImplementation.jl") +struct TestModel <: AbstractModelDefinition end +QEDbase.fundamental_interaction_type(::TestModel) = :test_interaction + +struct TestModel_FAIL <: AbstractModelDefinition end @testset "hard interface" begin - TESTMODEL = TestImplementation.TestModel() + TESTMODEL = TestModel() @test fundamental_interaction_type(TESTMODEL) == :test_interaction end @testset "interface fail" begin - TESTMODEL_FAIL = TestImplementation.TestModel_FAIL() + TESTMODEL_FAIL = TestModel_FAIL() @test_throws MethodError fundamental_interaction_type(TESTMODEL_FAIL) end @testset "broadcast" begin test_func(model) = model - TESTMODEL = TestImplementation.TestModel() + TESTMODEL = TestModel() @test test_func.(TESTMODEL) == TESTMODEL end diff --git a/test/test_implementation/test_model.jl b/test/test_implementation/test_model.jl index 12b961f..87f4f04 100644 --- a/test/test_implementation/test_model.jl +++ b/test/test_implementation/test_model.jl @@ -1,4 +1,3 @@ - struct TestModel <: AbstractModelDefinition end QEDbase.fundamental_interaction_type(::TestModel) = :test_interaction From 4df433d0b86de0e98447e7aee80b1422c3a9734c Mon Sep 17 00:00:00 2001 From: Anton Reinhard Date: Thu, 6 Jun 2024 17:06:35 +0200 Subject: [PATCH 07/21] Split Lorentz interface into files, split errors off --- src/QEDbase.jl | 14 +- src/errors.jl | 57 ++++++ src/interfaces/four_momentum.jl | 14 ++ .../arithmetic.jl} | 190 ------------------ src/interfaces/lorentz_vectors/fields.jl | 75 +++++++ src/interfaces/lorentz_vectors/registry.jl | 38 ++++ src/interfaces/lorentz_vectors/types.jl | 49 +++++ src/interfaces/lorentz_vectors/utility.jl | 50 +++++ src/interfaces/setup.jl | 24 --- 9 files changed, 295 insertions(+), 216 deletions(-) create mode 100644 src/errors.jl rename src/interfaces/{lorentz.jl => lorentz_vectors/arithmetic.jl} (82%) create mode 100644 src/interfaces/lorentz_vectors/fields.jl create mode 100644 src/interfaces/lorentz_vectors/registry.jl create mode 100644 src/interfaces/lorentz_vectors/types.jl create mode 100644 src/interfaces/lorentz_vectors/utility.jl diff --git a/src/QEDbase.jl b/src/QEDbase.jl index f81a694..8eb765c 100644 --- a/src/QEDbase.jl +++ b/src/QEDbase.jl @@ -72,12 +72,15 @@ export number_incoming_particles, number_outgoing_particles export particles, number_particles # Abstract setup interface -export AbstractComputationSetup, InvalidInputError, compute +export AbstractComputationSetup, compute export AbstractProcessSetup, scattering_process, physical_model # Abstract phase space interface export AbstractCoordinateSystem, AbstractFrameOfReference, AbstractPhasespaceDefinition +# errors +export InvalidInputError, RegistryError, OnshellError + using StaticArrays using LinearAlgebra using DocStringExtensions @@ -86,11 +89,18 @@ using SimpleTraits using ArgCheck using ConstructionBase +include("errors.jl") include("utils.jl") include("interfaces/dirac_tensors.jl") include("interfaces/gamma_matrices.jl") -include("interfaces/lorentz.jl") + +include("interfaces/lorentz_vectors/types.jl") +include("interfaces/lorentz_vectors/registry.jl") +include("interfaces/lorentz_vectors/arithmetic.jl") +include("interfaces/lorentz_vectors/fields.jl") +include("interfaces/lorentz_vectors/utility.jl") + include("interfaces/four_momentum.jl") include("interfaces/model.jl") diff --git a/src/errors.jl b/src/errors.jl new file mode 100644 index 0000000..d0e2442 --- /dev/null +++ b/src/errors.jl @@ -0,0 +1,57 @@ +struct OnshellError{M,T} <: Exception + mom::M + mass::T +end + +function Base.showerror(io::IO, e::OnshellError) + return print( + io, + "OnshellError: The momentum $(e.mom) is not onshell w.r.t. the mass $(e.mass).\n mom*mom = $(e.mom*e.mom)", + ) +end + +""" +$(TYPEDEF) + +Exception raised, if a certain type `target_type` can not be registed for a certain interface, since there needs the function `func` to be impleemnted. + +# Fields +$(TYPEDFIELDS) + +""" +struct RegistryError <: Exception + func::Function + target_type::Any +end + +function Base.showerror(io::IO, err::RegistryError) + return println( + io, + "RegistryError:", + " The type <$(err.target_type)> must implement <$(err.func)> to be registered.", + ) +end + +""" +Abstract base type for exceptions indicating invalid input. See [`InvalidInputError`](@ref) for a simple concrete implementation. +Concrete implementations should at least implement + +```Julia + +Base.showerror(io::IO, err::CustomInvalidError) where {CustomInvalidError<:AbstractInvalidInputException} + +``` +""" +abstract type AbstractInvalidInputException <: Exception end + +""" + InvalidInputError(msg::String) + +Exception which is thrown if a given input is invalid, e.g. passed to [`_assert_valid_input`](@ref). +""" +struct InvalidInputError <: AbstractInvalidInputException + msg::String +end +function Base.showerror(io::IO, err::InvalidInputError) + return println(io, "InvalidInputError: $(err.msg)") +end diff --git a/src/interfaces/four_momentum.jl b/src/interfaces/four_momentum.jl index ed05534..733baac 100644 --- a/src/interfaces/four_momentum.jl +++ b/src/interfaces/four_momentum.jl @@ -7,3 +7,17 @@ Also see: [`SFourMomentum`](@ref) """ abstract type AbstractFourMomentum <: AbstractLorentzVector{Float64} end + +function Base.getproperty(P::TM, sym::Symbol) where {TM<:AbstractFourMomentum} + if sym == :t + return P.E + elseif sym == :x + return P.px + elseif sym == :y + return P.py + elseif sym == :z + return P.pz + else + return getfield(P, sym) + end +end diff --git a/src/interfaces/lorentz.jl b/src/interfaces/lorentz_vectors/arithmetic.jl similarity index 82% rename from src/interfaces/lorentz.jl rename to src/interfaces/lorentz_vectors/arithmetic.jl index 6495621..cd7b7f8 100644 --- a/src/interfaces/lorentz.jl +++ b/src/interfaces/lorentz_vectors/arithmetic.jl @@ -1,190 +1,3 @@ -####### -# -# Abstract types -# -####### -""" -$(TYPEDEF) - -Abstract type to model generic Lorentz vectors, i.e. elements of a minkowski-like space, where the component space is arbitray. -""" -abstract type AbstractLorentzVector{T} <: FieldVector{4,T} end - -""" -## Definition of LorentzVector interface. - -It enables several functions related to Lorentz vectors for any custom type just by implementing the API for the respective type. -For instance, say we want to implement a type, which shall act like a Lorentz vector. Then, we start with our custom type: - -```julia -struct CustomType - t::Float64 - x::Float64 - y::Float64 - z::Float64 -end -``` -The first group of functions to be implemented for this `CustomType` in order to connect this type to the Lorentz vector interface are the getter for the cartesian components. - -```julia -QEDbase.getT(lv::CustomType) = lv.t -QEDbase.getX(lv::CustomType) = lv.x -QEDbase.getY(lv::CustomType) = lv.y -QEDbase.getZ(lv::CustomType) = lv.z -``` -Make sure that you dispatch on the getter from `QEDbase` by defining `QEDbase.getT`, etc. - -With this done, we can aleady register the `CustomType` as a `LorentzVectorLike` type using the `register_LorentzVectorLike` function -```julia -register_LorentzVectorLike(CustomType) -``` -If something goes wrong, this function call will raise an `RegistryError` indicating what is missing. With this done, `CustomType` is ready to be used as a `LorentzVectorLike` -```julia -L = CustomType(2.0,1.0,0.0,-1.0) - -getTheta(L) # -getRapidity(L) # -``` - -""" - -""" - -$(TYPEDEF) - -Exception raised, if a certain type `target_type` can not be registed for a certain interface, since there needs the function `func` to be impleemnted. - -# Fields -$(TYPEDFIELDS) - -""" -struct RegistryError <: Exception - func::Function - target_type::Any -end - -function Base.showerror(io::IO, err::RegistryError) - return println( - io, - "RegistryError:", - " The type <$(err.target_type)> must implement <$(err.func)> to be registered.", - ) -end - -""" -$(SIGNATURES) - -Wrapper around `Base.hasmethod` with a more meaningful error message in the context of function registration. -""" -function _hasmethod_registry(fun::Function, ::Type{T}) where {T} - @argcheck hasmethod(fun, Tuple{T}) RegistryError(fun, T) -end - -@traitdef IsLorentzVectorLike{T} -@traitdef IsMutableLorentzVectorLike{T} - -""" - getT(lv) - -Return the 0-component of a given `LorentzVectorLike`. - -!!! example - - If `(t,x,y,z)` is a `LorentzVectorLike`, this is equivalent to `t`. - -""" -function getT end - -""" - getX(lv) - -Return the 1-component of a given `LorentzVectorLike`. - -!!! example - - If `(t,x,y,z)` is a `LorentzVectorLike`, this is equivalent to `x`. - -""" -function getX end - -""" - getY(lv) - -Return the 2-component of a given `LorentzVectorLike`. - -!!! example - - If `(t,x,y,z)` is a `LorentzVectorLike`, this is equivalent to `y`. - -""" -function getY end - -""" - getZ(lv) - -Return the 3-component of a given `LorentzVectorLike`. - -!!! example - - If `(t,x,y,z)` is a `LorentzVectorLike`, this is equivalent to `z`. - -""" -function getZ end - -""" - setT!(lv,value) - -Sets the 0-component of a given `LorentzVectorLike` to the given `value`. -""" -function setT! end - -""" - setX!(lv,value) - -Sets the 1-component of a given `LorentzVectorLike` to the given `value`. -""" -function setX! end - -""" - setY!(lv,value) - -Sets the 2-component of a given `LorentzVectorLike` to the given `value`. -""" -function setY! end - -""" - setZ!(lv,value) - -Sets the 3-component of a given `LorentzVectorLike` to the given `value`. -""" -function setZ! end - -""" -$(SIGNATURES) - -Function to register a custom type as a LorentzVectorLike. - -Ensure the passed custom type has implemented at least the function `getT, getX, getY, getZ` -and enables getter functions of the lorentz vector library for the given type. -If additionally the functions `setT!, setX!, setY!, setZ!` are implemened for the passed custom type, -also the setter functions of the Lorentz vector interface are enabled. -""" -function register_LorentzVectorLike(T::Type) - _hasmethod_registry(getT, T) - _hasmethod_registry(getX, T) - _hasmethod_registry(getY, T) - _hasmethod_registry(getZ, T) - - @eval @traitimpl IsLorentzVectorLike{$T} - - if hasmethod(setT!, Tuple{T,<:Union{}}) && - hasmethod(setX!, Tuple{T,<:Union{}}) && - hasmethod(setY!, Tuple{T,<:Union{}}) && - hasmethod(setZ!, Tuple{T,<:Union{}}) - @eval @traitimpl IsMutableLorentzVectorLike{$T} - end - return nothing -end """ minkowski_dot(v1,v2) @@ -215,9 +28,6 @@ Function alias for [`minkowski_dot`](@ref). """ const mdot = minkowski_dot -######## -# getter -######## """ getMagnitude2(lv) diff --git a/src/interfaces/lorentz_vectors/fields.jl b/src/interfaces/lorentz_vectors/fields.jl new file mode 100644 index 0000000..1a7fc9d --- /dev/null +++ b/src/interfaces/lorentz_vectors/fields.jl @@ -0,0 +1,75 @@ +""" + getT(lv) + +Return the 0-component of a given `LorentzVectorLike`. + +!!! example + + If `(t,x,y,z)` is a `LorentzVectorLike`, this is equivalent to `t`. + +""" +function getT end + +""" + getX(lv) + +Return the 1-component of a given `LorentzVectorLike`. + +!!! example + + If `(t,x,y,z)` is a `LorentzVectorLike`, this is equivalent to `x`. + +""" +function getX end + +""" + getY(lv) + +Return the 2-component of a given `LorentzVectorLike`. + +!!! example + + If `(t,x,y,z)` is a `LorentzVectorLike`, this is equivalent to `y`. + +""" +function getY end + +""" + getZ(lv) + +Return the 3-component of a given `LorentzVectorLike`. + +!!! example + + If `(t,x,y,z)` is a `LorentzVectorLike`, this is equivalent to `z`. + +""" +function getZ end + +""" + setT!(lv,value) + +Sets the 0-component of a given `LorentzVectorLike` to the given `value`. +""" +function setT! end + +""" + setX!(lv,value) + +Sets the 1-component of a given `LorentzVectorLike` to the given `value`. +""" +function setX! end + +""" + setY!(lv,value) + +Sets the 2-component of a given `LorentzVectorLike` to the given `value`. +""" +function setY! end + +""" + setZ!(lv,value) + +Sets the 3-component of a given `LorentzVectorLike` to the given `value`. +""" +function setZ! end diff --git a/src/interfaces/lorentz_vectors/registry.jl b/src/interfaces/lorentz_vectors/registry.jl new file mode 100644 index 0000000..643e7da --- /dev/null +++ b/src/interfaces/lorentz_vectors/registry.jl @@ -0,0 +1,38 @@ +""" +$(SIGNATURES) + +Wrapper around `Base.hasmethod` with a more meaningful error message in the context of function registration. +""" +function _hasmethod_registry(fun::Function, ::Type{T}) where {T} + @argcheck hasmethod(fun, Tuple{T}) RegistryError(fun, T) +end + +@traitdef IsLorentzVectorLike{T} +@traitdef IsMutableLorentzVectorLike{T} + +""" +$(SIGNATURES) + +Function to register a custom type as a LorentzVectorLike. + +Ensure the passed custom type has implemented at least the function `getT, getX, getY, getZ` +and enables getter functions of the lorentz vector library for the given type. +If additionally the functions `setT!, setX!, setY!, setZ!` are implemened for the passed custom type, +also the setter functions of the Lorentz vector interface are enabled. +""" +function register_LorentzVectorLike(T::Type) + _hasmethod_registry(getT, T) + _hasmethod_registry(getX, T) + _hasmethod_registry(getY, T) + _hasmethod_registry(getZ, T) + + @eval @traitimpl IsLorentzVectorLike{$T} + + if hasmethod(setT!, Tuple{T,<:Union{}}) && + hasmethod(setX!, Tuple{T,<:Union{}}) && + hasmethod(setY!, Tuple{T,<:Union{}}) && + hasmethod(setZ!, Tuple{T,<:Union{}}) + @eval @traitimpl IsMutableLorentzVectorLike{$T} + end + return nothing +end diff --git a/src/interfaces/lorentz_vectors/types.jl b/src/interfaces/lorentz_vectors/types.jl new file mode 100644 index 0000000..456aa37 --- /dev/null +++ b/src/interfaces/lorentz_vectors/types.jl @@ -0,0 +1,49 @@ +####### +# +# Abstract types +# +####### +""" +$(TYPEDEF) + +Abstract type to model generic Lorentz vectors, i.e. elements of a minkowski-like space, where the component space is arbitray. +""" +abstract type AbstractLorentzVector{T} <: FieldVector{4,T} end + +""" +## Definition of LorentzVector interface. + +It enables several functions related to Lorentz vectors for any custom type just by implementing the API for the respective type. +For instance, say we want to implement a type, which shall act like a Lorentz vector. Then, we start with our custom type: + +```julia +struct CustomType + t::Float64 + x::Float64 + y::Float64 + z::Float64 +end +``` +The first group of functions to be implemented for this `CustomType` in order to connect this type to the Lorentz vector interface are the getter for the cartesian components. + +```julia +QEDbase.getT(lv::CustomType) = lv.t +QEDbase.getX(lv::CustomType) = lv.x +QEDbase.getY(lv::CustomType) = lv.y +QEDbase.getZ(lv::CustomType) = lv.z +``` +Make sure that you dispatch on the getter from `QEDbase` by defining `QEDbase.getT`, etc. + +With this done, we can aleady register the `CustomType` as a `LorentzVectorLike` type using the `register_LorentzVectorLike` function +```julia +register_LorentzVectorLike(CustomType) +``` +If something goes wrong, this function call will raise an `RegistryError` indicating, what is missing. With this done, `CustomType` is ready to be used as a `LorentzVectorLike` +```julia +L = CustomType(2.0,1.0,0.0,-1.0) + +getTheta(L) # +getRapidity(L) # +``` + +""" diff --git a/src/interfaces/lorentz_vectors/utility.jl b/src/interfaces/lorentz_vectors/utility.jl new file mode 100644 index 0000000..bc2a656 --- /dev/null +++ b/src/interfaces/lorentz_vectors/utility.jl @@ -0,0 +1,50 @@ + +####### +# +# Utility functions on FourMomenta +# +####### + +function isonshell(mom::QEDbase.AbstractLorentzVector{T}) where {T<:Real} + mag2 = getMag2(mom) + E = getE(mom) + return isapprox(E^2, mag2; rtol=eps(T)) +end + +""" +$(SIGNATURES) + +On-shell check of a given four-momentum `mom` w.r.t. a given mass `mass`. + +!!! note "Precision" + For `AbstactFourMomentum`, the element type is fixed to `Float64`, limiting the precision of comparisons between elements. + The current implementation has been tested within the boundaries for energy scales `E` with `1e-9 <= E <= 1e5`. + In those bounds, the mass error, which is correctly detected as off-shell, is `1e-4` times the mean value of the components, but has at most the value `0.01`, e.g. at the high energy end. + The energy scales correspond to `0.5 meV` for the lower bound and `50 GeV` for the upper bound. + + +!!! todo "FourMomenta with real entries" + * if `AbstractFourMomentum` is updated to elementtypes `T<:Real`, the `AbstractLorentzVector` should be updated with the `AbstractFourMomentum`. +""" +function isonshell(mom::QEDbase.AbstractLorentzVector{T}, mass::Real) where {T<:Real} + if iszero(mass) + return isonshell(mom) + end + mag2 = getMag2(mom) + E = getE(mom) + return isapprox(E^2, (mass)^2 + mag2; atol=2 * eps(T), rtol=eps(T)) +end + +""" +$(SIGNATURES) + +Assertion if a FourMomentum `mom` is on-shell w.r.t a given mass `mass`. + +!!! note "See also" + The precision of this functions is explained in [`isonshell`](@ref). + +""" +function assert_onshell(mom::QEDbase.AbstractLorentzVector, mass::Real) + isonshell(mom, mass) || throw(OnshellError(mom, mass)) + return nothing +end diff --git a/src/interfaces/setup.jl b/src/interfaces/setup.jl index e495549..7b0f8ce 100644 --- a/src/interfaces/setup.jl +++ b/src/interfaces/setup.jl @@ -60,30 +60,6 @@ abstract type AbstractComputationSetup end # convenience function to check if an object is a computation setup _is_computation_setup(::AbstractComputationSetup) = true -""" -Abstract base type for exceptions indicating invalid input. See [`InvalidInputError`](@ref) for a simple concrete implementation. -Concrete implementations should at least implement - -```Julia - -Base.showerror(io::IO, err::CustomInvalidError) where {CustomInvalidError<:AbstractInvalidInputException} - -``` -""" -abstract type AbstractInvalidInputException <: Exception end - -""" - InvalidInputError(msg::String) - -Exception which is thrown if a given input is invalid, e.g. passed to [`_assert_valid_input`](@ref). -""" -struct InvalidInputError <: AbstractInvalidInputException - msg::String -end -function Base.showerror(io::IO, err::InvalidInputError) - return println(io, "InvalidInputError: $(err.msg)") -end - """ _assert_valid_input(stp::AbstractComputationSetup, input::Any) From af20aad507a41ad6b60b04809d153bc288aecc4a Mon Sep 17 00:00:00 2001 From: Anton Reinhard Date: Thu, 6 Jun 2024 17:22:03 +0200 Subject: [PATCH 08/21] Move functions of abstract types back --- Project.toml | 3 +- src/QEDbase.jl | 2 + src/interfaces/lorentz_vectors/arithmetic.jl | 8 ++ .../lorentz_vectors/dirac_interaction.jl | 38 +++++++ src/interfaces/lorentz_vectors/utility.jl | 1 - test/lorentz_vector.jl | 98 ------------------- 6 files changed, 49 insertions(+), 101 deletions(-) create mode 100644 src/interfaces/lorentz_vectors/dirac_interaction.jl delete mode 100644 test/lorentz_vector.jl diff --git a/Project.toml b/Project.toml index 948200c..078455a 100644 --- a/Project.toml +++ b/Project.toml @@ -16,7 +16,6 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" [compat] ArgCheck = "2.3.0" -ConstructionBase = "1.3.0" DocStringExtensions = "0.8.5, 0.9" PhysicalConstants = "0.2.1" SimpleTraits = "0.9.4" @@ -25,8 +24,8 @@ julia = "1.6" [extras] SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" -Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] test = ["SafeTestsets", "Test", "Suppressor"] diff --git a/src/QEDbase.jl b/src/QEDbase.jl index 8eb765c..5f955d9 100644 --- a/src/QEDbase.jl +++ b/src/QEDbase.jl @@ -1,5 +1,6 @@ module QEDbase +import Base: * import StaticArrays: similar_type export minkowski_dot, mdot, register_LorentzVectorLike @@ -98,6 +99,7 @@ include("interfaces/gamma_matrices.jl") include("interfaces/lorentz_vectors/types.jl") include("interfaces/lorentz_vectors/registry.jl") include("interfaces/lorentz_vectors/arithmetic.jl") +include("interfaces/lorentz_vectors/dirac_interaction.jl") include("interfaces/lorentz_vectors/fields.jl") include("interfaces/lorentz_vectors/utility.jl") diff --git a/src/interfaces/lorentz_vectors/arithmetic.jl b/src/interfaces/lorentz_vectors/arithmetic.jl index cd7b7f8..2247eaa 100644 --- a/src/interfaces/lorentz_vectors/arithmetic.jl +++ b/src/interfaces/lorentz_vectors/arithmetic.jl @@ -1,3 +1,11 @@ +function dot(p1::T1, p2::T2) where {T1<:AbstractLorentzVector,T2<:AbstractLorentzVector} + return mdot(p1, p2) +end +@inline function *( + p1::T1, p2::T2 +) where {T1<:AbstractLorentzVector,T2<:AbstractLorentzVector} + return dot(p1, p2) +end """ minkowski_dot(v1,v2) diff --git a/src/interfaces/lorentz_vectors/dirac_interaction.jl b/src/interfaces/lorentz_vectors/dirac_interaction.jl new file mode 100644 index 0000000..530d8ff --- /dev/null +++ b/src/interfaces/lorentz_vectors/dirac_interaction.jl @@ -0,0 +1,38 @@ +""" +$(TYPEDSIGNATURES) + +Product of generic Lorentz vector with a Dirac tensor from the left. Basically, the multiplication is piped to the components from the Lorentz vector. + +!!! note "Multiplication operator" + This also overloads the `*` operator for this types. +""" +function mul( + DM::T, L::TL +) where {T<:Union{AbstractDiracMatrix,AbstractDiracVector},TL<:AbstractLorentzVector} + return constructorof(TL)(DM * L[1], DM * L[2], DM * L[3], DM * L[4]) +end +@inline function *( + DM::T, L::TL +) where {T<:Union{AbstractDiracMatrix,AbstractDiracVector},TL<:AbstractLorentzVector} + return mul(DM, L) +end + +""" +$(TYPEDSIGNATURES) + +Product of generic Lorentz vector with a Dirac tensor from the right. Basically, the multiplication is piped to the components from the Lorentz vector. + +!!! note "Multiplication operator" + This also overloads the `*` operator for this types. + +""" +function mul( + L::TL, DM::T +) where {TL<:AbstractLorentzVector,T<:Union{AbstractDiracMatrix,AbstractDiracVector}} + return constructorof(TL)(L[1] * DM, L[2] * DM, L[3] * DM, L[4] * DM) +end +@inline function *( + L::TL, DM::T +) where {TL<:AbstractLorentzVector,T<:Union{AbstractDiracMatrix,AbstractDiracVector}} + return mul(L, DM) +end diff --git a/src/interfaces/lorentz_vectors/utility.jl b/src/interfaces/lorentz_vectors/utility.jl index bc2a656..46dd554 100644 --- a/src/interfaces/lorentz_vectors/utility.jl +++ b/src/interfaces/lorentz_vectors/utility.jl @@ -1,4 +1,3 @@ - ####### # # Utility functions on FourMomenta diff --git a/test/lorentz_vector.jl b/test/lorentz_vector.jl deleted file mode 100644 index 3a154f0..0000000 --- a/test/lorentz_vector.jl +++ /dev/null @@ -1,98 +0,0 @@ -using QEDbase -using StaticArrays - -unary_methods = [-, +] - -@testset "LorentzVectorType" for LorentzVectorType in [SLorentzVector, MLorentzVector] - @testset "General Properties" begin - LV = LorentzVectorType(rand(Float64, 4)) - - @test size(LV) == (4,) - @test length(LV) == 4 - @test eltype(LV) == Float64 - - @test @inferred(LorentzVectorType(1, 2, 3, 4)) == LorentzVectorType([1, 2, 3, 4]) - - @test LorentzVectorType(1, 2, 3, 4) == LorentzVectorType(SVector{4}(1, 2, 3, 4)) - - M = rand(Float64, 4, 4) - LV_mat = LorentzVectorType(M, M, M, M) - - @test size(LV_mat) == (4,) - @test length(LV_mat) == 4 - @test eltype(LV_mat) == typeof(M) - - @test_throws DimensionMismatch( - "No precise constructor for $(LorentzVectorType) found. Length of input was 2." - ) LorentzVectorType(1, 2) - - @test LV.t == LV[1] == getT(LV) - @test LV.x == LV[2] == getX(LV) - @test LV.y == LV[3] == getY(LV) - @test LV.z == LV[4] == getZ(LV) - end # General Properties - - @testset "Arithmetics" begin - num = 2.0 - LV1 = LorentzVectorType(1, 2, 3, 4) - LV2 = LorentzVectorType(4, 3, 2, 1) - - @test @inferred(LV1 + LV2) == LorentzVectorType(5, 5, 5, 5) - @test @inferred(LV1 - LV2) == LorentzVectorType(-3, -1, 1, 3) - - @test LV1 * LV2 == -12.0 - - @test @inferred(num * LV1) == LorentzVectorType(2, 4, 6, 8) - @test @inferred(LV1 * num) == LorentzVectorType(2, 4, 6, 8) - @test @inferred(LV1 / num) == LorentzVectorType(0.5, 1.0, 1.5, 2.0) - - num_cmplx = 1.0 + 1.0im - - #maybe use @inferred for type stability check - @test typeof(num_cmplx * LV1) === LorentzVectorType{ComplexF64} # type casting - @test eltype(num_cmplx * LV1) === ComplexF64 # type casting - @test typeof(LV1 * num_cmplx) === LorentzVectorType{ComplexF64} # type casting - @test eltype(LV1 * num_cmplx) === ComplexF64 # type casting - @test typeof(LV1 / num_cmplx) === LorentzVectorType{ComplexF64} # type casting - @test eltype(LV1 / num_cmplx) === ComplexF64 # type casting - - LV = LorentzVectorType(rand(4)) - for ops in unary_methods - res = ops(LV) - @test typeof(res) == typeof(LV) - @test SArray(res) == ops(SArray(LV)) - end - end # Arithmetics - - @testset "interface dirac tensor" begin - DM_eye = one(DiracMatrix) - BS = BiSpinor(1, 2, 3, 4) - aBS = AdjointBiSpinor(1, 2, 3, 4) - - LV_mat = LorentzVectorType(DM_eye, DM_eye, DM_eye, DM_eye) - LV_BS = LorentzVectorType(BS, BS, BS, BS) - LV_aBS = LorentzVectorType(aBS, aBS, aBS, aBS) - - @test @inferred(LV_mat * BS) == LV_BS - @test @inferred(aBS * LV_mat) == LV_aBS - @test @inferred(LV_aBS * LV_BS) == -60.0 + 0.0im - end #interface dirac tensor - - @testset "utility functions" begin - LV = LorentzVectorType(4, 3, 2, 1) - - @test isapprox(@inferred(getMagnitude(LV)), sqrt(14)) - end # utility functions -end # LorentzVectorType - -@testset "different LorentzVectorTypes" begin - @testset "Artihmetics" begin - LV1 = SLorentzVector(1, 2, 3, 4) - LV2 = MLorentzVector(4, 3, 2, 1) - - @test @inferred(LV1 + LV2) === SLorentzVector(5, 5, 5, 5) - @test @inferred(LV1 - LV2) === SLorentzVector(-3, -1, 1, 3) - - @test LV1 * LV2 == -12.0 - end -end From c71bd8acc92b5ff8ea38714f013a667e6433d178 Mon Sep 17 00:00:00 2001 From: Anton Reinhard Date: Thu, 6 Jun 2024 18:11:25 +0200 Subject: [PATCH 09/21] Fix and move tests --- test/dirac_tensor.jl | 125 --------- .../lorentz.jl} | 0 test/interfaces/process.jl | 129 --------- test/interfaces/setup.jl | 188 ------------- test/particle_spinors.jl | 88 ------ test/particles.jl | 113 -------- test/runtests.jl | 22 +- .../test_implementation/TestImplementation.jl | 38 --- test/test_implementation/groundtruths.jl | 265 ------------------ test/test_implementation/random_momenta.jl | 83 ------ test/test_implementation/test_model.jl | 4 - test/test_implementation/test_process.jl | 131 --------- test/test_implementation/utils.jl | 44 --- test/utils.jl | 10 - 14 files changed, 1 insertion(+), 1239 deletions(-) delete mode 100644 test/dirac_tensor.jl rename test/{lorentz_interface.jl => interfaces/lorentz.jl} (100%) delete mode 100644 test/interfaces/process.jl delete mode 100644 test/interfaces/setup.jl delete mode 100644 test/particle_spinors.jl delete mode 100644 test/test_implementation/TestImplementation.jl delete mode 100644 test/test_implementation/groundtruths.jl delete mode 100644 test/test_implementation/random_momenta.jl delete mode 100644 test/test_implementation/test_model.jl delete mode 100644 test/test_implementation/test_process.jl delete mode 100644 test/test_implementation/utils.jl delete mode 100644 test/utils.jl diff --git a/test/dirac_tensor.jl b/test/dirac_tensor.jl deleted file mode 100644 index d153d17..0000000 --- a/test/dirac_tensor.jl +++ /dev/null @@ -1,125 +0,0 @@ -using QEDbase -using StaticArrays - -unary_methods = [-, +] -binary_array_methods = [+, -] -binary_float_methods = [*, /] - -allowed_muls = Dict([ - (AdjointBiSpinor, BiSpinor) => ComplexF64, - (BiSpinor, AdjointBiSpinor) => DiracMatrix, - (AdjointBiSpinor, DiracMatrix) => AdjointBiSpinor, - (DiracMatrix, BiSpinor) => BiSpinor, - (DiracMatrix, DiracMatrix) => DiracMatrix, -]) - -groundtruth_mul(a::AdjointBiSpinor, b::BiSpinor) = transpose(SArray(a)) * SArray(b) -function groundtruth_mul(a::BiSpinor, b::AdjointBiSpinor) - return DiracMatrix(SArray(a) * transpose(SArray(b))) -end -function groundtruth_mul(a::AdjointBiSpinor, b::DiracMatrix) - return AdjointBiSpinor(transpose(SArray(a)) * SArray(b)) -end -groundtruth_mul(a::DiracMatrix, b::BiSpinor) = BiSpinor(SArray(a) * SArray(b)) -groundtruth_mul(a::DiracMatrix, b::DiracMatrix) = DiracMatrix(SArray(a) * SArray(b)) - -@testset "DiracTensor" begin - dirac_tensors = [ - BiSpinor(rand(ComplexF64, 4)), - AdjointBiSpinor(rand(ComplexF64, 4)), - DiracMatrix(rand(ComplexF64, 4, 4)), - ] - - @testset "BiSpinor" begin - BS = BiSpinor(rand(4)) - - @test size(BS) == (4,) - @test length(BS) == 4 - @test eltype(BS) == ComplexF64 - - @test @inferred(BiSpinor(1, 2, 3, 4)) == @inferred(BiSpinor([1, 2, 3, 4])) - - BS1 = BiSpinor(1, 2, 3, 4) - BS2 = BiSpinor(4, 3, 2, 1) - - @test @inferred(BS1 + BS2) == BiSpinor(5, 5, 5, 5) - @test @inferred(BS1 - BS2) == BiSpinor(-3, -1, 1, 3) - - @test_throws DimensionMismatch( - "No precise constructor for BiSpinor found. Length of input was 2." - ) BiSpinor(1, 2) - end #BiSpinor - - @testset "AdjointBiSpinor" begin - aBS = AdjointBiSpinor(rand(4)) - - @test size(aBS) == (4,) - @test length(aBS) == 4 - @test eltype(aBS) == ComplexF64 - - @test @inferred(AdjointBiSpinor(1, 2, 3, 4)) == - @inferred(AdjointBiSpinor([1, 2, 3, 4])) - - aBS1 = AdjointBiSpinor(1, 2, 3, 4) - aBS2 = AdjointBiSpinor(4, 3, 2, 1) - - @test @inferred(aBS1 + aBS2) == AdjointBiSpinor(5, 5, 5, 5) - @test @inferred(aBS1 - aBS2) == AdjointBiSpinor(-3, -1, 1, 3) - - @test_throws DimensionMismatch( - "No precise constructor for AdjointBiSpinor found. Length of input was 2." - ) AdjointBiSpinor(1, 2) - end #AdjointBiSpinor - - @testset "DiracMatrix" begin - DM = DiracMatrix(rand(4, 4)) - - @test size(DM) == (4, 4) - @test length(DM) == 16 - @test eltype(DM) == ComplexF64 - - DM1 = DiracMatrix(SDiagonal(1, 2, 3, 4)) - DM2 = DiracMatrix(SDiagonal(4, 3, 2, 1)) - - @test @inferred(DM1 + DM2) == DiracMatrix(SDiagonal(5, 5, 5, 5)) - @test @inferred(DM1 - DM2) == DiracMatrix(SDiagonal(-3, -1, 1, 3)) - - @test_throws DimensionMismatch( - "No precise constructor for DiracMatrix found. Length of input was 2." - ) DiracMatrix(1, 2) - end #DiracMatrix - - @testset "General Arithmetics" begin - num = rand(ComplexF64) - - for ten in dirac_tensors - @testset "$ops($(typeof(ten)))" for ops in unary_methods - res = ops(ten) - @test typeof(res) == typeof(ten) - @test SArray(res) == ops(SArray(ten)) - end - - @testset "num*$(typeof(ten))" begin - #@test typeof(res_float_mul) == typeof(ten) - @test @inferred(ten * num) == @inferred(num * ten) - @test SArray(num * ten) == num * SArray(ten) - end - - @testset "$(typeof(ten))/num" begin - res_float_div = ten / num - #@test typeof(res_float_div) == typeof(ten) - @test SArray(@inferred(ten / num)) == SArray(ten) / num - end - - @testset "$(typeof(ten))*$(typeof(ten2))" for ten2 in dirac_tensors - mul_comb = (typeof(ten), typeof(ten2)) - if mul_comb in keys(allowed_muls) - res = ten * ten2 - @test typeof(res) == allowed_muls[mul_comb] - @test isapprox(res, groundtruth_mul(ten, ten2)) - #issue: test raise of method error - end - end - end - end #Arithmetics -end #"DiracTensor" diff --git a/test/lorentz_interface.jl b/test/interfaces/lorentz.jl similarity index 100% rename from test/lorentz_interface.jl rename to test/interfaces/lorentz.jl diff --git a/test/interfaces/process.jl b/test/interfaces/process.jl deleted file mode 100644 index 3ecc139..0000000 --- a/test/interfaces/process.jl +++ /dev/null @@ -1,129 +0,0 @@ -using Random -using QEDbase - -RNG = MersenneTwister(137137) -ATOL = 0.0 -RTOL = sqrt(eps()) - -include("../test_implementation/TestImplementation.jl") - -@testset "($N_INCOMING,$N_OUTGOING)" for (N_INCOMING, N_OUTGOING) in Iterators.product( - (1, rand(RNG, 2:8)), (1, rand(RNG, 2:8)) -) - INCOMING_PARTICLES = Tuple(rand(RNG, TestImplementation.PARTICLE_SET, N_INCOMING)) - OUTGOING_PARTICLES = Tuple(rand(RNG, TestImplementation.PARTICLE_SET, N_OUTGOING)) - - TESTPROC = TestImplementation.TestProcess(INCOMING_PARTICLES, OUTGOING_PARTICLES) - TESTMODEL = TestImplementation.TestModel() - TESTPSDEF = TestImplementation.TestPhasespaceDef() - IN_PS = TestImplementation._rand_momenta(RNG, N_INCOMING) - OUT_PS = TestImplementation._rand_momenta(RNG, N_OUTGOING) - PSP = PhaseSpacePoint(TESTPROC, TESTMODEL, TESTPSDEF, IN_PS, OUT_PS) - - @testset "failed interface" begin - TESTPROC_FAIL_ALL = TestImplementation.TestProcess_FAIL_ALL( - INCOMING_PARTICLES, OUTGOING_PARTICLES - ) - TESTPROC_FAIL_DIFFCS = TestImplementation.TestProcess_FAIL_DIFFCS( - INCOMING_PARTICLES, OUTGOING_PARTICLES - ) - TESTMODEL_FAIL = TestImplementation.TestModel_FAIL() - TESTPSDEF_FAIL = TestImplementation.TestPhasespaceDef_FAIL() - - @testset "failed process interface" begin - @test_throws MethodError incoming_particles(TESTPROC_FAIL_ALL) - @test_throws MethodError outgoing_particles(TESTPROC_FAIL_ALL) - end - - @testset "$PROC $MODEL" for (PROC, MODEL) in Iterators.product( - (TESTPROC, TESTPROC_FAIL_DIFFCS), (TESTMODEL, TESTMODEL_FAIL) - ) - if TestImplementation._any_fail(PROC, MODEL) - psp = PhaseSpacePoint(PROC, MODEL, TESTPSDEF, IN_PS, OUT_PS) - @test_throws MethodError QEDbase._incident_flux(psp) - @test_throws MethodError QEDbase._averaging_norm(psp) - @test_throws MethodError QEDbase._matrix_element(psp) - end - - for PS_DEF in (TESTPSDEF, TESTPSDEF_FAIL) - if TestImplementation._any_fail(PROC, MODEL, PS_DEF) - psp = PhaseSpacePoint(PROC, MODEL, PS_DEF, IN_PS, OUT_PS) - @test_throws MethodError QEDbase._phase_space_factor(psp) - end - end - end - end - - @testset "broadcast" begin - test_func(proc::AbstractProcessDefinition) = proc - @test test_func.(TESTPROC) == TESTPROC - - test_func(model::AbstractModelDefinition) = model - @test test_func.(TESTMODEL) == TESTMODEL - end - - @testset "incoming/outgoing particles" begin - @test incoming_particles(TESTPROC) == INCOMING_PARTICLES - @test outgoing_particles(TESTPROC) == OUTGOING_PARTICLES - @test number_incoming_particles(TESTPROC) == N_INCOMING - @test number_outgoing_particles(TESTPROC) == N_OUTGOING - end - - @testset "incident flux" begin - test_incident_flux = QEDbase._incident_flux( - InPhaseSpacePoint(TESTPROC, TESTMODEL, TESTPSDEF, IN_PS) - ) - groundtruth = TestImplementation._groundtruth_incident_flux(IN_PS) - @test isapprox(test_incident_flux, groundtruth, atol=ATOL, rtol=RTOL) - - test_incident_flux = QEDbase._incident_flux( - PhaseSpacePoint(TESTPROC, TESTMODEL, TESTPSDEF, IN_PS, OUT_PS) - ) - @test isapprox(test_incident_flux, groundtruth, atol=ATOL, rtol=RTOL) - - @test_throws MethodError QEDbase._incident_flux( - OutPhaseSpacePoint(TESTPROC, TESTMODEL, TESTPSDEF, OUT_PS) - ) - end - - @testset "averaging norm" begin - test_avg_norm = QEDbase._averaging_norm(TESTPROC) - groundtruth = TestImplementation._groundtruth_averaging_norm(TESTPROC) - @test isapprox(test_avg_norm, groundtruth, atol=ATOL, rtol=RTOL) - end - - @testset "matrix element" begin - test_matrix_element = QEDbase._matrix_element(PSP) - groundtruth = TestImplementation._groundtruth_matrix_element(IN_PS, OUT_PS) - @test length(test_matrix_element) == length(groundtruth) - for i in eachindex(test_matrix_element) - @test isapprox(test_matrix_element[i], groundtruth[i], atol=ATOL, rtol=RTOL) - end - end - - @testset "is in phasespace" begin - @test QEDbase._is_in_phasespace(PSP) - - IN_PS_unphysical = (zero(SFourMomentum), IN_PS[2:end]...) - OUT_PS_unphysical = (OUT_PS[1:(end - 1)]..., ones(SFourMomentum)) - PSP_unphysical_in_ps = PhaseSpacePoint( - TESTPROC, TESTMODEL, TESTPSDEF, IN_PS_unphysical, OUT_PS - ) - PSP_unphysical_out_ps = PhaseSpacePoint( - TESTPROC, TESTMODEL, TESTPSDEF, IN_PS, OUT_PS_unphysical - ) - PSP_unphysical = PhaseSpacePoint( - TESTPROC, TESTMODEL, TESTPSDEF, IN_PS_unphysical, OUT_PS_unphysical - ) - - @test !QEDbase._is_in_phasespace(PSP_unphysical_in_ps) - @test !QEDbase._is_in_phasespace(PSP_unphysical_out_ps) - @test !QEDbase._is_in_phasespace(PSP_unphysical) - end - - @testset "phase space factor" begin - test_phase_space_factor = QEDbase._phase_space_factor(PSP) - groundtruth = TestImplementation._groundtruth_phase_space_factor(IN_PS, OUT_PS) - @test isapprox(test_phase_space_factor, groundtruth, atol=ATOL, rtol=RTOL) - end -end diff --git a/test/interfaces/setup.jl b/test/interfaces/setup.jl deleted file mode 100644 index 44d5157..0000000 --- a/test/interfaces/setup.jl +++ /dev/null @@ -1,188 +0,0 @@ -using Random -using Suppressor -using QEDbase - -RNG = MersenneTwister(137137) -ATOL = 0.0 -RTOL = sqrt(eps()) - -_groundtruth_compute(x) = x -_groundtruth_input_validation(x) = (x > 0) -struct TestException <: QEDbase.AbstractInvalidInputException end -function _groundtruth_valid_input_assert(x) - _groundtruth_input_validation(x) || throw(TestException()) - return nothing -end -_transform_to_invalid(x) = -abs(x) -_groundtruth_post_processing(x, y) = x + y - -# setups for which the interface is implemented -abstract type AbstractTestSetup <: AbstractComputationSetup end -QEDbase._compute(stp::AbstractTestSetup, x) = _groundtruth_compute(x) - -# setup with default implementations -struct TestSetupDefault <: AbstractTestSetup end - -# setup with custom _assert_valid_input -struct TestSetupCustomAssertValidInput <: AbstractTestSetup end -function QEDbase._assert_valid_input(stp::TestSetupCustomAssertValidInput, x) - return _groundtruth_valid_input_assert(x) -end - -# setup with custom post processing -struct TestSetupCustomPostProcessing <: AbstractTestSetup end -function QEDbase._post_processing(::TestSetupCustomPostProcessing, x, y) - return _groundtruth_post_processing(x, y) -end - -# setup with custom input validation and post processing -struct TestSetupCustom <: AbstractTestSetup end -function QEDbase._assert_valid_input(stp::TestSetupCustom, x) - return _groundtruth_valid_input_assert(x) -end -QEDbase._post_processing(::TestSetupCustom, x, y) = _groundtruth_post_processing(x, y) - -# setup which fail on computation with default implementations -struct TestSetupFAIL <: AbstractComputationSetup end - -# setup which fail on computation with custom input validation, where the -# invalid input will be caught before the computation. -struct TestSetupCustomValidationFAIL <: AbstractComputationSetup end -function QEDbase._assert_valid_input(stp::TestSetupCustomValidationFAIL, x) - return _groundtruth_valid_input_assert(x) -end - -# setup which fail on computation with custom post processing -struct TestSetupCustomPostProcessingFAIL <: AbstractComputationSetup end -function QEDbase._post_processing(::TestSetupCustomPostProcessingFAIL, x, y) - return _groundtruth_post_processing(x, y) -end -@testset "general computation setup interface" begin - @testset "interface fail" begin - rnd_input = rand(RNG) - - @test_throws MethodError QEDbase._compute(TestSetupFAIL(), rnd_input) - @test_throws MethodError compute(TestSetupFAIL(), rnd_input) - - @test_throws MethodError QEDbase._compute( - TestSetupCustomValidationFAIL(), rnd_input - ) - @test_throws MethodError compute(TestSetupCustomValidationFAIL(), rnd_input) - # invalid input should be caught without throwing a MethodError - @test_throws TestException compute( - TestSetupCustomValidationFAIL(), _transform_to_invalid(rnd_input) - ) - - @test_throws MethodError QEDbase._compute( - TestSetupCustomPostProcessingFAIL(), rnd_input - ) - @test_throws MethodError compute(TestSetupCustomPostProcessingFAIL(), rnd_input) - end - - @testset "default interface" begin - stp = TestSetupDefault() - - rnd_input = rand(RNG) - rnd_output = rand(RNG) - @test QEDbase._post_processing(stp, rnd_input, rnd_output) == rnd_output - @test isapprox( - QEDbase._compute(stp, rnd_input), - _groundtruth_compute(rnd_input), - atol=ATOL, - rtol=RTOL, - ) - @test isapprox( - compute(stp, rnd_input), _groundtruth_compute(rnd_input), atol=ATOL, rtol=RTOL - ) - end - - @testset "custom input validation" begin - stp = TestSetupCustomAssertValidInput() - rnd_input = rand(RNG) - @test QEDbase._assert_valid_input(stp, rnd_input) == nothing - @test_throws TestException QEDbase._assert_valid_input( - stp, _transform_to_invalid(rnd_input) - ) - @test_throws TestException compute(stp, _transform_to_invalid(rnd_input)) - end - - @testset "custom post processing" begin - stp = TestSetupCustomPostProcessing() - rnd_input = rand(RNG) - rnd_output = rand(RNG) - @test isapprox( - QEDbase._post_processing(stp, rnd_input, rnd_output), - _groundtruth_post_processing(rnd_input, rnd_output), - ) - @test isapprox( - compute(stp, rnd_input), - _groundtruth_post_processing(rnd_input, _groundtruth_compute(rnd_input)), - ) - end - - @testset "custom input validation and post processing" begin - stp = TestSetupCustom() - rnd_input = rand(RNG) - rnd_output = rand(RNG) - - @test_throws TestException() compute(stp, _transform_to_invalid(rnd_input)) - @test isapprox( - QEDbase._post_processing(stp, rnd_input, rnd_output), - _groundtruth_post_processing(rnd_input, rnd_output), - ) - @test isapprox( - compute(stp, rnd_input), - _groundtruth_post_processing(rnd_input, _groundtruth_compute(rnd_input)), - ) - end -end -# process setup - -struct TestParticle1 <: AbstractParticle end -struct TestParticle2 <: AbstractParticle end -struct TestParticle3 <: AbstractParticle end -struct TestParticle4 <: AbstractParticle end - -PARTICLE_SET = [TestParticle1(), TestParticle2(), TestParticle3(), TestParticle4()] - -struct TestProcess <: AbstractProcessDefinition end -struct TestModel <: AbstractModelDefinition end - -struct TestProcessSetup <: AbstractProcessSetup end -QEDbase.scattering_process(::TestProcessSetup) = TestProcess() -QEDbase.physical_model(::TestProcessSetup) = TestModel() - -struct TestProcessSetupFAIL <: AbstractProcessSetup end - -@testset "process setup interface" begin - @testset "interface fail" begin - rnd_input = rand(RNG) - @test_throws MethodError scattering_process(TestProcessSetupFAIL()) - @test_throws MethodError physical_model(TestProcessSetupFAIL()) - @test_throws MethodError QEDbase._compute(TestProcessSetupFAIL(), rnd_input) - end - - @testset "hard interface" begin - stp = TestProcessSetup() - - @test QEDbase._is_computation_setup(stp) - @test scattering_process(stp) == TestProcess() - @test physical_model(stp) == TestModel() - end - - @testset "($N_INCOMING,$N_OUTGOING)" for (N_INCOMING, N_OUTGOING) in Iterators.product( - (1, rand(RNG, 2:8)), (1, rand(RNG, 2:8)) - ) - INCOMING_PARTICLES = rand(RNG, PARTICLE_SET, N_INCOMING) - OUTGOING_PARTICLES = rand(RNG, PARTICLE_SET, N_OUTGOING) - - @suppress QEDbase.incoming_particles(::TestProcess) = INCOMING_PARTICLES - @suppress QEDbase.outgoing_particles(::TestProcess) = OUTGOING_PARTICLES - - @testset "delegated functions" begin - stp = TestProcessSetup() - @test number_incoming_particles(stp) == N_INCOMING - @test number_outgoing_particles(stp) == N_OUTGOING - end - end -end diff --git a/test/particle_spinors.jl b/test/particle_spinors.jl deleted file mode 100644 index 26172df..0000000 --- a/test/particle_spinors.jl +++ /dev/null @@ -1,88 +0,0 @@ -using QEDbase -using Random - -const ATOL = 1e-15 - -const SPINS = (1, 2) - -@testset "particle spinors" for LorentzVectorType in [SFourMomentum, MFourMomentum] - rng = MersenneTwister(1234) - x, y, z = rand(rng, 3) - mass = rand(rng) - P = LorentzVectorType(sqrt(x^2 + y^2 + z^2 + mass^2), x, y, z) - - U = SpinorU(P, mass) - Ubar = SpinorUbar(P, mass) - V = SpinorV(P, mass) - Vbar = SpinorVbar(P, mass) - - @testset "construction" begin - @test U isa IncomingFermionSpinor - @test Ubar isa OutgoingFermionSpinor - @test V isa OutgoingAntiFermionSpinor - @test Vbar isa IncomingAntiFermionSpinor - - for spin in SPINS - @test isapprox(U[spin], U.booster * BASE_PARTICLE_SPINOR[spin]) - @test isapprox(V[spin], V.booster * BASE_ANTIPARTICLE_SPINOR[spin]) - - @test isapprox(Ubar[spin], AdjointBiSpinor(U[spin]) * GAMMA[1]) - @test isapprox(Vbar[spin], AdjointBiSpinor(V[spin]) * GAMMA[1]) - end - end # construction - - @testset "normatlisation" begin - for s1 in SPINS - for s2 in SPINS - @test isapprox(Ubar[s1] * U[s2], 2 * mass * (s1 == s2)) - @test isapprox(Vbar[s1] * V[s2], -2 * mass * (s1 == s2)) - @test isapprox(Ubar[s1] * V[s2], 0.0) - @test isapprox(Vbar[s1] * U[s2], 0.0) - end - end - end # normatlisation - - @testset "completeness" begin - sumU = zero(DiracMatrix) - sumV = zero(DiracMatrix) - for spin in SPINS - sumU += U(spin) * Ubar(spin) - sumV += V(spin) * Vbar(spin) - end - - @test isapprox(sumU, (slashed(P) + mass * one(DiracMatrix))) - @test isapprox(sumV, (slashed(P) - mass * one(DiracMatrix))) - end # completeness - - @testset "diracs equation" begin - for spin in SPINS - @test isapprox( - (slashed(P) - mass * one(DiracMatrix)) * U[spin], zero(BiSpinor), atol=ATOL - ) - @test isapprox( - (slashed(P) + mass * one(DiracMatrix)) * V[spin], zero(BiSpinor), atol=ATOL - ) - @test isapprox( - Ubar[spin] * (slashed(P) - mass * one(DiracMatrix)), - zero(AdjointBiSpinor), - atol=ATOL, - ) - @test isapprox( - Vbar[spin] * (slashed(P) + mass * one(DiracMatrix)), - zero(AdjointBiSpinor), - atol=ATOL, - ) - end - end #diracs equation - - @testset "sandwich" begin - for s1 in SPINS - for s2 in SPINS - @test isapprox( - LorentzVectorType(Ubar[s1] * (GAMMA * U[s2])) * (s1 == s2), - 2 * P * (s1 == s2), - ) - end - end - end #sandwich -end # particle spinors diff --git a/test/particles.jl b/test/particles.jl index 6129c1b..0370278 100644 --- a/test/particles.jl +++ b/test/particles.jl @@ -2,35 +2,6 @@ using QEDbase using StaticArrays using Random -include("utils.jl") - -FERMION_STATES_GROUNDTRUTH_FACTORY = Dict( - (Incoming, Electron) => IncomingFermionSpinor, - (Outgoing, Electron) => OutgoingFermionSpinor, - (Incoming, Positron) => IncomingAntiFermionSpinor, - (Outgoing, Positron) => OutgoingAntiFermionSpinor, -) - -RNG = MersenneTwister(708583836976) -ATOL = eps() -RTOL = 0.0 -PHOTON_ENERGIES = (0.0, rand(RNG), rand(RNG) * 10) -COS_THETAS = (-1.0, -rand(RNG), 0.0, rand(RNG), 1.0) -# check every quadrant -PHIS = ( - 0.0, - rand(RNG) * pi / 2, - pi / 2, - (1.0 + rand(RNG)) * pi / 2, - pi, - (2 + rand(RNG)) * pi / 2, - 3 * pi / 2, - (3 + rand(RNG)) * pi / 2, - 2 * pi, -) - -X, Y, Z = rand(RNG, 3) - # test function to test scalar broadcasting test_broadcast(x::AbstractParticle) = x test_broadcast(x::ParticleDirection) = x @@ -66,38 +37,6 @@ end end @testset "fermion likes" begin - @testset "fermion" begin - struct TestFermion <: Fermion end - @test is_fermion(TestFermion()) - @test is_particle(TestFermion()) - @test !is_anti_particle(TestFermion()) - @test test_broadcast.(TestFermion()) == TestFermion() - - @testset "$p $d" for (p, d) in - Iterators.product((Electron, Positron), (Incoming, Outgoing)) - mom = SFourMomentum(sqrt(mass(p()) + X^2 + Y^2 + Z^2), X, Y, Z) - particle_mass = mass(p()) - groundtruth_states = FERMION_STATES_GROUNDTRUTH_FACTORY[(d, p)]( - mom, particle_mass - ) - groundtruth_tuple = SVector(groundtruth_states(1), groundtruth_states(2)) - @test base_state(p(), d(), mom, AllSpin()) == groundtruth_tuple - @test base_state(p(), d(), mom, SpinUp()) == groundtruth_tuple[1] - @test base_state(p(), d(), mom, SpinDown()) == groundtruth_tuple[2] - - @test QEDbase._as_svec(base_state(p(), d(), mom, AllSpin())) isa SVector - @test QEDbase._as_svec(base_state(p(), d(), mom, SpinUp())) isa SVector - @test QEDbase._as_svec(base_state(p(), d(), mom, SpinDown())) isa SVector - - @test QEDbase._as_svec(base_state(p(), d(), mom, AllSpin())) == - groundtruth_tuple - @test QEDbase._as_svec(base_state(p(), d(), mom, SpinUp()))[1] == - groundtruth_tuple[1] - @test QEDbase._as_svec(base_state(p(), d(), mom, SpinDown()))[1] == - groundtruth_tuple[2] - end - end - @testset "antifermion" begin struct TestAntiFermion <: AntiFermion end @test is_fermion(TestAntiFermion()) @@ -161,55 +100,3 @@ end @test test_broadcast.(TestMajoranaBoson()) == TestMajoranaBoson() end end - -@testset "photon" begin - @test !is_fermion(Photon()) - @test is_boson(Photon()) - @test is_particle(Photon()) - @test is_anti_particle(Photon()) - @test charge(Photon()) == 0.0 - @test mass(Photon()) == 0.0 - @test test_broadcast.(Photon()) == Photon() - - @testset "$D" for D in [Incoming, Outgoing] - @testset "$om $cth $phi" for (om, cth, phi) in - Iterators.product(PHOTON_ENERGIES, COS_THETAS, PHIS) - #@testset "$x $y $z" for (x,y,z) in Iterators.product(X_arr,Y_arr,Z_arr) - - mom = SFourMomentum(_cartesian_coordinates(om, om, cth, phi)) - both_photon_states = base_state(Photon(), D(), mom, AllPolarization()) - - # property test the photon states - @test isapprox((both_photon_states[1] * mom), 0.0, atol=ATOL, rtol=RTOL) - @test isapprox((both_photon_states[2] * mom), 0.0, atol=ATOL, rtol=RTOL) - @test isapprox( - (both_photon_states[1] * both_photon_states[1]), -1.0, atol=ATOL, rtol=RTOL - ) - @test isapprox( - (both_photon_states[2] * both_photon_states[2]), -1.0, atol=ATOL, rtol=RTOL - ) - @test isapprox( - (both_photon_states[1] * both_photon_states[2]), 0.0, atol=ATOL, rtol=RTOL - ) - - # test the single polarization states - @test base_state(Photon(), D(), mom, PolarizationX()) == both_photon_states[1] - @test base_state(Photon(), D(), mom, PolarizationY()) == both_photon_states[2] - @test base_state(Photon(), D(), mom, PolX()) == both_photon_states[1] - @test base_state(Photon(), D(), mom, PolY()) == both_photon_states[2] - - @test QEDbase._as_svec(base_state(Photon(), D(), mom, PolX())) isa SVector - @test QEDbase._as_svec(base_state(Photon(), D(), mom, PolY())) isa SVector - @test QEDbase._as_svec(base_state(Photon(), D(), mom, AllPol())) isa SVector - - @test QEDbase._as_svec(base_state(Photon(), D(), mom, PolX()))[1] == - both_photon_states[1] - @test QEDbase._as_svec(base_state(Photon(), D(), mom, PolY()))[1] == - both_photon_states[2] - @test QEDbase._as_svec(base_state(Photon(), D(), mom, AllPol()))[1] == - both_photon_states[1] - @test QEDbase._as_svec(base_state(Photon(), D(), mom, AllPol()))[2] == - both_photon_states[2] - end - end -end diff --git a/test/runtests.jl b/test/runtests.jl index 67b8aa1..4e2da7f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -8,28 +8,8 @@ begin include("interfaces/model.jl") end - @time @safetestset "process interface" begin - include("interfaces/process.jl") - end - - @time @safetestset "computation setup interface" begin - include("interfaces/setup.jl") - end - @time @safetestset "Lorentz interface" begin - include("lorentz_interface.jl") - end - - @time @safetestset "Dirac tensors" begin - include("dirac_tensor.jl") - end - - @time @safetestset "Lorentz Vectors" begin - include("lorentz_vector.jl") - end - - @time @safetestset "particle spinors" begin - include("particle_spinors.jl") + include("interfaces/lorentz.jl") end @time @safetestset "particles" begin diff --git a/test/test_implementation/TestImplementation.jl b/test/test_implementation/TestImplementation.jl deleted file mode 100644 index 4cd4aac..0000000 --- a/test/test_implementation/TestImplementation.jl +++ /dev/null @@ -1,38 +0,0 @@ -""" -This module provides a full implementation of the model and process interface. Its purpose is only for testing and it does not reflect any -real-world physics. - -The module exports: - -``` -TestParticle1 # set of test particles without properties -TestParticle2 -TestParticle3 -TestParticle4 -TestModel # dummy compute model -TestModel_FAIL # failing compute model -TestProcess # dummy scattering process -TestProcess_FAIL # failing scattering process -TestPhasespaceDef # dummy phase space definition -TestPhasespaceDef_FAIL # failing phase space definition -``` -The respective groundtruth implementations for the interface functions are stored in `groundtruths.jl` -""" -module TestImplementation - -export TestParticle1, TestParticle2, TestParticle3, TestParticle4, PARTICLE_SET -export TestModel, TestModel_FAIL -export TestProcess, TestProcess_FAIL -export TestPhasespaceDef, TestPhasespaceDef_FAIL - -using Random -using QEDbase -using StaticArrays - -include("groundtruths.jl") -include("test_model.jl") -include("test_process.jl") -include("random_momenta.jl") -include("utils.jl") - -end diff --git a/test/test_implementation/groundtruths.jl b/test/test_implementation/groundtruths.jl deleted file mode 100644 index f05a010..0000000 --- a/test/test_implementation/groundtruths.jl +++ /dev/null @@ -1,265 +0,0 @@ -import QEDbase.AbstractFourMomentum - -""" - _groundtruth_incident_flux(in_ps) - -Test implementation of the incident flux. Return the Minkowski square of the sum of the incoming momenta: - -```math -\\begin{align} -I = \\left(\\sum p_i\\right)^2, -\\end{align} -``` -where \$p_i\\in\\mathrm{ps_in}\$. -""" -function _groundtruth_incident_flux(in_ps) - s = sum(in_ps) - return s * s -end - -""" - _groundtruth_matrix_element(in_ps, out_ps) - -Test implementation for a matrix elements. Returns a list of three complex numbers without any physical meaning. -""" -function _groundtruth_matrix_element(in_ps, out_ps) - s_in = sum(in_ps) - s_out = sum(out_ps) - res = s_in * s_in + 1im * (s_out * s_out) - return (res, 2 * res, 3 * res) -end - -""" - _groundtruth_averaging_norm(proc) - -Test implementation of the averaging norm. Returns the inverse of the sum of all external particles of the passed process. -""" -function _groundtruth_averaging_norm(proc) - return 1.0 / (number_incoming_particles(proc) + number_outgoing_particles(proc)) -end - -""" - _groundtruth_is_in_phasespace(in_ps, out_ps) - -Test implementation of the phase space check. Return `false` if either the momentum of the first incoming particle is exactly `zero(SFourMomentum)`, or if the momentum of the last outgoing momentum is exactly `ones(SFourMomentum)`. Otherwise, return true. -""" -function _groundtruth_is_in_phasespace(in_ps, out_ps) - if in_ps[1] == SFourMomentum(zeros(4)) - return false - end - if out_ps[end] == ones(SFourMomentum) - return false - end - return true -end - -""" - _groundtruth_phase_space_factor(in_ps, out_ps) - -Test implementation of the phase space factor. Return the inverse of the product of the energies of all external particles. -""" -function _groundtruth_phase_space_factor(in_ps, out_ps) - en_in = getE.(in_ps) - en_out = getE.(out_ps) - return 1 / (prod(en_in) * prod(en_out)) -end - -function _groundtruth_generate_momenta(ps_coords) - moms = _furl_moms(ps_coords) - return moms -end -""" - _groundtruth_unsafe_probability(proc, in_ps, out_ps) - -Test implementation of the unsafe differential probability. Uses the test implementations of `_groundtruth_matrix_element`,`_groundtruth_averaging_norm` and `_groundtruth_phase_space_factor`. -""" -function _groundtruth_unsafe_probability(proc, in_ps, out_ps) - mat_el = _groundtruth_matrix_element(in_ps, out_ps) - mat_el_sq = abs2.(mat_el) - normalization = _groundtruth_averaging_norm(proc) - ps_fac = _groundtruth_phase_space_factor(in_ps, out_ps) - return sum(mat_el_sq) * ps_fac * normalization -end - -function _groundtruth_unsafe_probability( - proc, in_ps::AbstractVector, out_ps::AbstractMatrix -) - res = Vector{Float64}(undef, size(out_ps, 2)) - for i in 1:size(out_ps, 2) - res[i] = _groundtruth_unsafe_probability(proc, in_ps, view(out_ps, :, i)) - end - return res -end - -function _groundtruth_unsafe_probability( - proc, in_ps::AbstractMatrix, out_ps::AbstractVector -) - res = Vector{Float64}(undef, size(in_ps, 2)) - for i in 1:size(in_ps, 2) - res[i] = _groundtruth_unsafe_probability(proc, view(in_ps, :, i), out_ps) - end - return res -end - -function _groundtruth_unsafe_probability( - proc, in_ps::AbstractMatrix, out_ps::AbstractMatrix -) - res = Matrix{Float64}(undef, size(in_ps, 2), size(out_ps, 2)) - for i in 1:size(in_ps, 2) - for j in 1:size(out_ps, 2) - res[i, j] = _groundtruth_unsafe_probability( - proc, view(in_ps, :, i), view(out_ps, :, j) - ) - end - end - return res -end - -""" - _groundtruth_safe_probability(proc, in_ps, out_ps) - -Test implementation of the safe differential probability. Uses the test implementations of `_groundtruth_is_in_phasespace` and `_groundtruth_unsafe_probability`. -""" -function _groundtruth_safe_probability(proc, in_ps, out_ps) - if !_groundtruth_is_in_phasespace(in_ps, out_ps) - return zero(Float64) - end - return _groundtruth_unsafe_probability(proc, in_ps, out_ps) -end - -function _groundtruth_safe_probability(proc, in_ps::AbstractVector, out_ps::AbstractMatrix) - res = Vector{Float64}(undef, size(out_ps, 2)) - for i in 1:size(out_ps, 2) - res[i] = _groundtruth_safe_probability(proc, in_ps, view(out_ps, :, i)) - end - return res -end - -function _groundtruth_safe_probability(proc, in_ps::AbstractMatrix, out_ps::AbstractVector) - res = Vector{Float64}(undef, size(in_ps, 2)) - for i in 1:size(in_ps, 2) - res[i] = _groundtruth_safe_probability(proc, view(in_ps, :, i), out_ps) - end - return res -end - -function _groundtruth_safe_probability(proc, in_ps::AbstractMatrix, out_ps::AbstractMatrix) - res = Matrix{Float64}(undef, size(in_ps, 2), size(out_ps, 2)) - for i in 1:size(in_ps, 2) - for j in 1:size(out_ps, 2) - res[i, j] = _groundtruth_safe_probability( - proc, view(in_ps, :, i), view(out_ps, :, j) - ) - end - end - return res -end - -""" - _groundtruth_unsafe_diffCS(proc, in_ps, out_ps) - -Test implementation of the unsafe differential cross section. Uses the test implementations of `_groundtruth_incident_flux` and `_groundtruth_unsafe_probability`. -""" -function _groundtruth_unsafe_diffCS(proc, in_ps, out_ps) - init_flux = _groundtruth_incident_flux(in_ps) - return _groundtruth_unsafe_probability(proc, in_ps, out_ps) / (4 * init_flux) -end - -function _groundtruth_unsafe_diffCS(proc, in_ps::AbstractVector, out_ps::AbstractMatrix) - res = Vector{Float64}(undef, size(out_ps, 2)) - for i in 1:size(out_ps, 2) - res[i] = _groundtruth_unsafe_diffCS(proc, in_ps, view(out_ps, :, i)) - end - return res -end - -function _groundtruth_unsafe_diffCS(proc, in_ps::AbstractMatrix, out_ps::AbstractVector) - res = Vector{Float64}(undef, size(in_ps, 2)) - for i in 1:size(in_ps, 2) - res[i] = _groundtruth_unsafe_diffCS(proc, view(in_ps, :, i), out_ps) - end - return res -end - -function _groundtruth_unsafe_diffCS(proc, in_ps::AbstractMatrix, out_ps::AbstractMatrix) - res = Matrix{Float64}(undef, size(in_ps, 2), size(out_ps, 2)) - for i in 1:size(in_ps, 2) - for j in 1:size(out_ps, 2) - res[i, j] = _groundtruth_unsafe_diffCS( - proc, view(in_ps, :, i), view(out_ps, :, j) - ) - end - end - return res -end - -""" - _groundtruth_safe_diffCS(proc, in_ps, out_ps) - -Test implementation of the safe differential cross section. Uses the test implementations of `_groundtruth_is_in_phasespace` and `_groundtruth_unsafe_diffCS`. -""" -function _groundtruth_safe_diffCS(proc, in_ps, out_ps) - if !_groundtruth_is_in_phasespace(in_ps, out_ps) - return zero(Float64) - end - return _groundtruth_unsafe_diffCS(proc, in_ps, out_ps) -end - -function _groundtruth_safe_diffCS(proc, in_ps::AbstractVector, out_ps::AbstractMatrix) - res = Vector{Float64}(undef, size(out_ps, 2)) - for i in 1:size(out_ps, 2) - res[i] = _groundtruth_safe_diffCS(proc, in_ps, view(out_ps, :, i)) - end - return res -end - -function _groundtruth_safe_diffCS(proc, in_ps::AbstractMatrix, out_ps::AbstractVector) - res = Vector{Float64}(undef, size(in_ps, 2)) - for i in 1:size(in_ps, 2) - res[i] = _groundtruth_safe_diffCS(proc, view(in_ps, :, i), out_ps) - end - return res -end - -function _groundtruth_safe_diffCS(proc, in_ps::AbstractMatrix, out_ps::AbstractMatrix) - res = Matrix{Float64}(undef, size(in_ps, 2), size(out_ps, 2)) - for i in 1:size(in_ps, 2) - for j in 1:size(out_ps, 2) - res[i, j] = _groundtruth_safe_diffCS( - proc, view(in_ps, :, i), view(out_ps, :, j) - ) - end - end - return res -end - -""" - _groundtruth_total_probability(in_ps::AbstractVector) - -Test implementation of the total cross section. Return the Minkowski square of the sum the momenta of all incoming particles. -""" -function _groundtruth_total_probability( - in_ps::NTuple{N,T} -) where {N,T<:AbstractFourMomentum} - Ptot = sum(in_ps) - return Ptot * Ptot -end - -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::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_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 deleted file mode 100644 index b631551..0000000 --- a/test/test_implementation/random_momenta.jl +++ /dev/null @@ -1,83 +0,0 @@ - -""" -Return a tuple of random four momenta, i.e. a random phase space point. -""" -function _rand_momenta(rng::AbstractRNG, n) - return NTuple{n,SFourMomentum}(SFourMomentum(rand(rng, 4)) for _ in 1:n) -end - -""" -Return a vector of tuples of random four momenta, i.e. a collection of phase space points. -n1 is the size of the phase space point, n2 is the number of points. -""" -function _rand_momenta(rng::AbstractRNG, n1, n2) - moms = Vector{NTuple{n1,SFourMomentum}}(undef, n2) - for i in 1:n2 - moms[i] = _rand_momenta(rng, n1) - end - return moms -end - -""" -Return a random phase space point that is failing the incoming phase space constraint, -i.e. the first entry of the phase space is the null momentum. -""" -function _rand_in_momenta_failing(rng::AbstractRNG, n) - return (zero(SFourMomentum), _rand_momenta(rng, n - 1)...) -end - -""" -Return a random phase space point that is failing the outgoing phase space constraint, -i.e. the last entry of the phase space is the unit momentum. -""" -function _rand_out_momenta_failing(rng::AbstractRNG, n) - return (_rand_momenta(rng, n - 1)..., ones(SFourMomentum)) -end - -""" -Return a collection of incoming phase space points, where the first point is failing the phase space constraint, -i.e. the first entry of the vector is invalid but the others pass. -n1 is the size of the phase space point, n2 is the number of points. -""" -function _rand_in_momenta_failing_mix(rng::AbstractRNG, n1, n2) - moms = _rand_momenta(rng, n1, n2) - moms[1] = _rand_in_momenta_failing(rng, n1) - return moms -end - -""" -Return a collection of incoming phase space points, where all points are failing the phase space constraint, -i.e. their first entries are null momenta. -n1 is the size of the phase space point, n2 is the number of points. -""" -function _rand_in_momenta_failing_all(rng::AbstractRNG, n1, n2) - moms = Vector{NTuple{n1,SFourMomentum}}(undef, n2) - for i in 1:n2 - moms[i] = _rand_in_momenta_failing(rng, n1) - end - return moms -end - -""" -Return a vector of outgoing phase space points, where the first point is failing the phase space constraint, -i.e. the last entry of the vector is invalid but the others pass. -n1 is the size of the phase space point, n2 is the number of points. -""" -function _rand_out_momenta_failing_mix(rng::AbstractRNG, n1, n2) - moms = _rand_momenta(rng, n1, n2) - moms[end] = _rand_out_momenta_failing(rng, n1) - return moms -end - -""" -Return a vector of outgoing phase space points, where all points are failing the phase space constraint, -i.e. their last entries are unit momenta. -n1 is the size of the phase space point, n2 is the number of points. -""" -function _rand_out_momenta_failing_all(rng::AbstractRNG, n1, n2) - moms = Vector{NTuple{n1,SFourMomentum}}(undef, n2) - for i in 1:n2 - moms[i] = _rand_out_momenta_failing(rng, n1) - end - return moms -end diff --git a/test/test_implementation/test_model.jl b/test/test_implementation/test_model.jl deleted file mode 100644 index 87f4f04..0000000 --- a/test/test_implementation/test_model.jl +++ /dev/null @@ -1,4 +0,0 @@ -struct TestModel <: AbstractModelDefinition end -QEDbase.fundamental_interaction_type(::TestModel) = :test_interaction - -struct TestModel_FAIL <: AbstractModelDefinition end diff --git a/test/test_implementation/test_process.jl b/test/test_implementation/test_process.jl deleted file mode 100644 index 5a81146..0000000 --- a/test/test_implementation/test_process.jl +++ /dev/null @@ -1,131 +0,0 @@ -# dummy particles -struct TestParticleFermion <: FermionLike end -struct TestParticleBoson <: BosonLike end - -const PARTICLE_SET = [TestParticleFermion(), TestParticleBoson()] - -""" - TestProcess(rng,incoming_particles,outgoing_particles) - -""" -struct TestProcess{IP<:Tuple,OP<:Tuple} <: AbstractProcessDefinition - incoming_particles::IP - outgoing_particles::OP -end - -function TestProcess(rng::AbstractRNG, N_in::Int, N_out::Int) - in_particles = rand(rng, PARTICLE_SET, N_in) - out_particles = rand(rng, PARTICLE_SET, N_out) - return TestProcess(in_particles, out_particles) -end - -QEDbase.incoming_particles(proc::TestProcess) = proc.incoming_particles -QEDbase.outgoing_particles(proc::TestProcess) = proc.outgoing_particles - -struct TestProcess_FAIL{IP<:Tuple,OP<:Tuple} <: AbstractProcessDefinition - incoming_particles::IP - outgoing_particles::OP -end - -function TestProcess_FAIL(rng::AbstractRNG, N_in::Int, N_out::Int) - in_particles = Tuple(rand(rng, PARTICLE_SET, N_in)) - out_particles = Tuple(rand(rng, PARTICLE_SET, N_out)) - return TestProcess_FAIL(in_particles, out_particles) -end - -function QEDbase.in_phase_space_dimension(proc::TestProcess, ::TestModel) - return number_incoming_particles(proc) * 4 -end -function QEDbase.out_phase_space_dimension(proc::TestProcess, ::TestModel) - return number_outgoing_particles(proc) * 4 -end - -""" -Test process with no implemented interface. Should fail every usage except construction. -""" -struct TestProcess_FAIL_ALL{IP<:Tuple,OP<:Tuple} <: AbstractProcessDefinition - incoming_particles::IP - outgoing_particles::OP -end - -function TestProcess_FAIL_ALL(rng::AbstractRNG, N_in::Int, N_out::Int) - in_particles = Tuple(rand(rng, PARTICLE_SET, N_in)) - out_particles = Tuple(rand(rng, PARTICLE_SET, N_out)) - return TestProcess_FAIL_ALL(in_particles, out_particles) -end - -""" -Test process with no implemented interface except the incoming and outgoing particles. -Should fail every usage except construction of itself and the respective phase space point for given four-momenta. -""" -struct TestProcess_FAIL_DIFFCS{IP<:Tuple,OP<:Tuple} <: AbstractProcessDefinition - incoming_particles::IP - outgoing_particles::OP -end - -function TestProcess_FAIL_DIFFCS(rng::AbstractRNG, N_in::Int, N_out::Int) - in_particles = Tuple(rand(rng, PARTICLE_SET, N_in)) - out_particles = Tuple(rand(rng, PARTICLE_SET, N_out)) - return TestProcess_FAIL_DIFFCS(in_particles, out_particles) -end - -QEDbase.incoming_particles(proc::TestProcess_FAIL_DIFFCS) = proc.incoming_particles -QEDbase.outgoing_particles(proc::TestProcess_FAIL_DIFFCS) = proc.outgoing_particles - -# dummy phase space definition + failing phase space definition -struct TestPhasespaceDef <: AbstractPhasespaceDefinition end -struct TestPhasespaceDef_FAIL <: AbstractPhasespaceDefinition end - -# dummy implementation of the process interface - -function QEDbase._incident_flux(in_psp::InPhaseSpacePoint{<:TestProcess,<:TestModel}) - return _groundtruth_incident_flux(momenta(in_psp, Incoming())) -end - -function QEDbase._averaging_norm(proc::TestProcess) - return _groundtruth_averaging_norm(proc) -end - -function QEDbase._matrix_element(psp::PhaseSpacePoint{<:TestProcess,TestModel}) - in_ps = momenta(psp, Incoming()) - out_ps = momenta(psp, Outgoing()) - return _groundtruth_matrix_element(in_ps, out_ps) -end - -function QEDbase._is_in_phasespace(psp::PhaseSpacePoint{<:TestProcess,TestModel}) - in_ps = momenta(psp, Incoming()) - out_ps = momenta(psp, Outgoing()) - return _groundtruth_is_in_phasespace(in_ps, out_ps) -end - -function QEDbase._phase_space_factor( - psp::PhaseSpacePoint{<:TestProcess,TestModel,TestPhasespaceDef} -) - in_ps = momenta(psp, Incoming()) - out_ps = momenta(psp, Outgoing()) - return _groundtruth_phase_space_factor(in_ps, out_ps) -end - -function QEDbase._generate_incoming_momenta( - proc::TestProcess, - model::TestModel, - phase_space_def::TestPhasespaceDef, - in_phase_space::NTuple{N,T}, -) where {N,T<:Real} - return _groundtruth_generate_momenta(in_phase_space) -end - -function QEDbase._generate_outgoing_momenta( - proc::TestProcess, - model::TestModel, - phase_space_def::TestPhasespaceDef, - out_phase_space::NTuple{N,T}, -) where {N,T<:Real} - return _groundtruth_generate_momenta(out_phase_space) -end - -function QEDbase._total_probability( - in_psp::InPhaseSpacePoint{<:TestProcess,<:TestModel,<:TestPhasespaceDef} -) - return _groundtruth_total_probability(momenta(in_psp, Incoming())) -end diff --git a/test/test_implementation/utils.jl b/test/test_implementation/utils.jl deleted file mode 100644 index 3a141d9..0000000 --- a/test/test_implementation/utils.jl +++ /dev/null @@ -1,44 +0,0 @@ - -# Check if any failed type is in the input -_any_fail(x...) = true -_any_fail(::TestProcess, ::TestModel) = false -_any_fail(::TestProcess, ::TestModel, ::TestPhasespaceDef) = false - -# unrolls all elements of a list of four-momenta into vector of coordinates -function _unroll_moms(ps_moms::AbstractVector{T}) where {T<:QEDbase.AbstractFourMomentum} - return collect(Iterators.flatten(ps_moms)) -end - -function _unroll_moms(ps_moms::AbstractMatrix{T}) where {T<:QEDbase.AbstractFourMomentum} - res = Matrix{eltype(T)}(undef, size(ps_moms, 1) * 4, size(ps_moms, 2)) - for i in 1:size(ps_moms, 2) - res[:, i] .= _unroll_moms(view(ps_moms, :, i)) - end - return res -end - -flat_components(moms::AbstractVecOrMat) = _unroll_moms(moms) -flat_components(moms::Tuple) = Tuple(_unroll_moms([moms...])) - -# collect components of four-momenta from a vector of coordinates -function __furl_moms(ps_coords::AbstractVector{T}) where {T<:Real} - return SFourMomentum.(eachcol(reshape(ps_coords, 4, :))) -end - -function _furl_moms(ps_coords::AbstractVector{T}) where {T<:Real} - @assert length(ps_coords) % 4 == 0 - return __furl_moms(ps_coords) -end - -function _furl_moms(ps_coords::AbstractMatrix{T}) where {T<:Real} - @assert size(ps_coords, 1) % 4 == 0 - res = Matrix{SFourMomentum}(undef, Int(size(ps_coords, 1)//4), size(ps_coords, 2)) - for i in 1:size(ps_coords, 2) - res[:, i] .= __furl_moms(view(ps_coords, :, i)) - end - return res -end - -function _furl_moms(moms::NTuple{N,Float64}) where {N} - return Tuple(_furl_moms(Vector{Float64}([moms...]))) -end diff --git a/test/utils.jl b/test/utils.jl deleted file mode 100644 index 522a9bd..0000000 --- a/test/utils.jl +++ /dev/null @@ -1,10 +0,0 @@ -""" -Returns the cartesian coordinates (E,px,py,pz) for given spherical coordiantes -(E, rho, cos_theta, phi), where rho denotes the length of the respective three-momentum, -theta is the polar- and phi the azimuthal angle. -""" -function _cartesian_coordinates(E, rho, cth, phi) - sth = sqrt(1 - cth^2) - sphi, cphi = sincos(phi) - return (E, rho * sth * cphi, rho * sth * sphi, rho * cth) -end From 62607f2052d27b41adf657dbf3b92e1b0a42d98d Mon Sep 17 00:00:00 2001 From: Anton Reinhard Date: Thu, 6 Jun 2024 18:20:38 +0200 Subject: [PATCH 10/21] Fix type casting issue --- src/QEDbase.jl | 2 +- src/interfaces/lorentz_vectors/arithmetic.jl | 13 +++++++++++++ 2 files changed, 14 insertions(+), 1 deletion(-) diff --git a/src/QEDbase.jl b/src/QEDbase.jl index 5f955d9..1bde44d 100644 --- a/src/QEDbase.jl +++ b/src/QEDbase.jl @@ -1,6 +1,6 @@ module QEDbase -import Base: * +import Base: *, / import StaticArrays: similar_type export minkowski_dot, mdot, register_LorentzVectorLike diff --git a/src/interfaces/lorentz_vectors/arithmetic.jl b/src/interfaces/lorentz_vectors/arithmetic.jl index 2247eaa..f7ef976 100644 --- a/src/interfaces/lorentz_vectors/arithmetic.jl +++ b/src/interfaces/lorentz_vectors/arithmetic.jl @@ -1,3 +1,16 @@ +# TODO: these overloads shouldn't be necessary, but are for some reason +@inline function *(L::TL, NUM::Number) where {TL<:AbstractLorentzVector} + return constructorof(TL)(L[1] * NUM, L[2] * NUM, L[3] * NUM, L[4] * NUM) +end + +@inline function *(NUM::Number, L::TL) where {TL<:AbstractLorentzVector} + return constructorof(TL)(L[1] * NUM, L[2] * NUM, L[3] * NUM, L[4] * NUM) +end + +@inline function /(L::TL, NUM::Number) where {TL<:AbstractLorentzVector} + return constructorof(TL)(L[1] / NUM, L[2] / NUM, L[3] / NUM, L[4] / NUM) +end + function dot(p1::T1, p2::T2) where {T1<:AbstractLorentzVector,T2<:AbstractLorentzVector} return mdot(p1, p2) end From d9c271a74c61b8e9c0f1a9bd4e92175876fff463 Mon Sep 17 00:00:00 2001 From: Anton Reinhard Date: Thu, 6 Jun 2024 18:35:34 +0200 Subject: [PATCH 11/21] Move SpinorConstructionError to errors.jl --- src/QEDbase.jl | 2 +- src/errors.jl | 8 ++++++++ 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/src/QEDbase.jl b/src/QEDbase.jl index 1bde44d..cf5c6ed 100644 --- a/src/QEDbase.jl +++ b/src/QEDbase.jl @@ -80,7 +80,7 @@ export AbstractProcessSetup, scattering_process, physical_model export AbstractCoordinateSystem, AbstractFrameOfReference, AbstractPhasespaceDefinition # errors -export InvalidInputError, RegistryError, OnshellError +export InvalidInputError, RegistryError, OnshellError, SpinorConstructionError using StaticArrays using LinearAlgebra diff --git a/src/errors.jl b/src/errors.jl index d0e2442..aacd618 100644 --- a/src/errors.jl +++ b/src/errors.jl @@ -55,3 +55,11 @@ end function Base.showerror(io::IO, err::InvalidInputError) return println(io, "InvalidInputError: $(err.msg)") end + +mutable struct SpinorConstructionError <: Exception + var::String +end + +function Base.showerror(io::IO, e::SpinorConstructionError) + return print(io, "SpinorConstructionError: ", e.var) +end From bc811344e751fd4ef1558bc44c9287dd7736abff Mon Sep 17 00:00:00 2001 From: Anton Reinhard Date: Wed, 12 Jun 2024 11:04:44 +0200 Subject: [PATCH 12/21] Add AbstractParticleSpinor to exports --- src/QEDbase.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/QEDbase.jl b/src/QEDbase.jl index cf5c6ed..0bb263a 100644 --- a/src/QEDbase.jl +++ b/src/QEDbase.jl @@ -40,7 +40,7 @@ export base_state export mass, charge # particle types -export AbstractParticleType +export AbstractParticleType, AbstractParticleSpinor export FermionLike, Fermion, AntiFermion, MajoranaFermion export BosonLike, Boson, AntiBoson, MajoranaBoson export Electron, Positron, Photon From abdfb4c6ec84f14a75a629146065626988b6d8d4 Mon Sep 17 00:00:00 2001 From: Anton Reinhard Date: Wed, 19 Jun 2024 15:39:29 +0200 Subject: [PATCH 13/21] Move concrete particle types to core --- ...{BuildDelopyDoc.yml => BuildDeployDoc.yml} | 0 docs/Project.toml | 1 + src/QEDbase.jl | 13 +- src/interfaces/particle.jl | 180 +++++++++++++ src/interfaces/particle_functions.jl | 145 ----------- src/interfaces/particle_types.jl | 243 ------------------ test/coverage/.gitkeep | 0 test/particle_properties.jl | 37 +++ test/particles.jl | 102 -------- test/runtests.jl | 4 +- 10 files changed, 223 insertions(+), 502 deletions(-) rename .github/workflows/{BuildDelopyDoc.yml => BuildDeployDoc.yml} (100%) delete mode 100644 src/interfaces/particle_functions.jl delete mode 100644 src/interfaces/particle_types.jl delete mode 100644 test/coverage/.gitkeep create mode 100644 test/particle_properties.jl delete mode 100644 test/particles.jl diff --git a/.github/workflows/BuildDelopyDoc.yml b/.github/workflows/BuildDeployDoc.yml similarity index 100% rename from .github/workflows/BuildDelopyDoc.yml rename to .github/workflows/BuildDeployDoc.yml diff --git a/docs/Project.toml b/docs/Project.toml index b927490..a82680f 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -4,3 +4,4 @@ Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DocumenterCitations = "daee34ce-89f3-4625-b898-19384cb65244" DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8" QEDbase = "10e22c08-3ccb-4172-bfcf-7d7aa3d04d93" +QEDcore = "35dc0263-cb5f-4c33-a114-1d7f54ab753e" diff --git a/src/QEDbase.jl b/src/QEDbase.jl index 0bb263a..6c41a00 100644 --- a/src/QEDbase.jl +++ b/src/QEDbase.jl @@ -34,19 +34,14 @@ export mul export AbstractGammaRepresentation # particle interface -export AbstractParticle +export AbstractParticle, AbstractParticleType export is_fermion, is_boson, is_particle, is_anti_particle -export base_state export mass, charge +export base_state, propagator -# particle types -export AbstractParticleType, AbstractParticleSpinor -export FermionLike, Fermion, AntiFermion, MajoranaFermion -export BosonLike, Boson, AntiBoson, MajoranaBoson -export Electron, Positron, Photon +# directions export ParticleDirection, Incoming, Outgoing export is_incoming, is_outgoing -export propagator # polarizations and spins export AbstractSpinOrPolarization, AbstractPolarization, AbstractSpin @@ -107,8 +102,6 @@ include("interfaces/four_momentum.jl") include("interfaces/model.jl") include("interfaces/particle.jl") -include("interfaces/particle_types.jl") -include("interfaces/particle_functions.jl") include("particles/direction.jl") include("particles/spin_pol.jl") diff --git a/src/interfaces/particle.jl b/src/interfaces/particle.jl index 494a310..b18a8b8 100644 --- a/src/interfaces/particle.jl +++ b/src/interfaces/particle.jl @@ -27,6 +27,16 @@ These functions must be implemented in order to have the subtype of `AbstractPar abstract type AbstractParticle end Base.broadcastable(part::AbstractParticle) = Ref(part) +""" + AbstractParticleType <: AbstractParticle + +This is the abstract base type for every species of particles. All functionalities defined on subtypes of `AbstractParticleType` should be static, i.e. known at compile time. +For adding runtime information, e.g. four-momenta or particle states, to a particle, consider implementing a concrete subtype of [`AbstractParticle`](@ref) instead, which may have a type parameter `P<:AbstractParticleType`. + +Concrete built-in subtypes of `AbstractParticleType` are available in `QEDcore.jl` and should always be singletons.. +""" +abstract type AbstractParticleType <: AbstractParticle end + """ $(TYPEDSIGNATURES) @@ -76,3 +86,173 @@ Interface function for particles. Return the electric charge of a particle (in u This needs to be implemented for each concrete subtype of [`AbstractParticle`](@ref). """ function charge end + +""" + propagator(particle::AbstractParticleType, mom::QEDbase.AbstractFourMomentum, [mass::Real]) + +Return the propagator of a particle for a given four-momentum. If `mass` is passed, the respective propagator for massive particles is used, if not, it is assumed the particle passed in is massless. + +!!! note "Convention" + + There are two types of implementations for propagators given in `QEDProcesses`: + For a `BosonLike` particle with four-momentum ``k`` and mass ``m``, the propagator is given as + + ```math + D(k) = \\frac{1}{k^2 - m^2}. + ``` + + For a `FermionLike` particle with four-momentum ``p`` and mass ``m``, the propagator is given as + + ```math + S(p) = \\frac{\\gamma^\\mu p_\\mu + mass}{p^2 - m^2}. + ``` + +!!! warning + + This function does not throw when the given particle is off-shell. If an off-shell particle is passed, the function `propagator` returns `Inf`. + +""" +function propagator end + +""" +```julia + base_state( + particle::AbstractParticleType, + direction::ParticleDirection, + momentum::QEDbase.AbstractFourMomentum, + [spin_or_pol::AbstractSpinOrPolarization] + ) +``` + +Return the base state of a directed on-shell particle with a given four-momentum. For internal usage only. + +# Input + +- `particle` -- the type of the particle, e.g. `QEDcore.Electron`, `QEDcore.Positron`, or `QEDcore.Photon`. +- `direction` -- the direction of the particle, i.e. [`Incoming`](@ref) or [`Outgoing`](@ref). +- `momentum` -- the four-momentum of the particle +- `[spin_or_pol]` -- if given, the spin or polarization of the particle, e.g. [`SpinUp`](@ref)/[`SpinDown`](@ref) or [`PolarizationX`](@ref)/[`PolarizationY`](@ref). + +# Output +The output type of `base_state` depends on wether the spin or polarization of the particle passed in is specified or not. + +If `spin_or_pol` is passed, the output of `base_state` is + +```julia +base_state(::Fermion, ::Incoming, mom, spin_or_pol) # -> BiSpinor +base_state(::AntiFermion, ::Incoming, mom, spin_or_pol) # -> AdjointBiSpinor +base_state(::Fermion, ::Outgoing, mom, spin_or_pol) # -> AdjointBiSpinor +base_state(::AntiFermion, ::Outgoing, mom, spin_or_pol) # -> BiSpinor +base_state(::Photon, ::Incoming, mom, spin_or_pol) # -> SLorentzVector{ComplexF64} +base_state(::Photon, ::Outgoing, mom, spin_or_pol) # -> SLorentzVector{ComplexF64} +``` + +If `spin_or_pol` is of type [`AllPolarization`](@ref) or [`AllSpin`](@ref), the output is an `SVector` with both spin/polarization alignments: + +```julia +base_state(::Fermion, ::Incoming, mom) # -> SVector{2,BiSpinor} +base_state(::AntiFermion, ::Incoming, mom) # -> SVector{2,AdjointBiSpinor} +base_state(::Fermion, ::Outgoing, mom) # -> SVector{2,AdjointBiSpinor} +base_state(::AntiFermion, ::Outgoing, mom) # -> SVector{2,BiSpinor} +base_state(::Photon, ::Incoming, mom) # -> SVector{2,SLorentzVector{ComplexF64}} +base_state(::Photon, ::Outgoing, mom) # -> SVector{2,SLorentzVector{ComplexF64}} +``` + +# Example + +```julia +using QEDbase + +mass = 1.0 # set electron mass to 1.0 +px,py,pz = rand(3) # generate random spatial components +E = sqrt(px^2 + py^2 + pz^2 + mass^2) # compute energy, i.e. the electron is on-shell +mom = SFourMomentum(E, px, py, pz) # initialize the four-momentum of the electron + +# compute the state of an incoming electron with spin = SpinUp +# note: base_state is not exported! +electron_state = base_state(QEDcore.Electron(), Incoming(), mom, SpinUp()) +``` + +```jldoctest +julia> using QEDbase; using QEDcore + +julia> mass = 1.0; px,py,pz = (0.1, 0.2, 0.3); E = sqrt(px^2 + py^2 + pz^2 + mass^2); mom = SFourMomentum(E, px, py, pz) +4-element SFourMomentum with indices SOneTo(4): + 1.0677078252031311 + 0.1 + 0.2 + 0.3 + +julia> electron_state = base_state(Electron(), Incoming(), mom, SpinUp()) +4-element BiSpinor with indices SOneTo(4): + 1.4379526505428235 + 0.0im + 0.0 + 0.0im + -0.20862995724285552 + 0.0im + -0.06954331908095185 - 0.1390866381619037im + +julia> electron_states = base_state(Electron(), Incoming(), mom, AllSpin()) +2-element StaticArraysCore.SVector{2, BiSpinor} with indices SOneTo(2): + [1.4379526505428235 + 0.0im, 0.0 + 0.0im, -0.20862995724285552 + 0.0im, -0.06954331908095185 - 0.1390866381619037im] + [0.0 + 0.0im, 1.4379526505428235 + 0.0im, -0.06954331908095185 + 0.1390866381619037im, 0.20862995724285552 + 0.0im] +``` + +!!! note "Iterator convenience" + The returned objects of `base_state` can be consistently wrapped in an `SVector` for iteration using [`_as_svec`](@ref). + + This way, a loop like the following becomes possible when `spin` may be definite or indefinite. + ```julia + for state in QEDbase._as_svec(base_state(Electron(), Incoming(), momentum, spin)) + # ... + end + ``` + +!!! note "Conventions" + + For an incoming fermion with momentum ``p``, we use the explicit formula: + + ```math + u_\\sigma(p) = \\frac{\\gamma^\\mu p_\\mu + m}{\\sqrt{\\vert p_0\\vert + m}} \\eta_\\sigma, + ``` + + where the elementary base spinors are given as + + ```math + \\eta_1 = (1, 0, 0, 0)^T\\\\ + \\eta_2 = (0, 1, 0, 0)^T + ``` + + For an outgoing anti-fermion with momentum ``p``, we use the explicit formula: + + ```math + v_\\sigma(p) = \\frac{-\\gamma^\\mu p_\\mu + m}{\\sqrt{\\vert p_0\\vert + m}} \\chi_\\sigma, + ``` + + where the elementary base spinors are given as + + ```math + \\chi_1 = (0, 0, 1, 0)^T\\\\ + \\chi_2 = (0, 0, 0, 1)^T + ``` + + For outgoing fermions and incoming anti-fermions with momentum ``p``, the base state is given as the Dirac-adjoint of the respective incoming fermion or outgoing anti-fermion state: + + ```math + \\overline{u}_\\sigma(p) = u_\\sigma^\\dagger \\gamma^0\\\\ + \\overline{v}_\\sigma(p) = v_\\sigma^\\dagger \\gamma^0 + ``` + + where ``v_\\sigma`` is the base state of the respective outgoing anti-fermion. + + For a photon with four-momentum ``k^\\mu = \\omega (1, \\cos\\phi \\sin\\theta, \\sin\\phi \\sin\\theta, \\cos\\theta)``, the two polarization vectors are given as + + ```math + \\begin{align*} + \\epsilon^\\mu_1 &= (0, \\cos\\theta \\cos\\phi, \\cos\\theta \\sin\\phi, -\\sin\\theta),\\\\ + \\epsilon^\\mu_2 &= (0, -\\sin\\phi, \\cos\\phi, 0). + \\end{align*} + ``` + +!!! warning + In the current implementation there are **no checks** built-in, which verify the passed arguments whether they describe on-shell particles, i.e. `p*p≈mass^2`. Using `base_state` with off-shell particles will cause unpredictable behavior. +""" +function base_state end diff --git a/src/interfaces/particle_functions.jl b/src/interfaces/particle_functions.jl deleted file mode 100644 index 9d89943..0000000 --- a/src/interfaces/particle_functions.jl +++ /dev/null @@ -1,145 +0,0 @@ - -""" -```julia - base_state( - particle::AbstractParticleType, - direction::ParticleDirection, - momentum::QEDbase.AbstractFourMomentum, - [spin_or_pol::AbstractSpinOrPolarization] - ) -``` - -Return the base state of a directed on-shell particle with a given four-momentum. For internal usage only. - -# Input - -- `particle` -- the type of the particle, e.g. [`Electron`](@ref), [`Positron`](@ref), or [`Photon`](@ref). -- `direction` -- the direction of the particle, i.e. [`Incoming`](@ref) or [`Outgoing`](@ref). -- `momentum` -- the four-momentum of the particle -- `[spin_or_pol]` -- if given, the spin or polarization of the particle, e.g. [`SpinUp`](@ref)/[`SpinDown`](@ref) or [`PolarizationX`](@ref)/[`PolarizationY`](@ref). - -# Output -The output type of `base_state` depends on wether the spin or polarization of the particle passed in is specified or not. - -If `spin_or_pol` is passed, the output of `base_state` is - -```julia -base_state(::Fermion, ::Incoming, mom, spin_or_pol) # -> BiSpinor -base_state(::AntiFermion, ::Incoming, mom, spin_or_pol) # -> AdjointBiSpinor -base_state(::Fermion, ::Outgoing, mom, spin_or_pol) # -> AdjointBiSpinor -base_state(::AntiFermion, ::Outgoing, mom, spin_or_pol) # -> BiSpinor -base_state(::Photon, ::Incoming, mom, spin_or_pol) # -> SLorentzVector{ComplexF64} -base_state(::Photon, ::Outgoing, mom, spin_or_pol) # -> SLorentzVector{ComplexF64} -``` - -If `spin_or_pol` is of type [`AllPolarization`](@ref) or [`AllSpin`](@ref), the output is an `SVector` with both spin/polarization alignments: - -```julia -base_state(::Fermion, ::Incoming, mom) # -> SVector{2,BiSpinor} -base_state(::AntiFermion, ::Incoming, mom) # -> SVector{2,AdjointBiSpinor} -base_state(::Fermion, ::Outgoing, mom) # -> SVector{2,AdjointBiSpinor} -base_state(::AntiFermion, ::Outgoing, mom) # -> SVector{2,BiSpinor} -base_state(::Photon, ::Incoming, mom) # -> SVector{2,SLorentzVector{ComplexF64}} -base_state(::Photon, ::Outgoing, mom) # -> SVector{2,SLorentzVector{ComplexF64}} -``` - -# Example - -```julia - -using QEDbase - -mass = 1.0 # set electron mass to 1.0 -px,py,pz = rand(3) # generate random spatial components -E = sqrt(px^2 + py^2 + pz^2 + mass^2) # compute energy, i.e. the electron is on-shell -mom = SFourMomentum(E, px, py, pz) # initialize the four-momentum of the electron - -# compute the state of an incoming electron with spin = SpinUp -# note: base_state is not exported! -electron_state = base_state(Electron(), Incoming(), mom, SpinUp()) - -``` - -```jldoctest -julia> using QEDbase; using QEDcore - -julia> mass = 1.0; px,py,pz = (0.1, 0.2, 0.3); E = sqrt(px^2 + py^2 + pz^2 + mass^2); mom = SFourMomentum(E, px, py, pz) -4-element SFourMomentum with indices SOneTo(4): - 1.0677078252031311 - 0.1 - 0.2 - 0.3 - -julia> electron_state = base_state(Electron(), Incoming(), mom, SpinUp()) -4-element BiSpinor with indices SOneTo(4): - 1.4379526505428235 + 0.0im - 0.0 + 0.0im - -0.20862995724285552 + 0.0im - -0.06954331908095185 - 0.1390866381619037im - -julia> electron_states = base_state(Electron(), Incoming(), mom, AllSpin()) -2-element StaticArraysCore.SVector{2, BiSpinor} with indices SOneTo(2): - [1.4379526505428235 + 0.0im, 0.0 + 0.0im, -0.20862995724285552 + 0.0im, -0.06954331908095185 - 0.1390866381619037im] - [0.0 + 0.0im, 1.4379526505428235 + 0.0im, -0.06954331908095185 + 0.1390866381619037im, 0.20862995724285552 + 0.0im] -``` - -!!! note "Iterator convenience" - The returned objects of `base_state` can be consistently wrapped in an `SVector` for iteration using [`_as_svec`](@ref). - - This way, a loop like the following becomes possible when `spin` may be definite or indefinite. - ```julia - for state in QEDbase._as_svec(base_state(Electron(), Incoming(), momentum, spin)) - # ... - end - ``` - -!!! note "Conventions" - - For an incoming fermion with momentum ``p``, we use the explicit formula: - - ```math - u_\\sigma(p) = \\frac{\\gamma^\\mu p_\\mu + m}{\\sqrt{\\vert p_0\\vert + m}} \\eta_\\sigma, - ``` - - where the elementary base spinors are given as - - ```math - \\eta_1 = (1, 0, 0, 0)^T\\\\ - \\eta_2 = (0, 1, 0, 0)^T - ``` - - For an outgoing anti-fermion with momentum ``p``, we use the explicit formula: - - ```math - v_\\sigma(p) = \\frac{-\\gamma^\\mu p_\\mu + m}{\\sqrt{\\vert p_0\\vert + m}} \\chi_\\sigma, - ``` - - where the elementary base spinors are given as - - ```math - \\chi_1 = (0, 0, 1, 0)^T\\\\ - \\chi_2 = (0, 0, 0, 1)^T - ``` - - For outgoing fermions and incoming anti-fermions with momentum ``p``, the base state is given as the Dirac-adjoint of the respective incoming fermion or outgoing anti-fermion state: - - ```math - \\overline{u}_\\sigma(p) = u_\\sigma^\\dagger \\gamma^0\\\\ - \\overline{v}_\\sigma(p) = v_\\sigma^\\dagger \\gamma^0 - ``` - - where ``v_\\sigma`` is the base state of the respective outgoing anti-fermion. - - For a photon with four-momentum ``k^\\mu = \\omega (1, \\cos\\phi \\sin\\theta, \\sin\\phi \\sin\\theta, \\cos\\theta)``, the two polarization vectors are given as - - ```math - \\begin{align*} - \\epsilon^\\mu_1 &= (0, \\cos\\theta \\cos\\phi, \\cos\\theta \\sin\\phi, -\\sin\\theta),\\\\ - \\epsilon^\\mu_2 &= (0, -\\sin\\phi, \\cos\\phi, 0). - \\end{align*} - ``` - -!!! warning - In the current implementation there are **no checks** built-in, which verify the passed arguments whether they describe on-shell particles, i.e. `p*p≈mass^2`. Using `base_state` with off-shell particles will cause unpredictable behavior. -""" -function base_state end diff --git a/src/interfaces/particle_types.jl b/src/interfaces/particle_types.jl deleted file mode 100644 index 121327d..0000000 --- a/src/interfaces/particle_types.jl +++ /dev/null @@ -1,243 +0,0 @@ -############### -# The particle types -# -# In this file, we define the types of particles, their states and directions, and -# implement the abstact particle interface accordingly. -############### - -""" - AbstractParticleSpinor - -TBW -""" -abstract type AbstractParticleSpinor end - -""" - AbstractParticleType <: AbstractParticle - -This is the abstract base type for every species of particles. All functionalities defined on subtypes of `AbstractParticleType` should be static, i.e. known at compile time. -For adding runtime information, e.g. four-momenta or particle states, to a particle, consider implementing a concrete subtype of [`AbstractParticle`](@ref) instead, which may have a type parameter `P<:AbstractParticleType`. - -Concrete built-in subtypes of `AbstractParticleType` are [`Electron`](@ref), [`Positron`](@ref) and [`Photon`](@ref). -""" -abstract type AbstractParticleType <: AbstractParticle end - -""" -Abstract base types for particle species that act like fermions in the sense of particle statistics. - -!!! note "particle interface" - Every concrete subtype of [`FermionLike`](@ref) has `is_fermion(::FermionLike) = true`. -""" -abstract type FermionLike <: AbstractParticleType end - -is_fermion(::FermionLike) = true - -""" -Abstract base type for fermions as distinct from [`AntiFermion`](@ref)s. - -!!! note "particle interface" - All subtypes of `Fermion` have - ```julia - is_fermion(::Fermion) = true - is_particle(::Fermion) = true - is_anti_particle(::Fermion) = false - ``` - -""" -abstract type Fermion <: FermionLike end - -is_particle(::Fermion) = true - -is_anti_particle(::Fermion) = false - -""" -Abstract base type for anti-fermions as distinct from its particle counterpart `Fermion`. - -!!! note "particle interface" - All subtypes of `AntiFermion` have - ```julia - is_fermion(::AntiFermion) = true - is_particle(::AntiFermion) = false - is_anti_particle(::AntiFermion) = true - ``` - -""" -abstract type AntiFermion <: FermionLike end - -is_particle(::AntiFermion) = false - -is_anti_particle(::AntiFermion) = true - -""" -Abstract base type for majorana-fermions, i.e. fermions which are their own anti-particles. - -!!! note "particle interface" - All subtypes of `MajoranaFermion` have - ```julia - is_fermion(::MajoranaFermion) = true - is_particle(::MajoranaFermion) = true - is_anti_particle(::MajoranaFermion) = true - ``` - -""" -abstract type MajoranaFermion <: FermionLike end - -is_particle(::MajoranaFermion) = true - -is_anti_particle(::MajoranaFermion) = true - -""" -Concrete type for *electrons* as a particle species. Mostly used for dispatch. - -```jldoctest -julia> using QEDbase - -julia> Electron() -electron -``` - -!!! note "particle interface" - Besides being a subtype of [`Fermion`](@ref), objects of type `Electron` have - - ```julia - mass(::Electron) = 1.0 - charge(::Electron) = -1.0 - ``` -""" -struct Electron <: Fermion end -mass(::Electron) = 1.0 -charge(::Electron) = -1.0 -Base.show(io::IO, ::Electron) = print(io, "electron") - -""" -Concrete type for *positrons* as a particle species. Mostly used for dispatch. - -```jldoctest -julia> using QEDbase - -julia> Positron() -positron -``` - -!!! note "particle interface" - Besides being a subtype of [`AntiFermion`](@ref), objects of type `Positron` have - - ```julia - mass(::Positron) = 1.0 - charge(::Positron) = 1.0 - ``` - -""" -struct Positron <: AntiFermion end -mass(::Positron) = 1.0 -charge(::Positron) = 1.0 -Base.show(io::IO, ::Positron) = print(io, "positron") - -""" -Abstract base types for particle species that act like bosons in the sense of particle statistics. - -!!! note "particle interface" - Every concrete subtype of `BosonLike` has `is_boson(::BosonLike) = true`. -""" -abstract type BosonLike <: AbstractParticleType end - -is_boson(::BosonLike) = true - -""" -Abstract base type for bosons as distinct from its anti-particle counterpart [`AntiBoson`](@ref). - -!!! note "particle interface" - All subtypes of `Boson` have - ```julia - is_boson(::Boson) = true - is_particle(::Boson) = true - is_anti_particle(::Boson) = false - ``` - -""" -abstract type Boson <: BosonLike end -is_particle(::Boson) = true -is_anti_particle(::Boson) = false - -""" -Abstract base type for anti-bosons as distinct from its particle counterpart [`Boson`](@ref). - -!!! note "particle interface" - All subtypes of `AntiBoson` have - ```julia - is_boson(::AntiBoson) = true - is_particle(::AntiBoson) = false - is_anti_particle(::AntiBoson) = true - ``` - -""" -abstract type AntiBoson <: BosonLike end -is_particle(::AntiBoson) = false -is_anti_particle(::AntiBoson) = true - -""" -Abstract base type for majorana-bosons, i.e. bosons which are their own anti-particles. - -!!! note "particle interface" - All subtypes of `MajoranaBoson` have - ```julia - is_boson(::MajoranaBoson) = true - is_particle(::MajoranaBoson) = true - is_anti_particle(::MajoranaBoson) = true - ``` - -""" -abstract type MajoranaBoson <: BosonLike end -is_particle(::MajoranaBoson) = true -is_anti_particle(::MajoranaBoson) = true - -""" -Concrete type for the *photons* as a particle species. Mostly used for dispatch. - -```jldoctest -julia> using QEDbase - -julia> Photon() -photon -``` - -!!! note "particle interface" - Besides being a subtype of `MajoranaBoson`, `Photon` has - - ```julia - mass(::Photon) = 0.0 - charge(::Photon) = 0.0 - ``` - -""" -struct Photon <: MajoranaBoson end -mass(::Photon) = 0.0 -charge(::Photon) = 0.0 -Base.show(io::IO, ::Photon) = print(io, "photon") - -""" - propagator(particle::AbstractParticleType, mom::QEDbase.AbstractFourMomentum, [mass::Real]) - -Return the propagator of a particle for a given four-momentum. If `mass` is passed, the respective propagator for massive particles is used, if not, it is assumed the particle passed in is massless. - -!!! note "Convention" - - There are two types of implementations for propagators given in `QEDProcesses`: - For a `BosonLike` particle with four-momentum ``k`` and mass ``m``, the propagator is given as - - ```math - D(k) = \\frac{1}{k^2 - m^2}. - ``` - - For a `FermionLike` particle with four-momentum ``p`` and mass ``m``, the propagator is given as - - ```math - S(p) = \\frac{\\gamma^\\mu p_\\mu + mass}{p^2 - m^2}. - ``` - -!!! warning - - This function does not throw when the given particle is off-shell. If an off-shell particle is passed, the function `propagator` returns `Inf`. - -""" -function propagator end diff --git a/test/coverage/.gitkeep b/test/coverage/.gitkeep deleted file mode 100644 index e69de29..0000000 diff --git a/test/particle_properties.jl b/test/particle_properties.jl new file mode 100644 index 0000000..a28be08 --- /dev/null +++ b/test/particle_properties.jl @@ -0,0 +1,37 @@ +using QEDbase +using StaticArrays +using Random + +# test function to test scalar broadcasting +test_broadcast(x::AbstractParticle) = x +test_broadcast(x::ParticleDirection) = x +test_broadcast(x::AbstractSpinOrPolarization) = x + +@testset "scalar broadcasting" begin + @testset "directions" begin + @testset "$dir" for dir in (Incoming(), Outgoing()) + @test test_broadcast.(dir) == dir + end + end + + @testset "spins and polarization" begin + @testset "$spin_or_pol" for spin_or_pol in ( + SpinUp(), SpinDown(), AllSpin(), PolX(), PolY(), AllPol() + ) + @test test_broadcast.(spin_or_pol) == spin_or_pol + end + end +end + +@testset "multiplicity of spins or pols" begin + @testset "single" begin + @testset "$spin_or_pol" for spin_or_pol in (SpinUp(), SpinDown(), PolX(), PolY()) + @test multiplicity(spin_or_pol) == 1 + end + end + @testset "multiple" begin + @testset "$spin_or_pol" for spin_or_pol in (AllSpin(), AllPol()) + @test multiplicity(spin_or_pol) == 2 + end + end +end diff --git a/test/particles.jl b/test/particles.jl deleted file mode 100644 index 0370278..0000000 --- a/test/particles.jl +++ /dev/null @@ -1,102 +0,0 @@ -using QEDbase -using StaticArrays -using Random - -# test function to test scalar broadcasting -test_broadcast(x::AbstractParticle) = x -test_broadcast(x::ParticleDirection) = x -test_broadcast(x::AbstractSpinOrPolarization) = x - -@testset "scalar broadcasting" begin - @testset "directions" begin - @testset "$dir" for dir in (Incoming(), Outgoing()) - @test test_broadcast.(dir) == dir - end - end - - @testset "spins and polarization" begin - @testset "$spin_or_pol" for spin_or_pol in ( - SpinUp(), SpinDown(), AllSpin(), PolX(), PolY(), AllPol() - ) - @test test_broadcast.(spin_or_pol) == spin_or_pol - end - end -end - -@testset "multiplicity of spins or pols" begin - @testset "single" begin - @testset "$spin_or_pol" for spin_or_pol in (SpinUp(), SpinDown(), PolX(), PolY()) - @test multiplicity(spin_or_pol) == 1 - end - end - @testset "multiple" begin - @testset "$spin_or_pol" for spin_or_pol in (AllSpin(), AllPol()) - @test multiplicity(spin_or_pol) == 2 - end - end -end - -@testset "fermion likes" begin - @testset "antifermion" begin - struct TestAntiFermion <: AntiFermion end - @test is_fermion(TestAntiFermion()) - @test !is_particle(TestAntiFermion()) - @test is_anti_particle(TestAntiFermion()) - @test test_broadcast.(TestAntiFermion()) == TestAntiFermion() - end - - @testset "majorana fermion" begin - struct TestMajoranaFermion <: MajoranaFermion end - @test is_fermion(TestMajoranaFermion()) - @test is_particle(TestMajoranaFermion()) - @test is_anti_particle(TestMajoranaFermion()) - @test test_broadcast.(TestMajoranaFermion()) == TestMajoranaFermion() - end - - @testset "electron" begin - @test is_fermion(Electron()) - @test is_particle(Electron()) - @test !is_anti_particle(Electron()) - @test mass(Electron()) == 1.0 - @test charge(Electron()) == -1.0 - @test test_broadcast.(Electron()) == Electron() - end - - @testset "positron" begin - @test is_fermion(Positron()) - @test !is_particle(Positron()) - @test is_anti_particle(Positron()) - @test mass(Positron()) == 1.0 - @test charge(Positron()) == 1.0 - @test test_broadcast.(Positron()) == Positron() - end -end - -@testset "boson likes" begin - @testset "boson" begin - struct TestBoson <: Boson end - @test !is_fermion(TestBoson()) - @test is_boson(TestBoson()) - @test is_particle(TestBoson()) - @test !is_anti_particle(TestBoson()) - @test test_broadcast.(TestBoson()) == TestBoson() - end - - @testset "antiboson" begin - struct TestAntiBoson <: AntiBoson end - @test !is_fermion(TestAntiBoson()) - @test is_boson(TestAntiBoson()) - @test !is_particle(TestAntiBoson()) - @test is_anti_particle(TestAntiBoson()) - @test test_broadcast.(TestAntiBoson()) == TestAntiBoson() - end - - @testset "majorana boson" begin - struct TestMajoranaBoson <: MajoranaBoson end - @test !is_fermion(TestMajoranaBoson()) - @test is_boson(TestMajoranaBoson()) - @test is_particle(TestMajoranaBoson()) - @test is_anti_particle(TestMajoranaBoson()) - @test test_broadcast.(TestMajoranaBoson()) == TestMajoranaBoson() - end -end diff --git a/test/runtests.jl b/test/runtests.jl index 4e2da7f..8c6062c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,7 +3,7 @@ using Test using SafeTestsets begin - # # Interfaces + # Interfaces @time @safetestset "model interface" begin include("interfaces/model.jl") end @@ -13,6 +13,6 @@ begin end @time @safetestset "particles" begin - include("particles.jl") + include("particle_properties.jl") end end From 38c0c8b7b731288fd4141bf2013f45b71eec4634 Mon Sep 17 00:00:00 2001 From: Anton Reinhard Date: Wed, 19 Jun 2024 18:46:34 +0200 Subject: [PATCH 14/21] Fix spin_pol doc references --- src/particles/spin_pol.jl | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/src/particles/spin_pol.jl b/src/particles/spin_pol.jl index 06f6e2f..8b677cf 100644 --- a/src/particles/spin_pol.jl +++ b/src/particles/spin_pol.jl @@ -1,30 +1,30 @@ """ -Abstract base type for the spin or polarization of [`FermionLike`](@ref) or [`BosonLike`](@ref) particles, respectively. +Abstract base type for the spin or polarization of [`is_fermion`](@ref) or [`is_boson`](@ref) particles, respectively. """ abstract type AbstractSpinOrPolarization end Base.broadcastable(spin_or_pol::AbstractSpinOrPolarization) = Ref(spin_or_pol) """ -Abstract base type for the spin of [`FermionLike`](@ref) particles. +Abstract base type for the spin of [`is_fermion`](@ref) particles. """ abstract type AbstractSpin <: AbstractSpinOrPolarization end """ -Abstract base type for definite spins of [`FermionLike`](@ref) particles. +Abstract base type for definite spins of [`is_fermion`](@ref) particles. Concrete types are [`SpinUp`](@ref) and [`SpinDown`](@ref). """ abstract type AbstractDefiniteSpin <: AbstractSpin end """ -Abstract base type for indefinite spins of [`FermionLike`](@ref) particles. +Abstract base type for indefinite spins of [`is_fermion`](@ref) particles. One concrete type is [`AllSpin`](@ref). """ abstract type AbstractIndefiniteSpin <: AbstractSpin end """ -Concrete type indicating that a [`FermionLike`](@ref) has spin-up. +Concrete type indicating that an [`is_fermion`](@ref) particle has spin-up. ```jldoctest julia> using QEDbase @@ -37,7 +37,7 @@ struct SpinUp <: AbstractDefiniteSpin end Base.show(io::IO, ::SpinUp) = print(io, "spin up") """ -Concrete type indicating that a [`FermionLike`](@ref) has spin-down. +Concrete type indicating that an [`is_fermion`](@ref) particle has spin-down. ```jldoctest julia> using QEDbase @@ -50,7 +50,7 @@ struct SpinDown <: AbstractDefiniteSpin end Base.show(io::IO, ::SpinDown) = print(io, "spin down") """ -Concrete type indicating that a [`FermionLike`](@ref) has an indefinite spin and the differential cross section calculation should average or sum over all spins, depending on the direction ([`Incoming`](@ref) or [`Outgoing`](@ref)) of the particle in question. +Concrete type indicating that a particle with [`is_fermion`](@ref) has an indefinite spin and the differential cross section calculation should average or sum over all spins, depending on the direction ([`Incoming`](@ref) or [`Outgoing`](@ref)) of the particle in question. ```jldoctest julia> using QEDbase @@ -77,26 +77,26 @@ _spin_index(::SpinUp) = 1 _spin_index(::SpinDown) = 2 """ -Abstract base type for the polarization of [`BosonLike`](@ref) particles. +Abstract base type for the polarization of [`is_boson`](@ref) particles. """ abstract type AbstractPolarization <: AbstractSpinOrPolarization end """ -Abstract base type for definite polarizations of [`BosonLike`](@ref) particles. +Abstract base type for definite polarizations of [`is_boson`](@ref) particles. Concrete types are [`PolarizationX`](@ref) and [`PolarizationY`](@ref). """ abstract type AbstractDefinitePolarization <: AbstractPolarization end """ -Abstract base type for indefinite polarizations of [`BosonLike`](@ref) particles. +Abstract base type for indefinite polarizations of [`is_boson`](@ref) particles. One concrete type is [`AllPolarization`](@ref). """ abstract type AbstractIndefinitePolarization <: AbstractPolarization end """ -Concrete type indicating that a [`BosonLike`](@ref) has an indefinite polarization and the differential cross section calculation should average or sum over all polarizations, depending on the direction ([`Incoming`](@ref) or [`Outgoing`](@ref)) of the particle in question. +Concrete type indicating that a particle with [`is_boson`](@ref) has an indefinite polarization and the differential cross section calculation should average or sum over all polarizations, depending on the direction ([`Incoming`](@ref) or [`Outgoing`](@ref)) of the particle in question. ```jldoctest julia> using QEDbase @@ -121,7 +121,7 @@ const AllPol = AllPolarization Base.show(io::IO, ::AllPol) = print(io, "all polarizations") """ -Concrete type which indicates, that a [`BosonLike`](@ref) has polarization in ``x``-direction. +Concrete type which indicates, that a particle with [`is_boson`](@ref) has polarization in ``x``-direction. ```jldoctest julia> using QEDbase @@ -133,7 +133,7 @@ x-polarized !!! note "Coordinate axes" The notion of axes, e.g. ``x``- and ``y``-direction is just to distinguish two orthogonal polarization directions. - However, if the three-momentum of the [`BosonLike`](@ref) is aligned to the ``z``-axis of a coordinate system, the polarization axes define the ``x``- or ``y``-axis, respectively. + However, if the three-momentum of the [`is_boson`](@ref) particle is aligned to the ``z``-axis of a coordinate system, the polarization axes define the ``x``- or ``y``-axis, respectively. !!! info "Alias" @@ -150,7 +150,7 @@ const PolX = PolarizationX Base.show(io::IO, ::PolX) = print(io, "x-polarized") """ -Concrete type which indicates, that a [`BosonLike`](@ref) has polarization in ``y``-direction. +Concrete type which indicates, that an [`is_boson`](@ref) particle has polarization in ``y``-direction. ```jldoctest julia> using QEDbase @@ -162,7 +162,7 @@ y-polarized !!! note "Coordinate axes" The notion of axes, e.g. ``x``- and ``y``-direction is just to distinguish two orthogonal polarization directions. - However, if the three-momentum of the [`BosonLike`](@ref) is aligned to the ``z``-axis of a coordinate system, the polarization axes define the ``x``- or ``y``-axis, respectively. + However, if the three-momentum of the [`is_boson`](@ref) particle is aligned to the ``z``-axis of a coordinate system, the polarization axes define the ``x``- or ``y``-axis, respectively. !!! info "Alias" From 8eac3199e8634190ef2f3f9beda6da2e4a05f303 Mon Sep 17 00:00:00 2001 From: Anton Reinhard Date: Wed, 19 Jun 2024 20:41:18 +0200 Subject: [PATCH 15/21] Remove unnecessary gitkeep --- .formatting/.gitkeep | 0 1 file changed, 0 insertions(+), 0 deletions(-) delete mode 100644 .formatting/.gitkeep diff --git a/.formatting/.gitkeep b/.formatting/.gitkeep deleted file mode 100644 index e69de29..0000000 From 0f652647780dc3923dd6a6809b5167a3cf102cc8 Mon Sep 17 00:00:00 2001 From: Anton Reinhard Date: Wed, 19 Jun 2024 20:49:06 +0200 Subject: [PATCH 16/21] Fix doctests --- docs/make.jl | 2 ++ src/interfaces/particle.jl | 3 ++- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/docs/make.jl b/docs/make.jl index 4e0e719..b663201 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -11,6 +11,8 @@ using QEDbase #DocMeta.setdocmeta!(QEDbase, :DocTestSetup, :(using QEDbase); recursive=true) +Pkg.add(; url="https://github.com/QEDjl-project/QEDcore.jl", rev="dev") + using DocumenterCitations bib = CitationBibliography(joinpath(@__DIR__, "Bibliography.bib")) diff --git a/src/interfaces/particle.jl b/src/interfaces/particle.jl index b18a8b8..834609c 100644 --- a/src/interfaces/particle.jl +++ b/src/interfaces/particle.jl @@ -114,6 +114,7 @@ Return the propagator of a particle for a given four-momentum. If `mass` is pass """ function propagator end +# TODO: Turn the doctest on again when QEDcore has QEDprocesses functionality """ ```julia base_state( @@ -173,7 +174,7 @@ mom = SFourMomentum(E, px, py, pz) # initialize the four-momentum of the el electron_state = base_state(QEDcore.Electron(), Incoming(), mom, SpinUp()) ``` -```jldoctest +```Julia julia> using QEDbase; using QEDcore julia> mass = 1.0; px,py,pz = (0.1, 0.2, 0.3); E = sqrt(px^2 + py^2 + pz^2 + mass^2); mom = SFourMomentum(E, px, py, pz) From 2704d70461c4e85c6ba9a5c7b36484d98548ad92 Mon Sep 17 00:00:00 2001 From: Anton Reinhard Date: Thu, 20 Jun 2024 15:52:11 +0200 Subject: [PATCH 17/21] Fix some of the docs and add rudimentary AbstractPhaseSpacePoint and AbstracParticleStateful --- docs/src/index.md | 4 +++ docs/src/library/api.md | 46 +---------------------------- src/QEDbase.jl | 2 ++ src/interfaces/four_momentum.jl | 1 - src/interfaces/particle_stateful.jl | 4 +++ src/interfaces/phase_space_point.jl | 31 +++++++++++++++++++ src/interfaces/process.jl | 15 ++++++++-- 7 files changed, 54 insertions(+), 49 deletions(-) create mode 100644 src/interfaces/particle_stateful.jl create mode 100644 src/interfaces/phase_space_point.jl diff --git a/docs/src/index.md b/docs/src/index.md index 65e77c5..13d7c86 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -5,3 +5,7 @@ CurrentModule = QEDbase # QEDbase Documentation for [QEDbase](https://gitlab.hzdr.de/QEDjl/QEDbase.jl). This library provides (more or less) simple implementations of Lorentz vectors, Dirac spinors and Dirac matrices. These are usually used in order to do calculations in particle physics, e.g. they are part of the set of Feynman rules for a given quantum field theory (see [Peskin:1995ev](@cite)). + +```@autodocs +Modules = [QEDbase] +``` diff --git a/docs/src/library/api.md b/docs/src/library/api.md index f00ed57..cc3ffff 100644 --- a/docs/src/library/api.md +++ b/docs/src/library/api.md @@ -1,47 +1,3 @@ # QEDbase -```@autodocs -Modules = [QEDbase] -Pages = ["QEDbase.jl"] -Order = [:module] -``` - -# Lorentz Vector Interface - -```@autodocs -Modules = [QEDbase] -Pages = ["lorentz_interface.jl"] -Order = [:type,:function] -``` - -# Lorentz Vectors - -```@autodocs -Modules = [QEDbase] -Pages = ["lorentz_vector.jl"] -Order = [:type,:function] -``` - -# Four Momentum - -```@autodocs -Modules = [QEDbase] -Pages = ["four_momentum.jl"] -Order = [:type,:function] -``` - -# Dirac Tensors - -```@autodocs -Modules = [QEDbase] -Pages = ["dirac_tensors.jl"] -Order = [:type,:function] -``` - -# Particles - -```@autodocs -Modules = [QEDbase] -Pages = ["particles/particle_direction.jl", "particles/particle_spin_pol.jl", "particles/particle_spinors.jl", "particles/particle_states.jl", "particles/particle_types.jl", "interfaces/particle_interface.jl"] -Order = [:type,:function] -``` +:warning: Under construction diff --git a/src/QEDbase.jl b/src/QEDbase.jl index 6c41a00..49bbaaa 100644 --- a/src/QEDbase.jl +++ b/src/QEDbase.jl @@ -108,6 +108,8 @@ include("particles/spin_pol.jl") include("interfaces/phase_space.jl") include("interfaces/process.jl") +include("interfaces/particle_stateful.jl") +include("interfaces/phase_space_point.jl") include("interfaces/setup.jl") end #QEDbase diff --git a/src/interfaces/four_momentum.jl b/src/interfaces/four_momentum.jl index 733baac..adffa1b 100644 --- a/src/interfaces/four_momentum.jl +++ b/src/interfaces/four_momentum.jl @@ -5,7 +5,6 @@ Abstract base type for four-momentas, representing one energy and three spacial Also see: [`SFourMomentum`](@ref) """ - abstract type AbstractFourMomentum <: AbstractLorentzVector{Float64} end function Base.getproperty(P::TM, sym::Symbol) where {TM<:AbstractFourMomentum} diff --git a/src/interfaces/particle_stateful.jl b/src/interfaces/particle_stateful.jl new file mode 100644 index 0000000..2938c8c --- /dev/null +++ b/src/interfaces/particle_stateful.jl @@ -0,0 +1,4 @@ + +abstract type AbstractParticleStateful{ + DIR<:ParticleDirection,SPECIES<:AbstractParticleType,ELEMENT<:AbstractFourMomentum +} <: AbstractParticle end diff --git a/src/interfaces/phase_space_point.jl b/src/interfaces/phase_space_point.jl new file mode 100644 index 0000000..79171fb --- /dev/null +++ b/src/interfaces/phase_space_point.jl @@ -0,0 +1,31 @@ + +abstract type AbstractPhaseSpacePoint{ + PROC<:AbstractProcessDefinition, + MODEL<:AbstractModelDefinition, + PSDEF<:AbstractPhasespaceDefinition, + IN_PARTICLES<:Tuple{Vararg{AbstractParticleStateful}}, + OUT_PARTICLES<:Tuple{Vararg{AbstractParticleStateful}}, + ELEMENT<:AbstractFourMomentum, +} end + +""" + AbstractInPhaseSpacePoint + +A partial type specialization on [`AbstractPhaseSpacePoint`](@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: [`AbstractOutPhaseSpacePoint`](@ref) +""" +AbstractInPhaseSpacePoint{P,M,D,IN,OUT,E} = AbstractPhaseSpacePoint{ + P,M,D,IN,OUT,E +} where {PS<:AbstractParticleStateful,IN<:Tuple{PS,Vararg},OUT<:Tuple{Vararg}} + +""" + AbstractOutPhaseSpacePoint + +A partial type specialization on [`AbstractPhaseSpacePoint`](@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: [`AbstractInPhaseSpacePoint`](@ref) +""" +AbstractOutPhaseSpacePoint{P,M,D,IN,OUT,E} = AbstractPhaseSpacePoint{ + P,M,D,IN,OUT,E +} where {PS<:AbstractParticleStateful,IN<:Tuple{Vararg},OUT<:Tuple{PS,Vararg}} diff --git a/src/interfaces/process.jl b/src/interfaces/process.jl index 469a1cd..d60bca0 100644 --- a/src/interfaces/process.jl +++ b/src/interfaces/process.jl @@ -70,8 +70,7 @@ function outgoing_particles end MODEL <: AbstractModelDefinition, } -Interface function which returns the incident flux of the given scattering process for a given [`InPhaseSpacePoint`](@ref). - +Interface function which returns the incident flux of the given scattering process for a given `InPhaseSpacePoint`. """ function _incident_flux end @@ -148,7 +147,7 @@ function out_phase_space_dimension end } 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 [`InPhaseSpacePoint`](@ref). +given process and model for a passed `InPhaseSpacePoint`. !!! note "total cross section" @@ -157,6 +156,16 @@ given process and model for a passed [`InPhaseSpacePoint`](@ref). """ function _total_probability end +""" + total_cross_section(in_psp::InPhaseSpacePoint) + +Return the total cross section for a given `InPhaseSpacePoint`. +""" +#function total_cross_section(in_psp::AbstractInPhaseSpacePoint) +# I = 1 / (4 * _incident_flux(in_psp)) +# return I * _total_probability(in_psp) +#end + ####################### # # utility functions From 6f35b7eed03b464dc2d2216d6918da733283dcc2 Mon Sep 17 00:00:00 2001 From: Anton Reinhard Date: Fri, 21 Jun 2024 11:53:35 +0200 Subject: [PATCH 18/21] Fix doc problems and complete PSP and PS interfaces and descriptions --- README.md | 4 +- src/QEDbase.jl | 10 ++++- src/interfaces/four_momentum.jl | 2 +- src/interfaces/particle_stateful.jl | 48 +++++++++++++++++++++ src/interfaces/phase_space_point.jl | 65 +++++++++++++++++++++++++++++ src/interfaces/process.jl | 10 ----- src/particles/direction.jl | 14 +++++++ src/total_cross_section.jl | 9 ++++ 8 files changed, 148 insertions(+), 14 deletions(-) create mode 100644 src/total_cross_section.jl diff --git a/README.md b/README.md index a2f4cc0..e310949 100644 --- a/README.md +++ b/README.md @@ -47,7 +47,7 @@ To install the locally downloaded package on Windows, change to the parent direc One can define a static four momentum component wise: ```julia -julia> using QEDbase +julia> using QEDbase; using QEDcore juila> mass = rand()*10 @@ -73,7 +73,7 @@ julia> @assert isapprox(getPlus(mom), 0.5*(E+pz)) julia> @assert isapprox(getPerp(mom), px^2 + py^2) ``` -and a lot more (see [here](www.docs-to-the-lorentz-interface-getter.jl) for a complete list). There is also a mutable version of a four vector in `QEDbase.jl`, where the Lorentz-vector interface provides setters to different properties as well (see [here](www.docs-to-the-lorentz-interface-setter.jl) for details). +and a lot more (see [here](www.docs-to-the-lorentz-interface-getter.jl) for a complete list). There is also a mutable version of a four vector in `QEDcore.jl`, where the Lorentz-vector interface provides setters to different properties as well (see [here](www.docs-to-the-lorentz-interface-setter.jl) for details). ## Testing diff --git a/src/QEDbase.jl b/src/QEDbase.jl index 49bbaaa..ed81be0 100644 --- a/src/QEDbase.jl +++ b/src/QEDbase.jl @@ -71,9 +71,15 @@ export particles, number_particles export AbstractComputationSetup, compute export AbstractProcessSetup, scattering_process, physical_model -# Abstract phase space interface +# Abstract phase space definition interface export AbstractCoordinateSystem, AbstractFrameOfReference, AbstractPhasespaceDefinition +# Abstract phase space point interface +export AbstractParticleStateful, AbstractPhaseSpacePoint +export particle_direction, particle_species, momentum +export process, model, phase_space_definition, momenta +export AbstractInPhaseSpacePoint, AbstractOutPhaseSpacePoint + # errors export InvalidInputError, RegistryError, OnshellError, SpinorConstructionError @@ -112,4 +118,6 @@ include("interfaces/particle_stateful.jl") include("interfaces/phase_space_point.jl") include("interfaces/setup.jl") +include("total_cross_section.jl") + end #QEDbase diff --git a/src/interfaces/four_momentum.jl b/src/interfaces/four_momentum.jl index adffa1b..f79de00 100644 --- a/src/interfaces/four_momentum.jl +++ b/src/interfaces/four_momentum.jl @@ -3,7 +3,7 @@ Abstract base type for four-momentas, representing one energy and three spacial components. -Also see: [`SFourMomentum`](@ref) +Also see: `QEDcore.SFourMomentum`, `QEDcore.MFourMomentum` """ abstract type AbstractFourMomentum <: AbstractLorentzVector{Float64} end diff --git a/src/interfaces/particle_stateful.jl b/src/interfaces/particle_stateful.jl index 2938c8c..9337c4e 100644 --- a/src/interfaces/particle_stateful.jl +++ b/src/interfaces/particle_stateful.jl @@ -1,4 +1,52 @@ +""" + AbstractParticleStateful <: QEDbase.AbstractParticle + +Abstract base type for the representation of a particle with a state. It requires the following interface functions to be provided: + +- [`particle_direction`](@ref): Returning the particle's direction. +- [`particle_species`](@ref): Returning the particle's species. +- [`momentum`](@ref): Returning the particle's momentum. + +Implementations for [`is_fermion`](@ref), [`is_boson`](@ref), [`is_particle`](@ref), [`is_anti_particle`](@ref), [`is_incoming`](@ref), [`is_outgoing`](@ref), [`mass`](@ref), and [`charge`](@ref) are automatically provided using the interface functions above to fulfill the [`QEDbase.AbstractParticle`](@ref) interface. + +""" abstract type AbstractParticleStateful{ DIR<:ParticleDirection,SPECIES<:AbstractParticleType,ELEMENT<:AbstractFourMomentum } <: AbstractParticle end + +""" + particle_direction(part::AbstractParticleStateful) + +Interface function that must return the particle's [`ParticleDirection`](@ref). +""" +function particle_direction end + +""" + particle_species(part::AbstractParticleStateful) + +Interface function that must return the particle's [`AbstractParticleType`](@ref), e.g. `QEDcore.Electron()`. +""" +function particle_species end + +""" + momentum(part::AbstractParticleStateful) + +Interface function that must return the particle's [`AbstractFourMomentum`](@ref). +""" +function momentum end + +# implement particle interface +@inline is_incoming(particle::AbstractParticleStateful) = + is_incoming(particle_direction(particle)) +@inline is_outgoing(particle::AbstractParticleStateful) = + is_outgoing(particle_direction(particle)) +@inline is_fermion(particle::AbstractParticleStateful) = + is_fermion(particle_species(particle)) +@inline is_boson(particle::AbstractParticleStateful) = is_boson(particle_species(particle)) +@inline is_particle(particle::AbstractParticleStateful) = + is_particle(particle_species(particle)) +@inline is_anti_particle(particle::AbstractParticleStateful) = + is_anti_particle(particle_species(particle)) +@inline mass(particle::AbstractParticleStateful) = mass(particle_species(particle)) +@inline charge(particle::AbstractParticleStateful) = charge(particle_species(particle)) diff --git a/src/interfaces/phase_space_point.jl b/src/interfaces/phase_space_point.jl index 79171fb..11b4a59 100644 --- a/src/interfaces/phase_space_point.jl +++ b/src/interfaces/phase_space_point.jl @@ -1,4 +1,31 @@ +""" + AbstractPhaseSpacePoint{PROC, MODEL, PSDEF, IN_PARTICLES, OUT_PARTICLES, ELEMENT} + +Representation of a point in the phase space of a process. It has several template arguments: +- `PROC <: `[`AbstractProcessDefinition`](@ref) +- `MODEL <: `[`AbstractModelDefinition`](@ref) +- `PSDEF <: `[`AbstractPhasespaceDefinition`](@ref) +- `IN_PARTICLES <: `Tuple{Vararg{AbstractParticleStateful}}`: The tuple type of all the incoming [`AbstractParticleStateful`](@ref)s. +- `OUT_PARTICLES <: `Tuple{Vararg{AbstractParticleStateful}}`: The tuple type of all the outgoing [`AbstractParticleStateful`](@ref)s. +- `ELEMENT <: `[`AbstractFourMomentum`](@ref): The type of the [`AbstractParticleStateful`](@ref)s and also the return type of the [`momentum`](@ref) interface function. + +The following interface functions must be provided: +- `Base.getindex(psp::AbstractPhaseSpacePoint, dir::ParticleDirection, n::Int)`: Return the nth [`AbstractParticleStateful`](@ref) of the given direction. Throw `BoundsError` for invalid indices. +- [`process`](@ref): Return the process. +- [`model`](@ref): Return the model. +- [`phase_space_definition`](@ref): Return the phase space definition. + +From this, the following functions are automatically derived: +- `momentum(psp::AbstractPhaseSpacePoint, dir::ParticleDirection, n::Int)::ELEMENT`: Return the momentum of the nth [`AbstractParticleStateful`](@ref) of the given direction. +- `momenta(psp::PhaseSpacePoint, ::ParticleDirection)::ELEMENT`: Return a `Tuple` of all the momenta of the given direction. +Furthermore, an implementation of an [`AbstractPhaseSpacePoint`](@ref) has to verify on construction that it is valid, i.e., the following conditions are fulfilled: +- `IN_PARTICLES` must match `incoming_particles(::PROC)` in length, order, and type **or** be an empty `Tuple`. +- `OUT_PARTICLES` must match the `outgoing_particles(::PROC)` in length, order, and type, **or** be an empty `Tuple`. +- `IN_PARTICLES` and `OUT_PARTICLES` may not both be empty. + +If `IN_PARTICLES` is non-empty, `AbstractPhaseSpacePoint <: AbstractInPhaseSpacePoint` is true. Likewise, if `OUT_PARTICLES` is non-empty, `AbstractPhaseSpacePoint <: AbstractOutPhaseSpacePoint` is true. Consequently, if both `IN_PARTICLES` and `OUT_PARTICLES` are non-empty, both `<:` statements are true. +""" abstract type AbstractPhaseSpacePoint{ PROC<:AbstractProcessDefinition, MODEL<:AbstractModelDefinition, @@ -8,6 +35,44 @@ abstract type AbstractPhaseSpacePoint{ ELEMENT<:AbstractFourMomentum, } end +""" + process(psp::AbstractPhaseSpacePoint) + +Return the phase space point's set process. +""" +function process end + +""" + model(psp::AbstractPhaseSpacePoint) + +Return the phase space point's set model. +""" +function model end + +""" + phase_space_definition(psp::AbstractPhaseSpacePoint) + +Return the phase space point's set phase space definition. +""" +function phase_space_definition end + +""" + momentum(psp::AbstractPhaseSpacePoint, dir::ParticleDirection, n::Int) + +Returns the momentum of the `n`th particle in the given [`AbstractPhaseSpacePoint`](@ref) which has direction `dir`. If `n` is outside the valid range for this phase space point, a `BoundsError` is thrown. +""" +function momentum(psp::AbstractPhaseSpacePoint, dir::ParticleDirection, n::Int) + return psp[dir, n].mom +end + +""" + momenta(psp::AbstractPhaseSpacePoint, ::ParticleDirection) + +Return a `Tuple` of all the particles' momenta for the given `ParticleDirection`. +""" +momenta(psp::AbstractPhaseSpacePoint, ::QEDbase.Incoming) = momentum.(psp.in_particles) +momenta(psp::AbstractPhaseSpacePoint, ::QEDbase.Outgoing) = momentum.(psp.out_particles) + """ AbstractInPhaseSpacePoint diff --git a/src/interfaces/process.jl b/src/interfaces/process.jl index d60bca0..5573037 100644 --- a/src/interfaces/process.jl +++ b/src/interfaces/process.jl @@ -156,16 +156,6 @@ given process and model for a passed `InPhaseSpacePoint`. """ function _total_probability end -""" - total_cross_section(in_psp::InPhaseSpacePoint) - -Return the total cross section for a given `InPhaseSpacePoint`. -""" -#function total_cross_section(in_psp::AbstractInPhaseSpacePoint) -# I = 1 / (4 * _incident_flux(in_psp)) -# return I * _total_probability(in_psp) -#end - ####################### # # utility functions diff --git a/src/particles/direction.jl b/src/particles/direction.jl index 9064007..45444b4 100644 --- a/src/particles/direction.jl +++ b/src/particles/direction.jl @@ -4,6 +4,20 @@ Abstract base type for the directions of particles in the context of processes, abstract type ParticleDirection end Base.broadcastable(dir::ParticleDirection) = Ref(dir) +""" + is_incoming(dir::ParticleDirection) + +Convenience function that returns true for [`Incoming`](@ref) and false otherwise. +""" +function is_incoming end + +""" + is_outgoing(dir::ParticleDirection) + +Convenience function that returns true for [`Outgoing`](@ref) and alse otherwise. +""" +function is_outgoing end + """ Incoming <: ParticleDirection diff --git a/src/total_cross_section.jl b/src/total_cross_section.jl new file mode 100644 index 0000000..84f55db --- /dev/null +++ b/src/total_cross_section.jl @@ -0,0 +1,9 @@ +""" + total_cross_section(in_psp::AbstractInPhaseSpacePoint) + +Return the total cross section for a given [`AbstractInPhaseSpacePoint`](@ref). +""" +function total_cross_section(in_psp::AbstractInPhaseSpacePoint) + I = 1 / (4 * _incident_flux(in_psp)) + return I * _total_probability(in_psp) +end From 4a30308596177430dde262de3eef596dc22d1ea5 Mon Sep 17 00:00:00 2001 From: Anton Reinhard Date: Sat, 22 Jun 2024 00:08:46 +0200 Subject: [PATCH 19/21] Apply suggestions from code review Co-authored-by: Uwe Hernandez Acosta --- docs/make.jl | 1 + src/interfaces/process.jl | 2 +- src/particles/spin_pol.jl | 24 ++++++++++++------------ 3 files changed, 14 insertions(+), 13 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index b663201..afc3b91 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -11,6 +11,7 @@ using QEDbase #DocMeta.setdocmeta!(QEDbase, :DocTestSetup, :(using QEDbase); recursive=true) +# TODO: remove before release Pkg.add(; url="https://github.com/QEDjl-project/QEDcore.jl", rev="dev") using DocumenterCitations diff --git a/src/interfaces/process.jl b/src/interfaces/process.jl index 5573037..f783de5 100644 --- a/src/interfaces/process.jl +++ b/src/interfaces/process.jl @@ -45,7 +45,7 @@ to enable the calculation of total probabilities and cross sections. """ abstract type AbstractProcessDefinition end -# broadcast every model as a scalar +# broadcast every process as a scalar Broadcast.broadcastable(proc::AbstractProcessDefinition) = Ref(proc) """ diff --git a/src/particles/spin_pol.jl b/src/particles/spin_pol.jl index 8b677cf..9695405 100644 --- a/src/particles/spin_pol.jl +++ b/src/particles/spin_pol.jl @@ -1,30 +1,30 @@ """ -Abstract base type for the spin or polarization of [`is_fermion`](@ref) or [`is_boson`](@ref) particles, respectively. +Abstract base type for the spin or polarization of particles with [`is_fermion`](@ref) or [`is_boson`](@ref), respectively. """ abstract type AbstractSpinOrPolarization end Base.broadcastable(spin_or_pol::AbstractSpinOrPolarization) = Ref(spin_or_pol) """ -Abstract base type for the spin of [`is_fermion`](@ref) particles. +Abstract base type for the spin of particles with [`is_fermion`](@ref). """ abstract type AbstractSpin <: AbstractSpinOrPolarization end """ -Abstract base type for definite spins of [`is_fermion`](@ref) particles. +Abstract base type for definite spins of particles with [`is_fermion`](@ref). Concrete types are [`SpinUp`](@ref) and [`SpinDown`](@ref). """ abstract type AbstractDefiniteSpin <: AbstractSpin end """ -Abstract base type for indefinite spins of [`is_fermion`](@ref) particles. +Abstract base type for indefinite spins of particles with [`is_fermion`](@ref). One concrete type is [`AllSpin`](@ref). """ abstract type AbstractIndefiniteSpin <: AbstractSpin end """ -Concrete type indicating that an [`is_fermion`](@ref) particle has spin-up. +Concrete type indicating that a particle with [`is_fermion`](@ref) has spin-up. ```jldoctest julia> using QEDbase @@ -37,7 +37,7 @@ struct SpinUp <: AbstractDefiniteSpin end Base.show(io::IO, ::SpinUp) = print(io, "spin up") """ -Concrete type indicating that an [`is_fermion`](@ref) particle has spin-down. +Concrete type indicating that a particle with [`is_fermion`](@ref) has spin-down. ```jldoctest julia> using QEDbase @@ -77,19 +77,19 @@ _spin_index(::SpinUp) = 1 _spin_index(::SpinDown) = 2 """ -Abstract base type for the polarization of [`is_boson`](@ref) particles. +Abstract base type for the polarization of particles with [`is_boson`](@ref). """ abstract type AbstractPolarization <: AbstractSpinOrPolarization end """ -Abstract base type for definite polarizations of [`is_boson`](@ref) particles. +Abstract base type for definite polarization of particles with [`is_boson`](@ref). Concrete types are [`PolarizationX`](@ref) and [`PolarizationY`](@ref). """ abstract type AbstractDefinitePolarization <: AbstractPolarization end """ -Abstract base type for indefinite polarizations of [`is_boson`](@ref) particles. +Abstract base type for indefinite polarization of particles with [`is_boson`](@ref). One concrete type is [`AllPolarization`](@ref). """ @@ -133,7 +133,7 @@ x-polarized !!! note "Coordinate axes" The notion of axes, e.g. ``x``- and ``y``-direction is just to distinguish two orthogonal polarization directions. - However, if the three-momentum of the [`is_boson`](@ref) particle is aligned to the ``z``-axis of a coordinate system, the polarization axes define the ``x``- or ``y``-axis, respectively. + However, if the three-momentum of the particle with [`is_boson`](@ref) is aligned to the ``z``-axis of a coordinate system, the polarization axes define the ``x``- or ``y``-axis, respectively. !!! info "Alias" @@ -150,7 +150,7 @@ const PolX = PolarizationX Base.show(io::IO, ::PolX) = print(io, "x-polarized") """ -Concrete type which indicates, that an [`is_boson`](@ref) particle has polarization in ``y``-direction. +Concrete type which indicates, that a particle with [`is_boson`](@ref) has polarization in ``y``-direction. ```jldoctest julia> using QEDbase @@ -162,7 +162,7 @@ y-polarized !!! note "Coordinate axes" The notion of axes, e.g. ``x``- and ``y``-direction is just to distinguish two orthogonal polarization directions. - However, if the three-momentum of the [`is_boson`](@ref) particle is aligned to the ``z``-axis of a coordinate system, the polarization axes define the ``x``- or ``y``-axis, respectively. + However, if the three-momentum of the particle with [`is_boson`](@ref) is aligned to the ``z``-axis of a coordinate system, the polarization axes define the ``x``- or ``y``-axis, respectively. !!! info "Alias" From 4427de191b4898082e053f3aab3a364963c3a0b8 Mon Sep 17 00:00:00 2001 From: AntonReinhard Date: Sat, 22 Jun 2024 00:35:04 +0200 Subject: [PATCH 20/21] Review fixes --- src/QEDbase.jl | 3 --- src/interfaces/lorentz_vectors/arithmetic.jl | 13 +------------ src/interfaces/lorentz_vectors/dirac_interaction.jl | 2 ++ src/interfaces/particle.jl | 3 ++- src/particles/direction.jl | 6 ++++-- 5 files changed, 9 insertions(+), 18 deletions(-) diff --git a/src/QEDbase.jl b/src/QEDbase.jl index ed81be0..07a7803 100644 --- a/src/QEDbase.jl +++ b/src/QEDbase.jl @@ -1,8 +1,5 @@ module QEDbase -import Base: *, / -import StaticArrays: similar_type - export minkowski_dot, mdot, register_LorentzVectorLike export getT, getX, getY, getZ export getMagnitude2, getMag2, getMagnitude, getMag diff --git a/src/interfaces/lorentz_vectors/arithmetic.jl b/src/interfaces/lorentz_vectors/arithmetic.jl index f7ef976..6f6e056 100644 --- a/src/interfaces/lorentz_vectors/arithmetic.jl +++ b/src/interfaces/lorentz_vectors/arithmetic.jl @@ -1,15 +1,4 @@ -# TODO: these overloads shouldn't be necessary, but are for some reason -@inline function *(L::TL, NUM::Number) where {TL<:AbstractLorentzVector} - return constructorof(TL)(L[1] * NUM, L[2] * NUM, L[3] * NUM, L[4] * NUM) -end - -@inline function *(NUM::Number, L::TL) where {TL<:AbstractLorentzVector} - return constructorof(TL)(L[1] * NUM, L[2] * NUM, L[3] * NUM, L[4] * NUM) -end - -@inline function /(L::TL, NUM::Number) where {TL<:AbstractLorentzVector} - return constructorof(TL)(L[1] / NUM, L[2] / NUM, L[3] / NUM, L[4] / NUM) -end +import Base: * function dot(p1::T1, p2::T2) where {T1<:AbstractLorentzVector,T2<:AbstractLorentzVector} return mdot(p1, p2) diff --git a/src/interfaces/lorentz_vectors/dirac_interaction.jl b/src/interfaces/lorentz_vectors/dirac_interaction.jl index 530d8ff..990d949 100644 --- a/src/interfaces/lorentz_vectors/dirac_interaction.jl +++ b/src/interfaces/lorentz_vectors/dirac_interaction.jl @@ -1,3 +1,5 @@ +import Base: * + """ $(TYPEDSIGNATURES) diff --git a/src/interfaces/particle.jl b/src/interfaces/particle.jl index 834609c..f71d693 100644 --- a/src/interfaces/particle.jl +++ b/src/interfaces/particle.jl @@ -129,7 +129,7 @@ Return the base state of a directed on-shell particle with a given four-momentum # Input -- `particle` -- the type of the particle, e.g. `QEDcore.Electron`, `QEDcore.Positron`, or `QEDcore.Photon`. +- `particle` -- the type of the particle, i.e., an instance of an [`AbstractParticleType`](@ref), e.g. `QEDcore.Electron`, `QEDcore.Positron`, or `QEDcore.Photon`. - `direction` -- the direction of the particle, i.e. [`Incoming`](@ref) or [`Outgoing`](@ref). - `momentum` -- the four-momentum of the particle - `[spin_or_pol]` -- if given, the spin or polarization of the particle, e.g. [`SpinUp`](@ref)/[`SpinDown`](@ref) or [`PolarizationX`](@ref)/[`PolarizationY`](@ref). @@ -174,6 +174,7 @@ mom = SFourMomentum(E, px, py, pz) # initialize the four-momentum of the el electron_state = base_state(QEDcore.Electron(), Incoming(), mom, SpinUp()) ``` +TODO: Reenable doctests ```Julia julia> using QEDbase; using QEDcore diff --git a/src/particles/direction.jl b/src/particles/direction.jl index 45444b4..66ac90f 100644 --- a/src/particles/direction.jl +++ b/src/particles/direction.jl @@ -6,15 +6,17 @@ Base.broadcastable(dir::ParticleDirection) = Ref(dir) """ is_incoming(dir::ParticleDirection) + is_incoming(particle::AbstractParticleStateful) -Convenience function that returns true for [`Incoming`](@ref) and false otherwise. +Convenience function that returns true for [`Incoming`](@ref) and incoming [`AbstractParticleStateful`](@ref) and false otherwise. """ function is_incoming end """ is_outgoing(dir::ParticleDirection) + is_outgoing(particle::AbstractParticleStateful) -Convenience function that returns true for [`Outgoing`](@ref) and alse otherwise. +Convenience function that returns true for [`Outgoing`](@ref) and outgoing [`AbstractParticleStateful`](@ref) and false otherwise. """ function is_outgoing end From 9f93aa5379147b5ae967d5af6761abca9f462524 Mon Sep 17 00:00:00 2001 From: AntonReinhard Date: Sat, 22 Jun 2024 00:49:34 +0200 Subject: [PATCH 21/21] Remove ELEMENT from AbstractPSP type parameters --- src/interfaces/phase_space_point.jl | 28 ++++++++++++++++------------ 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/src/interfaces/phase_space_point.jl b/src/interfaces/phase_space_point.jl index 11b4a59..7dc427a 100644 --- a/src/interfaces/phase_space_point.jl +++ b/src/interfaces/phase_space_point.jl @@ -1,5 +1,5 @@ """ - AbstractPhaseSpacePoint{PROC, MODEL, PSDEF, IN_PARTICLES, OUT_PARTICLES, ELEMENT} + AbstractPhaseSpacePoint{PROC, MODEL, PSDEF, IN_PARTICLES, OUT_PARTICLES} Representation of a point in the phase space of a process. It has several template arguments: - `PROC <: `[`AbstractProcessDefinition`](@ref) @@ -7,17 +7,17 @@ Representation of a point in the phase space of a process. It has several templa - `PSDEF <: `[`AbstractPhasespaceDefinition`](@ref) - `IN_PARTICLES <: `Tuple{Vararg{AbstractParticleStateful}}`: The tuple type of all the incoming [`AbstractParticleStateful`](@ref)s. - `OUT_PARTICLES <: `Tuple{Vararg{AbstractParticleStateful}}`: The tuple type of all the outgoing [`AbstractParticleStateful`](@ref)s. -- `ELEMENT <: `[`AbstractFourMomentum`](@ref): The type of the [`AbstractParticleStateful`](@ref)s and also the return type of the [`momentum`](@ref) interface function. The following interface functions must be provided: - `Base.getindex(psp::AbstractPhaseSpacePoint, dir::ParticleDirection, n::Int)`: Return the nth [`AbstractParticleStateful`](@ref) of the given direction. Throw `BoundsError` for invalid indices. +- `particles(psp::AbstractPhaseSpacePoint, dir::ParticleDirection)`: Return the particle tuple (type `IN_PARTICLES` or `OUT_PARTICLES` depending on `dir`) - [`process`](@ref): Return the process. - [`model`](@ref): Return the model. - [`phase_space_definition`](@ref): Return the phase space definition. From this, the following functions are automatically derived: -- `momentum(psp::AbstractPhaseSpacePoint, dir::ParticleDirection, n::Int)::ELEMENT`: Return the momentum of the nth [`AbstractParticleStateful`](@ref) of the given direction. -- `momenta(psp::PhaseSpacePoint, ::ParticleDirection)::ELEMENT`: Return a `Tuple` of all the momenta of the given direction. +- `momentum(psp::AbstractPhaseSpacePoint, dir::ParticleDirection, n::Int)`: Return the momentum of the nth [`AbstractParticleStateful`](@ref) of the given direction. +- `momenta(psp::PhaseSpacePoint, ::ParticleDirection)`: Return a `Tuple` of all the momenta of the given direction. Furthermore, an implementation of an [`AbstractPhaseSpacePoint`](@ref) has to verify on construction that it is valid, i.e., the following conditions are fulfilled: - `IN_PARTICLES` must match `incoming_particles(::PROC)` in length, order, and type **or** be an empty `Tuple`. @@ -32,7 +32,6 @@ abstract type AbstractPhaseSpacePoint{ PSDEF<:AbstractPhasespaceDefinition, IN_PARTICLES<:Tuple{Vararg{AbstractParticleStateful}}, OUT_PARTICLES<:Tuple{Vararg{AbstractParticleStateful}}, - ELEMENT<:AbstractFourMomentum, } end """ @@ -62,7 +61,7 @@ function phase_space_definition end Returns the momentum of the `n`th particle in the given [`AbstractPhaseSpacePoint`](@ref) which has direction `dir`. If `n` is outside the valid range for this phase space point, a `BoundsError` is thrown. """ function momentum(psp::AbstractPhaseSpacePoint, dir::ParticleDirection, n::Int) - return psp[dir, n].mom + return momentum(psp[dir, n]) end """ @@ -70,8 +69,13 @@ end Return a `Tuple` of all the particles' momenta for the given `ParticleDirection`. """ -momenta(psp::AbstractPhaseSpacePoint, ::QEDbase.Incoming) = momentum.(psp.in_particles) -momenta(psp::AbstractPhaseSpacePoint, ::QEDbase.Outgoing) = momentum.(psp.out_particles) +function momenta(psp::AbstractPhaseSpacePoint, ::QEDbase.Incoming) + return momentum.(particles(psp, QEDbase.Incoming())) +end + +function momenta(psp::AbstractPhaseSpacePoint, ::QEDbase.Outgoing) + return momentum.(particles(psp, QEDbase.Outgoing())) +end """ AbstractInPhaseSpacePoint @@ -80,8 +84,8 @@ A partial type specialization on [`AbstractPhaseSpacePoint`](@ref) which can be See also: [`AbstractOutPhaseSpacePoint`](@ref) """ -AbstractInPhaseSpacePoint{P,M,D,IN,OUT,E} = AbstractPhaseSpacePoint{ - P,M,D,IN,OUT,E +AbstractInPhaseSpacePoint{P,M,D,IN,OUT} = AbstractPhaseSpacePoint{ + P,M,D,IN,OUT } where {PS<:AbstractParticleStateful,IN<:Tuple{PS,Vararg},OUT<:Tuple{Vararg}} """ @@ -91,6 +95,6 @@ A partial type specialization on [`AbstractPhaseSpacePoint`](@ref) which can be See also: [`AbstractInPhaseSpacePoint`](@ref) """ -AbstractOutPhaseSpacePoint{P,M,D,IN,OUT,E} = AbstractPhaseSpacePoint{ - P,M,D,IN,OUT,E +AbstractOutPhaseSpacePoint{P,M,D,IN,OUT} = AbstractPhaseSpacePoint{ + P,M,D,IN,OUT } where {PS<:AbstractParticleStateful,IN<:Tuple{Vararg},OUT<:Tuple{PS,Vararg}}