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 some of the packages of the QEDjl-project for base functionality and the ScatteringProcess type.

using QEDcore
using QEDprocesses

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 = ScatteringProcess(
    (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
)

generic QED process
    incoming: electron (spin up), photon (x-polarized)
    outgoing: electron (spin up), electron (spin up), electron (spin up), positron (spin up), positron (spin up)

The feynman_diagrams function returns an iterator for all possible Feynman diagrams for this scattering process. With its length overload, we can check how many diagrams there are.

length(feynman_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: 619, QEDFeynmanDiagrams.ComputeTask_PropagatePairs: 42, QEDFeynmanDiagrams.ComputeTask_CollectTriples: 1, 
         QEDFeynmanDiagrams.ComputeTask_Triple: 81, QEDFeynmanDiagrams.ComputeTask_BaseState: 7, QEDFeynmanDiagrams.ComputeTask_CollectPairs: 42, 
         QEDFeynmanDiagrams.ComputeTask_Propagator: 42, QEDFeynmanDiagrams.ComputeTask_Pair: 69, QEDFeynmanDiagrams.ComputeTask_SpinPolCumulation: 1, 
         ComputableDAGs.DataTask: 334
  Edges: 950
  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__)

With the DAG, the process, and RuntimeGeneratedFunctions initalized, we can now generate the actual computable function:

func = get_compute_function(dag, proc, cpu_st(), @__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,
    PerturbativeQED(),
    PhasespaceDefinition(SphericalCoordinateSystem(), ElectronRestFrame()),
    tuple((rand(SFourMomentum) for _ in 1:number_incoming_particles(proc))...),
    tuple((rand(SFourMomentum) for _ in 1:number_outgoing_particles(proc))...),
)

PhaseSpacePoint:
    process: generic QED process "ek -> eeepp"
    model: perturbative QED
    phasespace definition: spherical coordinates in electron rest frame
    incoming particles:
     -> incoming electron: [0.039288920956634055, 0.1793051376783037, 0.040877702045857056, 0.8305673928496544]
     -> incoming photon: [0.0967942367870559, 0.4433805983076561, 0.04432540734441481, 0.8845963980802497]
    outgoing particles:
     -> outgoing electron: [0.8609764342654231, 0.3371047141795719, 0.2170733705691935, 0.935870102264508]
     -> outgoing electron: [0.42783534512226973, 0.3799614103866835, 0.18329084462688572, 0.3780260818471026]
     -> outgoing electron: [0.33378974015700114, 0.5651096608238709, 0.7021171552338386, 0.5473296494075333]
     -> outgoing positron: [0.48896102867270874, 0.9286437808226403, 0.4844472978351566, 0.3967672811366705]
     -> outgoing positron: [0.2500619643269336, 0.5349195902602311, 0.6356564476530698, 0.07552437474520801]

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

func(psp)

0.14374032150177865

If we want, we can benchmark the execution speed too:

using BenchmarkTools
@benchmark func($psp)

BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  13.235 μs …  32.421 μs  ┊ GC (min … max): 0.00% … 0.00%
 Time  (median):     13.396 μs               ┊ GC (median):    0.00%
 Time  (mean ± σ):   13.477 μs ± 748.855 ns  ┊ GC (mean ± σ):  0.00% ± 0.00%

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

 Memory estimate: 256 bytes, allocs estimate: 2.

This page was generated using Literate.jl.