Add QEDprocess functionality again

src/QEDcore.jl

@@ -28,6 +28,13 @@ export IncomingFermionSpinor,
 export SpinorU, SpinorUbar, SpinorV, SpinorVbar
 export @valid_spinor_input
 
+# phase space
+export SphericalCoordinateSystem
+export CenterOfMomentumFrame, ElectronRestFrame
+export PhasespaceDefinition
+export ParticleStateful, PhaseSpacePoint, InPhaseSpacePoint, OutPhaseSpacePoint
+export spin, polarization, particle_direction, particle_species, momentum, momenta, getindex
+
 using QEDbase: QEDbase
 using DocStringExtensions
 using StaticArrays
@@ -36,10 +43,18 @@ using SimpleTraits
 include("algebraic_objects/dirac_tensors/types.jl")
 include("algebraic_objects/dirac_tensors/multiplication.jl")
 
+include("phase_spaces/types.jl")
+include("phase_spaces/access.jl")
+include("phase_spaces/create.jl")
+include("phase_spaces/print.jl")
+include("phase_spaces/utility.jl")
+
 include("algebraic_objects/four_momentum.jl")
 include("algebraic_objects/lorentz_vector.jl")
 include("algebraic_objects/gamma_matrices.jl")
 
+include("particles/propagators.jl")
+include("particles/states.jl")
 include("particles/spinors.jl")
 include("particles/particle_types.jl")
 

src/algebraic_objects/lorentz_vector.jl

@@ -98,4 +98,4 @@ end
 # TODO: this breaks incremental compilation because it's trying to eval permanent changes in a different module
 #register_LorentzVectorLike(MLorentzVector)
 @traitimpl QEDbase.IsLorentzVectorLike{MLorentzVector}
-@traitimpl QEDbase.IsMutableLorentzVectorLike{MLorentzVector}
+@traitimpl QEDbase.IsMutableLorentzVectorLike{MLorentzVector}

src/particles/propagators.jl

@@ -0,0 +1,33 @@
+import QEDbase: propagator
+
+function _scalar_propagator(K::QEDbase.AbstractFourMomentum, mass::Real)
+    return one(mass) / (K * K - mass^2)
+end
+
+function _scalar_propagator(K::QEDbase.AbstractFourMomentum)
+    return one(getT(K)) / (K * K)
+end
+
+function _fermion_propagator(P::QEDbase.AbstractFourMomentum, mass::Real)
+    return (slashed(P) + mass * one(DiracMatrix)) * _scalar_propagator(P, mass)
+end
+
+function _fermion_propagator(P::QEDbase.AbstractFourMomentum)
+    return (slashed(P)) * _scalar_propagator(P)
+end
+
+function QEDbase.propagator(
+    particle_type::QEDbase.BosonLike, K::QEDbase.AbstractFourMomentum
+)
+    return _scalar_propagator(K, QEDbase.mass(particle_type))
+end
+
+function QEDbase.propagator(particle_type::QEDbase.Photon, K::QEDbase.AbstractFourMomentum)
+    return _scalar_propagator(K)
+end
+
+function QEDbase.propagator(
+    particle_type::QEDbase.FermionLike, P::QEDbase.AbstractFourMomentum
+)
+    return _fermion_propagator(P, QEDbase.mass(particle_type))
+end

src/particles/states.jl

