Reference

AbstractWignerRotation

Abstract type for representing Wigner rotations in three-body decays. Subtypes include TrivialWignerRotation and WignerRotation{N} for N = 0,2,3.

DalitzPlotPoint{I, S}

A structure representing a point in the Dalitz plot.

Fields

• σs::I: Mandelstam variables
• two_λs::S: Spin configuration

MandelstamTuple{T}

A named tuple representing Mandelstam variables (; σ1, σ2, σ3) for a three-body system.

MassTuple{T}

A named tuple representing the masses of a three-body system. Contains masses m₁, m₂, m₃ of the decay products and m₀ of the parent particle.

ParityTuple

A named tuple representing parities of a three-body system. Contains parities P₁, P₂, P₃ of the decay products and P₀ of the parent particle.

SpinParity

A structure representing spin and parity of a particle.

Fields

• two_j::Int: Twice the spin value
• p::Char: Parity ('+' or '-')

SpinTuple

A named tuple representing the spins of a three-body system. Contains twice the helicities of the particles (2h₁, 2h₂, 2h₃, 2h₀).

ThreeBodyDecay(descriptor)

Constructs a ThreeBodyDecay object using one argument, a descriptor. The descriptor is a list of pairs, names .=> zip(couplings, chains).

Examples

ThreeBodyDecay("K892" .=> zip([1.0, -1.0, 0.2im], [chain1, chain2, chain3]))

ThreeBodyDecay(; chains, couplings, names)

Constructs a ThreeBodyDecay object with the given parameters.

Arguments

• chains: An array of chains involved in the decay. The length of this array should match the lengths of couplings and names.
• couplings: An array of coupling constants for each chain in the decay. The length of this array should match the lengths of chains and names.
• names: An array of names for each chain, or names of resonances in the decay. The length of this array should match the lengths of chains and couplings.

Returns

• A ThreeBodyDecay object with the specified chains, couplings, and names.

Examples

ThreeBodyDecay(
chains=[chain1, chain2, chain3],
couplings=[1.0, -1.0, 0.2im],
names=["L1405", "L1405", "K892"])

ThreeBodySystem{T,K}

A structure representing a three-body system with masses and spins.

Fields

• ms::T: Masses of the system (MassTuple)
• two_js::K: Spins of the system (SpinTuple)

TrivialWignerRotation <: AbstractWignerRotation

Represents a trivial Wigner rotation (identity transformation). Used when the reference frame and system frame are the same.

Fields

• k::Int: Index of the spectator particle

VertexFunction{R<:Recoupling,F}

A struct that contains a recoupling and a form factor. There are two constructors:

VertexFunction(h::Recoupling)        # translates to a trivial form factor
VertexFunction(h::Recoupling, ff::F) # with a form factor
# amplitude(V::VertexFunction{<:Recoupling, yourF}, args...) where {yourF} needs to be defined

WignerRotation{N} <: AbstractWignerRotation

Represents a Wigner rotation in three-body decays. The type parameter N determines the kind of rotation:

• N=0: Rotation between total system frames
• N=2: Rotation involving two distinct indices
• N=3: Rotation involving three distinct indices

Fields

• k::Int: Index of the spectator particle
• ispositive::Bool: Direction of rotation
• iseven::Bool: Parity of the rotation

struct minusone
@x_str(s::String)

A utility type and macro for handling -1 powers in calculations. The macro creates a minusone instance that returns -1 when raised to odd powers and 1 when raised to even powers.

Example

x"-1"^3  # returns -1
x"-1"^2  # returns 1

Base.vcat(models::ThreeBodyDecay...)

Concatenates multiple ThreeBodyDecay objects into a single ThreeBodyDecay. Argument is variable number of ThreeBodyDecay objects.

An example

extended_model = vcat(model[2], model[2:3], model)

⊗(p1::Char, p2::Char)
⊗(jp1::SpinParity, jp2::SpinParity)

Compute the tensor product of two parities or spin-parity states.

Arguments

• p1, p2: Parity characters ('+' or '-')
• jp1, jp2: SpinParity objects

Returns

• For parities: The combined parity
• For spin-parity: Array of possible combined states

DecayChainLS(; k, Xlineshape, jp, Ps, tbs)

