Refactoring of differential cross section functionality

Project.toml

@@ -4,10 +4,8 @@ authors = ["Uwe Hernandez Acosta <[email protected]>", "Simeon Ehrig", "Klaus
 version = "0.1.0"
 
 [deps]
-DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 QEDbase = "10e22c08-3ccb-4172-bfcf-7d7aa3d04d93"
 
 [compat]
-DocStringExtensions = "0.9"
 QEDbase = "0.1"
 julia = "1.6"

src/QEDprocesses.jl

@@ -7,7 +7,14 @@ export AbstractModelDefinition, fundamental_interaction_type
 export AbstractProcessDefinition, incoming_particles, outgoing_particles
 export in_phasespace_dimension, out_phasespace_dimension
 export number_incoming_particles, number_outgoing_particles
-export differential_cross_section, total_cross_section
+
+# probabilities
+export differential_probability, unsafe_differential_probability
+export total_probability
+
+# differential cross section
+export differential_cross_section, unsafe_differential_cross_section
+export total_cross_section
 
 # Abstract setup interface
 export AbstractComputationSetup, InvalidInputError, compute
@@ -16,12 +23,19 @@ export AbstractProcessSetup, scattering_process, physical_model
 # propagator
 export propagator
 
-using DocStringExtensions
+# phase space 
+export AbstractCoordinateSystem, SphericalCoordinateSystem
+export AbstractFrameOfReference, CenterOfMomentumFrame, ElectronRestFrame
+export AbstractPhasespaceDefinition, PhasespaceDefinition
+
 using QEDbase
 
 include("utils.jl")
 include("interfaces/model_interface.jl")
 include("interfaces/process_interface.jl")
 include("interfaces/setup_interface.jl")
 include("propagators.jl")
+include("phase_spaces.jl")
+include("probabilities.jl")
+include("cross_sections.jl")
 end

src/cross_sections.jl

@@ -0,0 +1,304 @@
+########################
+# differential cross sections and probabilities.
+#
+# This file contains default implementations for differential and total cross
+# sections based on the scattering process interface
+########################
+
+############
+#
+# differential cross sections
+#
+############
+
+# differential cross sections without energy momentum conservation check
+# single in phase space point/ single out phase space point
+function _unsafe_differential_cross_section(
+    proc::AbstractProcessDefinition,
+    model::AbstractModelDefinition,
+    in_phase_space_def::AbstractPhasespaceDefinition,
+    in_phase_space::AbstractVector{T},
+    out_phase_space_def::AbstractPhasespaceDefinition,
+    out_phase_space::AbstractVector{T},
+) where {T<:QEDbase.AbstractFourMomentum}
+    I = 1 / (4 * _incident_flux(proc, model, in_phase_space))
+
+    return I * _unsafe_differential_probability(
+        proc,
+        model,
+        in_phase_space_def,
+        in_phase_space,
+        out_phase_space_def,
+        out_phase_space,
+    )
+end
+
+# differential cross sections without energy momentum conservation check
+# single in phase space point/ several out phase space points
+function _unsafe_differential_cross_section(
+    proc::AbstractProcessDefinition,
+    model::AbstractModelDefinition,
+    in_phase_space_def::AbstractPhasespaceDefinition,
+    in_phase_space::AbstractVector{T},
+    out_phase_space_def::AbstractPhasespaceDefinition,
+    out_phase_space::AbstractMatrix{T},
+) where {T<:QEDbase.AbstractFourMomentum}
+    res = Vector{eltype(T)}(undef, size(out_phase_space, 2))
+    for i in 1:size(out_phase_space, 2)
+        res[i] = _unsafe_differential_cross_section(
+            proc,
+            model,
+            in_phase_space_def,
+            in_phase_space,
+            out_phase_space_def,
+            view(out_phase_space, :, i),
+        )
+    end
+    return res
+end
+
+# differential cross sections without energy momentum conservation check
+# several in phase space points/ one or several out phase space points
+function _unsafe_differential_cross_section(
+    proc::AbstractProcessDefinition,
+    model::AbstractModelDefinition,
+    in_phase_space_def::AbstractPhasespaceDefinition,
+    in_phase_space::AbstractMatrix{T},
+    out_phase_space_def::AbstractPhasespaceDefinition,
+    out_phase_space::AbstractVecOrMat{T},
+) where {T<:QEDbase.AbstractFourMomentum}
+    res = Matrix{eltype(T)}(undef, size(in_phase_space, 2), size(out_phase_space, 2))
+    for i in 1:size(in_phase_space, 2)
+        res[i, :] .= _unsafe_differential_cross_section(
+            proc,
+            model,
+            in_phase_space_def,
+            view(in_phase_space, :, i),
+            out_phase_space_def,
+            out_phase_space,
+        )
+    end
+    return res
+end
+
+"""
+
+    function unsafe_differential_cross_section(
+        proc::AbstractProcessDefinition,
+        model::AbstractModelDefinition,
+        in_phase_space_def::AbstractPhasespaceDefinition,
+        in_phase_space::AbstractVecOrMat{T},
+        out_phase_space_def::AbstractPhasespaceDefinition,
+        out_phase_space::AbstractVecOrMat{T},
+    ) where {T<:QEDbase.AbstractFourMomentum}
+
+Return the differential cross section without checking if the given phase space(s) are physical.
+"""
+function unsafe_differential_cross_section(
+    proc::AbstractProcessDefinition,
+    model::AbstractModelDefinition,
+    in_phase_space_def::AbstractPhasespaceDefinition,
+    in_phase_space::AbstractVecOrMat{T},
+    out_phase_space_def::AbstractPhasespaceDefinition,
+    out_phase_space::AbstractVecOrMat{T},
+) where {T<:QEDbase.AbstractFourMomentum}
+    size(in_phase_space, 1) == number_incoming_particles(proc) || throw(
+        DimensionMismatch(
+            "The number of incoming particles <$(number_incoming_particles(proc))> is inconsistent with input size <$(size(in_phase_space,1))>",
+        ),
+    )
+
+    size(out_phase_space, 1) == number_outgoing_particles(proc) || throw(
+        DimensionMismatch(
+            "The number of outgoing particles <$(number_outgoing_particles(proc))> is inconsistent with input size <$(size(out_phase_space,1))>",
+        ),
+    )
+
+    return _unsafe_differential_cross_section(
+        proc,
+        model,
+        in_phase_space_def,
+        in_phase_space,
+        out_phase_space_def,
+        out_phase_space,
+    )
+end
+
+# differential cross sections with energy momentum conservation check
+# single in phase space point/ single out phase space point
+function _differential_cross_section(
+    proc::AbstractProcessDefinition,
+    model::AbstractModelDefinition,
+    in_phase_space_def::AbstractPhasespaceDefinition,
+    in_phase_space::AbstractVector{T},
+    out_phase_space_def::AbstractPhasespaceDefinition,
+    out_phase_space::AbstractVector{T},
+) where {T<:QEDbase.AbstractFourMomentum}
+    if !_is_in_phasespace(
+        proc,
+        model,
+        in_phase_space_def,
+        in_phase_space,
+        out_phase_space_def,
+        out_phase_space,
+    )
+        return zero(eltype(T))
+    end
+
+    return _unsafe_differential_cross_section(
+        proc,
+        model,
+        in_phase_space_def,
+        in_phase_space,
+        out_phase_space_def,
+        out_phase_space,
+    )
+end
+
+# differential cross sections with energy momentum conservation check
+# single in phase space point/ several out phase space points
+function _differential_cross_section(
+    proc::AbstractProcessDefinition,
+    model::AbstractModelDefinition,
+    in_phase_space_def::AbstractPhasespaceDefinition,
+    in_phase_space::AbstractVector{T},
+    out_phase_space_def::AbstractPhasespaceDefinition,
+    out_phase_space::AbstractMatrix{T},
+) where {T<:QEDbase.AbstractFourMomentum}
+    res = Vector{eltype(T)}(undef, size(out_phase_space, 2))
+    for i in 1:size(out_phase_space, 2)
+        res[i] = _differential_cross_section(
+            proc,
+            model,
+            in_phase_space_def,
+            in_phase_space,
+            out_phase_space_def,
+            view(out_phase_space, :, i),
+        )
+    end
+    return res
+end
+
+# differential cross sections with energy momentum conservation check
+# several in phase space points/ one or several out phase space points
+function _differential_cross_section(
+    proc::AbstractProcessDefinition,
+    model::AbstractModelDefinition,
+    in_phase_space_def::AbstractPhasespaceDefinition,
+    in_phase_space::AbstractMatrix{T},
+    out_phase_space_def::AbstractPhasespaceDefinition,
+    out_phase_space::AbstractVecOrMat{T},
+) where {T<:QEDbase.AbstractFourMomentum}
+    res = Matrix{eltype(T)}(undef, size(in_phase_space, 2), size(out_phase_space, 2))
+    for i in 1:size(in_phase_space, 2)
+        res[i, :] .= _differential_cross_section(
+            proc,
+            model,
+            in_phase_space_def,
+            view(in_phase_space, :, i),
+            out_phase_space_def,
+            out_phase_space,
+        )
+    end
+    return res
+end
+
+"""
+    differential_cross_section(
+        proc::AbstractProcessDefinition,
+        model::AbstractModelDefinition,
+        in_phase_space_def::AbstractPhasespaceDefinition,
+        in_phase_space::AbstractVecOrMat{T},
+        out_phase_space_def::AbstractPhasespaceDefinition,
+        out_phase_space::AbstractVecOrMat{T},
+    ) where {T<:QEDbase.AbstractFourMomentum}
+
+If the given phase spaces are physical, return differential cross section. Zero otherwise
+
+"""
+function differential_cross_section(
+    proc::AbstractProcessDefinition,
+    model::AbstractModelDefinition,
+    in_phase_space_def::AbstractPhasespaceDefinition,
+    in_phase_space::AbstractVecOrMat{T},
+    out_phase_space_def::AbstractPhasespaceDefinition,
+    out_phase_space::AbstractVecOrMat{T},
+) where {T<:QEDbase.AbstractFourMomentum}
+    size(in_phase_space, 1) == number_incoming_particles(proc) || throw(
+        DimensionMismatch(
+            "The number of incoming particles <$(number_incoming_particles(proc))> is inconsistent with input size <$(size(in_phase_space,1))>",
+        ),
+    )
+
+    size(out_phase_space, 1) == number_outgoing_particles(proc) || throw(
+        DimensionMismatch(
+            "The number of outgoing particles <$(number_outgoing_particles(proc))> is inconsistent with input size <$(size(out_phase_space,1))>",
+        ),
+    )
+
+    return _differential_cross_section(
+        proc,
+        model,
+        in_phase_space_def,
+        in_phase_space,
+        out_phase_space_def,
+        out_phase_space,
+    )
+end
+
+############
+# Total cross sections
+############
+
+# total cross section on single phase space point
+function _total_cross_section(
+    proc::AbstractProcessDefinition,
+    model::AbstractModelDefinition,
+    in_phase_space_def::AbstractPhasespaceDefinition,
+    in_phase_space::AbstractVector{T},
+) where {T<:QEDbase.AbstractFourMomentum}
+    I = 1 / (4 * _incident_flux(proc, model, in_phase_space))
+
+    return I * _total_probability(proc, model, in_phase_space_def, in_phase_space)
+end
+
+# total cross section on several phase space points
+function _total_cross_section(
+    proc::AbstractProcessDefinition,
+    model::AbstractModelDefinition,
+    in_phase_space_def::AbstractPhasespaceDefinition,
+    in_phase_space::AbstractMatrix{T},
+) where {T<:QEDbase.AbstractFourMomentum}
+    res = Vector{eltype(T)}(undef, size(in_phase_space, 2))
+    for i in 1:size(in_phase_space, 2)
+        res[i] = _total_cross_section(
+            proc, model, in_phase_space_def, view(in_phase_space, :, i)
+        )
+    end
+    return res
+end
+
+"""
+    total_cross_section(
+        proc::AbstractProcessDefinition,
+        model::AbstractModelDefinition,
+        in_phase_space_def::AbstractPhasespaceDefinition,
+        in_phase_space::AbstractVecOrMat{T},
+    ) where {T<:QEDbase.AbstractFourMomentum}
+
+Return the total cross section for a given combination of scattering process and compute model.
+"""
+function total_cross_section(
+    proc::AbstractProcessDefinition,
+    model::AbstractModelDefinition,
+    in_phase_space_def::AbstractPhasespaceDefinition,
+    in_phase_space::AbstractVecOrMat{T},
+) where {T<:QEDbase.AbstractFourMomentum}
+    size(in_phase_space, 1) == number_incoming_particles(proc) || throw(
+        DimensionMismatch(
+            "The number of incoming particles <$(number_incoming_particles(proc))> is inconsistent with input size <$(size(in_phase_space,1))>",
+        ),
+    )
+
+    return _total_cross_section(proc, model, in_phase_space_def, in_phase_space)
+end