n-Pair Trident Scattering Process

In this file, we set up an n-pair trident scattering process. A trident process looks like $k e^- \to e^- (e^- e^+)^n$.

You can download this file as a jupyter notebook.

using QEDFeynmanDiagrams

We need QEDcore of the QEDjl-project for base functionality and a process type, for which we can use the Mocks submodule from QEDbase for this tutorial. Downstream, a ScatteringProcess from QEDprocesses.jl could be used, for example.

using QEDcore
using QEDbase.Mocks

Let's decide how many pairs our trident should produce:

n = 2;

Now we setup the scattering process accordingly. We only consider a single spin/polarization combination here. For an example with many spin and polarization combinations, refer to the n-photon Compton example

proc = QEDbase.Mocks.MockProcessSP(
    (Electron(), Photon()),                                                         # incoming particles
    (Electron(), ntuple(_ -> Electron(), n)..., ntuple(_ -> Positron(), n)...),     # outgoing particles
    (SpinUp(), PolX()),                                                             # incoming particle spin/pols
    (SpinUp(), ntuple(_ -> SpinUp(), 2 * n)...),                                    # outgoing particle spin/pols
)

QEDbase.Mocks.MockProcessSP{Tuple{Electron, Photon}, Tuple{Electron, Electron, Electron, Positron, Positron}, Tuple{SpinUp, PolarizationX}, NTuple{5, SpinUp}}((electron, photon), (electron, electron, electron, positron, positron), (spin up, x-polarized), (spin up, spin up, spin up, spin up, spin up))

The number_of_diagrams function returns how many valid Feynman diagrams there are for a given process.

number_of_diagrams(proc)

252

Next, we can generate the DAG representing the computation for our scattering process' squared matrix element. This uses ComputableDAGs.jl.

dag = graph(proc)

Graph:
  Nodes: Total: 943, QEDFeynmanDiagrams.ComputeTask_Pair: 60, QEDFeynmanDiagrams.ComputeTask_BaseState: 7, 
         QEDFeynmanDiagrams.ComputeTask_SpinPolCumulation: 1, QEDFeynmanDiagrams.ComputeTask_Propagator: 60, DataTask: 505, 
         QEDFeynmanDiagrams.ComputeTask_PairNegated: 9, QEDFeynmanDiagrams.ComputeTask_TripleNegated: 150, QEDFeynmanDiagrams.ComputeTask_PropagatePairs: 60, 
         QEDFeynmanDiagrams.ComputeTask_CollectPairs: 60, QEDFeynmanDiagrams.ComputeTask_CollectTriples: 1, QEDFeynmanDiagrams.ComputeTask_Triple: 30
  Edges: 1553
  Total Compute Effort: 0.0
  Total Data Transfer: 0.0
  Total Compute Intensity: 0.0

To continue, we will need ComputableDAGs.jl. Since ComputableDAGs.jl uses RuntimeGeneratedFunctions as the return type of ComputableDAGs.get_compute_function, we need to initialize it in our current module.

using ComputableDAGs
using RuntimeGeneratedFunctions
RuntimeGeneratedFunctions.init(@__MODULE__)

Now we need an input for the function, which is a QEDcore.PhaseSpacePoint. For now, we generate random momenta for every particle. In the future, QEDevents will be able to generate physical PhaseSpacePoints.

psp = PhaseSpacePoint(
    proc,
    MockModel(),
    FlatPhaseSpaceLayout(TwoBodyRestSystem()),
    tuple((rand(SFourMomentum) for _ in 1:number_incoming_particles(proc))...),
    tuple((rand(SFourMomentum) for _ in 1:number_outgoing_particles(proc))...),
)

PhaseSpacePoint:
    process: QEDbase.Mocks.MockProcessSP{Tuple{Electron, Photon}, Tuple{Electron, Electron, Electron, Positron, Positron}, Tuple{SpinUp, PolarizationX}, NTuple{5, SpinUp}}((electron, photon), (electron, electron, electron, positron, positron), (spin up, x-polarized), (spin up, spin up, spin up, spin up, spin up))
    model: mock model
    phase space layout: FlatPhaseSpaceLayout{TwoBodyTargetSystem{Energy{2}}}(TwoBodyTargetSystem{Energy{2}}(Energy{2}()))
    incoming particles:
     -> incoming electron: [0.889885967374767, 0.7146083470091047, 0.3118198268944159, 0.6300075046341156]
     -> incoming photon: [0.004460596024637997, 0.9592628771528654, 0.09347519967629536, 0.13461144493888721]
    outgoing particles:
     -> outgoing electron: [0.03334292571676678, 0.22778306403086812, 0.2151782025115494, 0.08997890356381799]
     -> outgoing electron: [0.00204232593156628, 0.08730200210381789, 0.013681964605307462, 0.7031998120666759]
     -> outgoing electron: [0.3597173549181092, 0.8058116217905882, 0.4084957521877135, 0.661622674503487]
     -> outgoing positron: [0.7393014125724023, 0.28162302337975775, 0.4510910215969175, 0.6297054396315858]
     -> outgoing positron: [0.42333293511317915, 0.9792108559373601, 0.5439388426171936, 0.16987140644018894]

With the DAG, the process, RuntimeGeneratedFunctions initialized, and an input type to use we can now generate the actual computable function:

func = get_compute_function(
    dag, proc, cpu_st(), @__MODULE__; concrete_input_type = typeof(psp)
);

Finally, we can test that the function actually runs and computes something by simply calling it on the PhaseSpacePoint:

func(psp)

0.0075088249397024355

We can benchmark the execution speed too:

using BenchmarkTools
@benchmark func($psp)

BenchmarkTools.Trial: 10000 samples with 1 evaluation per sample.
 Range (min … max):  25.407 μs … 86.101 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     25.899 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   27.052 μs ±  4.157 μs  ┊ GC (mean ± σ):  0.00% ± 0.00%

  ▇█▅▃▂▁▁▁▁              ▁▁                                   ▂
  ██████████▇▇▇▇▇▆▆▇▆▆▆▆▇███▇▆▆▇▇▆▆▆▆▇▇▅▆▆▅▅▅▅▅▅▆▅▄▅▄▅▆▅▆▆▁▅▆ █
  25.4 μs      Histogram: log(frequency) by time      48.7 μs <

 Memory estimate: 256 bytes, allocs estimate: 2.

This page was generated using Literate.jl.