Constructs a decay chain with the smallest spin-orbit coupling.

Arguments

• k: Index of the spectator particle.
• Xlineshape: Lambda function for the lineshape of the resonance.
• jp: Spin-parity of the resonance (e.g., jp = "1/2-").
• Ps: Parities of the three-body system (e.g., Ps = ThreeBodyParities('+','+','+'; P0='+')).
• tbs: Three-body system structure.

Example

```julia DecayChainLS( k = 1, Xlineshape = x -> 1 / (x - 1im), jp = "1/2-", Ps = ThreeBodyParities('+', '+', '+'; P0 = '+'), tbs = ThreeBodySystem(1.0, 2.0, 3.0; m0 = 4.0) )

DecayChainsLS(; k, Xlineshape, jp, Ps, tbs)

Generate decay chains with all possible couplings based on the specified parameters.

Arguments

• k: index of spectator that specifies the decay chain.
• Xlineshape: Lambda function for the lineshape of the resonance.
• jp: Spin-parity quantum numbers of the resonance (e.g., jp = "1/2-").
• Ps: Parity quantum numbers of the three-body system (e.g., Ps = ThreeBodyParities('+','+','+'; P0='+')).
• tbs: Three-body-system structure that defines the involved particles and their properties.

Returns

• An array of DecayChain objects representing all possible couplings for the given decay configuration.

Notes

The function computes the possible coupling combinations (ls x LS). For each combination, a DecayChain object is created with the appropriate recoupling terms (Hij, HRk).

Example

```julia DecayChainsLS( k=3, Xlineshape=σ->1.0, jp="1/2-", Ps=ThreeBodyParities('+','+','+'; P0='+'), tbs=ThreeBodySystem( ThreeBodyMasses(1, 2, 3; m0=7.0), ThreeBodySpins(1, 2, 0; h0=3) ) )

Invariants(ms::MassTuple{T}; σ1, σ2, σ3)

Construct a tuple of Mandelstam invariants from given invariants and the mass tuple.

Arguments

• ms::MassTuple{T}: Masses of the system
• σ1, σ2, σ3: Optional invariants (two or three must be provided)

Returns

• MandelstamTuple{T}: The complete set of invariants

Throws

• ErrorException if not exactly two or three invariants are provided
• ErrorException if the invariants violate mass constraints

Kallen(x, y, z)

Calculate the Källén function λ(x,y,z) = x² + y² + z² - 2xy - 2yz - 2zx. This function appears frequently in relativistic kinematics calculations.

Arguments

• x, y, z: Real numbers or complex values

Returns

• The value of the Källén function

Examples

julia> Kallen(10.0, 1.0, 1.0)  # For a decay A -> B + C, calculate λ(mA², mB², mC²)
60.0

julia> Kallen(5.0, 1.0, 1.0)  # Another example with different masses
5.0

See also sqrtKallenFact, Kibble.

Kibble(σs, msq)

Calculate the Kibble function φ(σ₁,σ₂,σ₃) for three-body decays, which is defined in terms of Källén functions. The Kibble function determines the physical boundaries of the three-body phase space.

The function is calculated as: φ(σ₁,σ₂,σ₃) = λ(λ₁,λ₂,λ₃) where λᵢ = λ(M², mᵢ², σᵢ) are Källén functions of the total mass squared M², the i-th particle mass squared mᵢ², and the corresponding Mandelstam variable σᵢ.

Arguments

• σs: Tuple of Mandelstam variables (σ₁,σ₂,σ₃)
• msq: Tuple of squared masses (m₁²,m₂²,m₃²,M²)

Returns

• Value of the Kibble function

Examples

julia> msq = (1.0, 1.0, 1.0, 16.0);  # squared masses: m₁², m₂², m₃², M²

julia> σs = (2.0, 2.0, 2.0);         # Mandelstam variables

julia> Kibble(σs, msq)               # returns the Kibble function value
-77763.0

julia> Kibble(σs, msq) < 0           # physical configuration if negative
true

See also Kallen, isphysical.

ThreeBodyMasses(m1, m2, m3; m0)
ThreeBodyMasses(; m1, m2, m3, m0)

Construct a MassTuple for a three-body system.

Arguments

