diff --git a/notebooks/diagram_generation.ipynb b/notebooks/diagram_generation.ipynb
index 25d3e49..c4d164c 100644
--- a/notebooks/diagram_generation.ipynb
+++ b/notebooks/diagram_generation.ipynb
@@ -9,54 +9,25 @@
      "name": "stderr",
      "output_type": "stream",
      "text": [
-      "WARNING: Method definition number_particles(QEDbase.AbstractProcessDefinition, DIR, PT) where {DIR<:QEDbase.ParticleDirection, PT<:QEDbase.AbstractParticleType} in module QEDbase at /home/antonr/.julia/packages/QEDbase/gWHxL/src/interfaces/process.jl:208 overwritten in module FeynmanDiagramGenerator at /home/antonr/repos/FeynmanDiagramGenerator/src/QEDprocesses_patch.jl:23.\n",
+      "WARNING: Method definition (::Type{QEDcore.ParticleStateful{DIR, SPECIES, ELEMENT} where ELEMENT<:QEDbase.AbstractFourMomentum})(QEDbase.AbstractFourMomentum) where {DIR<:QEDbase.ParticleDirection, SPECIES<:QEDbase.AbstractParticleType} in module QEDcore at /home/antonr/.julia/packages/QEDcore/uVldP/src/phase_spaces/create.jl:7 overwritten in module MetagraphOptimization at /home/antonr/.julia/packages/MetagraphOptimization/mvCVq/src/QEDprocesses_patch.jl:15.\n",
       "ERROR: Method overwriting is not permitted during Module precompilation. Use `__precompile__(false)` to opt-out of precompilation.\n"
      ]
     }
    ],
    "source": [
+    "using Revise\n",
     "using FeynmanDiagramGenerator"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "254-element Vector{VirtualParticle{GenericQEDProcess{Tuple{Photon, Photon, Photon, Photon, Photon, Photon, Photon, Electron}, Tuple{Photon, Electron}}, PT, 8, 2} where PT<:AbstractParticleType}:\n",
-       " positron: \t00000000 | 11\n",
-       " electron: \t00000001 | 10\n",
-       " positron: \t00000010 | 01\n",
-       " electron: \t00000011 | 00\n",
-       " positron: \t00000100 | 01\n",
-       " electron: \t00000101 | 00\n",
-       " positron: \t00001000 | 01\n",
-       " electron: \t00001001 | 00\n",
-       " positron: \t00010000 | 01\n",
-       " electron: \t00010001 | 00\n",
-       " ⋮\n",
-       " electron: \t11100001 | 10\n",
-       " positron: \t11100010 | 01\n",
-       " electron: \t11100011 | 00\n",
-       " positron: \t11100100 | 01\n",
-       " electron: \t11100101 | 00\n",
-       " positron: \t11101000 | 01\n",
-       " electron: \t11101001 | 00\n",
-       " positron: \t11110000 | 01\n",
-       " electron: \t11110001 | 00"
-      ]
-     },
-     "metadata": {},
-     "output_type": "display_data"
-    }
-   ],
+   "outputs": [],
    "source": [
-    "proc = GenericQEDProcess(7, 1, 1, 1, 0, 0)\n",
+    "proc = GenericQEDProcess(3, 1, 1, 1, 0, 0)\n",
     "all_particles = Set()\n",
-    "for fd in feynman_diagrams(proc) \n",
+    "for fd in feynman_diagrams(proc)\n",
     "    push!(all_particles, virtual_particles(proc, fd)...)\n",
     "end\n",
     "all_particles = sort([all_particles...])\n",
@@ -65,125 +36,27 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "Dict{VirtualParticle, Vector{Tuple{VirtualParticle, VirtualParticle}}} with 254 entries:\n",
-       "  positron: \t01000010 | 11 => [(photon: \t00000000 | 10, positron: \t01000010 | 01),…\n",
-       "  positron: \t11010000 | 11 => [(photon: \t00000000 | 10, positron: \t11010000 | 01),…\n",
-       "  electron: \t10100011 | 00 => [(photon: \t00000010 | 00, electron: \t10100001 | 00),…\n",
-       "  positron: \t11000000 | 11 => [(photon: \t00000000 | 10, positron: \t11000000 | 01),…\n",
-       "  positron: \t00010100 | 01 => [(photon: \t00000100 | 00, positron: \t00010000 | 01),…\n",
-       "  positron: \t01001000 | 11 => [(photon: \t00000000 | 10, positron: \t01001000 | 01),…\n",
-       "  positron: \t11000110 | 01 => [(photon: \t00000010 | 00, positron: \t11000100 | 01),…\n",
-       "  electron: \t10001011 | 10 => [(photon: \t00000000 | 10, electron: \t10001011 | 00),…\n",
-       "  positron: \t00100000 | 11 => [(photon: \t00000000 | 10, positron: \t00100000 | 01),…\n",
-       "  positron: \t01010000 | 11 => [(photon: \t00000000 | 10, positron: \t01010000 | 01),…\n",
-       "  electron: \t11101001 | 00 => [(photon: \t00001000 | 00, electron: \t11100001 | 00),…\n",
-       "  electron: \t00100101 | 10 => [(photon: \t00000000 | 10, electron: \t00100101 | 00),…\n",
-       "  electron: \t11000101 | 00 => [(photon: \t00000100 | 00, electron: \t11000001 | 00),…\n",
-       "  electron: \t01101001 | 00 => [(photon: \t00001000 | 00, electron: \t01100001 | 00),…\n",
-       "  positron: \t00101000 | 01 => [(photon: \t00001000 | 00, positron: \t00100000 | 01),…\n",
-       "  positron: \t00011000 | 11 => [(photon: \t00000000 | 10, positron: \t00011000 | 01),…\n",
-       "  positron: \t00011010 | 01 => [(photon: \t00000010 | 00, positron: \t00011000 | 01),…\n",
-       "  electron: \t10000101 | 10 => [(photon: \t00000000 | 10, electron: \t10000101 | 00),…\n",
-       "  electron: \t10000111 | 00 => [(photon: \t00000010 | 00, electron: \t10000101 | 00),…\n",
-       "  ⋮                       => ⋮"
-      ]
-     },
-     "metadata": {},
-     "output_type": "display_data"
-    }
-   ],
+   "outputs": [],
    "source": [
     "pairs = FeynmanDiagramGenerator.particle_pairs(all_particles)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "280-element Vector{Tuple{VirtualParticle, VirtualParticle, VirtualParticle}}:\n",
-       " (photon: \t00000000 | 10, electron: \t01001011 | 00, positron: \t10110100 | 01)\n",
-       " (photon: \t00000000 | 10, electron: \t00010111 | 00, positron: \t11101000 | 01)\n",
-       " (photon: \t00000000 | 10, electron: \t00110011 | 00, positron: \t11001100 | 01)\n",
-       " (photon: \t00000000 | 10, electron: \t10111001 | 00, positron: \t01000110 | 01)\n",
-       " (photon: \t00000000 | 10, electron: \t10110101 | 00, positron: \t01001010 | 01)\n",
-       " (photon: \t00000000 | 10, electron: \t11110001 | 00, positron: \t00001110 | 01)\n",
-       " (photon: \t00000000 | 10, electron: \t11010011 | 00, positron: \t00101100 | 01)\n",
-       " (photon: \t00000000 | 10, electron: \t01001101 | 00, positron: \t10110010 | 01)\n",
-       " (photon: \t00000000 | 10, electron: \t11101001 | 00, positron: \t00010110 | 01)\n",
-       " (photon: \t00000000 | 10, electron: \t10011011 | 00, positron: \t01100100 | 01)\n",
-       " ⋮\n",
-       " (photon: \t01000000 | 00, electron: \t00111001 | 00, positron: \t10000110 | 11)\n",
-       " (photon: \t01000000 | 00, electron: \t00000111 | 10, positron: \t10111000 | 01)\n",
-       " (photon: \t01000000 | 00, electron: \t00011011 | 00, positron: \t10100100 | 11)\n",
-       " (photon: \t01000000 | 00, electron: \t10101011 | 00, positron: \t00010100 | 11)\n",
-       " (photon: \t01000000 | 00, electron: \t10000111 | 10, positron: \t00111000 | 01)\n",
-       " (photon: \t01000000 | 00, electron: \t10100101 | 10, positron: \t00011010 | 01)\n",
-       " (photon: \t01000000 | 00, electron: \t00100111 | 00, positron: \t10011000 | 11)\n",
-       " (photon: \t01000000 | 00, electron: \t10001101 | 10, positron: \t00110010 | 01)\n",
-       " (photon: \t01000000 | 00, electron: \t10011101 | 00, positron: \t00100010 | 11)"
-      ]
-     },
-     "metadata": {},
-     "output_type": "display_data"
-    }
-   ],
+   "outputs": [],
    "source": [
     "triples = FeynmanDiagramGenerator.total_particle_triples(all_particles)"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 9,
+   "execution_count": null,
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "s: 40320, should be: 40320\n",
-      "number of triples: 280\n"
-     ]
-    },
-    {
-     "data": {
-      "text/plain": [
-       "280-element Vector{Tuple{VirtualParticle, VirtualParticle, VirtualParticle}}:\n",
-       " (photon: \t00000000 | 10, electron: \t00001111 | 00, positron: \t11110000 | 01)\n",
-       " (photon: \t00000000 | 10, electron: \t00010111 | 00, positron: \t11101000 | 01)\n",
-       " (photon: \t00000000 | 10, electron: \t00011011 | 00, positron: \t11100100 | 01)\n",
-       " (photon: \t00000000 | 10, electron: \t00011101 | 00, positron: \t11100010 | 01)\n",
-       " (photon: \t00000000 | 10, electron: \t00100111 | 00, positron: \t11011000 | 01)\n",
-       " (photon: \t00000000 | 10, electron: \t00101011 | 00, positron: \t11010100 | 01)\n",
-       " (photon: \t00000000 | 10, electron: \t00101101 | 00, positron: \t11010010 | 01)\n",
-       " (photon: \t00000000 | 10, electron: \t00110011 | 00, positron: \t11001100 | 01)\n",
-       " (photon: \t00000000 | 10, electron: \t00110101 | 00, positron: \t11001010 | 01)\n",
-       " (photon: \t00000000 | 10, electron: \t00111001 | 00, positron: \t11000110 | 01)\n",
-       " ⋮\n",
-       " (photon: \t01000000 | 00, electron: \t10100101 | 10, positron: \t00011010 | 01)\n",
-       " (photon: \t01000000 | 00, electron: \t10100111 | 00, positron: \t00011000 | 11)\n",
-       " (photon: \t01000000 | 00, electron: \t10101001 | 10, positron: \t00010110 | 01)\n",
-       " (photon: \t01000000 | 00, electron: \t10101011 | 00, positron: \t00010100 | 11)\n",
-       " (photon: \t01000000 | 00, electron: \t10101101 | 00, positron: \t00010010 | 11)\n",
-       " (photon: \t01000000 | 00, electron: \t10110001 | 10, positron: \t00001110 | 01)\n",
-       " (photon: \t01000000 | 00, electron: \t10110011 | 00, positron: \t00001100 | 11)\n",
-       " (photon: \t01000000 | 00, electron: \t10110101 | 00, positron: \t00001010 | 11)\n",
-       " (photon: \t01000000 | 00, electron: \t10111001 | 00, positron: \t00000110 | 11)"
-      ]
-     },
-     "metadata": {},
-     "output_type": "display_data"
-    }
-   ],
+   "outputs": [],
    "source": [
     "function n(vp::VirtualParticle)\n",
     "    if !haskey(pairs, vp)\n",
@@ -206,6 +79,15 @@
     "\n",
     "sort(triples)"
    ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "generate_DAG(proc)"
+   ]
   }
  ],
  "metadata": {
diff --git a/src/FeynmanDiagramGenerator.jl b/src/FeynmanDiagramGenerator.jl
index 320c8b6..75ce75b 100644
--- a/src/FeynmanDiagramGenerator.jl
+++ b/src/FeynmanDiagramGenerator.jl
@@ -17,7 +17,7 @@ export FlatMatrix
 export GenericQEDProcess, isphysical
 
 export AbstractTreeLevelFeynmanDiagram, FeynmanVertex, FeynmanDiagram
-export external_particles, virtual_particles, process
+export external_particles, virtual_particles, process, generate_DAG
 
 export VirtualParticle, number_contributions
 export contributions, in_contributions, out_contributions
diff --git a/src/diagrams/diagrams.jl b/src/diagrams/diagrams.jl
index 77134e6..b966bb4 100644
--- a/src/diagrams/diagrams.jl
+++ b/src/diagrams/diagrams.jl
@@ -39,12 +39,12 @@ 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}
+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}
@@ -59,19 +59,26 @@ struct FeynmanDiagram{N,E,U,T,M,FM} <: AbstractTreeLevelFeynmanDiagram where {N,
         ::Val{E},
         ::Val{U},
         ::Val{T},
-        ::Val{M}
+        ::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))
+        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}
+struct VirtualParticle{
+    PROC<:QEDbase.AbstractProcessDefinition,PT<:QEDbase.AbstractParticleType,I,O
+}
     proc::PROC
     species::PT
     in_particle_contributions::NTuple{I,Bool}
