283 lines
9.2 KiB
Julia
283 lines
9.2 KiB
Julia
using Combinatorics
|
|
using QEDprocesses
|
|
using QEDbase
|
|
|
|
struct FeynmanDiagramDefinition
|
|
n::Int
|
|
end
|
|
|
|
struct FeynmanDiagramTopology
|
|
# largest subtree
|
|
subtree1::AbstractTree
|
|
# second largest subtree
|
|
subtree2::AbstractTree
|
|
# third largest subtree
|
|
subtree3::AbstractTree
|
|
end
|
|
|
|
struct TopologyPartition
|
|
leaves1::Int
|
|
leaves2::Int
|
|
leaves3::Int
|
|
end
|
|
|
|
function is_valid(partition::TopologyPartition)
|
|
if partition.leaves2 > partition.leaves1 || partition.leaves3 > partition.leaves2
|
|
return false
|
|
end
|
|
if partition.leaves1 > partition.leaves2 + partition.leaves3
|
|
return false
|
|
end
|
|
if partition.leaves1 <= 0 || partition.leaves2 <= 0 || partition.leaves3 <= 0
|
|
return false
|
|
end
|
|
return true
|
|
end
|
|
|
|
function leaves(partition::TopologyPartition)
|
|
return partition.leaves1 + partition.leaves2 + partition.leaves3
|
|
end
|
|
|
|
|
|
#
|
|
# Feynman Diagram, tree-level, QED
|
|
#
|
|
|
|
struct FeynmanDiagram{N,E,U,T,M,FM} <: AbstractTreeLevelFeynmanDiagram where {N,E,U,T,M,FM<:FlatMatrix}
|
|
diagram_structure::FM
|
|
|
|
electron_permutation::NTuple{E,Int}
|
|
muon_permutation::NTuple{U,Int}
|
|
tauon_permutation::NTuple{T,Int}
|
|
|
|
function FeynmanDiagram(
|
|
structure::Vector{Vector{Int}},
|
|
elec_perm::Vector{Int},
|
|
muon_perm::Vector{Int},
|
|
tauon_perm::Vector{Int},
|
|
::Val{E},
|
|
::Val{U},
|
|
::Val{T},
|
|
::Val{M}
|
|
) where {E,U,T,M}
|
|
@assert E == length(elec_perm)
|
|
@assert U == length(muon_perm)
|
|
@assert T == length(tauon_perm)
|
|
N = E + U + T
|
|
|
|
return new{N,E,U,T,M,FlatMatrix{Int64,2 * N - 2 + M,N}}(FlatMatrix(structure), NTuple{E,Int}(elec_perm), NTuple{U,Int}(muon_perm), NTuple{T,Int}(tauon_perm))
|
|
end
|
|
end
|
|
|
|
# representation of a virtual particle and the return type of the virtual_particles function
|
|
struct VirtualParticle{PROC<:QEDbase.AbstractProcessDefinition,PT<:QEDbase.AbstractParticleType,I,O}
|
|
proc::PROC
|
|
species::PT
|
|
in_particle_contributions::NTuple{I,Bool}
|
|
out_particle_contributions::NTuple{O,Bool}
|
|
end
|
|
|
|
import Base: +
|
|
|
|
# "addition" of the bool tuples
|
|
# realistically, there should never be "colliding" 1s. if there are there is probably an error and this should be asserted
|
|
function +(a::Tuple{NTuple{I,Bool},NTuple{O,Bool}}, b::Tuple{NTuple{I,Bool},NTuple{O,Bool}}) where {I,O}
|
|
return (ntuple(i -> a[1][i] || b[1][i], I), ntuple(i -> a[2][i] || b[2][i], O))
|
|
end
|
|
|
|
function virtual_particles(proc::QEDbase.AbstractProcessDefinition, diagram::FeynmanDiagram{N,E,U,T,M,FM}) where {N,E,U,T,M,FM}
|
|
I = number_incoming_particles(proc)
|
|
O = number_outgoing_particles(proc)
|
|
|
|
# map of all known particles' momentum composition
|
|
known_particles = Dict{Int64,Tuple{NTuple{I,Bool},NTuple{O,Bool}}}()
|
|
|
|
|
|
# 1: insert all the external ones (won't be returned), they all have exactly one 1 in their composition
|
|
# TODO
|
|
|
|
|
|
# 2: Loop:
|
|
# while there are incomplete fermion lines:
|
|
# take a fermion line where there is max=1 particle ∉ known_particles
|
|
# walk the fermion line, assign each virtual particle the momentum composition of the previous (or initial fermion if start) "+" the connected particle
|
|
# when/if the unknown particle is encountered, start walking from the other side
|
|
# when they meet at the unknown particle, assign the unknown particle Photon and left side - right side momentum contribution
|
|
# TODO
|
|
|
|
# 3: minimalize the contributions, i.e., if the number of contributing particles > half of all particles, invert both vectors
|
|
# if it's exactly half of all particles, think of some consistent way to break the symmetry, e.g. swap if the first particle is not contributing
|
|
# TODO
|
|
|
|
# 4: convert the known_particles Dict to an NTuple and remove the external particles (those with only 1 contributing momentum)
|
|
# TODO
|
|
|
|
return NTuple{?,VirtualParticle}()
|
|
end
|
|
|
|
function vertices(diagram::FeynmanDiagram{N,E,U,T,M,FM}) where {N,E,U,T,M,FM}
|
|
|
|
return NTuple{N,VertexType}()
|
|
end
|
|
|
|
|
|
#
|
|
# Generate Feynman Diagrams
|
|
#
|
|
|
|
mutable struct ExternalPhotonIterator
|
|
N::Int # number of fermion lines
|
|
M::Int # number of external photons
|
|
fermion_structures::Vector{Vector{Vector{Int}}}
|
|
fermion_structure_index::Int
|
|
photon_structures::Vector{Vector{Vector{Int}}}
|
|
photon_structure_index::Int
|
|
end
|
|
|
|
function _feynman_structures(n::Int, m::Int)
|
|
f = labelled_plane_trees(n)
|
|
return ExternalPhotonIterator(
|
|
n,
|
|
m,
|
|
f,
|
|
1,
|
|
_external_photon_multiplicity(f[1], n, m),
|
|
1
|
|
)
|
|
end
|
|
|
|
function Base.length(it::ExternalPhotonIterator)
|
|
return factorial(it.M + 3 * it.N - 3, 2 * it.N - 1)
|
|
end
|
|
|
|
function Base.iterate(iter::ExternalPhotonIterator)
|
|
return (iter.photon_structures[iter.photon_structure_index], nothing)
|
|
end
|
|
|
|
function Base.iterate(iter::ExternalPhotonIterator, n::Nothing)
|
|
if iter.photon_structure_index == length(iter.photon_structures)
|
|
if iter.fermion_structure_index == length(iter.fermion_structures)
|
|
return nothing
|
|
end
|
|
|
|
iter.fermion_structure_index += 1
|
|
iter.photon_structure_index = 1
|
|
iter.photon_structures = _external_photon_multiplicity(iter.fermion_structures[iter.fermion_structure_index], iter.N, iter.M)
|
|
else
|
|
iter.photon_structure_index += 1
|
|
end
|
|
|
|
return (iter.photon_structures[iter.photon_structure_index], nothing)
|
|
end
|
|
|
|
function _external_photon_multiplicity(fermion_structure::Vector{Vector{Int}}, n::Int, m::Int)
|
|
if m == 0
|
|
return [fermion_structure]
|
|
end
|
|
|
|
res = Vector{Vector{Vector{Int}}}()
|
|
|
|
# go through lines
|
|
for line_index in eachindex(fermion_structure)
|
|
|
|
# go through indices start to end
|
|
for index in 1:length(fermion_structure[line_index])+1
|
|
# copy to prevent muting
|
|
new_photon_structure = deepcopy(fermion_structure)
|
|
# add new photon
|
|
insert!(new_photon_structure[line_index], index, n + 1)
|
|
# recurse
|
|
append!(res, _external_photon_multiplicity(new_photon_structure, n + 1, m - 1))
|
|
end
|
|
end
|
|
|
|
return res
|
|
end
|
|
|
|
mutable struct FeynmanDiagramIterator{E,U,T,M}
|
|
e::Val{E} # number of electron lines (indices 1 - e)
|
|
e_perms::Vector{Vector{Int}} # list of all the possible permutations of the electrons
|
|
e_index::Int
|
|
u::Val{U} # number of muon lines (indices e+1 - e+u)
|
|
u_perms::Vector{Vector{Int}}
|
|
u_index::Int
|
|
t::Val{T} # number of tauon lines (indices e+u+1 - e+u+t)
|
|
t_perms::Vector{Vector{Int}}
|
|
t_index::Int
|
|
m::Val{M} # number of external photons
|
|
photon_structure_iter::ExternalPhotonIterator
|
|
photon_structure::Vector{Vector{Int}} # current structure that's being permuted
|
|
end
|
|
|
|
function Base.length(it::FeynmanDiagramIterator{E,U,T,M}) where {E,U,T,M}
|
|
N = E + U + T
|
|
return factorial(M + 3 * N - 3, 2 * N - 1) * factorial(E) * factorial(U) * factorial(T)
|
|
end
|
|
|
|
function Base.iterate(iter::FeynmanDiagramIterator{E,U,T,M}) where {E,U,T,M}
|
|
N = E + U + T
|
|
f = FeynmanDiagram(iter.photon_structure, iter.e_perms[iter.e_index], iter.u_perms[iter.u_index], iter.t_perms[iter.t_index], iter.e, iter.u, iter.t, iter.m)
|
|
return (
|
|
f,
|
|
nothing
|
|
)
|
|
end
|
|
|
|
function Base.iterate(iter::FeynmanDiagramIterator{E,U,T,M}, ::Nothing) where {E,U,T,M}
|
|
iter.t_index += 1
|
|
|
|
if iter.t_index > length(iter.t_perms)
|
|
iter.t_index = 1
|
|
iter.u_index += 1
|
|
end
|
|
|
|
if iter.u_index > length(iter.u_perms)
|
|
iter.u_index = 1
|
|
iter.e_index += 1
|
|
end
|
|
|
|
if iter.e_index > length(iter.e_perms)
|
|
iter.e_index = 1
|
|
photon_iter_result = iterate(iter.photon_structure_iter, nothing)
|
|
if isnothing(photon_iter_result)
|
|
return nothing
|
|
end
|
|
(iter.photon_structure, _) = photon_iter_result
|
|
end
|
|
|
|
N = E + U + T
|
|
f = FeynmanDiagram(iter.photon_structure, iter.e_perms[iter.e_index], iter.u_perms[iter.u_index], iter.t_perms[iter.t_index], iter.e, iter.u, iter.t, iter.m)
|
|
return (
|
|
f,
|
|
nothing
|
|
)
|
|
end
|
|
|
|
function feynman_diagrams(proc::PROC) where {PROC<:GenericQEDProcess}
|
|
return feynman_diagrams(incoming_particles(proc), outgoing_particles(proc))
|
|
end
|
|
|
|
function feynman_diagrams(in_particles::Tuple, out_particles::Tuple)
|
|
_assert_particle_type_tuple(in_particles)
|
|
_assert_particle_type_tuple(out_particles)
|
|
|
|
count(x -> x isa Electron, in_particles) + count(x -> x isa Positron, out_particles) ==
|
|
count(x -> x isa Positron, in_particles) + count(x -> x isa Electron, out_particles) ||
|
|
throw(InvalidInputError("the given particles do not make a physical process"))
|
|
|
|
# get the fermion line counts and external photon count
|
|
e = count(x -> x isa Electron, in_particles) + count(x -> x isa Positron, out_particles)
|
|
m = count(x -> x isa Photon, in_particles) + count(x -> x isa Photon, out_particles)
|
|
# TODO: do this the same way as for e when muons and tauons are a part of QED.jl
|
|
u = 0
|
|
t = 0
|
|
|
|
f_iter = _feynman_structures(e + u + t, m)
|
|
e_perms = collect(permutations(Int[1:e;]))
|
|
u_perms = collect(permutations(Int[e+1:e+u;]))
|
|
t_perms = collect(permutations(Int[e+u+1:e+u+t;]))
|
|
first_photon_structure, _ = iterate(f_iter)
|
|
|
|
return FeynmanDiagramIterator(Val(e), e_perms, 1, Val(u), u_perms, 1, Val(t), t_perms, 1, Val(m), f_iter, first_photon_structure)
|
|
end
|