• m1, m2, m3: Masses of the decay products
• m0: Mass of the parent particle

Returns

• MassTuple{T}: A named tuple containing the masses

Throws

• ErrorException if m₀ is less than the sum of m₁, m₂, and m₃

ThreeBodyParities(P1, P2, P3; P0)

Construct a ParityTuple for a three-body system.

Arguments

• P1, P2, P3: Parities of the decay products
• P0: Parity of the parent particle

Returns

• ParityTuple: A named tuple containing the parities

ThreeBodySpinParities(jp1, jp2, jp3; jp0)

Construct spin and parity information for a three-body system using shortcuts like "x/2±", or jp"x/2±".

Arguments

• jp1, jp2, jp3: SpinParity objects for decay products
• jp0: SpinParity object for parent particle

Returns

• Tuple{SpinTuple,ParityTuple}: Tuple of spin structure and parity structure

ThreeBodySpins(two_h1_or_h1, two_h2_or_h2, two_h3_or_h3; h0=nothing, two_h0=nothing)

Construct a SpinTuple for a three-body system. Depending on the key arguments, h0 or two_h0, the position arguments can be either twice the helicity or helicity itself.

Arguments

• two_h1_or_h1: Twice the helicity or helicity of first particle
• two_h2_or_h2: Twice the helicity or helicity of second particle
• two_h3_or_h3: Twice the helicity or helicity of third particle
• h0: Key argument for helicity of parent particle (optional)
• two_h0: Key argument for twice the helicity of parent particle (optional)

Returns

• SpinTuple: A named tuple containing the spins

Throws

• ErrorException if neither h0 nor two_h0 is provided
• ErrorException if baryon number is not conserved

aligned_four_vectors(σs,ms; k::Int)

Computes the four-momenta of the three particles in the center of momentum frame aligning the k-th particle with the -z-axis.

Arguments

• σs: Tuple of mandelstam variables,
• ms: Tuple of masses in the order m1, m2, m3, m0.
• k: Index of the particle to be aligned with the -z-axis.

Returns

A tuple of three four-momenta in the form of (px, py, pz, E).

Examples

julia> ms = ThreeBodyMasses(1.0, 1.0, 1.0; m0=4.0);

julia> σs = x2σs([0.5, 0.5], ms; k=2);

julia> p1, p2, p3 = aligned_four_vectors(σs, ms; k=1);

julia> p1[3] ≈ -breakup_Rk(σs[1], ms; k=1)  # first particle aligned with -z axis
true

julia> all(p -> abs(p[2]) < 1e-10, (p1, p2, p3))  # all momenta in x-z plane
true

amplitude(dc::DecayChain, orientation_angles::PlaneOrientation, σs::MandelstamTuple; kw...)

Compute the amplitude for a given decay chain with orientation angles applied to the rotation of the final state.

Arguments

• dc::DecayChain: The decay chain object containing the system's configuration (e.g., spin, parity, etc.).
• orientation_angles::PlaneOrientation: Named tuple representing the plane orientation angles (α, cosβ, γ) for the final state.
• σs::MandelstamTuple: Tuple representing the Mandelstam variables that describe the kinematic invariants of the process.
• kw...: Additional keyword arguments to be passed to the underlying amplitude calculation (e.g refζs reference kinematics).

Returns

• A four dimensional amplitude array

Details

The function first computes the array of aligned amplitudes. Then it contract it with a Wigner D-matrix.

amplitude(dc::DecayChain, σs::MandelstamTuple, two_λs; refζs = (1, 2, 3, 1))

Compute the total amplitude for a given decay chain and kinematic configuration.

Arguments

• dc::DecayChain: The decay-chain object.
• σs::MandelstamTuple: Tuple representing Mandelstam variables that describe the kinematic invariants of the process.
• two_λs: A collection of helicity values for the particles involved in the decay.
• refζs: Reference ζ indices for alignment rotations (default is (1, 2, 3, 1)).

Returns

• A four-dimensional array of amplitude values.

amplitude(dc::DecayChain, σs::MandelstamTuple, two_λs; refζs = (1, 2, 3, 1))

Compute the total amplitude for a given decay chain, kinematic configuration, and all possible helicity values.

Arguments

