diff --git a/examples/congruent_in_ph.jl b/examples/congruent_in_ph.jl
index a0639bf..1c91757 100644
--- a/examples/congruent_in_ph.jl
+++ b/examples/congruent_in_ph.jl
@@ -1,8 +1,10 @@
 using MetagraphOptimization
 using QEDbase
 using QEDcore
+using QEDprocesses
 using Random
 using UUIDs
+using CSV
 
 RNG = Random.MersenneTwister(123)
 
@@ -21,103 +23,79 @@ function mock_machine()
     )
 end
 
-function congruent_input(processDescription::QEDProcessDescription, omega::Number)
+function congruent_input_momenta(processDescription::GenericQEDProcess, omega::Number)
     # generate an input sample for given e + nk -> e' + k' process, where the nk are equal
     massSum = 0
     inputMasses = Vector{Float64}()
-    for (particle, n) in processDescription.inParticles
-        for _ in 1:n
-            massSum += mass(particle)
-            push!(inputMasses, mass(particle))
-        end
+    for particle in incoming_particles(processDescription)
+        massSum += mass(particle)
+        push!(inputMasses, mass(particle))
     end
     outputMasses = Vector{Float64}()
-    for (particle, n) in processDescription.outParticles
-        for _ in 1:n
-            massSum += mass(particle)
-            push!(outputMasses, mass(particle))
-        end
+    for particle in outgoing_particles(processDescription)
+        massSum += mass(particle)
+        push!(outputMasses, mass(particle))
     end
 
-    initialMomenta = [
-        i == 1 ? SFourMomentum(1, 0, 0, 0) : SFourMomentum(omega, 0, 0, omega) for
-        i in 1:length(inputMasses)
-    ]
+    initial_momenta =
+        [i == 1 ? SFourMomentum(1, 0, 0, 0) : SFourMomentum(omega, 0, 0, omega) for i in 1:length(inputMasses)]
 
     # add some extra random mass to allow for some momentum
-    ss = sqrt(sum(initialMomenta) * sum(initialMomenta))
+    ss = sqrt(sum(initial_momenta) * sum(initial_momenta))
+    final_momenta = MetagraphOptimization.generate_physical_massive_moms(RNG, ss, outputMasses)
 
-    result = Vector{QEDProcessInput}()
-    sizehint!(result, 16)
-
-    spin_pol_combinations = Iterators.product(
-        [SpinUp, SpinDown], [SpinUp, SpinDown], [PolX, PolY], [PolX, PolY]
-    )
-    for (in_spin, out_spin, in_pol, out_pol) in spin_pol_combinations
-
-        # get the electron first, then the n photons
-        particles = Vector{QEDParticle}()
-
-        for (particle, n) in processDescription.inParticles
-            if particle <: FermionStateful
-                mom = initialMomenta[1]
-                push!(particles, particle(mom, in_spin()))
-            elseif particle <: PhotonStateful
-                for i in 1:n
-                    mom = initialMomenta[i + 1]
-                    push!(particles, particle(mom, in_pol()))
-                end
-            else
-                @assert false
-            end
-        end
-
-        final_momenta = MetagraphOptimization.generate_physical_massive_moms(
-            RNG, ss, outputMasses
-        )
-        index = 1
-        for (particle, n) in processDescription.outParticles
-            for _ in 1:n
-                if particle <: FermionStateful
-                    push!(particles, particle(final_momenta[index], out_spin()))
-                elseif particle <: PhotonStateful
-                    push!(particles, particle(final_momenta[index], out_pol()))
-                end
-                index += 1
-            end
-        end
-
-        inFerms = MetagraphOptimization._svector_from_type(
-            processDescription, FermionStateful{Incoming,in_spin}, particles
-        )
-        outFerms = MetagraphOptimization._svector_from_type(
-            processDescription, FermionStateful{Outgoing,out_spin}, particles
-        )
-        inAntiferms = MetagraphOptimization._svector_from_type(
-            processDescription, AntiFermionStateful{Incoming,in_spin}, particles
-        )
-        outAntiferms = MetagraphOptimization._svector_from_type(
-            processDescription, AntiFermionStateful{Outgoing,out_spin}, particles
-        )
-        inPhotons = MetagraphOptimization._svector_from_type(
-            processDescription, PhotonStateful{Incoming,in_pol}, particles
-        )
-        outPhotons = MetagraphOptimization._svector_from_type(
-            processDescription, PhotonStateful{Outgoing,out_pol}, particles
-        )
-
-        processInput = QEDProcessInput(
-            processDescription,
-            inFerms,
-            outFerms,
-            inAntiferms,
-            outAntiferms,
-            inPhotons,
-            outPhotons,
-        )
-
-        push!(result, processInput)
-    end
-
-    return result
+    return (tuple(initial_momenta...), tuple(final_momenta...))
+end
+
+function build_psp(processDescription::GenericQEDProcess, momenta)
+    return PhaseSpacePoint(
+        processDescription,
+        PerturbativeQED(),
+        PhasespaceDefinition(SphericalCoordinateSystem(), ElectronRestFrame()),
+        momenta[1],
+        momenta[2],
+    )
+end
+
+n = 1000000
+photons = 4
+omega = 2e-3
+
+# n is the number of incoming photons
+# omega is the number
+
+println("Generating $n inputs for $photons photons, omega=$omega...")
+
+# temp process to generate momenta
+temp_process = parse_process("k"^photons * "e->ke", QEDModel(), PolX(), SpinUp(), PolX(), SpinUp())
+input_momenta = [congruent_input_momenta(temp_process, omega) for _ in 1:n]
+results = [0.0im for _ in 1:n]
+
+for (in_pol, in_spin, out_pol, out_spin) in
+    Iterators.product([PolX(), PolY()], [SpinUp(), SpinDown()], [PolX(), PolY()], [SpinUp(), SpinDown()])
+
+    print("Calculating for spin/pol config: $in_pol, $in_spin -> $out_pol, $out_spin... Preparing inputs... ")
+    process = parse_process("k"^photons * "e->ke", QEDModel(), in_pol, in_spin, out_pol, out_spin)
+    inputs = build_psp.(Ref(process), input_momenta)
+
+    print("Preparing graph... ")
+    # prepare function
+    graph = gen_graph(process)
+    optimize_to_fixpoint!(ReductionOptimizer(), graph)
+    print("Preparing function... ")
+    func = get_compute_function(graph, process, mock_machine())
+
+    print("Calculating... ")
+    for i in 1:n
+        results[i] += func(inputs[i])^2
+    end
+    println("Done.")
+end
+
+println("Writing results")
+open("$(photons)_congruent_photons_omega_$(omega).csv", "w") do f
+    println(f, "out_photon_momentum;out_electron_momentum;result")
+    for (momentum, result) in Iterators.zip(input_momenta, results)
+        println(f, "$(momentum[2][1]);$(momentum[2][2]);$(result)")
+    end
 end