GenericQEDProcess

.github/workflows/formatter.yaml

@@ -9,7 +9,7 @@ jobs:
       - name: install Julia
         uses: julia-actions/setup-julia@v1
         with:
-          version: 1.10
+          version: "1.10"
       - name: Install Julia requirements
         run: julia --project=${GITHUB_WORKSPACE}/.formatting -e 'import Pkg; Pkg.instantiate()'
       - name: Check code style

src/QEDprocesses.jl

@@ -4,14 +4,12 @@ module QEDprocesses
 export ALPHA,
     ALPHA_SQUARE, ELEMENTARY_CHARGE, ELEMENTARY_CHARGE_SQUARE, ELECTRONMASS, ONE_OVER_FOURPI
 
-# propagator
-export propagator
-
 # specific compute models
 export PerturbativeQED
 
 # specific scattering processes
 export Compton, omega_prime
+export ScatteringProcess, isphysical
 
 using QEDbase
 using QEDcore
@@ -23,6 +21,11 @@ include("utils.jl")
 
 include("models/models.jl")
 
+# generic qed process
+include("processes/generic_process/utility.jl")
+include("processes/generic_process/process.jl")
+include("processes/generic_process/perturbative/cross_section.jl")
+
 # one photon compton
 include("processes/one_photon_compton/process.jl")
 include("processes/one_photon_compton/perturbative/kinematics.jl")

src/models/perturbative_qed.jl

@@ -28,3 +28,17 @@ function Base.show(io::IO, ::PerturbativeQED)
     print(io, "perturbative QED")
     return nothing
 end
+
+"""
+    isphysical(proc::AbstractProcessDefinition, model::PerturbativeQED)
+
+A utility function that returns whether a given `AbstractProcessDefinition` conserves the number and charge of fermions and has at least 2 participating particles.
+"""
+function isphysical(proc::AbstractProcessDefinition, ::PerturbativeQED)
+    return (
+        number_particles(proc, Incoming(), Electron()) +
+        number_particles(proc, Outgoing(), Positron()) ==
+        number_particles(proc, Incoming(), Positron()) +
+        number_particles(proc, Outgoing(), Electron())
+    ) && number_particles(proc, Incoming()) + number_particles(proc, Outgoing()) >= 2
+end

src/processes/generic_process/perturbative/cross_section.jl

@@ -0,0 +1,20 @@
+# None of these functions are currently implemented.
+# They will be implemented using the QEDFeynmanDiagrams project when that is released.
+
+#=
+function QEDbase._incident_flux(psp::InPhaseSpacePoint{<:ScatteringProcess,PerturbativeQED}) end
+
+function QEDbase._matrix_element_squared(psp::PhaseSpacePoint{<:ScatteringProcess,PerturbativeQED})
+    proc = process(psp)
+    # by using FeynmanDiagramGenerator.jl
+    return proc.matrix_element_squared(psp)
+end
+
+function QEDbase._averaging_norm(proc::<:ScatteringProcess) end
+
+function QEDbase._is_in_phasespace(psp::PhaseSpacePoint{<:ScatteringProcess,PerturbativeQED}) end
+
+function QEDbase._phase_space_factor(
+    psp::PhaseSpacePoint{<:ScatteringProcess,PerturbativeQED,CustomPhasespaceDefinition}
+) end
+=#

src/processes/generic_process/process.jl