• dc::DecayChain: The decay-chain object.
• σs::MandelstamTuple: Tuple representing Mandelstam variables that describe the kinematic invariants of the process.
• refζs: Reference ζ indices for alignment rotations (default is (1, 2, 3, 1)).

Returns

• A four-dimensional array of amplitude values.

border(ms::MassTuple{T}; Nx::Int=300) where T

Calculate the border of the Dalitz plot.

Arguments

• ms::MassTuple{T}: Masses of the system
• Nx::Int: Number of points to generate

Returns

• Vector{MandelstamTuple{T}}: Points on the border

breakup(m, m1, m2)
breakup_Rk(σk, ms; k)
breakup_ij(σk, ms; k)

Calculate the breakup momentum for a two-body system.

• breakup: General form for any two-body decay
• breakup_Rk: For decay of the total system into a subsystem (ij), and a particle k
• breakup_ij: For decay of the subsystem (ij) into particles i and j

Arguments

• m: Mass of the decaying system
• m1, m2: Masses of the decay products
• σk: Invariant mass squared of the isobar
• ms: MassTuple containing all masses
• k: Spectator index

Returns

• Break-up momentum magnitude

Examples

julia> ms = ThreeBodyMasses(1.0, 1.0, 1.0; m0=4.0);

julia> p = breakup_Rk(2.5^2, ms; k=1)
0.897587913521567

julia> q = breakup_ij(2.5^2, ms; k=1)
0.75

julia> breakup(4.0, 2.5, 1.0)  # direct use of breakup function
0.897587913521567

cosζ(wr::AbstractWignerRotation, σs, msq)

Calculate the cosine of the Wigner rotation angle ζ for a given kinematic configuration. Different methods are implemented for different types of Wigner rotations.

Arguments

• wr: A Wigner rotation object
• σs: Tuple of Mandelstam variables
• msq: Tuple of squared masses

Returns

• Cosine of the Wigner rotation angle

Example

ms = ThreeBodyMasses(1.0, 1.0, 1.0; m0=4.0)
σs = x2σs([0.5, 0.5], ms; k=1)
w = wr(2, 1, 1)
cosζ(w, σs, ms^2)  # Get rotation angle cosine

cosθij(k,σs,msq)

Isobar decay angle for the chain-k, i.e. an angle of between vectors pi and -pk in the (ij) rest frame.

Explicit forms: cosθ23, cosθ31, cosθ12.

plot(intensity, ms;
    xbins = 100,
    ybins = 100,
    iσx = 1, iσy = 2,
    xlims = (:auto, :auto),
    ylims = (:auto, :auto)
)
# also
plot(ms, intensity)
dalitzplot(ms, intensity)
dalitzplot(intensity, ms)

A plotting recipe for a Dalitz plot.

This recipe generates a Dalitz plot as a heatmap, visualizing the intensity of a function over a specified range of invariants.

Parameters:

• intensity::Function: A real function of the invariants, (m23², m31², m12²), returning a value at a given kinematic point.
• ms::MassTuple: A tuple representing the masses of the particles involved in the system, used for determination of the borders and function dispatch.

Keyword Arguments:

• iσx: Index of the first invariant to use for the x-axis, 1->m23², 2->m31², and 3->m12². Defaults to 1.
• iσy: Index of the second invariant to use for the y-axis. Defaults to 2.
• xbins: Number of points for the x-axis grid.
• ybins: Number of points for the y-axis grid.
• xlims: Limits for the x-axis in terms of the invariant range. Defaults to lims(iσx, ms) (calculated automatically). Can be a tuple with :auto for automatic limits (e.g., (:auto, 4.4)).
• ylims: Limits for the y-axis in terms of the invariant range. Defaults to lims(iσy, ms) (calculated automatically). Can be a tuple with :auto for automatic limits (e.g., (2.2, :auto)).

Output:

The recipe generates:

1. A grid of invariant values for σx and σy axes.
2. A 2D array of values, where each element is either:
   • NaN if the corresponding kinematic region is forbidden by the Kibble condition.
   • The intensity calculated from the provided intensity function.

Usage:

Just call a plot command,

plot(intensity, ms)
plot(ms, intensity)
dalitzplot(ms, intensity)
dalitzplot(intensity, ms)