@@ -80,7 +87,10 @@ end
 
 function Base.show(io::IO, vp::VirtualParticle)
     pr = x -> x ? "1" : "0"
-    print(io, "$(vp.species): \t$(*(pr.(vp.in_particle_contributions)...)) | $(*(pr.(vp.out_particle_contributions)...))")
+    return print(
+        io,
+        "$(vp.species): \t$(*(pr.(vp.in_particle_contributions)...)) | $(*(pr.(vp.out_particle_contributions)...))",
+    )
 end
 
 function QEDbase.process(vp::VirtualParticle)
@@ -96,11 +106,14 @@ out_contributions(vp::VirtualParticle) = vp.out_particle_contributions
 contributions(vp::VirtualParticle) = ((in_contributions(vp), out_contributions(vp)))
 
 is_virtual(vp::VirtualParticle) = number_contributions(vp) > 1
+is_external(vp::VirtualParticle) = number_contributions(vp) == 1
 
 import Base: +
 
 # "addition" of the bool tuples
-function +(a::Tuple{NTuple{I,Bool},NTuple{O,Bool}}, b::Tuple{NTuple{I,Bool},NTuple{O,Bool}}) where {I,O}
+function +(
+    a::Tuple{NTuple{I,Bool},NTuple{O,Bool}}, b::Tuple{NTuple{I,Bool},NTuple{O,Bool}}
+) where {I,O}
     # realistically, there should never be "colliding" 1s. if there are there is probably an error and this should be asserted
     #= for (i, j) in zip(a[1], b[1]) @assert !(i && j) end
     for (i, j) in zip(a[2], b[2]) @assert !(i && j) end =#
