perturbative Compton

Project.toml

@@ -5,6 +5,7 @@ version = "0.1.0"
 
 [deps]
 QEDbase = "10e22c08-3ccb-4172-bfcf-7d7aa3d04d93"
+QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]

src/QEDprocesses.jl

@@ -46,6 +46,7 @@ export Compton, omega_prime
 
 using QEDbase
 using StaticArrays
+using QuadGK
 
 include("constants.jl")
 include("utils.jl")

src/processes/one_photon_compton/one_photon_compton.jl

@@ -8,3 +8,4 @@
 include("process.jl")
 include("perturbative/kinematics.jl")
 include("perturbative/cross_section.jl")
+include("perturbative/total_probability.jl")

src/processes/one_photon_compton/perturbative/kinematics.jl

@@ -1,4 +1,3 @@
-
 @inline function _pert_omega_prime(omega, cth; mass=1.0)
     return omega / (1 + omega / mass * (1 - cth))
 end
@@ -13,6 +12,20 @@ function generate_momenta(
     return _generate_momenta(proc, model, in_ps_def, in_ps, out_ps)
 end
 
+function _generate_incoming_momenta(
+    proc::Compton,
+    model::PerturbativeQED,
+    in_ps_def::PhasespaceDefinition{SphericalCoordinateSystem,ElectronRestFrame},
+    in_ps::AbstractVector{T},
+) where {T<:Real}
+    om = in_ps[1]
+
+    P = SFourMomentum(one(om), zero(om), zero(om), zero(om))
+    K = SFourMomentum(om, zero(om), zero(om), om)
+
+    return [P, K]
+end
+
 function _generate_momenta(
     proc::Compton,
     model::PerturbativeQED,

src/processes/one_photon_compton/perturbative/total_probability.jl

@@ -0,0 +1,20 @@
+
+function _total_probability(
+    proc::Compton,
+    model::PerturbativeQED,
+    in_phase_space_def::AbstractPhasespaceDefinition,
+    in_phase_space::AbstractVector{T},
+) where {T<:QEDbase.AbstractFourMomentum}
+    omega = getE(in_phase_space[2])
+
+    function func(x)
+        return unsafe_differential_probability(
+            proc, model, in_phase_space_def, [omega], [x, 0.0]
+        )
+    end
+
+    tot_prob, _ = quadgk(func, -1, 1; rtol=sqrt(eps(omega)))
+
+    tot_prob *= 2 * pi # phi integration is trivial
+    return tot_prob
+end

test/processes/one_photon_compton/perturbative.jl

@@ -5,19 +5,19 @@ using Random
 using StaticArrays
 using QuadGK
 
-RNG = MersenneTwister(77697185)
-ATOL = eps()
-RTOL = sqrt(eps())
+const RNG = MersenneTwister(77697185)
+const ATOL = eps()
+const RTOL = sqrt(eps())
 
 include("groundtruths.jl")
 
-MODEL = PerturbativeQED()
-PS_DEF = PhasespaceDefinition(SphericalCoordinateSystem(), ElectronRestFrame())
-OMEGAS = (1e-6 * rand(RNG), 1e-3 * rand(RNG), rand(RNG), 1e3 * rand(RNG))
+const MODEL = PerturbativeQED()
+const PS_DEF = PhasespaceDefinition(SphericalCoordinateSystem(), ElectronRestFrame())
+const OMEGAS = (1e-6 * rand(RNG), 1e-3 * rand(RNG), rand(RNG), 1e3 * rand(RNG))
 #OMEGAS = (rand(RNG),)
 
-COS_THETAS = [-1.0, 2 * rand(RNG) - 1, 0.0, 1.0]
-PHIS = [0, 2 * pi, rand(RNG) * 2 * pi]
+const COS_THETAS = [-1.0, 2 * rand(RNG) - 1, 0.0, 1.0]
+const PHIS = [0, 2 * pi, rand(RNG) * 2 * pi]
 
 @testset "perturbative kinematics" begin
     PROC = Compton()
@@ -39,7 +39,7 @@ PHIS = [0, 2 * pi, rand(RNG) * 2 * pi]
     end
 end
 
-@testset "perturbative cross section" begin
+@testset "perturbative" begin
     @testset "$omega" for omega in OMEGAS
         @testset "differential cross section" begin
             @testset "spin and pol summed" begin
@@ -93,26 +93,27 @@ end
                 end
             end
         end
+        @testset "total cross section" begin
+            @testset "spin and pol summed" begin
+                PROC = Compton()
+                # Klein-Nishina: total cross section
+                function klein_nishina_total_cross_section(in_ps)
+                    function func(x)
+                        return unsafe_differential_cross_section(
+                            Compton(), PerturbativeQED(), PS_DEF, in_ps, [x, 0.0]
+                        )
+                    end
+                    res, err = quadgk(func, -1, 1)
+
+                    # note: mul by 2pi instead of the phi-integration
+                    return 2 * pi * res
+                end
+
+                IN_COORDS = [omega]
+                groundtruth = klein_nishina_total_cross_section(IN_COORDS)
+                test_val = @inferred total_cross_section(PROC, MODEL, PS_DEF, IN_COORDS)
+                @test isapprox(test_val, groundtruth, atol=ATOL, rtol=RTOL)
+            end
+        end
     end
-    # @testset "total cross section" begin
-    #     @testset "spin and pol summed" begin
-    #         PROC = Compton()
-    #         # Klein-Nishina: total cross section
-    #         function klein_nishina_total_cross_section(om, mass)
-    #             function func(x)
-    #                 return differential_cross_section_on_coord(
-    #                     Compton(), PerturbativeQED(), om, [x, 0.0]
-    #                 )
-    #             end
-    #             res, err = quadgk(func, -1, 1)
-    #
-    #             # note: mul by 2pi instead of the phi-integration
-    #             return 2 * pi * res
-    #         end
-    #
-    #         groundtruth = klein_nishina_total_cross_section(OMEGA, MASS)
-    #         test_val = @inferred total_cross_section_on_coord(PROC, MODEL, OMEGA)
-    #         @test isapprox(test_val, groundtruth, atol=ATOL, rtol=RTOL)
-    #     end
-    # end
 end