231 lines
6.8 KiB
Julia
231 lines
6.8 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{PROC,N,E,U,T,M} <: AbstractTreeLevelFeynmanDiagram where {P<:AbstractProcessDefinition,N,E,U,T,M}
|
|
# E, U, T, and M can be inferred from the PROC, but not necessarily in 0 runtime
|
|
proc::PROC
|
|
|
|
diagram_structure::NTuple{N,Vector{Int}}
|
|
|
|
electron_permutation::NTuple{E,Int}
|
|
muon_permutation::NTuple{U,Int}
|
|
tauon_permutation::NTuple{T,Int}
|
|
|
|
function FeynmanDiagram(
|
|
proc::PROC,
|
|
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 {PROC<:AbstractProcessDefinition,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{PROC,N,E,U,T,M}(proc, NTuple{N,Vector{Int}}(structure), NTuple{E,Int}(elec_perm), NTuple{U,Int}(muon_perm), NTuple{T,Int}(tauon_perm))
|
|
end
|
|
end
|
|
|
|
process(diagram::FeynmanDiagram) = d.proc
|
|
|
|
function virtual_particles(diagram::FeynmanDiagram)
|
|
|
|
return NTuple{N,Tuple{QEDbase.AbstractParticleType,BitArray}}()
|
|
end
|
|
|
|
function vertices(::AbstractTreeLevelFeynmanDiagram)
|
|
|
|
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{PROC<:AbstractProcessDefinition,E,U,T,M}
|
|
process::PROC
|
|
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{<:AbstractProcessDefinition,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)
|
|
return (
|
|
FeynmanDiagram(iter.process, 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),
|
|
nothing
|
|
)
|
|
end
|
|
|
|
function Base.iterate(iter::FeynmanDiagramIterator, ::Nothing)
|
|
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
|
|
|
|
return (
|
|
FeynmanDiagram(iter.process, 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),
|
|
nothing
|
|
)
|
|
end
|
|
|
|
function feynman_diagrams(proc::PROC, e::Int, u::Int, t::Int, m::Int) where {PROC<:AbstractProcessDefinition}
|
|
if (e + u + t) * 2 + m < 4
|
|
throw(InvalidInputError("At least 4 particles must participate in a scattering process for a physical result!"))
|
|
end
|
|
|
|
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(proc, Val(e), e_perms, 1, Val(u), u_perms, 1, Val(t), t_perms, 1, Val(m), f_iter, first_photon_structure)
|
|
end
|