@@ -111,7 +124,9 @@ end
 invert(::Electron) = Positron()
 invert(::Positron) = Electron()
 invert(::Photon) = Photon()
-invert(::AbstractParticleType) = throw(InvalidInputError("unimplemented for this particle type"))
+function invert(::AbstractParticleType)
+    throw(InvalidInputError("unimplemented for this particle type"))
+end
 
 function invert(virtual_particle::VirtualParticle)
     I = length(virtual_particle.in_particle_contributions)
@@ -120,7 +135,8 @@ function invert(virtual_particle::VirtualParticle)
         virtual_particle.proc,
         invert(virtual_particle.species),
         ntuple(x -> !virtual_particle.in_particle_contributions[x], I),
-        ntuple(x -> !virtual_particle.out_particle_contributions[x], O))
+        ntuple(x -> !virtual_particle.out_particle_contributions[x], O),
+    )
 end
 
 # normalize the representation
@@ -138,7 +154,12 @@ function normalize(virtual_particle::VirtualParticle)
     end
 end
 
-function _momentum_contribution(proc::AbstractProcessDefinition, dir::ParticleDirection, species::AbstractParticleType, index::Int)
+function _momentum_contribution(
+    proc::AbstractProcessDefinition,
+    dir::ParticleDirection,
+    species::AbstractParticleType,
+    index::Int,
+)
     I = number_incoming_particles(proc)
     O = number_outgoing_particles(proc)
 
