rubydragon/metagraphoptimization.jl
1
0
Fork 0
Code Issues 6 Pull Requests Actions Packages Projects Releases Wiki Activity

Correct qed implementation and test compton

Browse Source
This commit is contained in:
Anton Reinhard 2023-11-29 16:43:47 +01:00
parent ba0c75c8dc
commit 6314539f2c
5 changed files with 40 additions and 38 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

19
src/models/qed/compute.jl
View File

@ -35,8 +35,6 @@ function compute(
p3 = QED_conserve_momentum(data1.p, data2.p)
P3 = interaction_result(P1, P2)
println("Virtual particle: $(P3(p3))")
state = QED_vertex()
if (typeof(data1.v) <: AdjointBiSpinor)
state = data1.v * state
@ -70,16 +68,11 @@ function compute(
P1 <: Union{AntiFermionStateful, FermionStateful},
P2 <: Union{AntiFermionStateful, FermionStateful},
}
println("S2: $(direction(data1.p)) $(data1.p) - $(direction(data2.p)) $(data2.p)")
# TODO: assert that data1 and data2 are opposites
#=@assert isapprox(abs(data1.p.momentum.E), abs(data2.p.momentum.E), rtol = 0.001, atol = sqrt(eps())) "E: $(data1.p.momentum.E) vs. $(data2.p.momentum.E)"
@assert isapprox(data1.p.momentum.px, -data2.p.momentum.px, rtol = 0.001, atol = sqrt(eps())) "px: $(data1.p.momentum.px) vs. $(data2.p.momentum.px)"
@assert isapprox(data1.p.momentum.py, -data2.p.momentum.py, rtol = 0.001, atol = sqrt(eps())) "py: $(data1.p.momentum.py) vs. $(data2.p.momentum.py)"
@assert isapprox(data1.p.momentum.pz, -data2.p.momentum.pz, rtol = 0.001, atol = sqrt(eps())) "pz: $(data1.p.momentum.pz) vs. $(data2.p.momentum.pz)"
=#
@assert isapprox(data1.p.momentum, data2.p.momentum, rtol = sqrt(eps()), atol = sqrt(eps())) "$(data1.p.momentum) vs. $(data2.p.momentum)"
inner = QED_inner_edge(data2.p)
# inner edge is just a scalar, data1 and data2 are bispinor/adjointbispinnor, need to keep correct order
inner = QED_inner_edge(propagation_result(P1)(data1.p))
# inner edge is just a "scalar", data1 and data2 are bispinor/adjointbispinnor, need to keep correct order
if typeof(data1.v) <: BiSpinor
return data2.v * inner * data1.v
else
@ -107,9 +100,9 @@ Compute inner edge (1 input particle, 1 output particle).
11 FLOP.
"""
function compute(::ComputeTaskQED_S1, data::QEDParticleValue{P})::QEDParticleValue where {P <: QEDParticle}
newP = propagation_result(P)
new_p = propagation_result(P)(data.p)
# inner edge is just a scalar, can multiply from either side
return QEDParticleValue{newP}(newP(data.p), data.v * QED_inner_edge(data.p))
return QEDParticleValue{newP}(new_p, data.v * QED_inner_edge(new_p))
end
"""

2
src/models/qed/create.jl
View File

@ -46,7 +46,7 @@ function gen_process_input(processDescription::QEDProcessDescription)
for (particle, n) in processDescription.outParticles
for _ in 1:n
mom = final_momenta[index]
push!(outputParticles, particle(SFourMomentum(-mom.E, mom.px, mom.py, mom.pz)))
push!(outputParticles, particle(SFourMomentum(mom.E, mom.px, mom.py, mom.pz)))
index += 1
end
end

23
src/models/qed/particle.jl
View File

@ -1,6 +1,9 @@
using QEDprocesses
import QEDbase.mass
# TODO check
const e = sqrt(4π / 137)
"""
QEDModel <: AbstractPhysicsModel
@ -60,7 +63,7 @@ A photon of the [`QEDModel`](@ref) with its state.
struct PhotonStateful{Direction <: ParticleDirection} <: QEDParticle{Direction}
momentum::SFourMomentum
# this will maybe change to the full polarization vector? or do i need both
polarization::AbstractPolarization
polarization::AbstractDefinitePolarization
end
PhotonStateful{Direction}(mom::SFourMomentum) where {Direction <: ParticleDirection} =
@ -76,7 +79,7 @@ A fermion of the [`QEDModel`](@ref) with its state.
"""
struct FermionStateful{Direction <: ParticleDirection} <: QEDParticle{Direction}
momentum::SFourMomentum
spin::AbstractSpin
spin::AbstractDefiniteSpin
# TODO: mass for electron/muon/tauon representation?
end
@ -93,7 +96,7 @@ An anti-fermion of the [`QEDModel`](@ref) with its state.
"""
struct AntiFermionStateful{Direction <: ParticleDirection} <: QEDParticle{Direction}
momentum::SFourMomentum
spin::AbstractSpin
spin::AbstractDefiniteSpin
# TODO: mass for electron/muon/tauon representation?
end
@ -210,6 +213,14 @@ end
@inline mass(::Type{<:AntiFermionStateful}) = 1.0
@inline mass(::Type{<:PhotonStateful}) = 0.0
@inline invert_momentum(p::FermionStateful{Dir}) where {Dir <: ParticleDirection} =
FermionStateful{Dir}(-p.momentum, p.spin)
@inline invert_momentum(p::AntiFermionStateful{Dir}) where {Dir <: ParticleDirection} =
AntiFermionStateful{Dir}(-p.momentum, p.spin)
@inline invert_momentum(k::PhotonStateful{Dir}) where {Dir <: ParticleDirection} =
PhotonStateful{Dir}(-k.momentum, k.polarization)
"""
caninteract(T1::Type{<:QEDParticle}, T2::Type{<:QEDParticle})
@ -272,11 +283,13 @@ end
Return the factor of a vertex in a QED feynman diagram.
"""
@inline function QED_vertex()::SLorentzVector{DiracMatrix}
return ComplexF64(0, 1) * gamma()
# Peskin-Schroeder notation
return -1im * e * gamma()
end
@inline function QED_inner_edge(p::QEDParticle)
return propagator(particle(p), p.momentum)
pos_mom = p.momentum
return propagator(particle(p), pos_mom)
end
"""

7
test/runtests.jl
View File

@ -1,5 +1,5 @@
using SafeTestsets
#=
@safetestset "Utility Unit Tests" begin
include("unit_tests_utility.jl")
end
@ -17,10 +17,10 @@ end
end
@safetestset "ABC-Model Unit Tests" begin
include("unit_tests_abcmodel.jl")
end=#
end
@safetestset "QED-Model Unit Tests" begin
include("unit_tests_qedmodel.jl")
end#=
end
@safetestset "Node Reduction Unit Tests" begin
include("node_reduction.jl")
end
@ -36,4 +36,3 @@ end
@safetestset "Known Graph Tests" begin
include("known_graphs.jl")
end
=#

27
test/unit_tests_qedmodel.jl
View File

@ -47,23 +47,22 @@ function compton_groundtruth(input::QEDProcessInput)
u_p1 = base_state(Electron(), Incoming(), p1.momentum, spin_or_pol(p1))
u_p2 = base_state(Electron(), Outgoing(), p2.momentum, spin_or_pol(p2))
eps_slashed_1 = base_state(Photon(), Incoming(), k1.momentum, spin_or_pol(k1))
eps_slashed_2 = base_state(Photon(), Outgoing(), k2.momentum, spin_or_pol(k2))
eps_1 = base_state(Photon(), Incoming(), k1.momentum, spin_or_pol(k1))
eps_2 = base_state(Photon(), Outgoing(), k2.momentum, spin_or_pol(k2))
virt1_mom = p2.momentum - k1.momentum
@test isapprox(p1.momentum - k2.momentum, virt1_mom)
virt2_mom = p1.momentum + k1.momentum
println("Groundtruth virtual particle (p2 - k1): $(virt1_mom)")
println("Groundtruth virtual particle (p1 + k1): $(virt2_mom)")
@test isapprox(p2.momentum + k2.momentum, virt2_mom)
s_p2_k1 = propagator(Electron(), virt1_mom)
s_p1_k1 = propagator(Electron(), virt2_mom)
diagram1 = u_p2 * eps_slashed_1 * QED_vertex() * s_p2_k1 * QED_vertex() * eps_slashed_2 * u_p1
diagram2 = u_p2 * eps_slashed_2 * QED_vertex() * s_p1_k1 * QED_vertex() * eps_slashed_1 * u_p1
diagram1 = u_p2 * (eps_1 * QED_vertex()) * s_p2_k1 * (eps_2 * QED_vertex()) * u_p1
diagram2 = u_p2 * (eps_2 * QED_vertex()) * s_p1_k1 * (eps_1 * QED_vertex()) * u_p1
return (diagram1 + diagram2)
return diagram1 + diagram2
end
@ -165,8 +164,8 @@ end
@test countmap(typeof.(input.outParticles)) == process.outParticles
@test isapprox(
sum(getfield.(input.inParticles, :momentum)) + sum(getfield.(input.outParticles, :momentum)),
SFourMomentum(0.0, 0.0, 0.0, 0.0);
sum(getfield.(input.inParticles, :momentum)),
sum(getfield.(input.outParticles, :momentum));
atol = sqrt(eps()),
)
end
@ -282,9 +281,7 @@ end
compton_function = get_compute_function(graph, process, machine)
input = gen_process_input(process)
input = [gen_process_input(process) for _ in 1:1000]
println("Function result: $(compton_function(input))")
println("Groundtruth: $(compton_groundtruth(input))")
@test isapprox(compton_function.(input), compton_groundtruth.(input))
end