dalitzprojection(function_of_σs, ms, integrator;
    k, # must be specified
    bins = 100,
    xlims = (:auto, :auto)
)
# also
dalitzprojection(ms, function_of_σs, integrator; k = 1)

A plotting recipe for Dalitz plot projections.

This recipe generates a 1D projection of the Dalitz plot functionofσs by integrating over one of the invariant mass coordinates, visualizing the functionofσs as a function of a single kinematic variable.

Parameters:

• ms::MassTuple: A tuple representing the masses of the particles involved in the system.
• function_of_σs::Function: A real function of the invariants, (m23², m31², m12²), returning a value at a given kinematic point.
• integrator: A numerical integration function (e.g., quadgk from QuadGK.jl). The integrator should accept (function, lower_limit, upper_limit) and return a tuple where the first element is the integral value.

Keyword Arguments:

• k: Index of the invariant to project onto (1, 2, or 3). The projection integrates over the other two coordinates.
• bins: Number of bins for the projection axis. Defaults to 100.
• xlims: Minimal and maximal values of the bin edges for the projection axis in terms of the invariant range. Defaults to lims(ms; k) (calculated automatically). Can be a tuple with :auto for automatic limits (e.g., (:auto, 4.4)).

Usage:

using QuadGK
dalitzprojection(ms, function_of_σs, quadgk; k = 1)
dalitzprojection(function_of_σs, ms, quadgk; k = 3, bins = 150)

Example:

using ThreeBodyDecays
using Plots
using QuadGK

ms = ThreeBodyMasses(0.938, 0.493, 0.0; m0 = 5.62)
function_of_σs = σs -> 1.0  # constant function_of_σs

# Project onto σ₁
dalitzprojection(ms, function_of_σs, quadgk; k = 1, bins = 100)

# Project onto σ₃ with custom limits
dalitzprojection(ms, function_of_σs, quadgk; k = 3, xlims = (20.0, 26.0))

ijk(k::Int)
ij_from_k(k::Int)

Return a tuple of indices (i,j,k) for a three-body system, where i,j are ordered cyclically from k. Used extensively in three-body decay calculations to determine the spectator particle (k) and the pair (i,j).

Arguments

• k::Int: Index of the spectator particle (1, 2, or 3)

Returns

• Tuple{Int,Int,Int}: A tuple of indices (i,j,k) where i,j are ordered cyclically from k

Example

ijk(1) # returns (2, 3, 1)
ijk(2) # returns (3, 1, 2)
ijk(3) # returns (1, 2, 3)

isphysical(σs::MandelstamTuple, ms::MassTuple) -> Bool
inphrange(σs::MandelstamTuple, ms::MassTuple) -> Bool
inrange(x, r) -> Bool

Check if a set of Mandelstam variables represents a physically valid configuration in the three-body phase space. A configuration is physical if:

1. The Kibble function is negative (φ(σ₁,σ₂,σ₃) < 0)
2. All Mandelstam variables are within their kinematically allowed ranges

Arguments

• σs::MandelstamTuple: Tuple of Mandelstam variables (σ₁,σ₂,σ₃)
• ms::MassTuple: Tuple of masses (m₁,m₂,m₃,M)
• x: Value to check (for inrange)
• r: Range tuple (min,max) to check against (for inrange)

Returns

• Bool: true if the configuration is physical, false otherwise

Examples

julia> ms = ThreeBodyMasses(1.0, 1.0, 1.0; m0=4.0);

julia> σs = (6.25, 6.25, 6.5);

julia> isphysical(σs, ms)  # checks if this configuration is physically possible
true

julia> σs_unphysical = (20.0, 20.0, 20.0);  # values outside allowed ranges

julia> isphysical(σs_unphysical, ms)  # this configuration is not physical
false

See also Kibble, lims.

ispositive(wr::AbstractWignerRotation)

Determine if a Wigner rotation is in the positive direction.

Arguments

• wr: A Wigner rotation object

Returns

• true for positive rotations, false otherwise

jls_coupling(two_j1, two_λ1, two_j2, two_λ2, two_j, two_l, two_s)

Calculate the LS-coupling coefficient for two particles coupling to total angular momentum j. Uses Clebsch-Gordan coefficients to compute the coupling of orbital angular momentum l and total spin s to total j, with specified helicities.