@@ -0,0 +1,112 @@
+"""
+    ScatteringProcess <: AbstractProcessDefinition
+
+Generic implementation for scattering processes of arbitrary particles. Currently, only calculations in combination with `PerturbativeQED` are supported. 
+However, this is supposed to describe scattering processes with any number of incoming and outgoing particles, and any combination of spins or polarizations for the particles.
+
+The [`isphysical`](@ref) function can be used to check whether the process is possible in perturbative QED.
+
+!!! warning  
+    The computation of cross sections and probabilities is currently unimplemented.
+
+## Constructors
+
+    ScatteringProcess(
+        in_particles::Tuple{AbstractParticleType},
+        out_particles::Tuple{AbstractParticleType},
+        [in_sp::Tuple{AbstractSpinOrPolarization},
+        out_sp::Tuple{AbstractSpinOrPolarization}]
+    )
+
+Constructor for a ScatteringProcess with the given incoming and outgoing particles and their respective spins and pols.
+The constructor asserts that the particles are compatible with their respective spins and polarizations. If the assertion fails, an 
+`InvalidInputError` is thrown.
+
+The `in_sp` and `out_sp` parameters can be omitted in which case all spins and polarizations will be set to `AllSpin` and `AllPol` for every fermion and boson, respectively.
+"""
+struct ScatteringProcess{INT,OUTT,INSP,OUTSP} <:
+       AbstractProcessDefinition where {INT<:Tuple,OUTT<:Tuple,INSP<:Tuple,OUTSP<:Tuple}
+    incoming_particles::INT
+    outgoing_particles::OUTT
+
+    incoming_spin_pols::INSP
+    outgoing_spin_pols::OUTSP
+
+    function ScatteringProcess(
+        in_particles::NTuple{I,AbstractParticleType},
+        out_particles::NTuple{O,AbstractParticleType},
+        in_spin_pols::NTuple{I,AbstractSpinOrPolarization},
+        out_spin_pols::NTuple{O,AbstractSpinOrPolarization},
+    ) where {I,O}
+        _assert_spin_pol_particle_compatability(in_particles, in_spin_pols)
+        _assert_spin_pol_particle_compatability(out_particles, out_spin_pols)
+
+        return new{
+            typeof(in_particles),
+            typeof(out_particles),
+            typeof(in_spin_pols),
+            typeof(out_spin_pols),
+        }(
+            in_particles, out_particles, in_spin_pols, out_spin_pols
+        )
+    end
+end
+
+function ScatteringProcess(
+    in_particles::NTuple{I,AbstractParticleType},
+    out_particles::NTuple{O,AbstractParticleType},
+) where {I,O}
+    in_spin_pols = ntuple(
+        x -> is_fermion(in_particles[x]) ? AllSpin() : AllPolarization(),
+        length(in_particles),
+    )
+    out_spin_pols = ntuple(
+        x -> is_fermion(out_particles[x]) ? AllSpin() : AllPolarization(),
+        length(out_particles),
+    )
+    return ScatteringProcess(in_particles, out_particles, in_spin_pols, out_spin_pols)
+end
+
+function QEDbase.incoming_particles(proc::ScatteringProcess)
+    return proc.incoming_particles
+end
+function QEDbase.outgoing_particles(proc::ScatteringProcess)
+    return proc.outgoing_particles
+end
+function QEDbase.incoming_spin_pols(proc::ScatteringProcess)
+    return proc.incoming_spin_pols
+end
+function QEDbase.outgoing_spin_pols(proc::ScatteringProcess)
+    return proc.outgoing_spin_pols
+end
+
+function Base.show(io::IO, proc::ScatteringProcess)
+    print(io, "generic QED process \"")
+    for p in incoming_particles(proc)
+        print(io, _particle_to_letter(p))
+    end
+    print(io, " -> ")
+    for p in outgoing_particles(proc)
+        print(io, _particle_to_letter(p))
+    end
+    print(io, "\"")
+    return nothing
+end
+
+function Base.show(io::IO, ::MIME"text/plain", proc::ScatteringProcess)
+    println(io, "generic QED process")
+    for dir in (Incoming(), Outgoing())
+        first = true
+        for (p, sp) in zip(particles(proc, dir), spin_pols(proc, dir))
+            if !first
+                print(io, ", ")
+            else
+                print(io, "    $(dir): ")
+                first = false
+            end
+            print(io, "$(p) ($(sp))")
+        end
+        println(io)
+    end
+    return nothing
+end

src/processes/generic_process/utility.jl