@@ -0,0 +1,138 @@
+import QEDbase: base_state
+
+function _booster_fermion(mom::QEDbase.AbstractFourMomentum, mass::Real)
+    return (slashed(mom) + mass * one(DiracMatrix)) / (sqrt(abs(getT(mom)) + mass))
+end
+
+function _booster_antifermion(mom::QEDbase.AbstractFourMomentum, mass::Real)
+    return (mass * one(DiracMatrix) - slashed(mom)) / (sqrt(abs(getT(mom)) + mass))
+end
+
+function QEDbase.base_state(
+    particle::QEDbase.Fermion,
+    ::QEDbase.Incoming,
+    mom::QEDbase.AbstractFourMomentum,
+    spin::QEDbase.AbstractDefiniteSpin,
+)
+    booster = _booster_fermion(mom, QEDbase.mass(particle))
+    return BiSpinor(booster[:, QEDbase._spin_index(spin)])
+end
+
+function QEDbase.base_state(
+    particle::QEDbase.Fermion,
+    ::QEDbase.Incoming,
+    mom::QEDbase.AbstractFourMomentum,
+    spin::QEDbase.AllSpin,
+)
+    booster = _booster_fermion(mom, QEDbase.mass(particle))
+    return SVector(BiSpinor(booster[:, 1]), BiSpinor(booster[:, 2]))
+end
+
+function QEDbase.base_state(
+    particle::QEDbase.AntiFermion,
+    ::QEDbase.Incoming,
+    mom::QEDbase.AbstractFourMomentum,
+    spin::QEDbase.AbstractDefiniteSpin,
+)
+    booster = _booster_antifermion(mom, QEDbase.mass(particle))
+    return AdjointBiSpinor(BiSpinor(booster[:, QEDbase._spin_index(spin) + 2])) * GAMMA[1]
+end
+
+function QEDbase.base_state(
+    particle::QEDbase.AntiFermion,
+    ::QEDbase.Incoming,
+    mom::QEDbase.AbstractFourMomentum,
+    spin::QEDbase.AllSpin,
+)
+    booster = _booster_antifermion(mom, QEDbase.mass(particle))
+    return SVector(
+        AdjointBiSpinor(BiSpinor(booster[:, 3])) * GAMMA[1],
+        AdjointBiSpinor(BiSpinor(booster[:, 4])) * GAMMA[1],
+    )
+end
+
+function QEDbase.base_state(
+    particle::QEDbase.Fermion,
+    ::QEDbase.Outgoing,
+    mom::QEDbase.AbstractFourMomentum,
+    spin::QEDbase.AbstractDefiniteSpin,
+)
+    booster = _booster_fermion(mom, QEDbase.mass(particle))
+    return AdjointBiSpinor(BiSpinor(booster[:, QEDbase._spin_index(spin)])) * GAMMA[1]
+end
+
+function QEDbase.base_state(
+    particle::QEDbase.Fermion,
+    ::QEDbase.Outgoing,
+    mom::QEDbase.AbstractFourMomentum,
+    spin::QEDbase.AllSpin,
+)
+    booster = _booster_fermion(mom, QEDbase.mass(particle))
+    return SVector(
+        AdjointBiSpinor(BiSpinor(booster[:, 1])) * GAMMA[1],
+        AdjointBiSpinor(BiSpinor(booster[:, 2])) * GAMMA[1],
+    )
+end
+
+function QEDbase.base_state(
+    particle::QEDbase.AntiFermion,
+    ::QEDbase.Outgoing,
+    mom::QEDbase.AbstractFourMomentum,
+    spin::QEDbase.AbstractDefiniteSpin,
+)
+    booster = _booster_antifermion(mom, QEDbase.mass(particle))
+    return BiSpinor(booster[:, QEDbase._spin_index(spin) + 2])
+end
+
+function QEDbase.base_state(
+    particle::QEDbase.AntiFermion,
+    ::QEDbase.Outgoing,
+    mom::QEDbase.AbstractFourMomentum,
+    spin::QEDbase.AllSpin,
+)
+    booster = _booster_antifermion(mom, QEDbase.mass(particle))
+    return SVector(BiSpinor(booster[:, 3]), BiSpinor(booster[:, 4]))
+end
+
+function _photon_state(pol::QEDbase.AllPolarization, mom::QEDbase.AbstractFourMomentum)
+    cth = QEDbase.getCosTheta(mom)
+    sth = sqrt(1 - cth^2)
+    cos_phi = QEDbase.getCosPhi(mom)
+    sin_phi = QEDbase.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::QEDbase.PolarizationX, mom::QEDbase.AbstractFourMomentum)
+    cth = QEDbase.getCosTheta(mom)
+    sth = sqrt(1 - cth^2)
+    cos_phi = QEDbase.getCosPhi(mom)
+    sin_phi = QEDbase.getSinPhi(mom)
+    return SLorentzVector{Float64}(0.0, cth * cos_phi, cth * sin_phi, -sth)
+end
+
+function _photon_state(pol::QEDbase.PolarizationY, mom::QEDbase.AbstractFourMomentum)
+    cos_phi = QEDbase.getCosPhi(mom)
+    sin_phi = QEDbase.getSinPhi(mom)
+    return SLorentzVector{Float64}(0.0, -sin_phi, cos_phi, 0.0)
+end
+
+@inline function QEDbase.base_state(
+    particle::QEDbase.Photon,
+    ::QEDbase.ParticleDirection,
+    mom::QEDbase.AbstractFourMomentum,
+    pol::QEDbase.AllPolarization,
+)
+    return _photon_state(pol, mom)
+end
+
+@inline function QEDbase.base_state(
+    particle::QEDbase.Photon,
+    ::QEDbase.ParticleDirection,
+    mom::QEDbase.AbstractFourMomentum,
+    pol::QEDbase.AbstractPolarization,
+)
+    return _photon_state(pol, mom)
+end