Arguments

• two_j1, two_j2: Twice the spins of the particles
• two_λ1, two_λ2: Twice the helicities
• two_j: Twice the total angular momentum
• two_l: Twice the orbital angular momentum
• two_s: Twice the total spin

Returns

• The coupling coefficient value

Example

# For j₁=1/2, λ₁=1/2, j₂=1/2, λ₂=-1/2, j=1, l=1, s=1
julia> ThreeBodyDecays.jls_coupling(1, 1, 1, -1, 2, 2, 2)
-0.7071067811865477

letterL(l::Int)
letterL(l::String)

Convert an angular momentum quantum number to its spectroscopic notation. For l = 0,1,2,3,4,5 returns S,P,D,F,G,H respectively. For l ≥ 6 returns the first character of the string representation.

Arguments

• l: Angular momentum quantum number (as Int or String)

Returns

• Char: The spectroscopic notation for the given angular momentum

Example

letterL(0)  # returns 'S'
letterL(1)  # returns 'P'
letterL("2")  # returns 'D'

lims(ms::MassTuple; k::Int)
lims(k::Int, ms::MassTuple)
lims1(ms::MassTuple)
lims2(ms::MassTuple)
lims3(ms::MassTuple)

Calculate the kinematic limits (boundaries) for the Mandelstam variables σᵢ in a three-body decay. For each pair of particles (i,j), the invariant mass squared σₖ = (pᵢ + pⱼ)² must lie within physical limits:

• Lower bound: (mᵢ + mⱼ)² (threshold for producing particles i and j)
• Upper bound: (M - mₖ)² (maximum energy available when particle k recoils)

Arguments

• ms: Tuple of masses (m₁,m₂,m₃,M)
• k: Index specifying which pair of particles (1,2,3)

Returns

• Tuple (min,max) of the allowed range for σₖ

Example

ms = ThreeBodyMasses(1.0, 1.0, 1.0; m0=4.0)
lims1(ms)  # limits for σ₁ = (p₂ + p₃)²
lims2(ms)  # limits for σ₂ = (p₃ + p₁)²
lims3(ms)  # limits for σ₃ = (p₁ + p₂)²

See also isphysical, Kibble.

Phase for wigner d-functions for clockwise rotations

phase_space_integrand(function_σs, ms; k::Int)

Calculate the phase space integrand for a given function function_σs, and mass tuple ms. The key argument k specifies the mapping: σk->[0,1], zk->[0,1]. It returns an integrand function of x, x ∈ [0,1]x[0,1] domain to pass to a numerical integrator.

Arguments

• function_σs: A function that takes a MandelstamTuple and returns a scalar.
• ms: A scalar representing the mass.
• k: An integer representing the mapping index.

Usage

integrand = phase_space_integrand(function_σs, ms; k)

See also

• x2σs
• projection_integrand

polardalitz2invariants(θ, expansion_point)

For given polar angle θ, returns an (σ1,σ2,σ3) Tuple of polynomials of radius r(θ) around the expansion point. The polynomial works as a function of the r coordinate.

Arguments

• θ: Polar angle
• expansion_point: Tuple of expansion point coordinates

Returns

• Tuple{Polynomial,Polynomial,Polynomial}: Polynomials for each invariant

possible_l_s_L_S(jp::SpinParity, two_js::SpinTuple, Ps::ParityTuple; k)

Convert the output of possible_lsLS to a more readable format with half-integer values.

Arguments

Same as possible_lsLS

Returns

• Vector of named tuples containing (l,s,L,S) as strings representing half-integer values

Example

ms = ThreeBodyMasses(1.0, 1.0, 1.0; m0=4.0)
two_js = ThreeBodySpins(1, 1, 0; h0=2)
Ps = ThreeBodyParities('+','+','+'; P0='-')
qn = possible_l_s_L_S(jp"1-", two_js, Ps; k=1)
# Each element has fields like qn[1].l == "0", qn[1].s == "1/2", etc.

possible_ls(jp1::SpinParity, jp2::SpinParity; jp::SpinParity)
possible_ls_ij(jp::SpinParity, two_js::SpinTuple, Ps::ParityTuple; k::Int)
possible_ls_Rk(jp::SpinParity, two_js::SpinTuple, Ps::ParityTuple; k::Int)