@@ -151,7 +172,10 @@ function _momentum_contribution(proc::AbstractProcessDefinition, dir::ParticleDi
             particles_seen += 1
         end
         if particles_seen == index
-            return (((is_incoming(dir) && x == c for x in 1:I)...,), ((is_outgoing(dir) && x == c for x in 1:O)...,))
+            return (
+                ((is_incoming(dir) && x == c for x in 1:I)...,),
+                ((is_outgoing(dir) && x == c for x in 1:O)...,),
+            )
         end
     end
 
@@ -159,17 +183,25 @@ function _momentum_contribution(proc::AbstractProcessDefinition, dir::ParticleDi
 end
 
 function _fermion_type(proc::AbstractProcessDefinition, n::Int)
-    E = number_particles(proc, Incoming(), Electron()) + number_particles(proc, Outgoing(), Positron())
+    E =
+        number_particles(proc, Incoming(), Electron()) +
+        number_particles(proc, Outgoing(), Positron())
     U = 0 # TODO add muons
     T = 0 # TODO add tauons
-    M = number_particles(proc, Incoming(), Photon()) + number_particles(proc, Outgoing(), Photon())
+    M =
+        number_particles(proc, Incoming(), Photon()) +
+        number_particles(proc, Outgoing(), Photon())
     N = E + U + T
 
     # from the fermion index, get (Direction, Species, n) tuple, where n means it's the nth particle of that dir and species
     if (n > 0 && n <= E)
         electron_n = n
         if electron_n > number_particles(proc, Incoming(), Electron())
-            return (Outgoing(), Positron(), electron_n - number_particles(proc, Incoming(), Electron()))
+            return (
+                Outgoing(),
+                Positron(),
+                electron_n - number_particles(proc, Incoming(), Electron()),
+            )
         else
             return (Incoming(), Electron(), electron_n)
         end
@@ -185,7 +217,11 @@ function _fermion_type(proc::AbstractProcessDefinition, n::Int)
         # photon
         photon_n = n - N
         if photon_n > number_particles(proc, Incoming(), Photon())
-            return (Outgoing(), Photon(), photon_n - number_particles(proc, Incoming(), Photon()))
+            return (
+                Outgoing(),
+                Photon(),
+                photon_n - number_particles(proc, Incoming(), Photon()),
+            )
         else
             return (Incoming(), Photon(), photon_n)
         end
@@ -194,7 +230,11 @@ function _fermion_type(proc::AbstractProcessDefinition, n::Int)
         electron_n = n - N - M
         if electron_n > number_particles(proc, Outgoing(), Electron())
             # incoming positron
-            return (Incoming(), Positron(), electron_n - number_particles(proc, Outgoing(), Electron()))
+            return (
+                Incoming(),
+                Positron(),
+                electron_n - number_particles(proc, Outgoing(), Electron()),
+            )
         else
             # outgoing electron
             return (Outgoing(), Electron(), electron_n)
@@ -245,7 +285,8 @@ end
 Return true if the momenta contributions of `a` and `b` are disjunct.
 """
 function disjunct(a::VirtualParticle, b::VirtualParticle)
-    for (a_contrib, b_contrib) in Iterators.zip(Iterators.flatten.(contributions.((a, b)))...)
+    for (a_contrib, b_contrib) in
+        Iterators.zip(Iterators.flatten.(contributions.((a, b)))...)
         if b_contrib && a_contrib
             return false
         end
@@ -260,7 +301,8 @@ end
 Returns true if the set of particles contributing to `a` are contains the set of particles contributing to `b`.
 """
 function contains(a::VirtualParticle, b::VirtualParticle)
-    for (a_contrib, b_contrib) in Iterators.zip(Iterators.flatten.(contributions.((a, b)))...)
+    for (a_contrib, b_contrib) in
+        Iterators.zip(Iterators.flatten.(contributions.((a, b)))...)
         if b_contrib && !a_contrib
             return false
         end
@@ -279,7 +321,8 @@ function make_up(a::VirtualParticle, b::VirtualParticle, c::VirtualParticle)
         return false
     end
     # it should be unnecessary to check here that a and b can actually react. if a + b = c they must be able to if a, b and c all exist in the diagram.
-    for (a_contrib, b_contrib, c_contrib) in Iterators.zip(Iterators.flatten.(contributions.((a, b, c)))...)
+    for (a_contrib, b_contrib, c_contrib) in
+        Iterators.zip(Iterators.flatten.(contributions.((a, b, c)))...)
         if c_contrib != a_contrib + b_contrib
             return false
         end
@@ -294,7 +337,8 @@ end
 Return true if a, b and c combined contain all external particles exactly once.
 """
 function are_total(a::VirtualParticle, b::VirtualParticle, c::VirtualParticle)
-    for (a_contrib, b_contrib, c_contrib) in Iterators.zip(Iterators.flatten.(contributions.((a, b, c)))...)
+    for (a_contrib, b_contrib, c_contrib) in
+        Iterators.zip(Iterators.flatten.(contributions.((a, b, c)))...)
         if a_contrib + b_contrib + c_contrib != 1
             return false
         end
@@ -311,14 +355,17 @@ function particle_pairs(particles::Vector)
     all_particles = sort([_pseudo_virtual_particles(proc)..., particles...])
 
     # find pairs for every particle after the external ones (those can't have pairs)
-    for p_i in number_incoming_particles(proc)+number_outgoing_particles(proc)+1:length(all_particles)
+    for p_i in
+        (number_incoming_particles(proc) + number_outgoing_particles(proc) + 1):length(
+        all_particles
+    )
         p = all_particles[p_i]
         pairs[p] = Vector{Tuple{VirtualParticle,VirtualParticle}}()
 
         # only need to consider external particles and virtual particles that come before p_i
-        for p_a_i in 1:p_i-2
+        for p_a_i in 1:(p_i - 2)
             # and only partners between a and p_i
-            for p_b_i in p_a_i+1:p_i-1
+            for p_b_i in (p_a_i + 1):(p_i - 1)
                 p_a = all_particles[p_a_i]
                 p_b = all_particles[p_b_i]
 
@@ -343,11 +390,20 @@ function total_particle_triples(particles::Vector)
     photons = filter(p -> is_boson(particle_species(p)), working_set)
 
     # make electrons a set for fast deletion
-    electrons = Set(filter(p -> is_fermion(particle_species(p)) && is_particle(particle_species(p)), working_set))
+    electrons = Set(
+        filter(
+            p -> is_fermion(particle_species(p)) && is_particle(particle_species(p)),
+            working_set,
+        ),
+    )
 
     # make positrons a set for fast lookup
-    positrons = Set(filter(p -> is_fermion(particle_species(p)) && is_anti_particle(particle_species(p)), working_set))
-
+    positrons = Set(
+        filter(
+            p -> is_fermion(particle_species(p)) && is_anti_particle(particle_species(p)),
+            working_set,
+        ),
+    )
 
     # no participant can have more than half the external particles, so every possible particle is contained here
     # every photon has exactly one electron and positron partner
@@ -358,7 +414,11 @@ function total_particle_triples(particles::Vector)
             end
 
             # create the only partner the ph and e could have together, then look for it in the actual positrons
-            expected_p = invert(VirtualParticle(proc, particle_species(e), (contributions(ph) + contributions(e))...))
+            expected_p = invert(
+                VirtualParticle(
+                    proc, particle_species(e), (contributions(ph) + contributions(e))...
+                ),
+            )
 
             if expected_p in positrons
                 @assert are_total(ph, e, expected_p)
@@ -376,10 +436,16 @@ end
 Return a vector of `VirtualParticle` for each external particle. These are not actually virtual particles, but can be helpful as entry points.
 """
 function _pseudo_virtual_particles(proc::AbstractProcessDefinition)
-    return sort(_external_particle.(proc, [1:number_incoming_particles(proc)+number_outgoing_particles(proc);]))
+    return sort(
+        _external_particle.(
+            proc, [1:(number_incoming_particles(proc) + number_outgoing_particles(proc));]
+        ),
+    )
 end
 
-function virtual_particles(proc::AbstractProcessDefinition, diagram::FeynmanDiagram{N,E,U,T,M,FM}) where {N,E,U,T,M,FM}
+function virtual_particles(
+    proc::AbstractProcessDefinition, diagram::FeynmanDiagram{N,E,U,T,M,FM}
+) where {N,E,U,T,M,FM}
     fermion_lines = PriorityQueue{Int64,Int64}()
 
     # count number of internal photons in each fermion line and make a priority queue for fermion line => number of internal photons
@@ -452,11 +518,17 @@ function virtual_particles(proc::AbstractProcessDefinition, diagram::FeynmanDiag
                 else # if the binding particle is unknown
                     # we have arrived at the "middle" of the line
                     # this line will be the unknown particle for the other lines
-                    internal_photon_contributions[current_line] = cumulative_mom + unknown_photon_momentum
+                    internal_photon_contributions[current_line] =
+                        cumulative_mom + unknown_photon_momentum
                     # now we know that the fermion line that binding_particle binds to on the other end has one fewer unknown internal photons
                     fermion_lines[binding_particle] -= 1
                     # add the internal photon virtual particle
-                    push!(result, VirtualParticle(proc, Photon(), (cumulative_mom + unknown_photon_momentum)...))
+                    push!(
+                        result,
+                        VirtualParticle(
+                            proc, Photon(), (cumulative_mom + unknown_photon_momentum)...
+                        ),
+                    )
                     break
                 end
             else # binding_particle is an external photon
@@ -466,7 +538,7 @@ function virtual_particles(proc::AbstractProcessDefinition, diagram::FeynmanDiag
         end
     end
 
-    return normalize.(result)[1:end-1]
+    return normalize.(result)[1:(end - 1)]
 end
 
 #
@@ -484,14 +556,7 @@ 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
-    )
+    return ExternalPhotonIterator(n, m, f, 1, _external_photon_multiplicity(f[1], n, m), 1)
 end
 
 function Base.length(it::ExternalPhotonIterator)
@@ -510,7 +575,9 @@ function Base.iterate(iter::ExternalPhotonIterator, n::Nothing)
 
         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)
+        iter.photon_structures = _external_photon_multiplicity(
+            iter.fermion_structures[iter.fermion_structure_index], iter.N, iter.M
+        )
     else
         iter.photon_structure_index += 1
     end
@@ -518,7 +585,9 @@ function Base.iterate(iter::ExternalPhotonIterator, n::Nothing)
     return (iter.photon_structures[iter.photon_structure_index], nothing)
 end
 
-function _external_photon_multiplicity(fermion_structure::Vector{Vector{Int}}, n::Int, m::Int)
+function _external_photon_multiplicity(
+    fermion_structure::Vector{Vector{Int}}, n::Int, m::Int
+)
     if m == 0
         return [fermion_structure]
     end
@@ -529,7 +598,7 @@ function _external_photon_multiplicity(fermion_structure::Vector{Vector{Int}}, n
     for line_index in eachindex(fermion_structure)
 
         # go through indices start to end
-        for index in 1:length(fermion_structure[line_index])+1
+        for index in 1:(length(fermion_structure[line_index]) + 1)
             # copy to prevent muting
             new_photon_structure = deepcopy(fermion_structure)
             # add new photon
@@ -564,11 +633,17 @@ 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
+    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}
@@ -594,11 +669,17 @@ function Base.iterate(iter::FeynmanDiagramIterator{E,U,T,M}, ::Nothing) where {E
     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
+    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}
@@ -625,10 +706,39 @@ function feynman_diagrams(in_particles::Tuple, out_particles::Tuple)
     # left electrons -> left muons -> left tauons -> left photons -> right photons -> right electrons -> right muons -> right tauons
     # a "left" fermion is simply an incoming fermion or outgoing antifermion of the type, while a "left" photon is an incoming photon, and the reverse for the right ones
     f_iter = _feynman_structures(e + u + t, m)
-    e_perms = collect(permutations(Int[n+m+1:n+m+e;]))
-    u_perms = collect(permutations(Int[n+m+e+1:n+m+e+u;]))
-    t_perms = collect(permutations(Int[n+m+e+u+1:n+m+e+u+t;]))
+    e_perms = collect(permutations(Int[(n + m + 1):(n + m + e);]))
+    u_perms = collect(permutations(Int[(n + m + e + 1):(n + m + e + u);]))
+    t_perms = collect(permutations(Int[(n + m + e + u + 1):(n + m + 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)
+    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
+
+function virtual_particles(proc::GenericQEDProcess)
+    if !isempty(proc.virtual_particles_cache)
+        return proc.virtual_particles_cache
+    end
+
+    # use a set for deduplication
+    particles = Set{VirtualParticle}()
+    for fd in feynman_diagrams(proc)
+        push!(particles, virtual_particles(proc, fd)...)
+    end
+
+    # convert to vector
+    proc.virtual_particles_cache = [particles...]
+    return return proc.virtual_particles_cache
 end
diff --git a/src/generic_process_def.jl b/src/generic_process_def.jl
index 7690888..dd0db08 100644
--- a/src/generic_process_def.jl
+++ b/src/generic_process_def.jl
@@ -1,41 +1,115 @@
 _assert_particle_type_tuple(::Tuple{}) = nothing
-_assert_particle_type_tuple(t::Tuple{AbstractParticleType,Vararg}) = _assert_particle_type_tuple(t[2:end])
-_assert_particle_type_tuple(t::Any) = throw(InvalidInputError("invalid input, provide a tuple of AbstractParticleTypes to construct a GenericQEDProcess"))
+function _assert_particle_type_tuple(t::Tuple{AbstractParticleType,Vararg})
+    return _assert_particle_type_tuple(t[2:end])
+end
+function _assert_particle_type_tuple(t::Any)
+    throw(
+        InvalidInputError(
+            "invalid input, provide a tuple of AbstractParticleTypes to construct a GenericQEDProcess",
+        ),
+    )
+end
 
-struct GenericQEDProcess{INT,OUTT} <: QEDbase.AbstractProcessDefinition where {INT<:Tuple,OUTT<:Tuple}
+mutable struct GenericQEDProcess{INT,OUTT,INSP,OUTSP} <: AbstractProcessDefinition where {
+    INT<:Tuple,OUTT<:Tuple,INSP<:Tuple,OUTSP<:Tuple
+}
     incoming_particles::INT
     outgoing_particles::OUTT
 
-    function GenericQEDProcess(in_particles::INT, out_particles::OUTT) where {INT<:Tuple,OUTT<:Tuple}
+    incoming_spins_pols::INSP
+    outgoing_spins_pols::OUTSP
+
+    virtual_particles_cache::Vector
+
+    function GenericQEDProcess(
+        in_particles::INT, out_particles::OUTT, in_sp::INSP, out_sp::OUTSP
+    ) where {INT<:Tuple,OUTT<:Tuple,INSP<:Tuple,OUTSP<:Tuple}
         _assert_particle_type_tuple(in_particles)
         _assert_particle_type_tuple(out_particles)
 
-        return new{INT,OUTT}(in_particles, out_particles)
+        return new{INT,OUTT,INSP,OUTSP}(in_particles, out_particles, in_sp, out_sp, [])
     end
 
     """
         GenericQEDProcess(in_ph::Int, out_ph::Int, in_el::Int, out_el::Int, in_po::Int, out_po::Int)
 
     Convenience constructor from numbers of input/output photons, electrons and positrons.
+    Uses `AllSpin()` and `AllPol()` for every particle's spin/pol by default.
     """
-    function GenericQEDProcess(in_ph::Int, out_ph::Int, in_el::Int, out_el::Int, in_po::Int, out_po::Int)
-        in_p = ntuple(i -> i <= in_ph ? Photon() : i <= in_ph + in_el ? Electron() : Positron(), in_ph + in_el + in_po)
-        out_p = ntuple(i -> i <= out_ph ? Photon() : i <= out_ph + out_el ? Electron() : Positron(), out_ph + out_el + out_po)
-        return GenericQEDProcess(in_p, out_p)
+    function GenericQEDProcess(
+        in_ph::Int, out_ph::Int, in_el::Int, out_el::Int, in_po::Int, out_po::Int
+    )
+        in_p = ntuple(i -> if i <= in_ph
+            Photon()
+        elseif i <= in_ph + in_el
+            Electron()
+        else
+            Positron()
+        end, in_ph + in_el + in_po)
+        out_p = ntuple(i -> if i <= out_ph
+            Photon()
+        elseif i <= out_ph + out_el
+            Electron()
+        else
+            Positron()
+        end, out_ph + out_el + out_po)
+        in_sp = tuple([i <= in_ph ? AllPol() : AllSpin() for i in 1:length(in_p)]...)
+        out_sp = tuple([i <= out_ph ? AllPol() : AllSpin() for i in 1:length(out_p)]...)
+        return GenericQEDProcess(in_p, out_p, in_sp, out_sp)
     end
 end
 
-QEDprocesses.incoming_particles(proc::GenericQEDProcess{INT,OUTT}) where {INT,OUTT} = proc.incoming_particles
-QEDprocesses.outgoing_particles(proc::GenericQEDProcess{INT,OUTT}) where {INT,OUTT} = proc.outgoing_particles
+function spin_or_pol(
+    process::GenericQEDProcess,
+    dir::ParticleDirection,
+    species::AbstractParticleType,
+    n::Int,
+)
+    i = 0
+    c = n
+    for p in particles(process, dir)
+        i += 1
+        if p == species
+            c -= 1
+        end
+        if c == 0
+            break
+        end
+    end
+
+    if c != 0 || n <= 0
+        throw(
+            InvalidInputError(
+                "could not get $n-th spin/pol of $dir $species, does not exist"
+            ),
+        )
+    end
+
+    if is_incoming(dir)
+        return process.incoming_spins_pols[i]
+    elseif is_outgoing(dir)
+        return process.outgoing_spins_pols[i]
+    else
+        throw(InvalidInputError("unknown direction $(dir) given"))
+    end
+end
+
+function QEDprocesses.incoming_particles(proc::GenericQEDProcess{INT,OUTT}) where {INT,OUTT}
+    return proc.incoming_particles
+end
+function QEDprocesses.outgoing_particles(proc::GenericQEDProcess{INT,OUTT}) where {INT,OUTT}
+    return proc.outgoing_particles
+end
 
 function isphysical(proc::GenericQEDProcess)
-    return (number_particles(proc, Incoming(), Electron()) + number_particles(proc, Outgoing(), Positron()) ==
-            number_particles(proc, Incoming(), Positron()) + number_particles(proc, Outgoing(), Electron())) &&
-           number_particles(proc, Incoming()) + number_particles(proc, Outgoing()) >= 2
+    return (
+        number_particles(proc, Incoming(), Electron()) +
+        number_particles(proc, Outgoing(), Positron()) ==
+        number_particles(proc, Incoming(), Positron()) +
+        number_particles(proc, Outgoing(), Electron())
+    ) && number_particles(proc, Incoming()) + number_particles(proc, Outgoing()) >= 2
 end
 
 function matrix_element(proc::GenericQEDProcess, psp::PhaseSpacePoint)
-
-
     return nothing
 end
diff --git a/src/metagraph_impl/compute.jl b/src/metagraph_impl/compute.jl
index 5cb0590..8700670 100644
--- a/src/metagraph_impl/compute.jl
+++ b/src/metagraph_impl/compute.jl
@@ -7,7 +7,15 @@ struct ComputeTask_Triple <: AbstractComputeTask end            # from a triple
 struct ComputeTask_CollectTriples <: AbstractComputeTask end    # sum over triples results and 
 
 # import compute so we don't have to repeat it all the time
-import MetagraphOptimization.compute
+import MetagraphOptimization: compute, compute_effort
+
+compute_effort(::ComputeTask_BaseState) = 0
+compute_effort(::ComputeTask_Propagator) = 0
+compute_effort(::ComputeTask_Pair) = 0
+compute_effort(::ComputeTask_CollectPairs) = 0
+compute_effort(::ComputeTask_PropagatePairs) = 0
+compute_effort(::ComputeTask_Triple) = 0
+compute_effort(::ComputeTask_CollectTriples) = 0
 
 struct BaseStateInput{PS_T<:AbstractParticleStateful,SPIN_POL_T<:AbstractSpinOrPolarization}
     particle::PS_T
diff --git a/src/metagraph_impl/generation.jl b/src/metagraph_impl/generation.jl
index 0c4027a..77be83e 100644
--- a/src/metagraph_impl/generation.jl
+++ b/src/metagraph_impl/generation.jl
@@ -1,19 +1,264 @@
+using MetagraphOptimization: insert_node!, insert_edge!, make_node, make_edge
 
-function generate_DAG(proc::GenericQEDProcess)
-    # use a set for deduplication
-    particles = Set()
-    for fd in feynman_diagrams(proc)
-        push!(all_particles, virtual_particles(proc, fd)...)
+_construction_string(::Incoming) = "Incoming()"
+_construction_string(::Outgoing) = "Outgoing()"
+
+_construction_string(::Electron) = "Electron()"
+_construction_string(::Positron) = "Positron()"
+_construction_string(::Photon) = "Photon()"
+
+_construction_string(::PolX) = "PolX()"
+_construction_string(::PolY) = "PolY()"
+_construction_string(::SpinUp) = "SpinUp()"
+_construction_string(::SpinDown) = "SpinDown()"
+
+function _parse_particle(name::String)
+    local dir
+    if startswith(name, "inc_")
+        dir = Incoming()
+    elseif startswith(name, "out_")
+        dir = Outgoing()
+    else
+        throw(InvalidInputError("failed to parse particle direction from \"$name\""))
     end
 
-    # convert to vector
-    particles = [particles...]
+    name = name[4:end]
 
-    external_particles = _pseudo_virtual_particles(proc)
-    pairs = particle_pairs(particles)
-    triples = total_particle_triples(particles)
+    local species
+    if startswith(name, "el")
+        species = Electron()
+    elseif startswith(name, "ph")
+        species = Photon()
+    elseif startswith(name, "po")
+        species = Positron()
+    else
+        throw(InvalidInputError("failed to parse particle species from name \"$name\""))
+    end
+
+    name = name[3:end]
+
+    local spin_pol
+    if startswith(name, "su")
+        spin_pol = SpinUp()
+    elseif startswith(name, "sd")
+        spin_pol = SpinDown()
+    elseif startswith(name, "px")
+        spin_pol = PolX()
+    elseif startswith(name, "py")
+        spin_pol = PolY()
+    else
+        throw(
+            InvalidInputError(
+                "failed to parse particle spin or polarization from \"$name\""
+            ),
+        )
+    end
+
+    name = name[3:end]
+
+    index = parse(Int, name)
+    return (dir, species, spin_pol, index)
+end
+
+function input_expr(instance::GenericQEDProcess, name::String, psp_symbol::Symbol)
+    if startswith(name, "bs_")
+        (dir, species, spin_pol, index) = _parse_particle(name[4:end])
+        dir_str = _construction_string(dir)
+        species_str = _construction_string(species)
+        sp_str = _construction_string(spin_pol)
+
+        return Meta.parse(
+            "BaseStateInput(
+                ParticleStateful($dir_str, $species_str, momentum($psp_symbol, $dir_str, $species_str, $index)),
+                $sp_str,
+            )",
+        )
+
+    elseif startswith(name, "pr_")
+        index = parse(Int, name[4:end]) # get index of the virtual particle in the process
+
+        vp = virtual_particles(proc)[index]
+        return Meta.parse("PropagatorInput(
+                              VirtualParticle(
+                                process($psp_symbol),
+                                $species_str,
+                                $(vp.in_particle_contributions),
+                                $(vp.out_particle_contributions)
+                              ),
+                              $psp_symbol
+                          )")
+    else
+        throw(InvalidInputError("failed to parse node name \"$name\""))
+    end
+end
+
+_species_str(::Photon) = "ph"
+_species_str(::Electron) = "el"
+_species_str(::Positron) = "po"
+
+_spin_pol_str(::SpinUp) = "su"
+_spin_pol_str(::SpinDown) = "sd"
+_spin_pol_str(::PolX) = "px"
+_spin_pol_str(::PolY) = "py"
+
+_dir_str(::Incoming) = "inc"
+_dir_str(::Outgoing) = "out"
+
+_spin_pols(::AllSpin) = (SpinUp(), SpinDown())
+_spin_pols(::SpinUp) = (SpinUp(),)
+_spin_pols(::SpinDown) = (SpinDown(),)
+_spin_pols(::AllPol) = (PolX(), PolY())
+_spin_pols(::PolX) = (PolX(),)
+_spin_pols(::PolY) = (PolY(),)
+
+_is_external(p::VirtualParticle) = number_contributions(p) == 1
+
+function _base_state_name(p::VirtualParticle)
+    proc = process(p)
+
+    dir = sum(in_contributions(p)) == 1 ? Incoming() : Outgoing()
+
+    # find particle in the contributions
+    index = 0
+    contribs = is_incoming(dir) ? in_contributions(p) : out_contributions(p)
+    for p in particles(proc, dir)
+        index += 1
+        if contribs[index]
+            break
+        end
+    end
+
+    species = particles(proc, dir)[index]
+
+    # find particle index of *this species*
+    species_index = 0
+    for i in 1:index
+        if particles(proc, dir)[i] == species
+            species_index += 1
+        end
+    end
+
+    spin_pol = spin_or_pol(proc, dir, species, species_index)
+
+    return string.(
+        "bs_$(_dir_str(dir))_$(_species_str(species))_",
+        _spin_pol_str.(_spin_pols(spin_pol)),
+        "_$(species_index)",
+    )
+end
+
+function generate_DAG(proc::GenericQEDProcess)
+    external_particles = _pseudo_virtual_particles(proc) # external particles that will be input to base_state tasks
+    particles = virtual_particles(proc)                  # virtual particles that will be input to propagator tasks
+    pairs = particle_pairs(particles)                    # pairs to generate the pair tasks
+    triples = total_particle_triples(particles)          # triples to generate the triple tasks
 
     graph = DAG()
 
+    # -- Base State Tasks --
+    base_state_task_outputs = Dict()
+    for dir in (Incoming(), Outgoing())
+        for species in (Electron(), Positron(), Photon())
+            for index in 1:number_particles(proc, dir, species)
+                for spin_pol in _spin_pols(spin_or_pol(proc, dir, species, index))
+                    # gen entry nodes
+                    # names are "bs_<dir>_<species>_<spin/pol>_<index>"
+                    data_node_name = "bs_$(_dir_str(dir))_$(_species_str(species))_$(_spin_pol_str(spin_pol))_$(index)"
 
+                    data_in = insert_node!(
+                        graph,
+                        make_node(DataTask(0), data_node_name);
+                        track=false,
+                        invalidate_cache=false,
+                    )
+
+                    # generate initial base_state tasks
+                    compute_base_state = insert_node!(
+                        graph,
+                        make_node(ComputeTask_BaseState());
+                        track=false,
+                        invalidate_cache=false,
+                    )
+
+                    data_out = insert_node!(
+                        graph, make_node(DataTask(0)); track=false, invalidate_cache=false
+                    )
+
+                    insert_edge!(graph, data_in, compute_base_state)
+                    insert_edge!(graph, compute_base_state, data_out)
+
+                    base_state_task_outputs[data_node_name] = data_out
+                end
+            end
+        end
+    end
+
+    # -- Propagator Tasks --
+    propagator_task_outputs = Dict()
+    vp_index = 0
+    for vp in virtual_particles(proc)
+        vp_index += 1
+
+        data_node_name = "pr_$vp_index"
+
+        data_in = insert_node!(
+            graph, make_node(DataTask(0)); track=false, invalidate_cache=false
+        )
+        compute_vp_propagator = insert_node!(
+            graph, make_node(ComputeTask_Propagator()); track=false, invalidate_cache=false
+        )
+        data_out = insert_node!(
+            graph, make_node(DataTask(0)); track=false, invalidate_cache=false
+        )
+
+        insert_edge!(graph, data_in, compute_vp_propagator)
+        insert_edge!(graph, compute_vp_propagator, data_out)
+
+        propagator_task_outputs[data_node_name] = data_out
+    end
+
+    # -- Pair Tasks --
+    pair_task_outputs = Dict()
+    for (product_particle, input_particle_vector) in pairs
+        # for all spins/pols of particles in product_particles do ...
+
+        pair_task_outputs[product_particle] = Vector()
+
+        for input_particles in input_particle_vector
+            particles_data_out_nodes = (Vector(), Vector())
+            c = 0
+            for p in input_particles
+                c += 1
+                if (is_external(p))
+                    # grab from base_states (broadcast over _base_state_name because it is a tuple for different spin_pols)
+                    push!.(
+                        Ref(particles_data_out_nodes[c]),
+                        getindex.(Ref(base_state_task_outputs), _base_state_name(p)),
+                    )
+                else
+                    # grab from propagated particles
+                    push!(particles_date_out_nodes[c], pair_task_outputs[p])
+                end
+            end
+
+            for in_nodes in Iterators.product(input_particles...)
+                # make the compute pair nodes for every combination of the found input_particle_nodes to get all spin/pol combinations
+
+                #insert_node!(graph, )
+
+            end
+            # make the collect pair and propagate nodes
+
+        end
+
+        data_out_propagated = insert_node!(
+            graph, make_node(DataTask(0)); track=false, invalidate_caches=false
+        )
+
+        pair_task_outputs[p] = data_out_propagated
+
+        # ... end do
+    end
+
+    return graph
 end