diff --git a/Project.toml b/Project.toml
index ae82e04..73ce957 100644
--- a/Project.toml
+++ b/Project.toml
@@ -19,11 +19,17 @@ QEDprocesses = "46de9c38-1bb3-4547-a1ec-da24d767fdad"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
 
 [extras]
 CUDA_Runtime_jll = "76a88914-d11a-5bdc-97e0-2f5a05c973a2"
+QEDbase = "10e22c08-3ccb-4172-bfcf-7d7aa3d04d93"
+QEDcore = "35dc0263-cb5f-4c33-a114-1d7f54ab753e"
+QEDprocesses = "46de9c38-1bb3-4547-a1ec-da24d767fdad"
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test"]
+test = ["SafeTestsets", "Test", "QEDbase", "QEDcore", "QEDprocesses"]
diff --git a/notebooks/abc_model_large.ipynb b/notebooks/abc_model_large.ipynb
index e98871c..ef75af5 100644
--- a/notebooks/abc_model_large.ipynb
+++ b/notebooks/abc_model_large.ipynb
@@ -413,7 +413,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Julia 1.10.2",
+   "display_name": "Julia 1.10.4",
    "language": "julia",
    "name": "julia-1.10"
   },
@@ -421,7 +421,7 @@
    "file_extension": ".jl",
    "mimetype": "application/julia",
    "name": "julia",
-   "version": "1.10.2"
+   "version": "1.10.4"
   }
  },
  "nbformat": 4,
diff --git a/src/MetagraphOptimization.jl b/src/MetagraphOptimization.jl
index 429b7d4..77a85eb 100644
--- a/src/MetagraphOptimization.jl
+++ b/src/MetagraphOptimization.jl
@@ -6,6 +6,8 @@ A module containing tools to work on DAGs.
 module MetagraphOptimization
 
 using QEDbase
+using QEDcore
+using QEDprocesses
 
 # graph types
 export DAG
@@ -64,8 +66,7 @@ export ComputeTaskABC_Sum
 
 # QED model
 export FeynmanDiagram, FeynmanVertex, FeynmanTie, FeynmanParticle
-export PhotonStateful, FermionStateful, AntiFermionStateful
-export QEDParticle, QEDProcessDescription, QEDProcessInput, QEDModel
+export GenericQEDProcess, QEDModel
 export ComputeTaskQED_P
 export ComputeTaskQED_S1
 export ComputeTaskQED_S2
@@ -114,6 +115,7 @@ import Base.delete!
 import Base.insert!
 import Base.collect
 
+include("QEDprocesses_patch.jl")
 
 include("devices/interface.jl")
 include("task/type.jl")
@@ -175,6 +177,7 @@ include("models/print.jl")
 
 include("models/physics_models/interface.jl")
 
+#=
 include("models/physics_models/abc/types.jl")
 include("models/physics_models/abc/particle.jl")
 include("models/physics_models/abc/compute.jl")
@@ -182,6 +185,7 @@ include("models/physics_models/abc/create.jl")
 include("models/physics_models/abc/properties.jl")
 include("models/physics_models/abc/parse.jl")
 include("models/physics_models/abc/print.jl")
+=#
 
 include("models/physics_models/qed/types.jl")
 include("models/physics_models/qed/particle.jl")
@@ -199,7 +203,7 @@ include("devices/impl.jl")
 include("devices/numa/impl.jl")
 include("devices/cuda/impl.jl")
 include("devices/rocm/impl.jl")
-include("devices/oneapi/impl.jl")
+#include("devices/oneapi/impl.jl")
 
 include("scheduler/interface.jl")
 include("scheduler/greedy.jl")
diff --git a/src/QEDprocesses_patch.jl b/src/QEDprocesses_patch.jl
new file mode 100644
index 0000000..1eba98c
--- /dev/null
+++ b/src/QEDprocesses_patch.jl
@@ -0,0 +1,71 @@
+# patch QEDprocesses
+# see issue https://github.com/QEDjl-project/QEDprocesses.jl/issues/77
+@inline function QEDprocesses.number_particles(
+    proc_def::QEDbase.AbstractProcessDefinition,
+    dir::DIR,
+    ::PT,
+) where {DIR <: QEDbase.ParticleDirection, PT <: QEDbase.AbstractParticleType}
+    return count(x -> x isa PT, particles(proc_def, dir))
+end
+
+@inline function QEDprocesses.number_particles(
+    proc_def::QEDbase.AbstractProcessDefinition,
+    ::PS,
+) where {
+    DIR <: QEDbase.ParticleDirection,
+    PT <: QEDbase.AbstractParticleType,
+    EL <: AbstractFourMomentum,
+    PS <: ParticleStateful{DIR, PT, EL},
+}
+    return QEDprocesses.number_particles(proc_def, DIR(), PT())
+end
+
+@inline function QEDprocesses.number_particles(
+    proc_def::QEDbase.AbstractProcessDefinition,
+    ::Type{PS},
+) where {
+    DIR <: QEDbase.ParticleDirection,
+    PT <: QEDbase.AbstractParticleType,
+    EL <: AbstractFourMomentum,
+    PS <: ParticleStateful{DIR, PT, EL},
+}
+    return QEDprocesses.number_particles(proc_def, DIR(), PT())
+end
+
+@inline function QEDcore.ParticleStateful{DIR, SPECIES}(
+    mom::AbstractFourMomentum,
+) where {DIR <: ParticleDirection, SPECIES <: AbstractParticleType}
+    return ParticleStateful(DIR(), SPECIES(), mom)
+end
+
+@inline function QEDcore.ParticleStateful{DIR, SPECIES, EL}(
+    mom::EL,
+) where {DIR <: ParticleDirection, SPECIES <: AbstractParticleType, EL <: AbstractFourMomentum}
+    return ParticleStateful(DIR(), SPECIES(), mom)
+end
+
+@inline function QEDbase.momentum(
+    psp::AbstractPhaseSpacePoint{MODEL, PROC, PS_DEF, INT, OUTT},
+    dir::ParticleDirection,
+    species::AbstractParticleType,
+    n::Int,
+) where {MODEL, PROC, PS_DEF, INT, OUTT}
+    # TODO: can be done through fancy template recursion too with 0 overhead
+    i = 0
+    c = n
+    for p in particles(psp, dir)
+        i += 1
+        if particle_species(p) isa typeof(species)
+            c -= 1
+        end
+        if c == 0
+            break
+        end
+    end
+
+    if c != 0 || n <= 0
+        throw(InvalidInputError("could not get $n-th momentum of $dir $species, does not exist"))
+    end
+
+    return momenta(psp, dir)[i]
+end
diff --git a/src/code_gen/function.jl b/src/code_gen/function.jl
index 588b350..fc9c2ba 100644
--- a/src/code_gen/function.jl
+++ b/src/code_gen/function.jl
@@ -1,10 +1,10 @@
 """
-    get_compute_function(graph::DAG, process::AbstractProcessDescription, machine::Machine)
+    get_compute_function(graph::DAG, instance, machine::Machine)
 
-Return a function of signature `compute_<id>(input::AbstractProcessInput)`, which will return the result of the DAG computation on the given input.
+Return a function of signature `compute_<id>(input::input_type(instance))`, which will return the result of the DAG computation on the given input.
 """
-function get_compute_function(graph::DAG, process::AbstractProcessDescription, machine::Machine)
-    tape = gen_tape(graph, process, machine)
+function get_compute_function(graph::DAG, instance, machine::Machine)
+    tape = gen_tape(graph, instance, machine)
 
     initCaches = Expr(:block, tape.initCachesCode...)
     assignInputs = Expr(:block, expr_from_fc.(tape.inputAssignCode)...)
