updated psp usage

Project.toml

@@ -1,12 +1,6 @@
 name = "QEDprocesses"
 uuid = "46de9c38-1bb3-4547-a1ec-da24d767fdad"
-authors = [
-    "Uwe Hernandez Acosta <[email protected]>",
-    "Simeon Ehrig",
-    "Klaus Steiniger",
-    "Tom Jungnickel",
-    "Anton Reinhard",
-]
+authors = ["Uwe Hernandez Acosta <[email protected]>", "Simeon Ehrig", "Klaus Steiniger", "Tom Jungnickel", "Anton Reinhard"]
 version = "0.3.0"
 
 [deps]
@@ -16,8 +10,8 @@ QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
-QEDbase = "0.3"
-QEDcore = "0.2"
+QEDbase = "0.3.0"
+QEDcore = "0.2.0"
 QuadGK = "2"
 StaticArrays = "1"
 julia = "1.10"

src/QEDprocesses.jl

@@ -9,6 +9,10 @@ export PerturbativeQED
 
 # specific scattering processes
 export Compton, omega_prime
+export ComptonRestSystem
+export ComptonSphericalLayout
+
+# generic scattering process
 export ScatteringProcess, isphysical
 
 using QEDbase

src/models/perturbative_qed.jl

@@ -1,16 +1,16 @@
 