@@ -0,0 +1,30 @@
+# recursion success, all particles are compatible with their spin/pol
+_assert_spin_pol_particle_compatability(::Tuple{}, ::Tuple{}) = nothing
+
+# recursion base case: check first particle against first spin/pol, then recurse
+# note: the length of the tuples is expected to be the same, the constructor ensures this by using NTuples
+function _assert_spin_pol_particle_compatability(
+    particles::Tuple{AbstractParticleType,Vararg},
+    spin_pols::Tuple{AbstractSpinOrPolarization,Vararg},
+)
+    if is_fermion(particles[1]) && !(spin_pols[1] isa AbstractSpin)
+        throw(
+            InvalidInputError(
+                "particle \"$(particles[1])\" is a fermion and should have a spin, but has \"$(spin_pols[1])\"",
+            ),
+        )
+    end
+    if is_boson(particles[1]) && !(spin_pols[1] isa AbstractPolarization)
+        throw(
+            InvalidInputError(
+                "particle \"$(particles[1])\" is a boson and should have a polarization, but has \"$(spin_pols[1])\"",
+            ),
+        )
+    end
+    return _assert_spin_pol_particle_compatability(particles[2:end], spin_pols[2:end])
+end
+
+# this should move to QEDbase as part of the interface, see https://github.com/QEDjl-project/QEDbase.jl/issues/114
+_particle_to_letter(::Electron) = "e"
+_particle_to_letter(::Positron) = "p"
+_particle_to_letter(::Photon) = "k"

src/processes/one_photon_compton/perturbative/cross_section.jl

@@ -3,7 +3,7 @@
 # Implementation of the cross section interface
 #####
 
-function QEDbase._incident_flux(in_psp::InPhaseSpacePoint{<:Compton,<:PerturbativeQED})
+function QEDbase._incident_flux(in_psp::InPhaseSpacePoint{<:Compton,PerturbativeQED})
     return momentum(in_psp, Incoming(), 1) * momentum(in_psp, Incoming(), 2)
 end
 
@@ -39,9 +39,7 @@ end
         )
 end
 
-@inline function QEDbase._is_in_phasespace(
-    psp::PhaseSpacePoint{<:Compton,<:PerturbativeQED}
-)
+@inline function QEDbase._is_in_phasespace(psp::PhaseSpacePoint{<:Compton,PerturbativeQED})
     @inbounds if (
         !isapprox(
             momentum(psp, Incoming(), 1) + momentum(psp, Incoming(), 2),

src/processes/one_photon_compton/perturbative/total_probability.jl

@@ -1,4 +1,4 @@
-function QEDbase._total_probability(in_psp::InPhaseSpacePoint{<:Compton,<:PerturbativeQED})
+function QEDbase._total_probability(in_psp::InPhaseSpacePoint{<:Compton,PerturbativeQED})
     omega = getE(momentum(in_psp[Incoming(), 2]))
 
     function func(x)

test/processes/generic_process/groundtruths.jl

@@ -0,0 +1,25 @@
+POLS = [PolX(), PolY(), AllPol()]
+SPINS = [SpinUp(), SpinDown(), AllSpin()]
+
+function _groundtruth_is_physical(proc::ScatteringProcess, ::PerturbativeQED)
+    incoming_electrons = number_particles(proc, Incoming(), Electron())
+    incoming_positrons = number_particles(proc, Incoming(), Positron())
+    outgoing_electrons = number_particles(proc, Outgoing(), Electron())
+    outgoing_positrons = number_particles(proc, Outgoing(), Positron())
+
+    return incoming_electrons + outgoing_positrons ==
+           outgoing_electrons + incoming_positrons
+end
+
+function _groundtruth_spin_pols(particles)
+    return ntuple(
+        x -> is_fermion(particles[x]) ? AllSpin() : AllPolarization(), length(particles)
+    )
+end
+
+function _random_spin_pols(RNG, particles)
+    return ntuple(
+        x -> is_fermion(particles[x]) ? rand(RNG, SPINS) : rand(RNG, POLS),
+        length(particles),
+    )
+end