diff --git a/Project.toml b/Project.toml index 156915a..d2a9c16 100644 --- a/Project.toml +++ b/Project.toml @@ -1,12 +1,6 @@ name = "QEDprocesses" uuid = "46de9c38-1bb3-4547-a1ec-da24d767fdad" -authors = [ - "Uwe Hernandez Acosta ", - "Simeon Ehrig", - "Klaus Steiniger", - "Tom Jungnickel", - "Anton Reinhard", -] +authors = ["Uwe Hernandez Acosta ", "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" diff --git a/src/QEDprocesses.jl b/src/QEDprocesses.jl index 268d5ab..819461f 100644 --- a/src/QEDprocesses.jl +++ b/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 diff --git a/src/models/perturbative_qed.jl b/src/models/perturbative_qed.jl index b045619..add5b9f 100644 --- a/src/models/perturbative_qed.jl +++ b/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 diff --git a/src/patch_QEDbase.jl b/src/patch_QEDbase.jl index 3da67e1..cbe50a6 100644 --- a/src/patch_QEDbase.jl +++ b/src/patch_QEDbase.jl @@ -1,3 +1,4 @@ +# TODO: remove this! ############# # Patches for `QEDbase.jl` # remove if this went into `QEDbase.jl` diff --git a/src/processes/generic_process/process.jl b/src/processes/generic_process/process.jl index c642fce..be14ea7 100644 --- a/src/processes/generic_process/process.jl +++ b/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. diff --git a/src/processes/one_photon_compton/perturbative/cross_section.jl b/src/processes/one_photon_compton/perturbative/cross_section.jl index 5464339..f8a89a9 100644 --- a/src/processes/one_photon_compton/perturbative/cross_section.jl +++ b/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())) diff --git a/src/processes/one_photon_compton/perturbative/kinematics.jl b/src/processes/one_photon_compton/perturbative/kinematics.jl index b28ac1e..ab72f65 100644 --- a/src/processes/one_photon_compton/perturbative/kinematics.jl +++ b/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 diff --git a/src/processes/one_photon_compton/perturbative/total_probability.jl b/src/processes/one_photon_compton/perturbative/total_probability.jl index 838a730..a61e68f 100644 --- a/src/processes/one_photon_compton/perturbative/total_probability.jl +++ b/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 diff --git a/src/utils.jl b/src/utils.jl index 97eb6a2..4ffdf8f 100644 --- a/src/utils.jl +++ b/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 diff --git a/test/processes/one_photon_compton/perturbative.jl b/test/processes/one_photon_compton/perturbative/cross_section.jl similarity index 61% rename from test/processes/one_photon_compton/perturbative.jl rename to test/processes/one_photon_compton/perturbative/cross_section.jl index c7abf1f..4c16a88 100644 --- a/test/processes/one_photon_compton/perturbative.jl +++ b/test/processes/one_photon_compton/perturbative/cross_section.jl @@ -1,4 +1,3 @@ - using QEDbase using QEDcore using QEDprocesses @@ -13,43 +12,13 @@ const RTOL = sqrt(eps()) include("groundtruths.jl") const MODEL = PerturbativeQED() -const PS_DEF = PhasespaceDefinition(SphericalCoordinateSystem(), ElectronRestFrame()) -const OMEGAS = (1e-6 * rand(RNG), 1e-3 * rand(RNG), rand(RNG), 1e3 * rand(RNG)) +const IN_PSL= ComptonRestSystem(Energy(2)) +const OUT_PSL= ComptonSphericalLayout(IN_PSL) +const OMEGAS = (1e-6 * rand(RNG), 1e-3 * rand(RNG), rand(RNG), 1e3 * rand(RNG)) const COS_THETAS = [-1.0, 2 * rand(RNG) - 1, 0.0, 1.0] const PHIS = [0, 2 * pi, rand(RNG) * 2 * pi] -@testset "pretty-printing" begin - buf = IOBuffer() - print(buf, MODEL) - @test String(take!(buf)) == "perturbative QED" - - show(buf, MIME"text/plain"(), MODEL) - @test String(take!(buf)) == "perturbative QED" -end - -@testset "perturbative kinematics" begin - PROC = Compton() - @testset "momentum generation" begin - @testset "$om, $cth, $phi" for (om, cth, phi) in - Iterators.product(OMEGAS, COS_THETAS, PHIS) - IN_COORDS = (om,) - OUT_COORDS = (cth, phi) - IN_PS, OUT_PS = QEDbase._generate_momenta( - PROC, MODEL, PS_DEF, IN_COORDS, OUT_COORDS - ) - in_mom_square = getMass2.(IN_PS) - out_mom_square = getMass2.(OUT_PS) - in_masses = mass.(incoming_particles(PROC)) .^ 2 - out_masses = mass.(outgoing_particles(PROC)) .^ 2 - - # we need a larger ATOL than eps() here because the error is accumulated over several additions - @test all(isapprox.(in_mom_square, in_masses, atol=4 * ATOL, rtol=RTOL)) - @test all(isapprox.(out_mom_square, out_masses, atol=4 * ATOL, rtol=RTOL)) - end - end -end - @testset "perturbative" begin @testset "$omega" for omega in OMEGAS @testset "differential cross section" begin @@ -60,16 +29,14 @@ end Iterators.product(COS_THETAS, PHIS) IN_COORDS = (omega,) OUT_COORDS = (cos_theta, phi) - IN_PS, OUT_PS = QEDbase._generate_momenta( - PROC, MODEL, PS_DEF, IN_COORDS, OUT_COORDS - ) - groundtruth = _groundtruth_pert_compton_diffCS_spinsum_polsum( + PSP = PhaseSpacePoint(PROC, MODEL, OUT_PSL, IN_COORDS, OUT_COORDS) + test_val = unsafe_differential_cross_section(PSP) + + groundtruth = _groundtruth_pert_compton_diffCS_spinsum_polsum_elab_sph( omega, cos_theta, 1.0 ) - PSP = PhaseSpacePoint(PROC, MODEL, PS_DEF, IN_PS, OUT_PS) - test_val = unsafe_differential_cross_section(PSP) @test isapprox(test_val, groundtruth, atol=ATOL, rtol=RTOL) end @@ -82,17 +49,13 @@ end Iterators.product(COS_THETAS, PHIS) IN_COORDS = (omega,) OUT_COORDS = (cos_theta, phi) - IN_PS, OUT_PS = QEDbase._generate_momenta( - PROC, MODEL, PS_DEF, IN_COORDS, OUT_COORDS - ) + PSP = PhaseSpacePoint(PROC, MODEL, OUT_PSL, IN_COORDS, OUT_COORDS) + test_val = unsafe_differential_cross_section(PSP) - groundtruth = _groundtruth_pert_compton_diffCS_spinsum_xpol( + groundtruth = _groundtruth_pert_compton_diffCS_spinsum_xpol_elab_sph( omega, cos_theta, phi, 1.0 ) - PSP = PhaseSpacePoint(PROC, MODEL, PS_DEF, IN_PS, OUT_PS) - test_val = unsafe_differential_cross_section(PSP) - @test isapprox(test_val, groundtruth, atol=ATOL, rtol=RTOL) end end @@ -104,17 +67,13 @@ end Iterators.product(COS_THETAS, PHIS) IN_COORDS = (omega,) OUT_COORDS = (cos_theta, phi) - IN_PS, OUT_PS = QEDbase._generate_momenta( - PROC, MODEL, PS_DEF, IN_COORDS, OUT_COORDS - ) + PSP = PhaseSpacePoint(PROC, MODEL, OUT_PSL, IN_COORDS, OUT_COORDS) + test_val = unsafe_differential_cross_section(PSP) - groundtruth = _groundtruth_pert_compton_diffCS_spinsum_ypol( + groundtruth = _groundtruth_pert_compton_diffCS_spinsum_ypol_elab_sph( omega, cos_theta, phi, 1.0 ) - PSP = PhaseSpacePoint(PROC, MODEL, PS_DEF, IN_PS, OUT_PS) - test_val = unsafe_differential_cross_section(PSP) - @test isapprox(test_val, groundtruth, atol=ATOL, rtol=RTOL) end end @@ -126,9 +85,7 @@ end function klein_nishina_total_cross_section(in_ps) function func(x) return unsafe_differential_cross_section( - PhaseSpacePoint( - Compton(), PerturbativeQED(), PS_DEF, in_ps, (x, 0.0) - ), + PhaseSpacePoint(Compton(), MODEL, OUT_PSL, in_ps, (x, 0.0)) ) end res, err = quadgk(func, -1, 1) @@ -140,8 +97,9 @@ end IN_COORDS = (omega,) groundtruth = klein_nishina_total_cross_section(IN_COORDS) test_val = @inferred QEDprocesses.total_cross_section( - InPhaseSpacePoint(PROC, MODEL, PS_DEF, IN_COORDS) + InPhaseSpacePoint(PROC, MODEL, IN_PSL, IN_COORDS) ) + @test isapprox(test_val, groundtruth, atol=ATOL, rtol=RTOL) @testset "$cos_theta $phi" for (cos_theta, phi) in @@ -149,16 +107,16 @@ end OUT_COORDS = (cos_theta, phi) test_val = @inferred QEDprocesses.total_cross_section( - PhaseSpacePoint(PROC, MODEL, PS_DEF, IN_COORDS, OUT_COORDS) + PhaseSpacePoint(PROC, MODEL, OUT_PSL, IN_COORDS, OUT_COORDS) ) @test isapprox(test_val, groundtruth, atol=ATOL, rtol=RTOL) out_moms = momenta( - PhaseSpacePoint(PROC, MODEL, PS_DEF, IN_COORDS, OUT_COORDS), + PhaseSpacePoint(PROC, MODEL, OUT_PSL, IN_COORDS, OUT_COORDS), Outgoing(), ) @test_throws MethodError QEDprocesses.total_cross_section( - OutPhaseSpacePoint(PROC, MODEL, PS_DEF, out_moms) + OutPhaseSpacePoint(PROC, MODEL, OUT_PSL, out_moms) ) end end diff --git a/test/processes/one_photon_compton/groundtruths.jl b/test/processes/one_photon_compton/perturbative/groundtruths.jl similarity index 50% rename from test/processes/one_photon_compton/groundtruths.jl rename to test/processes/one_photon_compton/perturbative/groundtruths.jl index 3505733..decee93 100644 --- a/test/processes/one_photon_compton/groundtruths.jl +++ b/test/processes/one_photon_compton/perturbative/groundtruths.jl @@ -1,4 +1,8 @@ +function _pert_omega_prime_elab_sph(om,cth) + om/(1+om*(1-cth)) +end + """ Klein-Nishina: differential cross section @@ -6,21 +10,21 @@ source: Peskin, Schroeder. "Quantum field theory." (1995). eq: 5.91 note: we compute d sigma/ d Omega, and *not* d sigma/ d cos theta (the difference is a factor 2 pi) """ -function _groundtruth_pert_compton_diffCS_spinsum_polsum(om, cth, mass) +function _groundtruth_pert_compton_diffCS_spinsum_polsum_elab_sph(om, cth, mass) prefac = ALPHA_SQUARE / 2 - omp = QEDprocesses._pert_omega_prime(om, cth) + omp = _pert_omega_prime_elab_sph(om, cth) sth_sq = 1 - cth^2 return prefac * (omp / om)^2 * (omp / om + om / omp - sth_sq) end -function _groundtruth_pert_compton_diffCS_spinsum_xpol(omega, ctheta, phi, mass) - om_prime = QEDprocesses._pert_omega_prime(omega, ctheta) +function _groundtruth_pert_compton_diffCS_spinsum_xpol_elab_sph(omega, ctheta, phi, mass) + om_prime = _pert_omega_prime_elab_sph(omega, ctheta) om_prime_over_om = om_prime / omega return 0.5 * ALPHA_SQUARE / mass^2 * om_prime_over_om^2 * (om_prime_over_om + 1.0 / om_prime_over_om - 2 * (1 - ctheta^2) * cos(phi)^2) end -function _groundtruth_pert_compton_diffCS_spinsum_ypol(omega, ctheta, phi, mass) - return _groundtruth_pert_compton_diffCS_spinsum_xpol(omega, ctheta, phi + pi / 2, mass) +function _groundtruth_pert_compton_diffCS_spinsum_ypol_elab_sph(omega, ctheta, phi, mass) + return _groundtruth_pert_compton_diffCS_spinsum_xpol_elab_sph(omega, ctheta, phi + pi / 2, mass) end diff --git a/test/processes/one_photon_compton/perturbative/kinematics.jl b/test/processes/one_photon_compton/perturbative/kinematics.jl new file mode 100644 index 0000000..2f928de --- /dev/null +++ b/test/processes/one_photon_compton/perturbative/kinematics.jl @@ -0,0 +1,125 @@ + +using QEDbase +using QEDcore +using QEDprocesses +using Random + +const RNG = MersenneTwister(77697185) +const ATOL = eps() +const RTOL = sqrt(eps()) + +const PROC = Compton() +const MODEL = PerturbativeQED() + +const OMEGAS = (1e-6 * rand(RNG), 1e-3 * rand(RNG), rand(RNG), 1e3 * rand(RNG)) +const SQRT_S = (1.0, 1 + rand(RNG)) + +const COS_THETAS = [-1.0, 2 * rand(RNG) - 1, 0.0, 1.0] +const PHIS = [0, 2 * pi, rand(RNG) * 2 * pi] + +@testset "in-phase-space layout" begin + @testset "Compton rest frame" begin + @testset "default" begin + @test ComptonRestSystem() == ComptonRestSystem(Energy(2)) + end + + @testset "wrong rest system" begin + @testset "$coord_type" for coord_type in (Energy, SpatialMagnitude) + coord = coord_type(1) + @test_throws ArgumentError ComptonRestSystem(coord) + end + + @testset "no photon rapidity" begin + coord = Rapidity(2) + @test_throws ArgumentError ComptonRestSystem(coord) + end + end + + @testset "photon energy" begin + in_psl = ComptonRestSystem(Energy(2)) + + @testset "$om" for om in OMEGAS + in_coords = (om,) + in_psp = InPhaseSpacePoint(PROC, MODEL, in_psl, in_coords) + + test_P, test_K = momenta(in_psp, Incoming()) + + @test isapprox(getMass2(test_K), mass(Photon())^2, atol=ATOL, rtol=RTOL) + @test isapprox(getMass2(test_P), mass(Electron())^2, atol=ATOL, rtol=RTOL) + @test isapprox(getE(test_K), om, atol=ATOL, rtol=RTOL) + @test isapprox(getE(test_P), mass(Electron()), atol=ATOL, rtol=RTOL) + end + end + + @testset "photon magnitude" begin + in_psl = ComptonRestSystem(SpatialMagnitude(2)) + + # for photons: E === rho + @testset "$om" for om in OMEGAS + in_coords = (om,) + in_psp = InPhaseSpacePoint(PROC, MODEL, in_psl, in_coords) + + test_P, test_K = momenta(in_psp, Incoming()) + + @test isapprox(getMass2(test_K), mass(Photon())^2, atol=ATOL, rtol=RTOL) + @test isapprox(getMass2(test_P), mass(Electron())^2, atol=ATOL, rtol=RTOL) + @test isapprox(getE(test_K), om, atol=ATOL, rtol=RTOL) + @test isapprox(getE(test_P), mass(Electron()), atol=ATOL, rtol=RTOL) + end + end + + @testset "cms energy" begin + in_psl = ComptonRestSystem(CMSEnergy()) + + @testset "$ss" for ss in SQRT_S + in_coords = (ss,) + in_psp = InPhaseSpacePoint(PROC, MODEL, in_psl, in_coords) + + test_P, test_K = momenta(in_psp, Incoming()) + + @test isapprox(getMass2(test_K), mass(Photon())^2, atol=ATOL, rtol=RTOL) + @test isapprox(getMass2(test_P), mass(Electron())^2, atol=ATOL, rtol=RTOL) + @test isapprox(getMass(test_K + test_P), ss, atol=ATOL, rtol=RTOL) + @test isapprox(getE(test_P), mass(Electron()), atol=ATOL, rtol=RTOL) + end + end + end +end + +@testset "out-phase-space-layout" begin + @testset "Compton spherical system" begin + in_psl = ComptonRestSystem(Energy(2)) + out_psl = ComptonSphericalLayout(in_psl) + + @testset "wrong in_psl" begin + in_psl_fail = TwoBodyBeamSystem() + @test_throws ArgumentError ComptonSphericalLayout(in_psl_fail) + end + + @testset "$om, $cth, $phi" for (om, cth, phi) in + Iterators.product(OMEGAS, COS_THETAS, PHIS) + IN_COORDS = (om,) + OUT_COORDS = (cth, phi) + + test_psp = PhaseSpacePoint(PROC, MODEL, out_psl, IN_COORDS, OUT_COORDS) + @testset "mass shell" begin + @testset "$dir $part" for (dir, part) in Iterators.product( + (Incoming(), Outgoing()), (Electron(), Photon()) + ) + mom = momentum(test_psp, dir, part) + mom_square = mom * mom + mass_square = mass(part)^2 + + @test isapprox(mom_square, mass_square, atol=1e2 * ATOL, rtol=RTOL) + end + end + + @testset "energy-momentum conservation" begin + in_moms = momenta(test_psp, Incoming()) + out_moms = momenta(test_psp, Outgoing()) + + @test isapprox(sum(in_moms), sum(out_moms), atol=ATOL, rtol=RTOL) + end + end + end +end diff --git a/test/processes/run_process_test.jl b/test/processes/run_process_test.jl index d8bf4a3..179a5fb 100644 --- a/test/processes/run_process_test.jl +++ b/test/processes/run_process_test.jl @@ -6,6 +6,12 @@ end include("one_photon_compton/process.jl") end -@time @safetestset "perturbative one photon compton" begin - include("one_photon_compton/perturbative.jl") +@testset "perturbative one photon compton" begin + @safetestset "kinematics" begin + include("one_photon_compton/perturbative/kinematics.jl") + end + + @safetestset "cross section" begin + include("one_photon_compton/perturbative/cross_section.jl") + end end