-struct PerturbativeQED <: AbstractModelDefinition end
+struct PerturbativeQED <: AbstractPerturbativeModel end
 
 QEDbase.fundamental_interaction_type(::PerturbativeQED) = :electromagnetic
 
 """
     in_phase_space_dimension(proc::AbstractProcessDefinition, ::PerturbativeQED)
 
-Return the number of degrees of freedom to determine the incoming phase space for processes in PerturbativeQED. 
+Return the number of degrees of freedom to determine the incoming phase space for processes in PerturbativeQED.
 
 !!! note "Convention"
 
-    The current implementation only supports the case where two of the incoming particles collide head-on. 
+    The current implementation only supports the case where two of the incoming particles collide head-on.
 """
 function QEDbase.in_phase_space_dimension(
     proc::AbstractProcessDefinition, ::PerturbativeQED

src/patch_QEDbase.jl

@@ -1,3 +1,4 @@
+# TODO: remove this!
 #############
 # Patches for `QEDbase.jl`
 # remove if this went into `QEDbase.jl`

src/processes/generic_process/process.jl

@@ -1,12 +1,12 @@
 """
     ScatteringProcess <: AbstractProcessDefinition
 
-Generic implementation for scattering processes of arbitrary particles. Currently, only calculations in combination with `PerturbativeQED` are supported. 
+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  
+!!! warning
     The computation of cross sections and probabilities is currently unimplemented.
 
 ## Constructors
@@ -19,7 +19,7 @@ The [`isphysical`](@ref) function can be used to check whether the process is po
     )
 
 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 
+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.

src/processes/one_photon_compton/perturbative/cross_section.jl

@@ -56,7 +56,7 @@ end
 )
     in_ps = momenta(psp, Incoming())
     out_ps = momenta(psp, Outgoing())
-    return _pert_compton_ps_fac(psp.ps_def, in_ps[2], out_ps[2])
+    return _pert_compton_ps_fac(psp.psl, in_ps[2], out_ps[2])
 end
 
 #######
@@ -163,9 +163,8 @@ end
 #######
 
 function _pert_compton_ps_fac(
-    in_ps_def::PhasespaceDefinition{inCS,ElectronRestFrame}, in_photon_mom, out_photon_mom
-) where {inCS}
-    # TODO
+    in_psl::ComptonSphericalLayout{<:ComptonRestSystem}, in_photon_mom, out_photon_mom
+)
     omega = getE(in_photon_mom)
     omega_prime = getE(out_photon_mom)
     return omega_prime^2 / (16 * pi^2 * omega * mass(Electron()))

src/processes/one_photon_compton/perturbative/kinematics.jl

@@ -1,54 +1,84 @@
-@inline function _pert_omega_prime(omega, cth; mass=1.0)
-    return omega / (1 + omega / mass * (1 - cth))
+
+# incoming phase space layout
+# alias for electron rest system (based on QEDcore.TwoBodyRestSystem)
+
+const ComptonRestSystem{COORD} = TwoBodyRestSystem{1,COORD} where {COORD}
+ComptonRestSystem(coord) = ComptonRestSystem(Val(particle_index(coord)), coord)
+function ComptonRestSystem(
+    run_idx::Val{1}, coord::COORD
+) where {COORD<:QEDcore.AbstractSingleParticleCoordinate}
+    throw(
+        ArgumentError(
+            "the first incoming particle of Compton is an electron, which has no degrees of freedom in the Compton rest frame.",
+        ),
+    )
+end
+function ComptonRestSystem(
+    run_idx::Val{2}, coord::COORD
+) where {COORD<:QEDcore.AbstractSingleParticleCoordinate{2}}
+    return TwoBodyRestSystem{1}(coord)
+end
+ComptonRestSystem(coord::CMSEnergy) = TwoBodyRestSystem{1}(coord)
+function ComptonRestSystem(::Rapidity)
+    throw(
+        ArgumentError(
+            "there is no finite rapidity for a photon in the electron rest system"
+        ),
+    )
 end
+ComptonRestSystem() = ComptonRestSystem(Energy(2))
 
-function generate_momenta(
-    proc::Compton,
-    model::PerturbativeQED,
-    in_ps_def::PhasespaceDefinition{SphericalCoordinateSystem,ElectronRestFrame},
-    in_ps::NTuple{N,T},
-    out_ps::NTuple{M,T},
-) where {N,M,T<:Real}
-    return QEDbase._generate_momenta(proc, model, in_ps_def, in_ps, out_ps)
+# outgoing phase space layout
+# spherical coordinates in electron rest frame
+
+@inline function _pert_Compton_omega_prime(pt, cth)
+    Et = getE(pt)
+    rho2_t = getMag2(pt)
+    s = Et^2 - rho2_t
+
+    return (s - 1) / (2 * (Et - sqrt(rho2_t) * cth))
 end
 
-function QEDbase._generate_incoming_momenta(
-    proc::Compton,
-    model::PerturbativeQED,
-    in_ps_def::PhasespaceDefinition{SphericalCoordinateSystem,ElectronRestFrame},
-    in_ps::NTuple{N,T},
-) where {N,T<:Real}
-    om = in_ps[1]
+struct ComptonSphericalLayout{INPSL<:AbstractTwoBodyInPhaseSpaceLayout} <:
+       AbstractOutPhaseSpaceLayout{INPSL}
+    in_psl::INPSL
+end
 
-    P = SFourMomentum(one(om), zero(om), zero(om), zero(om))
-    K = SFourMomentum(om, zero(om), zero(om), om)
+function ComptonSphericalLayout(::TwoBodyRestSystem{2})
+    throw(
+        ArgumentError(
+            "the first incoming particle of Compton is an electron, which has no degrees of freedom in the Compton rest frame.",
+        ),
+    )
+end
 
-    return P, K
+function QEDbase.phase_space_dimension(
+    proc::Compton, model::PerturbativeQED, psl::ComptonSphericalLayout
+)
+    # cth, phi
+    return 2
 end
 
-function QEDbase._generate_momenta(
+QEDbase.in_phase_space_layout(psl::ComptonSphericalLayout) = psl.in_psl
+
+function QEDbase._build_momenta(
     proc::Compton,
     model::PerturbativeQED,
-    in_ps_def::PhasespaceDefinition{SphericalCoordinateSystem,ElectronRestFrame},
-    in_ps::NTuple{N,T},
-    out_ps::NTuple{M,T},
-) where {N,M,T<:Real}
-    omega = in_ps[1]
-    cth = out_ps[1]
-    phi = out_ps[2]
-    P, K, Pp, Kp = _generate_momenta_elab_sph(omega, cth, phi) # TODO: do this coord and frame dependent
-    in_moms = (P, K)
-    out_moms = (Pp, Kp)
-    return in_moms, out_moms
-end
-
-function _generate_momenta_elab_sph(om, cth, phi, m=1.0)
-    P = SFourMomentum(m, zero(m), zero(m), zero(m))
-    K = SFourMomentum(om, zero(om), zero(om), om)
-    omp = _pert_omega_prime(om, cth)
+    psl::ComptonSphericalLayout,
+    in_coords::NTuple{1,T},
+    out_coords::NTuple{2,T},
+) where {T<:Real}
+    P, K = QEDbase._build_momenta(proc, model, in_phase_space_layout(psl), in_coords)
+    Pt = P + K
+    cth, phi = @inbounds out_coords
+    omega_prime = _pert_Compton_omega_prime(Pt, cth)
     sth = sqrt(1 - cth^2)
     sphi, cphi = sincos(phi)
-    Kp = SFourMomentum(omp, omp * sth * cphi, omp * sth * sphi, omp * cth)
-    Pp = P + K - Kp
-    return P, K, Pp, Kp
+
+    Kp = SFourMomentum(
+        omega_prime, omega_prime * sth * cphi, omega_prime * sth * sphi, omega_prime * cth
+    )
+    Pp = Pt - Kp
+
+    return (P, K), (Pp, Kp)
 end

src/processes/one_photon_compton/perturbative/total_probability.jl

@@ -1,9 +1,14 @@
+
+_build_sph_out_psl(psl::ComptonSphericalLayout) = psl
+_build_sph_out_psl(psl::AbstractTwoBodyInPhaseSpaceLayout) = ComptonSphericalLayout(psl)
+
+
 function QEDbase._total_probability(in_psp::InPhaseSpacePoint{<:Compton,PerturbativeQED})
     omega = getE(momentum(in_psp[Incoming(), 2]))
 
     function func(x)
         return unsafe_differential_probability(
-            PhaseSpacePoint(in_psp.proc, in_psp.model, in_psp.ps_def, (omega,), (x, 0.0))
+            PhaseSpacePoint(in_psp.proc, in_psp.model, _build_sph_out_psl(in_psp.psl), (omega,), (x, 0.0))
         )
     end
 

src/utils.jl

@@ -8,7 +8,7 @@
 
     _base_component_type(array_of_lv::AbstractArray{LV}) where {LV<:AbstractLorentzVector}
 
-Return the type of the components of given Lorentz vectors, which are by themself elements of an 
+Return the type of the components of given Lorentz vectors, which are by themself elements of an
 `AbstractArray`.
 
 # Examples