Calculate possible orbital angular momentum (l) and spin (s) combinations for a two-body system.

• possible_ls: For any two particles with given spin-parities
• possible_ls_ij: For the (i,j) pair in a three-body decay
• possible_ls_Rk: For the resonance-spectator system

Arguments

• jp1, jp2: Spin-parity quantum numbers of the particles
• jp: Total spin-parity of the system
• two_js: Spins of all particles in the system
• Ps: Parities of all particles
• k: Spectator index

Returns

• Vector of tuples (l,s) containing allowed combinations

Example

# For a decay 1⁻ → 1/2⁺ + 1/2⁺
ls = possible_ls(jp"1/2+", jp"1/2+"; jp=jp"1-")
# For a three-body system
ms = ThreeBodyMasses(1.0, 1.0, 1.0; m0=4.0)
two_js = ThreeBodySpins(1, 1, 0; h0=2)
Ps = ThreeBodyParities('+','+','+'; P0='-')
ls_ij = possible_ls_ij(jp"1-", two_js, Ps; k=1)

possible_lsLS(jp::SpinParity, two_js::SpinTuple, Ps::ParityTuple; k::Int)

Calculate all possible combinations of orbital angular momenta and spins for both the isobar decay (l,s) and the total system decay (L,S).

Arguments

• jp: Spin-parity of the resonance
• two_js: Spins of all particles
• Ps: Parities of all particles
• k: Spectator index

Returns

• Vector of named tuples containing with all possible (twols, twoLS) combinations

Example

ms = ThreeBodyMasses(1.0, 1.0, 1.0; m0=4.0)
two_js = ThreeBodySpins(1, 1, 0; h0=2)
Ps = ThreeBodyParities('+','+','+'; P0='-')
lsLS = possible_lsLS(jp"1-", two_js, Ps; k=1)

project_cosθij_intergand(fs, ms, z; k)

Calculate the integrand for projecting onto cos(θij) for a given function fs, mass tuple ms, and cosine value z, with k specifying the spectator. Returns an integrand function of x ∈ [0,1] to pass to a numerical integrator.

Arguments

• fs: A function that takes a MandelstamTuple and returns a scalar.
• ms: A MassTuple containing the masses.
• z: The cosine value to project onto.
• k: The spectator index.

Returns

• An integrand function for numerical integration.

Example

let Nb = 100
    [sum(range(-1, 1, Nb)) do z
        integrand = project_cosθij_intergand(I, ms, z; k)
        quadgk(integrand, 0, 1)[1] * 2/Nb
    end for k in 1:3]
end

projectionintegrand(functionσs, ms, σk; k)

Calculate the projection integrand for a given function function_σs, mass tuple ms, and Mandelstam variable σk, with k specified by a keyword argument. It returns an integrand function of x, x ∈ [0,1] to pass to a numerical integrator.

Arguments

• function_σs: A function that takes a MandelstamTuple and returns a scalar.
• ms: A scalar representing the mass.
• σk: A scalar representing the Mandelstam variable.
• k: A scalar representing the momentum transfer (optional).

Usage

plot(4.2, 4.6) do e1
I = Base.Fix1(unpolarized_intensity, model)
integrand = projection_integrand(I, masses(model), e1^2; k = 3)
e1 * quadgk(integrand, 0, 1)[1]
end

randomPoint(ms::MassTuple)
randomPoint(two_js::SpinTuple)
randomPoint(tbs::ThreeBodySystem)

Generate a random point in the phase space.

Arguments

• ms::MassTuple: Masses of the system
• two_js::SpinTuple: Spins of the system
• tbs::ThreeBodySystem: Complete three-body system

Returns

• For masses: Random Mandelstam variables
• For spins: Random spin configuration
• For system: Random DalitzPlotPoint

shift_by_half(v)

Shift a vector of values by half the step size between elements. Used in plotting and binning operations to center values between grid points.

Arguments

• v: Vector of values

Returns

• Vector of values shifted by half the step size

Example

v = [1, 2, 3, 4]
shift_by_half(v)  # returns [1.5, 2.5, 3.5]

