diff --git a/.formatting/.gitkeep b/.formatting/.gitkeep deleted file mode 100644 index e69de29..0000000 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/Project.toml b/Project.toml index 75dc58d..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,7 +24,8 @@ julia = "1.6" [extras] SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" +Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["SafeTestsets","Test"] +test = ["SafeTestsets", "Test", "Suppressor"] 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/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/docs/make.jl b/docs/make.jl index 4e0e719..afc3b91 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -11,6 +11,9 @@ 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 bib = CitationBibliography(joinpath(@__DIR__, "Bibliography.bib")) 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 8439461..07a7803 100644 --- a/src/QEDbase.jl +++ b/src/QEDbase.jl @@ -1,12 +1,5 @@ module QEDbase -using SimpleTraits -using ArgCheck -using ConstructionBase - -import Base: * -import StaticArrays: similar_type - export minkowski_dot, mdot, register_LorentzVectorLike export getT, getX, getY, getZ export getMagnitude2, getMag2, getMagnitude, getMag @@ -28,32 +21,22 @@ 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 +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 -export FermionLike, Fermion, AntiFermion, MajoranaFermion -export BosonLike, Boson, AntiBoson, MajoranaBoson -export Electron, Positron, Photon +# directions export ParticleDirection, Incoming, Outgoing export is_incoming, is_outgoing @@ -65,22 +48,73 @@ export AbstractDefiniteSpin, AbstractIndefiniteSpin export SpinUp, SpinDown, AllSpin export multiplicity +# 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, compute +export AbstractProcessSetup, scattering_process, physical_model + +# 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 + using StaticArrays using LinearAlgebra using DocStringExtensions -include("dirac_tensors.jl") -include("lorentz_interface.jl") -include("lorentz_vector.jl") -include("gamma_matrices.jl") +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_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") + +include("interfaces/four_momentum.jl") +include("interfaces/model.jl") + +include("interfaces/particle.jl") + +include("particles/direction.jl") +include("particles/spin_pol.jl") -include("four_momentum.jl") # maybe go to a kinematics module!! +include("interfaces/phase_space.jl") +include("interfaces/process.jl") +include("interfaces/particle_stateful.jl") +include("interfaces/phase_space_point.jl") +include("interfaces/setup.jl") -include("interfaces/particle_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("total_cross_section.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/errors.jl b/src/errors.jl new file mode 100644 index 0000000..aacd618 --- /dev/null +++ b/src/errors.jl @@ -0,0 +1,65 @@ +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 + +mutable struct SpinorConstructionError <: Exception + var::String +end + +function Base.showerror(io::IO, e::SpinorConstructionError) + return print(io, "SpinorConstructionError: ", e.var) +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/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/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..f79de00 --- /dev/null +++ b/src/interfaces/four_momentum.jl @@ -0,0 +1,22 @@ +""" + AbstractFourMomentum + +Abstract base type for four-momentas, representing one energy and three spacial components. + +Also see: `QEDcore.SFourMomentum`, `QEDcore.MFourMomentum` +""" +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/gamma_matrices.jl b/src/interfaces/gamma_matrices.jl new file mode 100644 index 0000000..2bdc5e5 --- /dev/null +++ b/src/interfaces/gamma_matrices.jl @@ -0,0 +1,5 @@ +#### +# generic definition of the gamma matrices +#### + +abstract type AbstractGammaRepresentation end diff --git a/src/lorentz_interface.jl b/src/interfaces/lorentz_vectors/arithmetic.jl similarity index 82% rename from src/lorentz_interface.jl rename to src/interfaces/lorentz_vectors/arithmetic.jl index 8123e00..6f6e056 100644 --- a/src/lorentz_interface.jl +++ b/src/interfaces/lorentz_vectors/arithmetic.jl @@ -1,189 +1,15 @@ -""" -## 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 +import Base: * -function Base.showerror(io::IO, err::RegistryError) - return println( - io, - "RegistryError:", - " The type <$(err.target_type)> must implement <$(err.func)> to be registered.", - ) +function dot(p1::T1, p2::T2) where {T1<:AbstractLorentzVector,T2<:AbstractLorentzVector} + return mdot(p1, p2) 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) +@inline function *( + p1::T1, p2::T2 +) where {T1<:AbstractLorentzVector,T2<:AbstractLorentzVector} + return dot(p1, p2) 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) Return the Minkowski dot product of two `LorentzVectorLike`. @@ -212,11 +38,7 @@ Function alias for [`minkowski_dot`](@ref). """ 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 +63,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 +87,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 +105,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 +138,6 @@ const getMass = getInvariantMass ########################## """ - getE(lv) Return the energy component of a given `LorentzVectorLike`, i.e. its 0-component. @@ -335,7 +153,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 +165,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 +177,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 +189,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 +213,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 +230,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 +253,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 +277,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 +290,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 +298,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 +329,6 @@ end const getMt = getTransverseMass """ - getRapidity(lv) Return the [rapidity](https://en.wikipedia.org/wiki/Rapidity) for a given `LorentzVectorLike`. @@ -550,7 +357,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 +381,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 +396,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 +415,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 +430,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 +448,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 +470,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 +498,6 @@ end #### """ - setE!(lv,value) Sets the energy component of a given `LorentzVectorLike` to a given `value`. @@ -715,7 +514,6 @@ end const setEnergy! = setE! """ - setPx!(lv,value) Sets the 1-component of a given `LorentzVectorLike` to a given `value`. @@ -732,7 +530,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 +546,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 +564,6 @@ end # setter spherical coordinates """ - setTheta!(lv,value) Sets the theta angle of a `LorentzVectorLike` to a given `value`. @@ -791,7 +586,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 +610,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 +633,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 +656,6 @@ end # setter light cone coordinates """ - setPlus!(lv,value) Sets the plus component of a `LorentzVectorLike` to a given `value`. @@ -882,7 +673,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 +691,6 @@ end # transverse coordinates """ - setTransverseMomentum!(lv,value) Sets the transverse momentum of a `LorentzVectorLike` to a given `value`. @@ -928,7 +717,6 @@ const setPerp! = setTransverseMomentum! const setPt! = setTransverseMomentum! """ - setTransverseMass!(lv,value) Sets the transverse mass of a `LorentzVectorLike` to a given `value`. @@ -953,7 +741,6 @@ end const setMt! = setTransverseMass! """ - setRapidity!(lv,value) Sets the rapidity of a `LorentzVectorLike` to a given `value`. diff --git a/src/interfaces/lorentz_vectors/dirac_interaction.jl b/src/interfaces/lorentz_vectors/dirac_interaction.jl new file mode 100644 index 0000000..990d949 --- /dev/null +++ b/src/interfaces/lorentz_vectors/dirac_interaction.jl @@ -0,0 +1,40 @@ +import Base: * + +""" +$(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/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..46dd554 --- /dev/null +++ b/src/interfaces/lorentz_vectors/utility.jl @@ -0,0 +1,49 @@ +####### +# +# 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/model.jl b/src/interfaces/model.jl new file mode 100644 index 0000000..7b98325 --- /dev/null +++ b/src/interfaces/model.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/particles/particle_states.jl b/src/interfaces/particle.jl similarity index 50% rename from src/particles/particle_states.jl rename to src/interfaces/particle.jl index 4729939..f71d693 100644 --- a/src/particles/particle_states.jl +++ b/src/interfaces/particle.jl @@ -1,11 +1,120 @@ -function _booster_fermion(mom::QEDbase.AbstractFourMomentum, mass::Real) - return (slashed(mom) + mass * one(DiracMatrix)) / (sqrt(abs(getT(mom)) + mass)) -end +############### +# The particle interface +# +# In this file, we define the interface for working with particles in a general +# sense. +############### -function _booster_antifermion(mom::QEDbase.AbstractFourMomentum, mass::Real) - return (mass * one(DiracMatrix) - slashed(mom)) / (sqrt(abs(getT(mom)) + mass)) -end +""" +Abstract base type for every type which might be considered a *particle* in the context of `QED.jl`. For every (concrete) subtype of `AbstractParticle`, there are two kinds of interface functions implemented: static functions and property functions. +The static functions provide information on what kind of particle it is (defaults are written in square brackets) + +```julia + is_fermion(::AbstractParticle)::Bool [= false] + is_boson(::AbstractParticle)::Bool [= false] + is_particle(::AbstractParticle)::Bool [= true] + is_anti_particle(::AbstractParticle)::Bool [= false] +``` +If the output of those functions differ from the defaults for a subtype of `AbstractParticle`, these functions need to be overwritten. +The second type of functions define a hard interface for `AbstractParticle`: + +```julia + mass(::AbstractParticle)::Real + charge(::AbstractParticle)::Real +``` +These functions must be implemented in order to have the subtype of `AbstractParticle` work with the functionalities of `QEDprocesses.jl`. +""" +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) + +Interface function for particles. Return `true` if the passed subtype of [`AbstractParticle`](@ref) can be considered a *fermion* in the sense of particle statistics, and `false` otherwise. +The default implementation of `is_fermion` for every subtype of [`AbstractParticle`](@ref) will always return `false`. +""" +is_fermion(::AbstractParticle) = false + +""" + $(TYPEDSIGNATURES) + +Interface function for particles. Return `true` if the passed subtype of [`AbstractParticle`](@ref) can be considered a *boson* in the sense of particle statistics, and `false` otherwise. +The default implementation of `is_boson` for every subtype of [`AbstractParticle`](@ref) will always return `false`. +""" +is_boson(::AbstractParticle) = false + +""" + $(TYPEDSIGNATURES) + +Interface function for particles. Return `true` if the passed subtype of [`AbstractParticle`](@ref) can be considered a *particle* as distinct from anti-particles, and `false` otherwise. +The default implementation of `is_particle` for every subtype of [`AbstractParticle`](@ref) will always return `true`. +""" +is_particle(::AbstractParticle) = true + +""" + $(TYPEDSIGNATURES) +Interface function for particles. Return true if the passed subtype of [`AbstractParticle`](@ref) can be considered an *anti particle* as distinct from their particle counterpart, and `false` otherwise. +The default implementation of `is_anti_particle` for every subtype of [`AbstractParticle`](@ref) will always return `false`. +""" +is_anti_particle(::AbstractParticle) = false + +""" + mass(particle::AbstractParticle)::Real + +Interface function for particles. Return the rest mass of a particle (in units of the electron mass). + +This needs to be implemented for each concrete subtype of [`AbstractParticle`](@ref). +""" +function mass end + +""" + charge(::AbstractParticle)::Real + +Interface function for particles. Return the electric charge of a particle (in units of the elementary electric charge). + +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 + +# TODO: Turn the doctest on again when QEDcore has QEDprocesses functionality """ ```julia base_state( @@ -20,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. [`Electron`](@ref), [`Positron`](@ref), or [`Photon`](@ref). +- `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). @@ -53,7 +162,6 @@ base_state(::Photon, ::Outgoing, mom) # -> SVector{2,SLorentzVector{Complex # Example ```julia - using QEDbase mass = 1.0 # set electron mass to 1.0 @@ -63,12 +171,12 @@ mom = SFourMomentum(E, px, py, pz) # initialize the four-momentum of the el # compute the state of an incoming electron with spin = SpinUp # note: base_state is not exported! -electron_state = base_state(Electron(), Incoming(), mom, SpinUp()) - +electron_state = base_state(QEDcore.Electron(), Incoming(), mom, SpinUp()) ``` -```jldoctest -julia> using QEDbase +TODO: Reenable doctests +```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) 4-element SFourMomentum with indices SOneTo(4): @@ -150,133 +258,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/interfaces/particle_interface.jl b/src/interfaces/particle_interface.jl deleted file mode 100644 index 494a310..0000000 --- a/src/interfaces/particle_interface.jl +++ /dev/null @@ -1,78 +0,0 @@ -############### -# The particle interface -# -# In this file, we define the interface for working with particles in a general -# sense. -############### - -""" -Abstract base type for every type which might be considered a *particle* in the context of `QED.jl`. For every (concrete) subtype of `AbstractParticle`, there are two kinds of interface functions implemented: static functions and property functions. -The static functions provide information on what kind of particle it is (defaults are written in square brackets) - -```julia - is_fermion(::AbstractParticle)::Bool [= false] - is_boson(::AbstractParticle)::Bool [= false] - is_particle(::AbstractParticle)::Bool [= true] - is_anti_particle(::AbstractParticle)::Bool [= false] -``` -If the output of those functions differ from the defaults for a subtype of `AbstractParticle`, these functions need to be overwritten. -The second type of functions define a hard interface for `AbstractParticle`: - -```julia - mass(::AbstractParticle)::Real - charge(::AbstractParticle)::Real -``` -These functions must be implemented in order to have the subtype of `AbstractParticle` work with the functionalities of `QEDprocesses.jl`. -""" -abstract type AbstractParticle end -Base.broadcastable(part::AbstractParticle) = Ref(part) - -""" - $(TYPEDSIGNATURES) - -Interface function for particles. Return `true` if the passed subtype of [`AbstractParticle`](@ref) can be considered a *fermion* in the sense of particle statistics, and `false` otherwise. -The default implementation of `is_fermion` for every subtype of [`AbstractParticle`](@ref) will always return `false`. -""" -is_fermion(::AbstractParticle) = false - -""" - $(TYPEDSIGNATURES) - -Interface function for particles. Return `true` if the passed subtype of [`AbstractParticle`](@ref) can be considered a *boson* in the sense of particle statistics, and `false` otherwise. -The default implementation of `is_boson` for every subtype of [`AbstractParticle`](@ref) will always return `false`. -""" -is_boson(::AbstractParticle) = false - -""" - $(TYPEDSIGNATURES) - -Interface function for particles. Return `true` if the passed subtype of [`AbstractParticle`](@ref) can be considered a *particle* as distinct from anti-particles, and `false` otherwise. -The default implementation of `is_particle` for every subtype of [`AbstractParticle`](@ref) will always return `true`. -""" -is_particle(::AbstractParticle) = true - -""" - $(TYPEDSIGNATURES) - -Interface function for particles. Return true if the passed subtype of [`AbstractParticle`](@ref) can be considered an *anti particle* as distinct from their particle counterpart, and `false` otherwise. -The default implementation of `is_anti_particle` for every subtype of [`AbstractParticle`](@ref) will always return `false`. -""" -is_anti_particle(::AbstractParticle) = false - -""" - mass(particle::AbstractParticle)::Real - -Interface function for particles. Return the rest mass of a particle (in units of the electron mass). - -This needs to be implemented for each concrete subtype of [`AbstractParticle`](@ref). -""" -function mass end - -""" - charge(::AbstractParticle)::Real - -Interface function for particles. Return the electric charge of a particle (in units of the elementary electric charge). - -This needs to be implemented for each concrete subtype of [`AbstractParticle`](@ref). -""" -function charge end diff --git a/src/interfaces/particle_stateful.jl b/src/interfaces/particle_stateful.jl new file mode 100644 index 0000000..9337c4e --- /dev/null +++ b/src/interfaces/particle_stateful.jl @@ -0,0 +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.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/phase_space_point.jl b/src/interfaces/phase_space_point.jl new file mode 100644 index 0000000..7dc427a --- /dev/null +++ b/src/interfaces/phase_space_point.jl @@ -0,0 +1,100 @@ +""" + 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) +- `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. + +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)`: 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`. +- `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, + PSDEF<:AbstractPhasespaceDefinition, + IN_PARTICLES<:Tuple{Vararg{AbstractParticleStateful}}, + OUT_PARTICLES<:Tuple{Vararg{AbstractParticleStateful}}, +} 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 momentum(psp[dir, n]) +end + +""" + momenta(psp::AbstractPhaseSpacePoint, ::ParticleDirection) + +Return a `Tuple` of all the particles' momenta for the given `ParticleDirection`. +""" +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 + +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} = AbstractPhaseSpacePoint{ + P,M,D,IN,OUT +} 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} = AbstractPhaseSpacePoint{ + P,M,D,IN,OUT +} where {PS<:AbstractParticleStateful,IN<:Tuple{Vararg},OUT<:Tuple{PS,Vararg}} diff --git a/src/interfaces/process.jl b/src/interfaces/process.jl new file mode 100644 index 0000000..f783de5 --- /dev/null +++ b/src/interfaces/process.jl @@ -0,0 +1,256 @@ +############### +# 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 process 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`. +""" +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`. + +!!! 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) + +##### +# 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/setup.jl b/src/interfaces/setup.jl new file mode 100644 index 0000000..7b0f8ce --- /dev/null +++ b/src/interfaces/setup.jl @@ -0,0 +1,164 @@ +############### +# 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 + +""" + _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/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/particle_direction.jl b/src/particles/direction.jl similarity index 74% rename from src/particles/particle_direction.jl rename to src/particles/direction.jl index 9064007..66ac90f 100644 --- a/src/particles/particle_direction.jl +++ b/src/particles/direction.jl @@ -4,6 +4,22 @@ 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) + is_incoming(particle::AbstractParticleStateful) + +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 outgoing [`AbstractParticleStateful`](@ref) and false otherwise. +""" +function is_outgoing end + """ Incoming <: ParticleDirection diff --git a/src/particles/particle_spinors.jl b/src/particles/particle_spinors.jl deleted file mode 100644 index 56fb1c3..0000000 --- a/src/particles/particle_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/particles/particle_types.jl b/src/particles/particle_types.jl deleted file mode 100644 index 16282e1..0000000 --- a/src/particles/particle_types.jl +++ /dev/null @@ -1,209 +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. -############### - -""" - 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") diff --git a/src/particles/particle_spin_pol.jl b/src/particles/spin_pol.jl similarity index 64% rename from src/particles/particle_spin_pol.jl rename to src/particles/spin_pol.jl index 06f6e2f..9695405 100644 --- a/src/particles/particle_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 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 [`FermionLike`](@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 [`FermionLike`](@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 [`FermionLike`](@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 a [`FermionLike`](@ref) 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 a [`FermionLike`](@ref) has spin-down. +Concrete type indicating that a particle with [`is_fermion`](@ref) 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 particles with [`is_boson`](@ref). """ abstract type AbstractPolarization <: AbstractSpinOrPolarization end """ -Abstract base type for definite polarizations of [`BosonLike`](@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 [`BosonLike`](@ref) particles. +Abstract base type for indefinite polarization of particles with [`is_boson`](@ref). 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 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 a [`BosonLike`](@ref) 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 [`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 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" 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 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/coverage/.gitkeep b/test/coverage/.gitkeep deleted file mode 100644 index e69de29..0000000 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/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/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/model.jl b/test/interfaces/model.jl new file mode 100644 index 0000000..42ec606 --- /dev/null +++ b/test/interfaces/model.jl @@ -0,0 +1,22 @@ +using QEDbase + +struct TestModel <: AbstractModelDefinition end +QEDbase.fundamental_interaction_type(::TestModel) = :test_interaction + +struct TestModel_FAIL <: AbstractModelDefinition end + +@testset "hard interface" begin + TESTMODEL = TestModel() + @test fundamental_interaction_type(TESTMODEL) == :test_interaction +end + +@testset "interface fail" begin + TESTMODEL_FAIL = TestModel_FAIL() + @test_throws MethodError fundamental_interaction_type(TESTMODEL_FAIL) +end + +@testset "broadcast" begin + test_func(model) = model + TESTMODEL = TestModel() + @test test_func.(TESTMODEL) == TESTMODEL +end 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 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/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 deleted file mode 100644 index 6129c1b..0000000 --- a/test/particles.jl +++ /dev/null @@ -1,215 +0,0 @@ -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 -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 "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()) - @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 - -@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 3860a86..8c6062c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,30 +3,16 @@ using Test using SafeTestsets begin - @time @safetestset "Dirac tensors" begin - include("dirac_tensor.jl") - end - - @time @safetestset "Lorentz Vectors" begin - include("lorentz_vector.jl") + # Interfaces + @time @safetestset "model interface" begin + include("interfaces/model.jl") end @time @safetestset "Lorentz interface" begin - include("lorentz_interface.jl") - end - @time @safetestset "Gamma matrices" begin - include("gamma_matrices.jl") - end - - @time @safetestset "particle spinors" begin - include("particle_spinors.jl") - end - - @time @safetestset "four momentum" begin - include("four_momentum.jl") + include("interfaces/lorentz.jl") end @time @safetestset "particles" begin - include("particles.jl") + include("particle_properties.jl") end end diff --git a/test/utils.jl b/test/utils.jl deleted file mode 100644 index b4d4f8e..0000000 --- a/test/utils.jl +++ /dev/null @@ -1,11 +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