@@ -13,7 +13,7 @@ function get_compute_function(graph::DAG, process::AbstractProcessDescription, m
     functionId = to_var_name(UUIDs.uuid1(rng[1]))
     resSym = eval(gen_access_expr(entry_device(tape.machine), tape.outputSymbol))
     expr = Meta.parse(
-        "function compute_$(functionId)(data_input::AbstractProcessInput) $(initCaches); $(assignInputs); $code; return $resSym; end",
+        "function compute_$(functionId)(data_input::$(input_type(instance))) $(initCaches); $(assignInputs); $code; return $resSym; end",
     )
 
     func = eval(expr)
@@ -22,12 +22,12 @@ function get_compute_function(graph::DAG, process::AbstractProcessDescription, m
 end
 
 """
-    get_cuda_kernel(graph::DAG, process::AbstractProcessDescription, machine::Machine)
+    get_cuda_kernel(graph::DAG, instance, machine::Machine)
 
 Return a function of signature `compute_<id>(input::CuVector, output::CuVector, n::Int64)`, which will return the result of the DAG computation of the input on the given output variable.
 """
-function get_cuda_kernel(graph::DAG, process::AbstractProcessDescription, machine::Machine)
-    tape = gen_tape(graph, process, machine)
+function get_cuda_kernel(graph::DAG, instance, machine::Machine)
+    tape = gen_tape(graph, instance, machine)
 
     initCaches = Expr(:block, tape.initCachesCode...)
     assignInputs = Expr(:block, expr_from_fc.(tape.inputAssignCode)...)
@@ -54,19 +54,19 @@ function get_cuda_kernel(graph::DAG, process::AbstractProcessDescription, machin
 end
 
 """
-    execute(graph::DAG, process::AbstractProcessDescription, machine::Machine, input::AbstractProcessInput)
+    execute(graph::DAG, instance, machine::Machine, input)
 
 Execute the code of the given `graph` on the given input values.
 
 This is essentially shorthand for
 ```julia
-tape = gen_tape(graph, process, machine)
+tape = gen_tape(graph, instance, machine)
 return execute_tape(tape, input)
 ```
 
 See also: [`parse_dag`](@ref), [`parse_process`](@ref), [`gen_process_input`](@ref)
 """
-function execute(graph::DAG, process::AbstractProcessDescription, machine::Machine, input::AbstractProcessInput)
-    tape = gen_tape(graph, process, machine)
+function execute(graph::DAG, instance, machine::Machine, input)
+    tape = gen_tape(graph, instance, machine)
     return execute_tape(tape, input)
 end
diff --git a/src/code_gen/tape_machine.jl b/src/code_gen/tape_machine.jl
index e986ad4..db8ade7 100644
--- a/src/code_gen/tape_machine.jl
+++ b/src/code_gen/tape_machine.jl
@@ -5,7 +5,7 @@ function call_fc(fc::FunctionCall{VectorT, 0}, cache::Dict{Symbol, Any}) where {
 end
 
 function call_fc(fc::FunctionCall{VectorT, 1}, cache::Dict{Symbol, Any}) where {VectorT <: SVector{1}}
-    cache[fc.return_symbol] = fc.func(fc.additional_arguments[1], cache[fc.arguments[1]])
+    cache[fc.return_symbol] = fc.func(fc.value_arguments[1], cache[fc.arguments[1]])
     return nothing
 end
 
@@ -15,12 +15,12 @@ function call_fc(fc::FunctionCall{VectorT, 0}, cache::Dict{Symbol, Any}) where {
 end
 
 function call_fc(fc::FunctionCall{VectorT, 1}, cache::Dict{Symbol, Any}) where {VectorT <: SVector{2}}
-    cache[fc.return_symbol] = fc.func(fc.additional_arguments[1], cache[fc.arguments[1]], cache[fc.arguments[2]])
+    cache[fc.return_symbol] = fc.func(fc.value_arguments[1], cache[fc.arguments[1]], cache[fc.arguments[2]])
     return nothing
 end
 
 function call_fc(fc::FunctionCall{VectorT, 1}, cache::Dict{Symbol, Any}) where {VectorT}
-    cache[fc.return_symbol] = fc.func(fc.additional_arguments[1], getindex.(Ref(cache), fc.arguments)...)
+    cache[fc.return_symbol] = fc.func(fc.value_arguments[1], getindex.(Ref(cache), fc.arguments)...)
     return nothing
 end
 
@@ -32,7 +32,7 @@ Execute the given [`FunctionCall`](@ref) on the dictionary.
 Several more specialized versions of this function exist to reduce vector unrolling work for common cases.
 """
 function call_fc(fc::FunctionCall{VectorT, M}, cache::Dict{Symbol, Any}) where {VectorT, M}
-    cache[fc.return_symbol] = fc.func(fc.additional_arguments..., getindex.(Ref(cache), fc.arguments)...)
+    cache[fc.return_symbol] = fc.func(fc.value_arguments..., getindex.(Ref(cache), fc.arguments)...)
     return nothing
 end
 
@@ -48,12 +48,8 @@ end
 For a given function call, return an expression evaluating it.
 """
 function expr_from_fc(fc::FunctionCall{VectorT, M}) where {VectorT, M}
-    func_call = Expr(
-        :call,
-        Symbol(fc.func),
-        fc.additional_arguments...,
-        eval.(gen_access_expr.(Ref(fc.device), fc.arguments))...,
-    )
+    func_call =
+        Expr(:call, Symbol(fc.func), fc.value_arguments..., eval.(gen_access_expr.(Ref(fc.device), fc.arguments))...)
 
     expr = :($(eval(gen_access_expr(fc.device, fc.return_symbol))) = $func_call)
     return expr
@@ -86,27 +82,19 @@ Return a `Vector{Expr}` doing the input assignments from the given `problemInput
 """
 function gen_input_assignment_code(
     inputSymbols::Dict{String, Vector{Symbol}},
-    instance::AbstractProblemInstance,
+    instance,
     machine::Machine,
     problemInputSymbol::Symbol = :input,
 )
     assignInputs = Vector{FunctionCall}()
     for (name, symbols) in inputSymbols
-        (type, index) = type_index_from_name(model(instance), name)
         # make a function for this, since we can't use anonymous functions in the FunctionCall
 
         for symbol in symbols
             device = entry_device(machine)
             push!(
                 assignInputs,
-                FunctionCall(
-                    # x is the process input
-                    part_from_x,
-                    SVector{2, Any}(type, index),
-                    SVector{1, Symbol}(problemInputSymbol),
-                    symbol,
-                    device,
-                ),
+                FunctionCall(input_expr, SVector{1, Any}(name), SVector{1, Symbol}(problemInputSymbol), symbol, device),
             )
         end
     end
@@ -121,12 +109,7 @@ Generate the code for a given graph. The return value is a [`Tape`](@ref).
 
 See also: [`execute`](@ref), [`execute_tape`](@ref)
 """
-function gen_tape(
-    graph::DAG,
-    instance::AbstractProblemInstance,
-    machine::Machine,
-    scheduler::AbstractScheduler = GreedyScheduler(),
-)
+function gen_tape(graph::DAG, instance, machine::Machine, scheduler::AbstractScheduler = GreedyScheduler())
     schedule = schedule_dag(scheduler, graph, machine)
 
     # get inSymbols
@@ -145,7 +128,7 @@ function gen_tape(
     initCaches = gen_cache_init_code(machine)
     assignInputs = gen_input_assignment_code(inputSyms, instance, machine, :input)
 
-    return Tape(initCaches, assignInputs, schedule, inputSyms, outSym, Dict(), instance, machine)
+    return Tape{input_type(instance)}(initCaches, assignInputs, schedule, inputSyms, outSym, Dict(), instance, machine)
 end
 
 """
diff --git a/src/code_gen/type.jl b/src/code_gen/type.jl
index 46e9fd1..db420fd 100644
--- a/src/code_gen/type.jl
+++ b/src/code_gen/type.jl
@@ -3,6 +3,8 @@
     Tape{INPUT}
 
 TODO: update docs
+- `INPUT` the input type of the problem instance
+
 - `code::Vector{Expr}`: The julia expression containing the code for the whole graph.
 - `inputSymbols::Dict{String, Vector{Symbol}}`: A dictionary of symbols mapping the names of the input nodes of the graph to the symbols their inputs should be provided on.
 - `outputSymbol::Symbol`: The symbol of the final calculated value
@@ -14,6 +16,6 @@ struct Tape{INPUT}
     inputSymbols::Dict{String, Vector{Symbol}}
     outputSymbol::Symbol
     cache::Dict{Symbol, Any}
-    instance::AbstractProblemInstance
+    instance::Any
     machine::Machine
 end
diff --git a/src/models/interface.jl b/src/models/interface.jl
index 3ed1b34..8bf7514 100644
--- a/src/models/interface.jl
+++ b/src/models/interface.jl
@@ -39,8 +39,8 @@ Generate the [`DAG`](@ref) for the given [`AbstractProblemInstance`](@ref). Ever
 function graph end
 
 """
-    input_expr(::AbstractProblemInstance, input_sym::Symbol, name::String)
+    input_expr(::AbstractProblemInstance, name::String, input)
 
-For a given [`AbstractProblemInstance`](@ref), a `Symbol` on which the problem input is available (see [`input_type`](@ref)), and entry node name, return an `Expr` getting that specific input value from the 
+For a given [`AbstractProblemInstance`](@ref), the problem input (of type `input_type(...)`) and an entry node name, return a that specific input value from the input symbol.
 """
 function input_expr end
diff --git a/src/models/physics_models/interface.jl b/src/models/physics_models/interface.jl
index 17b90fd..5654735 100644
--- a/src/models/physics_models/interface.jl
+++ b/src/models/physics_models/interface.jl
@@ -29,6 +29,7 @@ See also: [`parse_process`](@ref), [`AbstractProblemInstance`](@ref)
 abstract type AbstractProcessDescription end
 
 
+#TODO: i don't think giving this a base type is a good idea, the input type should just be returned of some function, allowing anything as an input type
 """
     AbstractProcessInput
 
diff --git a/src/models/physics_models/qed/compute.jl b/src/models/physics_models/qed/compute.jl
index b30994f..ce94d16 100644
--- a/src/models/physics_models/qed/compute.jl
+++ b/src/models/physics_models/qed/compute.jl
@@ -1,13 +1,9 @@
 using StaticArrays
 
-"""
-    compute(::ComputeTaskQED_P, data::QEDParticleValue)
+function input_expr(name::String, psp::PhaseSpacePoint)
+    (type, index) = type_index_from_name(QEDModel(), name)
 
-Return the particle as is and initialize the Value.
-"""
-function compute(::ComputeTaskQED_P, data::QEDParticleValue{P}) where {P <: QEDParticle}
-    # TODO do we actually need this for anything?
-    return ParticleValue{P, DiracMatrix}(data.p, one(DiracMatrix))
+    return ParticleValue(type(momentum(psp, particle_direction(type), particle_species(type), index)), 0.0im)
 end
 
 """
@@ -15,9 +11,9 @@ end
 
 Compute an outer edge. Return the particle value with the same particle and the value multiplied by an outer_edge factor.
 """
-function compute(::ComputeTaskQED_U, data::PV) where {P <: QEDParticle, PV <: QEDParticleValue{P}}
+function compute(::ComputeTaskQED_U, data::ParticleValue{P, V}) where {P <: ParticleStateful, V <: ValueType}
     part::P = data.p
-    state = base_state(particle(part), direction(part), momentum(part), spin_or_pol(part))
+    state = base_state(particle_species(part), particle_direction(part), momentum(part), spin_or_pol(part))
     return ParticleValue{P, typeof(state)}(
         data.p,
         state, # will return a SLorentzVector{ComplexF64}, BiSpinor or AdjointBiSpinor
@@ -31,9 +27,9 @@ Compute a vertex. Preserve momentum and particle types (e + gamma->p etc.) to cr
 """
 function compute(
     ::ComputeTaskQED_V,
-    data1::PV1,
-    data2::PV2,
-) where {P1 <: QEDParticle, P2 <: QEDParticle, PV1 <: QEDParticleValue{P1}, PV2 <: QEDParticleValue{P2}}
+    data1::ParticleValue{P1, V1},
+    data2::ParticleValue{P2, V2},
+) where {P1 <: ParticleStateful, P2 <: ParticleStateful, V1 <: ValueType, V2 <: ValueType}
     p3 = QED_conserve_momentum(data1.p, data2.p)
     P3 = interaction_result(P1, P2)
     state = QED_vertex()
@@ -63,9 +59,16 @@ For valid inputs, both input particles should have the same momenta at this poin
 """
 function compute(
     ::ComputeTaskQED_S2,
-    data1::ParticleValue{P1},
-    data2::ParticleValue{P2},
-) where {P1 <: Union{AntiFermionStateful, FermionStateful}, P2 <: Union{AntiFermionStateful, FermionStateful}}
+    data1::ParticleValue{ParticleStateful{D1, S1}, V1},
+    data2::ParticleValue{ParticleStateful{D2, S2}, V2},
+) where {
+    D1 <: ParticleDirection,
+    D2 <: ParticleDirection,
+    S1 <: Union{Electron, Positron},
+    S2 <: Union{Electron, Positron},
+    V1 <: ValueType,
+    V2 <: ValueType,
+}
     #@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(propagation_result(P1)(momentum(data1.p)))
@@ -80,9 +83,9 @@ end
 
 function compute(
     ::ComputeTaskQED_S2,
-    data1::ParticleValue{P1},
-    data2::ParticleValue{P2},
-) where {P1 <: PhotonStateful, P2 <: PhotonStateful}
+    data1::ParticleValue{ParticleStateful{D1, Photon}, V1},
+    data2::ParticleValue{ParticleStateful{D2, Photon}, V2},
+) where {D1 <: ParticleDirection, D2 <: ParticleDirection, V1 <: ValueType, V2 <: ValueType}
     # TODO: assert that data1 and data2 are opposites
     inner = QED_inner_edge(data1.p)
     # inner edge is just a scalar, data1 and data2 are photon states that are just Complex numbers here
@@ -90,11 +93,11 @@ function compute(
 end
 
 """
-    compute(::ComputeTaskQED_S1, data::QEDParticleValue)
+    compute(::ComputeTaskQED_S1, data::ParticleValue)
 
 Compute inner edge (1 input particle, 1 output particle).
 """
-function compute(::ComputeTaskQED_S1, data::QEDParticleValue{P}) where {P <: QEDParticle}
+function compute(::ComputeTaskQED_S1, data::ParticleValue{P, V}) where {P <: ParticleStateful, V <: ValueType}
     newP = propagation_result(P)
     new_p = newP(momentum(data.p))
     # inner edge is just a scalar, can multiply from either side
diff --git a/src/models/physics_models/qed/create.jl b/src/models/physics_models/qed/create.jl
index 85ee3aa..b8a8a50 100644
--- a/src/models/physics_models/qed/create.jl
+++ b/src/models/physics_models/qed/create.jl
@@ -1,7 +1,7 @@
 
 ComputeTaskQED_Sum() = ComputeTaskQED_Sum(0)
 
-function _svector_from_type(processDescription::QEDProcessDescription, type, particles)
+function _svector_from_type(processDescription::GenericQEDProcess, type, particles)
     if haskey(in_particles(processDescription), type)
         return SVector{in_particles(processDescription)[type], type}(filter(x -> typeof(x) <: type, particles))
     end
@@ -12,72 +12,48 @@ function _svector_from_type(processDescription::QEDProcessDescription, type, par
 end
 
 """
-    gen_process_input(processDescription::QEDProcessDescription)
+    gen_process_input(processDescription::GenericQEDProcess)
 
-Return a ProcessInput of randomly generated [`QEDParticle`](@ref)s from a [`QEDProcessDescription`](@ref). The process description can be created manually or parsed from a string using [`parse_process`](@ref).
+Return a ProcessInput of randomly generated [`QEDParticle`](@ref)s from a [`GenericQEDProcess`](@ref). The process description can be created manually or parsed from a string using [`parse_process`](@ref).
 
 Note: This uses RAMBO to create a valid process with conservation of momentum and energy.
 """
-function gen_process_input(processDescription::QEDProcessDescription)
+function gen_process_input(processDescription::GenericQEDProcess)
     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
 
     # add some extra random mass to allow for some momentum
     massSum += rand(rng[threadid()]) * (length(inputMasses) + length(outputMasses))
 
-
-    particles = Vector{QEDParticle}()
-    initialMomenta = generate_initial_moms(massSum, inputMasses)
-    index = 1
-    for (particle, n) in processDescription.inParticles
-        for _ in 1:n
-            mom = initialMomenta[index]
-            push!(particles, particle(mom))
-            index += 1
-        end
-    end
-
+    initial_momenta = generate_initial_moms(massSum, inputMasses)
     final_momenta = generate_physical_massive_moms(rng[threadid()], massSum, outputMasses)
-    index = 1
-    for (particle, n) in processDescription.outParticles
-        for _ in 1:n
-            push!(particles, particle(final_momenta[index]))
-            index += 1
-        end
-    end
 
-    inFerms = _svector_from_type(processDescription, FermionStateful{Incoming, SpinUp}, particles)
-    outFerms = _svector_from_type(processDescription, FermionStateful{Outgoing, SpinUp}, particles)
-    inAntiferms = _svector_from_type(processDescription, AntiFermionStateful{Incoming, SpinUp}, particles)
-    outAntiferms = _svector_from_type(processDescription, AntiFermionStateful{Outgoing, SpinUp}, particles)
-    inPhotons = _svector_from_type(processDescription, PhotonStateful{Incoming, PolX}, particles)
-    outPhotons = _svector_from_type(processDescription, PhotonStateful{Outgoing, PolX}, particles)
-
-    processInput =
-        QEDProcessInput(processDescription, inFerms, outFerms, inAntiferms, outAntiferms, inPhotons, outPhotons)
+    processInput = PhaseSpacePoint(
+        processDescription,
+        PerturbativeQED(),
+        PhasespaceDefinition(SphericalCoordinateSystem(), ElectronRestFrame()),
+        tuple(initial_momenta...),
+        tuple(final_momenta...),
+    )
 
     return processInput
 end
 
 """
-    gen_graph(process_description::QEDProcessDescription)
+    gen_graph(process_description::GenericQEDProcess)
 
-For a given [`QEDProcessDescription`](@ref), return the [`DAG`](@ref) that computes it.
+For a given [`GenericQEDProcess`](@ref), return the [`DAG`](@ref) that computes it.
 """
-function gen_graph(process_description::QEDProcessDescription)
+function gen_graph(process_description::GenericQEDProcess)
     initial_diagram = FeynmanDiagram(process_description)
     diagrams = gen_diagrams(initial_diagram)
 
diff --git a/src/models/physics_models/qed/diagrams.jl b/src/models/physics_models/qed/diagrams.jl
index 9477eb6..e3f94a6 100644
--- a/src/models/physics_models/qed/diagrams.jl
+++ b/src/models/physics_models/qed/diagrams.jl
@@ -8,10 +8,10 @@ import Base.show
 """
     FeynmanParticle
 
-Representation of a particle for use in [`FeynmanDiagram`](@ref)s. Consist of the [`QEDParticle`](@ref) type and an id.
+Representation of a particle for use in [`FeynmanDiagram`](@ref)s. Consist of the `ParticleStateful` type and an id.
 """
 struct FeynmanParticle
-    particle::Type{<:QEDParticle}
+    particle::Type{<:ParticleStateful}
     id::Int
 end
 
@@ -51,31 +51,21 @@ struct FeynmanDiagram
 end
 
 """
-    FeynmanDiagram(pd::QEDProcessDescription)
+    FeynmanDiagram(pd::GenericQEDProcess)
 
 Create an initial [`FeynmanDiagram`](@ref) with only its initial particles set and no vertices or ties.
 
 Use [`gen_diagrams`](@ref) to generate all possible diagrams from this one.
 """
-function FeynmanDiagram(pd::QEDProcessDescription)
+function FeynmanDiagram(pd::GenericQEDProcess)
     parts = Vector{FeynmanParticle}()
-    for (type, n) in pd.inParticles
-        for i in 1:n
-            push!(parts, FeynmanParticle(type, i))
-        end
-    end
-    for (type, n) in pd.outParticles
-        for i in 1:n
-            push!(parts, FeynmanParticle(type, i))
-        end
-    end
+
     ids = Dict{Type, Int64}()
-    for t in types(QEDModel())
-        if (isincoming(t))
-            ids[t] = get(pd.inParticles, t, 0)
-        else
-            ids[t] = get(pd.outParticles, t, 0)
+    for type in types(model(pd))
+        for i in 1:number_particles(pd, type)
+            push!(parts, FeynmanParticle(type, i))
         end
+        ids[type] = number_particles(pd, type)
     end
 
     return FeynmanDiagram([], missing, parts, ids)
@@ -83,7 +73,7 @@ end
 
 function particle_after_tie(p::FeynmanParticle, t::FeynmanTie)
     if p == t.in1 || p == t.in2
-        return FeynmanParticle(FermionStateful{Incoming, SpinUp}, -1) # placeholder particle and id for tied particles
+        return FeynmanParticle(ParticleStateful{Incoming, Electron, SFourMomentum}, -1) # placeholder particle and id for tied particles
     end
     return p
 end
@@ -114,7 +104,7 @@ end
 Return a string representation of the [`FeynmanParticle`](@ref) in a format that is readable by [`type_index_from_name`](@ref).
 """
 function String(p::FeynmanParticle)
-    return "$(String(p.particle))$(String(direction(p.particle)))$(p.id)"
+    return "$(String(p.particle))$(String(particle_direction(p.particle)))$(p.id)"
 end
 
 function hash(v::FeynmanVertex)
@@ -166,11 +156,11 @@ copy(fd::FeynmanDiagram) =
     FeynmanDiagram(deepcopy(fd.vertices), copy(fd.tie[]), deepcopy(fd.particles), copy(fd.type_ids))
 
 """
-    id_for_type(d::FeynmanDiagram, t::Type{<:QEDParticle})
+    id_for_type(d::FeynmanDiagram, t::Type{<:ParticleStateful})
 
 Return the highest id of any particle of the given type in the diagram + 1.
 """
-function id_for_type(d::FeynmanDiagram, t::Type{<:QEDParticle})
+function id_for_type(d::FeynmanDiagram, t::Type{<:ParticleStateful})
     return d.type_ids[t] + 1
 end
 
@@ -439,18 +429,19 @@ function remove_duplicates(compare_set::Set{FeynmanDiagram})
     return result
 end
 
+
 """
     is_compton(fd::FeynmanDiagram)
 
 Returns true iff the given feynman diagram is an (empty) diagram of a compton process like ke->k^ne
 """
 function is_compton(fd::FeynmanDiagram)
-    return fd.type_ids[FermionStateful{Incoming, SpinUp}] == 1 &&
-           fd.type_ids[FermionStateful{Outgoing, SpinUp}] == 1 &&
-           fd.type_ids[AntiFermionStateful{Incoming, SpinUp}] == 0 &&
-           fd.type_ids[AntiFermionStateful{Outgoing, SpinUp}] == 0 &&
-           fd.type_ids[PhotonStateful{Incoming, PolX}] >= 1 &&
-           fd.type_ids[PhotonStateful{Outgoing, PolX}] >= 1
+    return fd.type_ids[ParticleStateful{Incoming, Electron, SFourMomentum}] == 1 &&
+           fd.type_ids[ParticleStateful{Outgoing, Electron, SFourMomentum}] == 1 &&
+           fd.type_ids[ParticleStateful{Incoming, Positron, SFourMomentum}] == 0 &&
+           fd.type_ids[ParticleStateful{Outgoing, Positron, SFourMomentum}] == 0 &&
+           fd.type_ids[ParticleStateful{Incoming, Photon, SFourMomentum}] >= 1 &&
+           fd.type_ids[ParticleStateful{Outgoing, Photon, SFourMomentum}] >= 1
 end
 
 """
@@ -460,8 +451,8 @@ Helper function for [`gen_compton_diagrams`](@Ref). Generates a single diagram f
 """
 function gen_compton_diagram_from_order(order::Vector{Int}, inFerm, outFerm, n::Int, m::Int)
     photons = vcat(
-        [FeynmanParticle(PhotonStateful{Incoming, PolX}, i) for i in 1:n],
-        [FeynmanParticle(PhotonStateful{Outgoing, PolX}, i) for i in 1:m],
+        [FeynmanParticle(ParticleStateful{Incoming, Photon, SFourMomentum}, i) for i in 1:n],
+        [FeynmanParticle(ParticleStateful{Outgoing, Photon, SFourMomentum}, i) for i in 1:m],
     )
 
     new_diagram = FeynmanDiagram(
@@ -469,10 +460,10 @@ function gen_compton_diagram_from_order(order::Vector{Int}, inFerm, outFerm, n::
         missing,
         [inFerm, outFerm, photons...],
         Dict{Type, Int64}(
-            FermionStateful{Incoming, SpinUp} => 1,
-            FermionStateful{Outgoing, SpinUp} => 1,
-            PhotonStateful{Incoming, PolX} => n,
-            PhotonStateful{Outgoing, PolX} => m,
+            ParticleStateful{Incoming, Electron, SFourMomentum} => 1,
+            ParticleStateful{Outgoing, Electron, SFourMomentum} => 1,
+            ParticleStateful{Incoming, Photon, SFourMomentum} => n,
+            ParticleStateful{Outgoing, Photon, SFourMomentum} => m,
         ),
     )
 
@@ -484,9 +475,9 @@ function gen_compton_diagram_from_order(order::Vector{Int}, inFerm, outFerm, n::
     while left_index <= right_index
         # left side
         v_left = FeynmanVertex(
-            FeynmanParticle(FermionStateful{Incoming, SpinUp}, iterations),
+            FeynmanParticle(ParticleStateful{Incoming, Electron, SFourMomentum}, iterations),
             photons[order[left_index]],
-            FeynmanParticle(FermionStateful{Incoming, SpinUp}, iterations + 1),
+            FeynmanParticle(ParticleStateful{Incoming, Electron, SFourMomentum}, iterations + 1),
         )
         left_index += 1
         add_vertex!(new_diagram, v_left)
@@ -497,9 +488,9 @@ function gen_compton_diagram_from_order(order::Vector{Int}, inFerm, outFerm, n::
 
         # right side
         v_right = FeynmanVertex(
-            FeynmanParticle(FermionStateful{Outgoing, SpinUp}, iterations),
+            FeynmanParticle(ParticleStateful{Outgoing, Electron, SFourMomentum}, iterations),
             photons[order[right_index]],
-            FeynmanParticle(FermionStateful{Outgoing, SpinUp}, iterations + 1),
+            FeynmanParticle(ParticleStateful{Outgoing, Electron, SFourMomentum}, iterations + 1),
         )
         right_index -= 1
         add_vertex!(new_diagram, v_right)
@@ -512,7 +503,7 @@ function gen_compton_diagram_from_order(order::Vector{Int}, inFerm, outFerm, n::
     return new_diagram
 end
 
-
+#=
 """
     gen_compton_diagram_from_order_one_side(order::Vector{Int}, inFerm, outFerm, n::Int, m::Int)
 
@@ -520,8 +511,8 @@ Helper function for [`gen_compton_diagrams`](@Ref). Generates a single diagram f
 """
 function gen_compton_diagram_from_order_one_side(order::Vector{Int}, inFerm, outFerm, n::Int, m::Int)
     photons = vcat(
-        [FeynmanParticle(PhotonStateful{Incoming, PolX}, i) for i in 1:n],
-        [FeynmanParticle(PhotonStateful{Outgoing, PolX}, i) for i in 1:m],
+        [FeynmanParticle(ParticleStateful{Incoming, Photon}, i) for i in 1:n],
+        [FeynmanParticle(ParticleStateful{Outgoing, Photon}, i) for i in 1:m],
     )
 
     new_diagram = FeynmanDiagram(
@@ -529,10 +520,10 @@ function gen_compton_diagram_from_order_one_side(order::Vector{Int}, inFerm, out
         missing,
         [inFerm, outFerm, photons...],
         Dict{Type, Int64}(
-            FermionStateful{Incoming, SpinUp} => 1,
-            FermionStateful{Outgoing, SpinUp} => 1,
-            PhotonStateful{Incoming, PolX} => n,
-            PhotonStateful{Outgoing, PolX} => m,
+            ParticleStateful{Incoming, Electron} => 1,
+            ParticleStateful{Outgoing, Electron} => 1,
+            ParticleStateful{Incoming, Photon} => n,
+            ParticleStateful{Outgoing, Photon} => m,
         ),
     )
 
@@ -544,9 +535,9 @@ function gen_compton_diagram_from_order_one_side(order::Vector{Int}, inFerm, out
     while left_index <= right_index
         # left side
         v_left = FeynmanVertex(
-            FeynmanParticle(FermionStateful{Incoming, SpinUp}, iterations),
+            FeynmanParticle(ParticleStateful{Incoming, Electron}, iterations),
             photons[order[left_index]],
-            FeynmanParticle(FermionStateful{Incoming, SpinUp}, iterations + 1),
+            FeynmanParticle(ParticleStateful{Incoming, Electron}, iterations + 1),
         )
         left_index += 1
         add_vertex!(new_diagram, v_left)
@@ -559,9 +550,9 @@ function gen_compton_diagram_from_order_one_side(order::Vector{Int}, inFerm, out
         if (iterations == 1)
             # right side
             v_right = FeynmanVertex(
-                FeynmanParticle(FermionStateful{Outgoing, SpinUp}, iterations),
+                FeynmanParticle(ParticleStateful{Outgoing, Electron}, iterations),
                 photons[order[right_index]],
-                FeynmanParticle(FermionStateful{Outgoing, SpinUp}, iterations + 1),
+                FeynmanParticle(ParticleStateful{Outgoing, Electron}, iterations + 1),
             )
             right_index -= 1
             add_vertex!(new_diagram, v_right)
@@ -575,7 +566,7 @@ function gen_compton_diagram_from_order_one_side(order::Vector{Int}, inFerm, out
     add_tie!(new_diagram, FeynmanTie(ps[1], ps[2]))
     return new_diagram
 end
-
+=#
 
 """
     gen_compton_diagrams(n::Int, m::Int)
@@ -583,8 +574,8 @@ end
 Special case diagram generation for Compton processes, i.e., processes of the form k^ne->k^me
 """
 function gen_compton_diagrams(n::Int, m::Int)
-    inFerm = FeynmanParticle(FermionStateful{Incoming, SpinUp}, 1)
-    outFerm = FeynmanParticle(FermionStateful{Outgoing, SpinUp}, 1)
+    inFerm = FeynmanParticle(ParticleStateful{Incoming, Electron, SFourMomentum}, 1)
+    outFerm = FeynmanParticle(ParticleStateful{Outgoing, Electron, SFourMomentum}, 1)
 
     perms = [permutations([i for i in 1:(n + m)])...]
 
@@ -603,8 +594,8 @@ end
 Special case diagram generation for Compton processes, i.e., processes of the form k^ne->k^me, but generating from one end, yielding larger diagrams
 """
 function gen_compton_diagrams_one_side(n::Int, m::Int)
-    inFerm = FeynmanParticle(FermionStateful{Incoming, SpinUp}, 1)
-    outFerm = FeynmanParticle(FermionStateful{Outgoing, SpinUp}, 1)
+    inFerm = FeynmanParticle(ParticleStateful{Incoming, Electron, SFourMomentum}, 1)
+    outFerm = FeynmanParticle(ParticleStateful{Outgoing, Electron, SFourMomentum}, 1)
 
     perms = [permutations([i for i in 1:(n + m)])...]
 
@@ -623,12 +614,15 @@ From a given feynman diagram in its initial state, e.g. when created through the
 """
 function gen_diagrams(fd::FeynmanDiagram)
     if is_compton(fd)
-        return gen_compton_diagrams_one_side(
-            fd.type_ids[PhotonStateful{Incoming, PolX}],
-            fd.type_ids[PhotonStateful{Outgoing, PolX}],
+        return gen_compton_diagrams(
+            fd.type_ids[ParticleStateful{Incoming, Photon, SFourMomentum}],
+            fd.type_ids[ParticleStateful{Outgoing, Photon, SFourMomentum}],
         )
     end
 
+    throw(error("Unimplemented for non-compton!"))
+
+    #=
     working = Set{FeynmanDiagram}()
     results = Set{FeynmanDiagram}()
 
@@ -667,4 +661,5 @@ function gen_diagrams(fd::FeynmanDiagram)
     end
 
     return remove_duplicates(results)
+    =#
 end
diff --git a/src/models/physics_models/qed/parse.jl b/src/models/physics_models/qed/parse.jl
index ffcac74..31f2e85 100644
--- a/src/models/physics_models/qed/parse.jl
+++ b/src/models/physics_models/qed/parse.jl
@@ -5,8 +5,16 @@
 Parse a string representation of a process, such as "ke->ke" into the corresponding [`QEDProcessDescription`](@ref).
 """
 function parse_process(str::AbstractString, model::QEDModel)
-    inParticles = Dict{Type, Int}()
-    outParticles = Dict{Type, Int}()
+    inParticles = Dict(
+        ParticleStateful{Incoming, Photon, SFourMomentum} => 0,
+        ParticleStateful{Incoming, Electron, SFourMomentum} => 0,
+        ParticleStateful{Incoming, Positron, SFourMomentum} => 0,
+    )
+    outParticles = Dict(
+        ParticleStateful{Outgoing, Photon, SFourMomentum} => 0,
+        ParticleStateful{Outgoing, Electron, SFourMomentum} => 0,
+        ParticleStateful{Outgoing, Positron, SFourMomentum} => 0,
+    )
 
     if !(contains(str, "->"))
         throw("Did not find -> while parsing process \"$str\"")
@@ -19,14 +27,14 @@ function parse_process(str::AbstractString, model::QEDModel)
     end
 
     for t in types(model)
-        if (isincoming(t))
+        if (is_incoming(t))
             inCount = count(x -> x == String(t)[1], inStr)
 
             if inCount != 0
                 inParticles[t] = inCount
             end
         end
-        if (isoutgoing(t))
+        if (is_outgoing(t))
             outCount = count(x -> x == String(t)[1], outStr)
             if outCount != 0
                 outParticles[t] = outCount
@@ -40,5 +48,13 @@ function parse_process(str::AbstractString, model::QEDModel)
         throw("Encountered unknown characters in the output part of process \"$str\"")
     end
 
-    return QEDProcessDescription(inParticles, outParticles)
+
+    return GenericQEDProcess(
+        inParticles[ParticleStateful{Incoming, Photon, SFourMomentum}],
+        outParticles[ParticleStateful{Outgoing, Photon, SFourMomentum}],
+        inParticles[ParticleStateful{Incoming, Electron, SFourMomentum}],
+        outParticles[ParticleStateful{Outgoing, Electron, SFourMomentum}],
+        inParticles[ParticleStateful{Incoming, Positron, SFourMomentum}],
+        outParticles[ParticleStateful{Outgoing, Positron, SFourMomentum}],
+    )
 end
diff --git a/src/models/physics_models/qed/particle.jl b/src/models/physics_models/qed/particle.jl
index 130dc2e..f41fa8b 100644
--- a/src/models/physics_models/qed/particle.jl
+++ b/src/models/physics_models/qed/particle.jl
@@ -10,25 +10,70 @@ Singleton definition for identification of the QED-Model.
 """
 struct QEDModel <: AbstractPhysicsModel end
 
-"""
-    QEDProcessDescription <: AbstractProcessDescription
+QEDbase.is_incoming(::Type{<:ParticleStateful{Incoming}}) = true
+QEDbase.is_outgoing(::Type{<:ParticleStateful{Outgoing}}) = true
+QEDbase.is_incoming(::Type{<:ParticleStateful{Outgoing}}) = false
+QEDbase.is_outgoing(::Type{<:ParticleStateful{Incoming}}) = false
 
-A description of a process in the QED-Model. Contains the input and output particles.
+QEDbase.particle_direction(::Type{<:ParticleStateful{DIR}}) where {DIR <: ParticleDirection} = DIR()
+QEDbase.particle_species(
+    ::Type{<:ParticleStateful{DIR, SPECIES}},
+) where {DIR <: ParticleDirection, SPECIES <: AbstractParticleType} = SPECIES()
 
-See also: [`in_particles`](@ref), [`out_particles`](@ref), [`parse_process`](@ref)
-"""
-struct QEDProcessDescription <: AbstractProcessDescription
-    inParticles::Dict{Type{<:QEDParticle{Incoming}}, Int}
-    outParticles::Dict{Type{<:QEDParticle{Outgoing}}, Int}
+_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"))
+
+struct GenericQEDProcess{INT, OUTT} <: AbstractProcessDefinition where {INT <: Tuple, OUTT <: Tuple}
+    incoming_particles::INT
+    outgoing_particles::OUTT
+
+    function GenericQEDProcess(in_particles::INT, out_particles::OUTT) where {INT <: Tuple, OUTT <: Tuple}
+        _assert_particle_type_tuple(in_particles)
+        _assert_particle_type_tuple(out_particles)
+
+        return new{INT, OUTT}(in_particles, out_particles)
+    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.
+    """
+    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)
+    end
 end
 
-QEDParticleValue{ParticleType <: ParticleStateful} = Union{
-    ParticleValue{ParticleType, BiSpinor},
-    ParticleValue{ParticleType, AdjointBiSpinor},
-    ParticleValue{ParticleType, DiracMatrix},
-    ParticleValue{ParticleType, SLorentzVector{Float64}},
-    ParticleValue{ParticleType, ComplexF64},
-}
+QEDprocesses.incoming_particles(proc::GenericQEDProcess) = proc.incoming_particles
+QEDprocesses.outgoing_particles(proc::GenericQEDProcess) = proc.outgoing_particles
+
+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
+end
+
+function input_type(p::GenericQEDProcess)
+    in_t = QEDcore._assemble_tuple_type(incoming_particles(p), Incoming(), SFourMomentum)
+    out_t = QEDcore._assemble_tuple_type(outgoing_particles(p), Outgoing(), SFourMomentum)
+    return PhaseSpacePoint{
+        typeof(p),
+        PerturbativeQED,
+        PhasespaceDefinition{SphericalCoordinateSystem, ElectronRestFrame},
+        Tuple{in_t...},
+        Tuple{out_t...},
+    }
+end
+
+ValueType = Union{BiSpinor, AdjointBiSpinor, DiracMatrix, SLorentzVector{Float64}, ComplexF64}
 
 """
     interaction_result(t1::Type{T1}, t2::Type{T2}) where {T1 <: QEDParticle, T2 <: QEDParticle}
@@ -39,39 +84,53 @@ function interaction_result(t1::Type{T1}, t2::Type{T2}) where {T1 <: ParticleSta
     @assert false "Invalid interaction between particles of types $t1 and $t2"
 end
 
-interaction_result(::Type{ParticleStateful{Incoming, Electron}}, ::Type{ParticleStateful{Outgoing, Electron}}) =
-    ParticleStateful{Incoming, Photon}
-interaction_result(::Type{ParticleStateful{Incoming, Electron}}, ::Type{ParticleStateful{Incoming, Positron}}) =
-    ParticleStateful{Incoming, Photon}
 interaction_result(
-    ::Type{ParticleStateful{Incoming, Electron}},
-    ::Type{ParticleStateful{<:ParticleDirection, Photon}},
-) = ParticleStateful{Outgoing, Electron}
+    ::Type{ParticleStateful{Incoming, Electron, EL}},
+    ::Type{ParticleStateful{Outgoing, Electron, EL}},
+) where {EL <: AbstractFourMomentum} = ParticleStateful{Incoming, Photon, EL}
+interaction_result(
+    ::Type{ParticleStateful{Incoming, Electron, EL}},
+    ::Type{ParticleStateful{Incoming, Positron, EL}},
+) where {EL <: AbstractFourMomentum} = ParticleStateful{Incoming, Photon, EL}
+interaction_result(
+    ::Type{ParticleStateful{Incoming, Electron, EL}},
+    ::Type{ParticleStateful{DIR, Photon, EL}},
+) where {EL <: AbstractFourMomentum, DIR <: ParticleDirection} = ParticleStateful{Outgoing, Electron, EL}
 
-interaction_result(::Type{ParticleStateful{Outgoing, Electron}}, ::Type{ParticleStateful{Incoming, Electron}}) =
-    ParticleStateful{Incoming, Photon}
-interaction_result(::Type{ParticleStateful{Outgoing, Electron}}, ::Type{ParticleStateful{Outgoing, Positron}}) =
-    ParticleStateful{Incoming, Photon}
 interaction_result(
-    ::Type{ParticleStateful{Outgoing, Electron}},
-    ::Type{ParticleStateful{<:ParticleDirection, Photon}},
-) = ParticleStateful{Incoming, Electron}
+    ::Type{ParticleStateful{Outgoing, Electron, EL}},
+    ::Type{ParticleStateful{Incoming, Electron, EL}},
+) where {EL <: AbstractFourMomentum} = ParticleStateful{Incoming, Photon, EL}
+interaction_result(
+    ::Type{ParticleStateful{Outgoing, Electron, EL}},
+    ::Type{ParticleStateful{Outgoing, Positron, EL}},
+) where {EL <: AbstractFourMomentum} = ParticleStateful{Incoming, Photon, EL}
+interaction_result(
+    ::Type{ParticleStateful{Outgoing, Electron, EL}},
+    ::Type{ParticleStateful{DIR, Photon, EL}},
+) where {EL <: AbstractFourMomentum, DIR <: ParticleDirection} = ParticleStateful{Incoming, Electron, EL}
 
 # antifermion mirror
-interaction_result(::Type{ParticleStateful{Incoming, Positron}}, t2::Type{<:ParticleStateful}) =
-    interaction_result(ParticleStateful{Outgoing, Electron}, t2)
-interaction_result(::Type{ParticleStateful{Outgoing, Positron}}, t2::Type{<:ParticleStateful}) =
-    interaction_result(ParticleStateful{Incoming, Electron}, t2)
+interaction_result(
+    ::Type{ParticleStateful{Incoming, Positron, EL}},
+    t2::Type{<:ParticleStateful},
+) where {EL <: AbstractFourMomentum} = interaction_result(ParticleStateful{Outgoing, Electron, EL}, t2)
+interaction_result(
+    ::Type{ParticleStateful{Outgoing, Positron, EL}},
+    t2::Type{<:ParticleStateful},
+) where {EL <: AbstractFourMomentum} = interaction_result(ParticleStateful{Incoming, Electron, EL}, t2)
 
 # photon commutativity
-interaction_result(t1::Type{<:ParticleStateful{<:ParticleDirection, Photon}}, t2::Type{<:ParticleStateful}) =
-    interaction_result(t2, t1)
+interaction_result(
+    t1::Type{ParticleStateful{DIR, Photon, EL}},
+    t2::Type{<:ParticleStateful},
+) where {EL <: AbstractFourMomentum, DIR <: ParticleDirection} = interaction_result(t2, t1)
 
 # but prevent stack overflow
 function interaction_result(
-    t1::Type{ParticleStateful{<:ParticleDirection, Photon}},
-    t2::Type{ParticleStateful{<:ParticleDirection, Photon}},
-)
+    t1::Type{ParticleStateful{DIR1, Photon, EL}},
+    t2::Type{ParticleStateful{DIR2, Photon, EL}},
+) where {DIR1 <: ParticleDirection, DIR2 <: ParticleDirection, EL <: AbstractFourMomentum}
     @assert false "Invalid interaction between particles of types $t1 and $t2"
 end
 
@@ -80,12 +139,18 @@ end
 
 Return the type of the inverted direction. E.g.
 """
-propagation_result(::Type{ParticleStateful{Incoming, Electron}}) = ParticleStateful{Outgoing, Electron}
-propagation_result(::Type{ParticleStateful{Outgoing, Electron}}) = ParticleStateful{Incoming, Electron}
-propagation_result(::Type{ParticleStateful{Incoming, Positron}}) = ParticleStateful{Outgoing, Positron}
-propagation_result(::Type{ParticleStateful{Outgoing, Positron}}) = ParticleStateful{Incoming, Positron}
-propagation_result(::Type{ParticleStateful{Incoming, Photon}}) = ParticleStateful{Outgoing, Photon}
-propagation_result(::Type{ParticleStateful{Outgoing, Photon}}) = ParticleStateful{Incoming, Photon}
+propagation_result(::Type{ParticleStateful{Incoming, Electron, EL}}) where {EL <: AbstractFourMomentum} =
+    ParticleStateful{Outgoing, Electron, EL}
+propagation_result(::Type{ParticleStateful{Outgoing, Electron, EL}}) where {EL <: AbstractFourMomentum} =
+    ParticleStateful{Incoming, Electron, EL}
+propagation_result(::Type{ParticleStateful{Incoming, Positron, EL}}) where {EL <: AbstractFourMomentum} =
+    ParticleStateful{Outgoing, Positron, EL}
+propagation_result(::Type{ParticleStateful{Outgoing, Positron, EL}}) where {EL <: AbstractFourMomentum} =
+    ParticleStateful{Incoming, Positron, EL}
+propagation_result(::Type{ParticleStateful{Incoming, Photon, EL}}) where {EL <: AbstractFourMomentum} =
+    ParticleStateful{Outgoing, Photon, EL}
+propagation_result(::Type{ParticleStateful{Outgoing, Photon, EL}}) where {EL <: AbstractFourMomentum} =
+    ParticleStateful{Incoming, Photon, EL}
 
 """
     types(::QEDModel)
@@ -94,12 +159,12 @@ Return a Vector of the possible types of particle in the [`QEDModel`](@ref).
 """
 function types(::QEDModel)
     return [
-        ParticleStateful{Incoming, Photon},
-        ParticleStateful{Outgoing, Photon},
-        ParticleStateful{Incoming, Electron},
-        ParticleStateful{Outgoing, Electron},
-        ParticleStateful{Incoming, Positron},
-        ParticleStateful{Outgoing, Positron},
+        ParticleStateful{Incoming, Photon, SFourMomentum},
+        ParticleStateful{Outgoing, Photon, SFourMomentum},
+        ParticleStateful{Incoming, Electron, SFourMomentum},
+        ParticleStateful{Outgoing, Electron, SFourMomentum},
+        ParticleStateful{Incoming, Positron, SFourMomentum},
+        ParticleStateful{Outgoing, Positron, SFourMomentum},
     ]
 end
 
@@ -116,94 +181,39 @@ String(::Type{SpinDown}) = "spindown"
 String(::Incoming) = "i"
 String(::Outgoing) = "o"
 
-function String(::Type{<:PhotonStateful})
+function String(::Type{<:ParticleStateful{DIR, Photon}}) where {DIR <: ParticleDirection}
     return "k"
 end
-function String(::Type{<:FermionStateful})
+function String(::Type{<:ParticleStateful{DIR, Electron}}) where {DIR <: ParticleDirection}
     return "e"
 end
-function String(::Type{<:AntiFermionStateful})
+function String(::Type{<:ParticleStateful{DIR, Positron}}) where {DIR <: ParticleDirection}
     return "p"
 end
 
-function unique_name(::Type{PhotonStateful{Dir, Pol}}) where {Dir, Pol}
-    return String(PhotonStateful) * String(Dir) * String(Pol)
-end
-function unique_name(::Type{FermionStateful{Dir, Spin}}) where {Dir, Spin}
-    return String(FermionStateful) * String(Dir) * String(Spin)
-end
-function unique_name(::Type{AntiFermionStateful{Dir, Spin}}) where {Dir, Spin}
-    return String(AntiFermionStateful) * String(Dir) * String(Spin)
-end
-
-@inline particle(::PhotonStateful) = Photon()
-@inline particle(::FermionStateful) = Electron()
-@inline particle(::AntiFermionStateful) = Positron()
-
-@inline momentum(p::PhotonStateful)::SFourMomentum = p.momentum
-@inline momentum(p::FermionStateful)::SFourMomentum = p.momentum
-@inline momentum(p::AntiFermionStateful)::SFourMomentum = p.momentum
-
-@inline spin_or_pol(p::PhotonStateful{Dir, Pol}) where {Dir, Pol <: AbstractDefinitePolarization} = Pol()
-@inline spin_or_pol(p::FermionStateful{Dir, Spin}) where {Dir, Spin <: AbstractDefiniteSpin} = Spin()
-@inline spin_or_pol(p::AntiFermionStateful{Dir, Spin}) where {Dir, Spin <: AbstractDefiniteSpin} = Spin()
-
-@inline direction(
-    ::Type{P},
-) where {P <: Union{FermionStateful{Incoming}, AntiFermionStateful{Incoming}, PhotonStateful{Incoming}}} = Incoming()
-@inline direction(
-    ::Type{P},
-) where {P <: Union{FermionStateful{Outgoing}, AntiFermionStateful{Outgoing}, PhotonStateful{Outgoing}}} = Outgoing()
-
-@inline direction(
-    ::P,
-) where {P <: Union{FermionStateful{Incoming}, AntiFermionStateful{Incoming}, PhotonStateful{Incoming}}} = Incoming()
-@inline direction(
-    ::P,
-) where {P <: Union{FermionStateful{Outgoing}, AntiFermionStateful{Outgoing}, PhotonStateful{Outgoing}}} = Outgoing()
-
-@inline isincoming(::QEDParticle{Incoming}) = true
-@inline isincoming(::QEDParticle{Outgoing}) = false
-@inline isoutgoing(::QEDParticle{Incoming}) = false
-@inline isoutgoing(::QEDParticle{Outgoing}) = true
-
-@inline isincoming(::Type{<:QEDParticle{Incoming}}) = true
-@inline isincoming(::Type{<:QEDParticle{Outgoing}}) = false
-@inline isoutgoing(::Type{<:QEDParticle{Incoming}}) = false
-@inline isoutgoing(::Type{<:QEDParticle{Outgoing}}) = true
-
-@inline QEDbase.mass(::Type{<:FermionStateful}) = 1.0
-@inline QEDbase.mass(::Type{<:AntiFermionStateful}) = 1.0
-@inline QEDbase.mass(::Type{<:PhotonStateful}) = 0.0
-
-@inline invert_momentum(p::FermionStateful{Dir, Spin}) where {Dir, Spin} =
-    FermionStateful{Dir, Spin}(-p.momentum, p.spin)
-@inline invert_momentum(p::AntiFermionStateful{Dir, Spin}) where {Dir, Spin} =
-    AntiFermionStateful{Dir, Spin}(-p.momentum, p.spin)
-@inline invert_momentum(k::PhotonStateful{Dir, Spin}) where {Dir, Spin} =
-    PhotonStateful{Dir, Spin}(-k.momentum, k.polarization)
-
-
 """
-    caninteract(T1::Type{<:QEDParticle}, T2::Type{<:QEDParticle})
+    caninteract(T1::Type{<:ParticleStateful}, T2::Type{<:ParticleStateful})
 
-For two given [`QEDParticle`](@ref) types, return whether they can interact at a vertex. This is equivalent to `!issame(T1, T2)`.
+For two given `ParticleStateful` types, return whether they can interact at a vertex. This is equivalent to `!issame(T1, T2)`.
 
 See also: [`issame`](@ref) and [`interaction_result`](@ref)
 """
-function caninteract(T1::Type{<:QEDParticle}, T2::Type{<:QEDParticle})
+function caninteract(
+    T1::Type{<:ParticleStateful{D1, S1}},
+    T2::Type{<:ParticleStateful{D2, S2}},
+) where {D1 <: ParticleDirection, S1 <: AbstractParticleType, D2 <: ParticleDirection, S2 <: AbstractParticleType}
     if (T1 == T2)
         return false
     end
-    if (T1 <: PhotonStateful && T2 <: PhotonStateful)
+    if (S1 == Photon && S2 == Photon)
         return false
     end
 
     for (P1, P2) in [(T1, T2), (T2, T1)]
-        if (P1 <: FermionStateful{Incoming} && P2 <: AntiFermionStateful{Outgoing})
+        if (P1 <: ParticleStateful{Incoming, Electron} && P2 <: ParticleStateful{Outgoing, Positron})
             return false
         end
-        if (P1 <: FermionStateful{Outgoing} && P2 <: AntiFermionStateful{Incoming})
+        if (P1 <: ParticleStateful{Outgoing, Electron} && P2 <: ParticleStateful{Incoming, Positron})
             return false
         end
     end
@@ -213,30 +223,30 @@ end
 
 function type_index_from_name(::QEDModel, name::String)
     if startswith(name, "ki")
-        return (PhotonStateful{Incoming, PolX}, parse(Int, name[3:end]))
+        return (ParticleStateful{Incoming, Photon, SFourMomentum}, parse(Int, name[3:end]))
     elseif startswith(name, "ko")
-        return (PhotonStateful{Outgoing, PolX}, parse(Int, name[3:end]))
+        return (ParticleStateful{Outgoing, Photon, SFourMomentum}, parse(Int, name[3:end]))
     elseif startswith(name, "ei")
-        return (FermionStateful{Incoming, SpinUp}, parse(Int, name[3:end]))
+        return (ParticleStateful{Incoming, Electron, SFourMomentum}, parse(Int, name[3:end]))
     elseif startswith(name, "eo")
-        return (FermionStateful{Outgoing, SpinUp}, parse(Int, name[3:end]))
+        return (ParticleStateful{Outgoing, Electron, SFourMomentum}, parse(Int, name[3:end]))
     elseif startswith(name, "pi")
-        return (AntiFermionStateful{Incoming, SpinUp}, parse(Int, name[3:end]))
+        return (ParticleStateful{Incoming, Positron, SFourMomentum}, parse(Int, name[3:end]))
     elseif startswith(name, "po")
-        return (AntiFermionStateful{Outgoing, SpinUp}, parse(Int, name[3:end]))
+        return (ParticleStateful{Outgoing, Positron, SFourMomentum}, parse(Int, name[3:end]))
     else
         throw("Invalid name for a particle in the QED model")
     end
 end
 
 """
-    issame(T1::Type{<:QEDParticle}, T2::Type{<:QEDParticle})
+    issame(T1::Type{<:ParticleStateful}, T2::Type{<:ParticleStateful})
 
-For two given [`QEDParticle`](@ref) types, return whether they are equivalent for the purpose of a Feynman Diagram. That means e.g. an `Incoming` `AntiFermion` is the same as an `Outgoing` `Fermion`. This is equivalent to `!caninteract(T1, T2)`.
+For two given `ParticleStateful` types, return whether they are equivalent for the purpose of a Feynman Diagram. That means e.g. an `Incoming` `AntiFermion` is the same as an `Outgoing` `Fermion`. This is equivalent to `!caninteract(T1, T2)`.
 
 See also: [`caninteract`](@ref) and [`interaction_result`](@ref)
 """
-function issame(T1::Type{<:QEDParticle}, T2::Type{<:QEDParticle})
+function issame(T1::Type{<:ParticleStateful}, T2::Type{<:ParticleStateful})
     return !caninteract(T1, T2)
 end
 
@@ -250,12 +260,12 @@ Return the factor of a vertex in a QED feynman diagram.
     return -1im * e * gamma()
 end
 
-@inline function QED_inner_edge(p::QEDParticle)
-    return propagator(particle(p), p.momentum)
+@inline function QED_inner_edge(p::ParticleStateful)
+    return propagator(particle_species(p), momentum(p))
 end
 
 """
-    QED_conserve_momentum(p1::QEDParticle, p2::QEDParticle)
+    QED_conserve_momentum(p1::ParticleStateful, p2::ParticleStateful)
 
 Calculate and return a new particle from two given interacting ones at a vertex.
 """
@@ -265,43 +275,26 @@ function QED_conserve_momentum(
 ) where {
     Dir1 <: ParticleDirection,
     Dir2 <: ParticleDirection,
-    SpinPol1 <: AbstractSpinOrPolarization,
-    SpinPol2 <: AbstractSpinOrPolarization,
-    P1 <: Union{FermionStateful{Dir1, SpinPol1}, AntiFermionStateful{Dir1, SpinPol1}, PhotonStateful{Dir1, SpinPol1}},
-    P2 <: Union{FermionStateful{Dir2, SpinPol2}, AntiFermionStateful{Dir2, SpinPol2}, PhotonStateful{Dir2, SpinPol2}},
+    Species1 <: AbstractParticleType,
+    Species2 <: AbstractParticleType,
+    P1 <: ParticleStateful{Dir1, Species1},
+    P2 <: ParticleStateful{Dir2, Species2},
 }
     P3 = interaction_result(P1, P2)
-    p1_mom = p1.momentum
+    p1_mom = momentum(p1)
     if (Dir1 <: Outgoing)
         p1_mom *= -1
     end
-    p2_mom = p2.momentum
+    p2_mom = momentum(p2)
     if (Dir2 <: Outgoing)
         p2_mom *= -1
     end
 
     p3_mom = p1_mom + p2_mom
-    if (typeof(direction(P3)) <: Incoming)
-        return P3(-p3_mom)
+    if (particle_direction(P3) isa Incoming)
+        return ParticleStateful(particle_direction(P3), particle_species(P3), -p3_mom)
     end
-    return P3(p3_mom)
-end
-
-"""
-    QEDProcessInput <: AbstractProcessInput
-
-Input for a QED Process. Contains the [`QEDProcessDescription`](@ref) of the process it is an input for, and the values of the in and out particles.
-
-See also: [`gen_process_input`](@ref)
-"""
-struct QEDProcessInput{N1, N2, N3, N4, N5, N6} <: AbstractProcessInput
-    process::QEDProcessDescription
-    inFerms::SVector{N1, FermionStateful{Incoming, SpinUp}}
-    outFerms::SVector{N2, FermionStateful{Outgoing, SpinUp}}
-    inAntiferms::SVector{N3, AntiFermionStateful{Incoming, SpinUp}}
-    outAntiferms::SVector{N4, AntiFermionStateful{Outgoing, SpinUp}}
-    inPhotons::SVector{N5, PhotonStateful{Incoming, PolX}}
-    outPhotons::SVector{N6, PhotonStateful{Outgoing, PolX}}
+    return ParticleStateful(particle_direction(P3), particle_species(P3), p3_mom)
 end
 
 """
@@ -309,37 +302,22 @@ end
 
 Return the model of this process description.
 """
-model(::QEDProcessDescription) = QEDModel()
-model(::QEDProcessInput) = QEDModel()
+model(::GenericQEDProcess) = QEDModel()
+model(::PhaseSpacePoint) = QEDModel()
 
-function copy(process::QEDProcessDescription)
-    return QEDProcessDescription(copy(process.inParticles), copy(process.outParticles))
-end
-
-==(p1::QEDProcessDescription, p2::QEDProcessDescription) =
-    p1.inParticles == p2.inParticles && p1.outParticles == p2.outParticles
-
-function in_particles(process::QEDProcessDescription)
-    return process.inParticles
-end
-
-function out_particles(process::QEDProcessDescription)
-    return process.outParticles
-end
-
-function get_particle(input::QEDProcessInput, t::Type{Particle}, n::Int)::Particle where {Particle}
-    if (t <: FermionStateful{Incoming})
-        return input.inFerms[n]
-    elseif (t <: FermionStateful{Outgoing})
-        return input.outFerms[n]
-    elseif (t <: AntiFermionStateful{Incoming})
-        return input.inAntiferms[n]
-    elseif (t <: AntiFermionStateful{Outgoing})
-        return input.outAntiferms[n]
-    elseif (t <: PhotonStateful{Incoming})
-        return input.inPhotons[n]
-    elseif (t <: PhotonStateful{Outgoing})
-        return input.outPhotons[n]
+function get_particle(
+    input::PhaseSpacePoint,
+    t::Type{ParticleStateful{DIR, SPECIES}},
+    n::Int,
+) where {DIR <: ParticleDirection, SPECIES <: AbstractParticleType}
+    i = 0
+    for p in particles(input, DIR())
+        if p isa t
+            i += 1
+            if i == n
+                return p
+            end
+        end
     end
     @assert false "Invalid type given"
 end
diff --git a/src/models/physics_models/qed/print.jl b/src/models/physics_models/qed/print.jl
index c5cb3f5..2ac8ac4 100644
--- a/src/models/physics_models/qed/print.jl
+++ b/src/models/physics_models/qed/print.jl
@@ -1,8 +1,9 @@
 
+#=
 """
-    show(io::IO, process::QEDProcessDescription)
+    show(io::IO, process::GenericQEDProcess)
 
-Pretty print an [`QEDProcessDescription`](@ref) (no newlines).
+Pretty print an [`GenericQEDProcess`](@ref) (no newlines).
 
 ```jldoctest
 julia> using MetagraphOptimization
@@ -14,7 +15,7 @@ julia> print(parse_process("kk->ep", QEDModel()))
 QED Process: 'kk->ep'
 ```
 """
-function show(io::IO, process::QEDProcessDescription)
+function show(io::IO, process::GenericQEDProcess)
     # types() gives the types in order (QED) instead of random like keys() would
     print(io, "QED Process: \'")
     for type in types(QEDModel())
@@ -34,7 +35,7 @@ end
 
 
 """
-    String(process::QEDProcessDescription)
+    String(process::GenericQEDProcess)
 
 Create a short string suitable as a filename or similar, describing the given process.
 
@@ -64,7 +65,9 @@ function String(process::QEDProcessDescription)
     end
     return str
 end
+=#
 
+#=
 """
     show(io::IO, processInput::QEDProcessInput)
 
@@ -92,7 +95,9 @@ function show(io::IO, processInput::QEDProcessInput)
     end
     return nothing
 end
+=#
 
+#=
 """
     show(io::IO, particle::T) where {T <: QEDParticle}
 
@@ -102,13 +107,14 @@ function show(io::IO, particle::T) where {T <: QEDParticle}
     print(io, "$(String(typeof(particle))): $(particle.momentum)")
     return nothing
 end
+=#
 
 """
     show(io::IO, particle::FeynmanParticle)
 
 Pretty print a [`FeynmanParticle`](@ref) (no newlines).
 """
-show(io::IO, p::FeynmanParticle) = print(io, "$(String(p.particle))_$(String(direction(p.particle)))_$(p.id)")
+show(io::IO, p::FeynmanParticle) = print(io, "$(String(p.particle))_$(String(particle_direction(p.particle)))_$(p.id)")
 
 """
     show(io::IO, particle::FeynmanVertex)
diff --git a/src/task/properties.jl b/src/task/properties.jl
index c608843..2c4b320 100644
--- a/src/task/properties.jl
+++ b/src/task/properties.jl
@@ -3,28 +3,21 @@
 
 Fallback implementation of the compute function of a compute task, throwing an error.
 """
-function compute(t::AbstractTask, data...)
-    return error("Need to implement compute()")
-end
+function compute end
 
 """
     compute_effort(t::AbstractTask)
 
 Fallback implementation of the compute effort of a task, throwing an error.
 """
-function compute_effort(t::AbstractTask)::Float64
-    # default implementation using compute
-    return error("Need to implement compute_effort()")
-end
+function compute_effort end
 
 """
     data(t::AbstractTask)
 
 Fallback implementation of the data of a task, throwing an error.
 """
-function data(t::AbstractTask)::Float64
-    return error("Need to implement data()")
-end
+function data end
 
 """
     compute_effort(t::AbstractDataTask)
diff --git a/src/utility.jl b/src/utility.jl
index ac02bbf..a610663 100644
--- a/src/utility.jl
+++ b/src/utility.jl
@@ -1,3 +1,6 @@
+using Roots
+using ForwardDiff
+
 """
     noop()
 
@@ -249,7 +252,6 @@ end
 
 first_derivative(func) = x -> ForwardDiff.derivative(func, float(x))
 
-
 function generate_physical_massive_moms(rng, ss, masses; x0 = 0.1)
     n = length(masses)
     massless_moms = generate_physical_massless_moms(rng, ss, n)
diff --git a/test/Project.toml b/test/Project.toml
deleted file mode 100644
index 41c02e3..0000000
--- a/test/Project.toml
+++ /dev/null
@@ -1,10 +0,0 @@
-[deps]
-AccurateArithmetic = "22286c92-06ac-501d-9306-4abd417d9753"
-QEDbase = "10e22c08-3ccb-4172-bfcf-7d7aa3d04d93"
-QEDprocesses = "46de9c38-1bb3-4547-a1ec-da24d767fdad"
-Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
-StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
-StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-UUIDs = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
diff --git a/test/runtests.jl b/test/runtests.jl
index fa37b10..277f5de 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,5 +1,6 @@
 using SafeTestsets
 
+#=
 @safetestset "Utility Unit Tests                  " begin
     include("unit_tests_utility.jl")
 end
@@ -18,9 +19,11 @@ end
 @safetestset "ABC-Model Unit Tests                " begin
     include("unit_tests_abcmodel.jl")
 end
+=#
 @safetestset "QED-Model Unit Tests                " begin
     include("unit_tests_qedmodel.jl")
 end
+#=
 @safetestset "QED Feynman Diagram Generation Tests" begin
     include("unit_tests_qed_diagrams.jl")
 end
@@ -39,3 +42,4 @@ end
 @safetestset "Known Graph Tests                   " begin
     include("known_graphs.jl")
 end
+=#
\ No newline at end of file
diff --git a/test/unit_tests_qed_diagrams.jl b/test/unit_tests_qed_diagrams.jl
index 802fb50..cd3486b 100644
--- a/test/unit_tests_qed_diagrams.jl
+++ b/test/unit_tests_qed_diagrams.jl
@@ -22,7 +22,7 @@ trident = ("Trident", parse_process("ke->epe", model), 8)
 
         n_particles = 0
         for type in types(model)
-            if (isincoming(type))
+            if (is_incoming(type))
                 n_particles += get(process.inParticles, type, 0)
             else
                 n_particles += get(process.outParticles, type, 0)
diff --git a/test/unit_tests_qedmodel.jl b/test/unit_tests_qedmodel.jl
index 8126605..6a9defe 100644
--- a/test/unit_tests_qedmodel.jl
+++ b/test/unit_tests_qedmodel.jl
@@ -1,5 +1,6 @@
 using MetagraphOptimization
 using QEDbase
+using QEDcore
 using QEDprocesses
 using StatsBase # for countmap
 using Random
@@ -9,8 +10,6 @@ import MetagraphOptimization.caninteract
 import MetagraphOptimization.issame
 import MetagraphOptimization.interaction_result
 import MetagraphOptimization.propagation_result
-import MetagraphOptimization.direction
-import MetagraphOptimization.spin_or_pol
 import MetagraphOptimization.QED_vertex
 
 def_momentum = SFourMomentum(1.0, 0.0, 0.0, 0.0)
@@ -18,32 +17,32 @@ def_momentum = SFourMomentum(1.0, 0.0, 0.0, 0.0)
 RNG = Random.MersenneTwister(0)
 
 testparticleTypes = [
-    PhotonStateful{Incoming, PolX},
-    PhotonStateful{Outgoing, PolX},
-    FermionStateful{Incoming, SpinUp},
-    FermionStateful{Outgoing, SpinUp},
-    AntiFermionStateful{Incoming, SpinUp},
-    AntiFermionStateful{Outgoing, SpinUp},
+    ParticleStateful{Incoming, Photon, SFourMomentum},
+    ParticleStateful{Outgoing, Photon, SFourMomentum},
+    ParticleStateful{Incoming, Electron, SFourMomentum},
+    ParticleStateful{Outgoing, Electron, SFourMomentum},
+    ParticleStateful{Incoming, Positron, SFourMomentum},
+    ParticleStateful{Outgoing, Positron, SFourMomentum},
 ]
 
 testparticleTypesPropagated = [
-    PhotonStateful{Outgoing, PolX},
-    PhotonStateful{Incoming, PolX},
-    FermionStateful{Outgoing, SpinUp},
-    FermionStateful{Incoming, SpinUp},
-    AntiFermionStateful{Outgoing, SpinUp},
-    AntiFermionStateful{Incoming, SpinUp},
+    ParticleStateful{Outgoing, Photon, SFourMomentum},
+    ParticleStateful{Incoming, Photon, SFourMomentum},
+    ParticleStateful{Outgoing, Electron, SFourMomentum},
+    ParticleStateful{Incoming, Electron, SFourMomentum},
+    ParticleStateful{Outgoing, Positron, SFourMomentum},
+    ParticleStateful{Incoming, Positron, SFourMomentum},
 ]
 
-function compton_groundtruth(input::QEDProcessInput)
+function compton_groundtruth(input::PhaseSpacePoint)
     # p1k1 -> p2k2
     # formula: −(ie)^2 (u(p2) slashed(ε1) S(p2 − k1) slashed(ε2) u(p1) + u(p2) slashed(ε2) S(p1 + k1) slashed(ε1) u(p1))
 
-    p1 = input.inFerms[1]
-    p2 = input.outFerms[1]
+    p1 = momentum(psp, Incoming(), 2)
+    p2 = momentum(psp, Outgoing(), 2)
 
-    k1 = input.inPhotons[1]
-    k2 = input.outPhotons[1]
+    k1 = momentum(psp, Incoming(), 1)
+    k2 = momentum(psp, Outgoing(), 1)
 
     u_p1 = base_state(Electron(), Incoming(), p1.momentum, spin_or_pol(p1))
     u_p2 = base_state(Electron(), Outgoing(), p2.momentum, spin_or_pol(p2))
@@ -88,8 +87,8 @@ end
         @test issame(typeof(resultParticle), interaction_result(p1, p2))
 
         totalMom = zero(SFourMomentum)
-        for (p, mom) in [(p1, testParticle1.momentum), (p2, testParticle2.momentum), (p3, resultParticle.momentum)]
-            if (typeof(direction(p)) <: Incoming)
+        for (p, mom) in [(p1, momentum(testParticle1)), (p2, momentum(testParticle2)), (p3, momentum(resultParticle))]
+            if (typeof(particle_direction(p)) <: Incoming)
                 totalMom += mom
             else
                 totalMom -= mom
@@ -103,7 +102,7 @@ end
 @testset "Propagation Result" begin
     for (p, propResult) in zip(testparticleTypes, testparticleTypesPropagated)
         @test issame(propagation_result(p), propResult)
-        @test direction(propagation_result(p)(def_momentum)) != direction(p(def_momentum))
+        @test particle_direction(propagation_result(p)(def_momentum)) != particle_direction(p(def_momentum))
     end
 end
 
@@ -117,38 +116,23 @@ end
     end
 
     @testset "Known processes" begin
-        compton_process = QEDProcessDescription(
-            Dict{Type, Int}(PhotonStateful{Incoming, PolX} => 1, FermionStateful{Incoming, SpinUp} => 1),
-            Dict{Type, Int}(PhotonStateful{Outgoing, PolX} => 1, FermionStateful{Outgoing, SpinUp} => 1),
-        )
+        compton_process = GenericQEDProcess(1, 1, 1, 1, 0, 0)
 
         @test parse_process("ke->ke", QEDModel()) == compton_process
 
-        positron_compton_process = QEDProcessDescription(
-            Dict{Type, Int}(PhotonStateful{Incoming, PolX} => 1, AntiFermionStateful{Incoming, SpinUp} => 1),
-            Dict{Type, Int}(PhotonStateful{Outgoing, PolX} => 1, AntiFermionStateful{Outgoing, SpinUp} => 1),
-        )
+        positron_compton_process = GenericQEDProcess(1, 1, 0, 0, 1, 1)
 
         @test parse_process("kp->kp", QEDModel()) == positron_compton_process
 
-        trident_process = QEDProcessDescription(
-            Dict{Type, Int}(PhotonStateful{Incoming, PolX} => 1, FermionStateful{Incoming, SpinUp} => 1),
-            Dict{Type, Int}(FermionStateful{Outgoing, SpinUp} => 2, AntiFermionStateful{Outgoing, SpinUp} => 1),
-        )
+        trident_process = GenericQEDProcess(1, 0, 1, 2, 0, 1)
 
         @test parse_process("ke->eep", QEDModel()) == trident_process
 
-        pair_production_process = QEDProcessDescription(
-            Dict{Type, Int}(PhotonStateful{Incoming, PolX} => 2),
-            Dict{Type, Int}(FermionStateful{Outgoing, SpinUp} => 1, AntiFermionStateful{Outgoing, SpinUp} => 1),
-        )
+        pair_production_process = GenericQEDProcess(2, 0, 0, 1, 0, 1)
 
         @test parse_process("kk->pe", QEDModel()) == pair_production_process
 
-        pair_annihilation_process = QEDProcessDescription(
-            Dict{Type, Int}(FermionStateful{Incoming, SpinUp} => 1, AntiFermionStateful{Incoming, SpinUp} => 1),
-            Dict{Type, Int}(PhotonStateful{Outgoing, PolX} => 2),
-        )
+        pair_annihilation_process = GenericQEDProcess(0, 2, 1, 0, 1, 0)
 
         @test parse_process("pe->kk", QEDModel()) == pair_annihilation_process
     end
@@ -161,12 +145,18 @@ end
 
         for i in 1:100
             input = gen_process_input(process)
-            @test length(input.inFerms) == get(process.inParticles, FermionStateful{Incoming, SpinUp}, 0)
-            @test length(input.inAntiferms) == get(process.inParticles, AntiFermionStateful{Incoming, SpinUp}, 0)
-            @test length(input.inPhotons) == get(process.inParticles, PhotonStateful{Incoming, PolX}, 0)
-            @test length(input.outFerms) == get(process.outParticles, FermionStateful{Outgoing, SpinUp}, 0)
-            @test length(input.outAntiferms) == get(process.outParticles, AntiFermionStateful{Outgoing, SpinUp}, 0)
-            @test length(input.outPhotons) == get(process.outParticles, PhotonStateful{Outgoing, PolX}, 0)
+            @test length(input.inFerms) ==
+                  get(process.inParticles, ParticleStateful{Incoming, Electron, SFourMomentum}, 0)
+            @test length(input.inAntiferms) ==
+                  get(process.inParticles, ParticleStateful{Incoming, Positron, SFourMomentum}, 0)
+            @test length(input.inPhotons) ==
+                  get(process.inParticles, ParticleStateful{Incoming, Photon, SFourMomentum}, 0)
+            @test length(input.outFerms) ==
+                  get(process.outParticles, ParticleStateful{Outgoing, Electron, SFourMomentum}, 0)
+            @test length(input.outAntiferms) ==
+                  get(process.outParticles, ParticleStateful{Outgoing, Positron, SFourMomentum}, 0)
+            @test length(input.outPhotons) ==
+                  get(process.outParticles, ParticleStateful{Outgoing, Photon, SFourMomentum}, 0)
 
             @test isapprox(
                 sum([
@@ -185,6 +175,7 @@ end
     end
 end
 
+#=
 @testset "Compton" begin
     import MetagraphOptimization.insert_node!
     import MetagraphOptimization.insert_edge!
@@ -347,3 +338,4 @@ end
         @test isapprox(compute_function.(input), reduced_compute_function.(input))
     end
 end
+=#
\ No newline at end of file