sqrtKallenFact(a, b, c)

Calculate the square root of the Källén function λ(a,b,c) in factorized form. The Källén function is defined as λ(a,b,c) = a² + b² + c² - 2ab - 2bc - 2ca. This function returns √λ(a,b,c) in a factorized form to improve numerical stability: √(a-(b+c)) √(a-(b-c)) √(a+(b+c)) √(a+(b-c)).

Arguments

• a, b, c: Input values for the Källén function

Returns

• Square root of the Källén function in factorized form

Examples

julia> sqrtKallenFact(10.0, 1.0, 1.0)  # √λ(10,1,1) = √9600
97.97958971132714

julia> sqrtKallenFact(4.0, 1.0, 1.0)  # √λ(4,1,1) = √192
13.856406460551018

See also Kallen, Kibble.

str2jp(pin::AbstractString)

Convert a string representation of spin-parity to a SpinParity object.

Arguments

• pin::AbstractString: String in format "x±" or "x/2±"

Returns

• SpinParity: The corresponding spin-parity object

Throws

• ErrorException if the string format is invalid

unpolarized_intensity(model::ThreeBodyDecay, σs; kw...)

Computes squared amplitude summed over spin projections.

wr(system_a, reference_b, particle_c=0)

Create a WignerRotation object for transforming between different reference frames in a three-body system.

Arguments

• system_a: Index of the isobar system being considered (1,2,3)
• reference_b: Index of the reference system (1,2,3)
• particle_c: Index of the particle (0 for parent particle, 1,2,3 for daughters)

Returns

• A WignerRotation object of the appropriate type

Example

# Rotation from system 2 to system 1 for particle 1 => zeta_21_for1
w = wr(2, 1, 1)
# Rotation between total system frames => zeta_12_for0
w0 = wr(1, 2, 0)

x2σs(x, ms::MassTuple; k::Int = last(findmin(Tuple(ms))))

Maps a pair of variables x to a physical set of mandelstam invariants σs using a linear transformation. x[1] is mapped to σ[k], and x[2] is mapped into cosθij, that is used to compute σ[i] and σ[j]. Uniform distribution of x does not lead to the phase space distribution of σs.

Arguments

• x : a pair of numbers
• ms : masses of the system as a MassTuple, see ThreeBodyMasses.
• k : the index for which the variable is not generated. By default, the function picks the coordinates where the Dalitz plot has the closest shape to the squared fitting box.

Returns

• an instance of MandelstamTuple with the squared masses.

Example

The phase space sample with 100 points can be generated as follows:

σs0 = x2σs([0.3, 0.2], ms)

See also y2σs.

y2σs(y, ms::MassTuple; k::Int = last(findmin(Tuple(ms))))

Maps a pair of variables to the plane of squared masses using a linear transformation, by shifting and scaling y[1] to σ[i], and y[2] to σ[j]. The indices of variables are determined by the argument k. The mapping does not guarantee that values are physical, however, physical values of σs are phase-space distributed for uniform y.

Arguments

• y : a pair of numbers
• ms : masses of the system as a MassTuple, see ThreeBodyMasses.
• k : the index for which the variable is not generated. By default, the function picks the coordinates

where the Dalitz plot has the closest shape to the squared fitting box.

Returns

• an instance of MandelstamTuple with the squared masses.

Example

The phase space sample with 100 points can be generated as follows:

data = let
    N = 100
    # map random variables to dalitz
    _data = mapslices(rand(N, 2); dims=2) do xy
        y2σs(xy, ms)
    end[:, 1]
    # select physical
    filter!(_data) do σs
        isphysical(σs, ms)
    end
    _data
end

See also x2σs.

σiofk(k,z,σj,msq)

Computes invariant σi = (p0 - pi)² from the scattering angle z=cosθij in the rest from of (i,j), given the mass of the system m(i,j)² = σk

Explicit forms: σ3of2, σ1of3, σ2of1.

See also σjofk

σjofk(z,σi,msq; k::Int)

Computes invariant σj = (p0-pj)² from the scattering angle z=cosθij in the rest from of (i,j), given the mass of the system m(i,j)² = σk

Explicit forms: σ3of1, σ1of2, σ2of3.

See also σiofk