src/phase_spaces/access.jl

@@ -0,0 +1,62 @@
+import QEDbase:
+    is_particle,
+    is_anti_particle,
+    is_fermion,
+    is_boson,
+    is_incoming,
+    is_outgoing,
+    mass,
+    charge
+
+# particle interface
+@inline is_incoming(particle::ParticleStateful) = is_incoming(particle.dir)
+@inline is_outgoing(particle::ParticleStateful) = is_outgoing(particle.dir)
+@inline QEDbase.is_fermion(particle::ParticleStateful) =
+    QEDbase.is_fermion(particle.species)
+@inline QEDbase.is_boson(particle::ParticleStateful) = QEDbase.is_boson(particle.species)
+@inline QEDbase.is_particle(particle::ParticleStateful) =
+    QEDbase.is_particle(particle.species)
+@inline QEDbase.is_anti_particle(particle::ParticleStateful) =
+    QEDbase.is_anti_particle(particle.species)
+@inline QEDbase.mass(particle::ParticleStateful) = QEDbase.mass(particle.species)
+@inline QEDbase.charge(particle::ParticleStateful) = QEDbase.charge(particle.species)
+
+# accessors
+particle_direction(part::ParticleStateful) = part.dir
+particle_species(part::ParticleStateful) = part.species
+momentum(part::ParticleStateful) = part.mom
+
+"""
+    momenta(psp::PhaseSpacePoint, ::ParticleDirection)
+
+Return a `Tuple` of all the particles' momenta for the given `ParticleDirection`.
+"""
+momenta(psp::PhaseSpacePoint, ::QEDbase.Incoming) = momentum.(psp.in_particles)
+momenta(psp::PhaseSpacePoint, ::QEDbase.Outgoing) = momentum.(psp.out_particles)
+
+"""
+    Base.getindex(psp::PhaseSpacePoint, dir::Incoming, n::Int)
+
+Overload for the array indexing operator `[]`. Returns the nth incoming particle in this phase space point.
+"""
+function Base.getindex(psp::PhaseSpacePoint, ::QEDbase.Incoming, n::Int)
+    return psp.in_particles[n]
+end
+
+"""
+    Base.getindex(psp::PhaseSpacePoint, dir::Outgoing, n::Int)
+
+Overload for the array indexing operator `[]`. Returns the nth outgoing particle in this phase space point.
+"""
+function Base.getindex(psp::PhaseSpacePoint, ::QEDbase.Outgoing, n::Int)
+    return psp.out_particles[n]
+end
+
+"""
+    momentum(psp::PhaseSpacePoint, dir::ParticleDirection, n::Int)
+
+Returns the momentum of the `n`th particle in the given [`PhaseSpacePoint`](@ref) which has direction `dir`. If `n` is outside the valid range for this phase space point, a `BoundsError` is thrown.
+"""
+function momentum(psp::PhaseSpacePoint, dir::QEDbase.ParticleDirection, n::Int)
+    return psp[dir, n].mom
+end