Welcome to AmpForm-DPD!#
This Python package is a (temporary) extension of AmpForm and provides a symbolic implementation of Dalitz-plot decomposition (10.1103/PhysRevD.101.034033) with SymPy. It has been extracted from the ComPWA/polarimetry repository, which is not yet public.
Installation#
The fastest way of installing this package is through PyPI (ampform-dpd):
python3 -m pip install ampform-dpd
This installs the latest version that you can find on the stable branch. The latest version on the main branch
can be installed as follows:
python3 -m pip install git+https://github.com/ComPWA/ampform@main
You can substitute main in the above command with any of the tags listed under the releases or any other commit hashes. However, we highly recommend using the more dynamic, ‘editable installation’ instead. This goes as follows:
Get the source code (see the Pro Git Book):
git clone https://github.com/ComPWA/ampform-dpd cd ampform-dpd
[Recommended] Create a virtual environment (see here or the tip below).
Install the project in ‘editable installation’ with additional dependencies for the developer:
python3 -m pip install -e . --group dev
That’s all! Have a look at the Welcome to AmpForm-DPD! page to try out the package, and see Help developing for tips on how to work with this ‘editable’ developer setup!
Tip
It’s easiest to install the project in a virtual environment with uv. In that case, to install in editable mode, just run:
uv sync
source .venv/bin/activate
This way of installing is also safer, because it pins all dependencies.
Physics#
Dalitz-plot decomposition allows us to separate variables that affect the angular distribution from variables that describe the dynamics. It allows rewriting a transition amplitude \(T\) as
Here, \(\Lambda\) and \(\nu\) indicate the allowed spin projections of the initial state, \(\{\lambda\}\) are the allowed spin projections of the final state (e.g. \(\{\lambda\}=\lambda_1,\lambda_3,\lambda_3\) for a three-body decay). The Euler angles \(\alpha,\beta,\gamma\) are obtained by choosing a specific aligned center-of-momentum frame (“aligned CM”), see Fig. 2 in Ref [1], which gives us an “aligned” transition amplitude \(O^\nu_{\{\lambda\}}\) that only depends on dynamic variables \(\{\sigma\}\) (in the case of a three-body decay, the three Mandelstam variables \(\sigma_1,\sigma_2,\sigma_3\)).
These aligned transition amplitudes are then combined into an observable differential cross section (intensity distribution), using a spin density matrix \(\rho_{_{\Lambda,\Lambda'}}\) for the spin projections \(\Lambda\) of the initial state,
Given the right alignment, the aligned transition amplitude can be written as
Notice the general structure:
Summations: The outer sum is taken over the three decay chain combinations \((ij)k \in \left\{(23)1, (31)2, (12)3\right\}\). Next, we sum over the spin magnitudes \(s\) of all resonances[1], the corresponding allowed helicities \(\tau\), and allowed spin projections \(\{\lambda'\}\) of the final state.
Dynamics: The function \(X_s\) only depends on a single Mandelstam variable and carries all the dynamical information about the decay chain. Typically, these are your \(K\)-matrix or Breit-Wigner lineshape functions.
Isobars: There is a Wigner \(d\)-function and a helicity coupling \(H\) for each isobar in the three-body decay chain: the \(0\to(ij),k\) production node and the \((ij)\to i,j\) decay node. The argument of these Wigner \(d\)-functions are the polar angles. The factors \(\eta_J=\sqrt{2S+1}\) and \(\eta_s=\sqrt{2s+1}\) are normalization factors. The phase \((-1)^{j-\lambda}\) is added to both helicity couplings to convert to the Jacob-Wick particle-2 convention.[2] The convention treats the first and the second particle unequally, however, it enables the simple relation of the helicity couplings to the \(LS\) couplings explained below.
Wigner rotations: The last three Wigner \(d\)-functions represent Wigner rotations that appear when rotating the boosted frames of the production and decay isobar amplitudes back to the space-fixed CM frame.
If \(k=1\), we have \(\hat\theta_{k(1)}=0\), so the Wigner \(d\) function for the production isobar reduces to a Kronecker delta, \(d^J_{\nu,\tau-\lambda'_k}\!\left(\hat\theta_{k(1)}\right) = \delta_{\nu,\tau-\lambda'_k}\).
Equation (1) is written in terms of helicity couplings, but can be rewritten in terms of \(LS\) couplings, using
The dynamics function is dependent on the \(LS\) values and we write \(X_s^{LS;l's'}\) instead of \(X_s\).
Usage#
The core class of this package is the DalitzPlotDecompositionBuilder. It requires a ThreeBodyDecay object to formulate amplitude models. You can create such a decay object by hand, but it is easier to generate it with QRules. Here is an example for the decay \(J/\psi \to p \bar{\varSigma}^- \bar{K}^0\). In the last step, we “normalize” the state IDs of the incoming and outgoing states so that they match the index naming scheme of the Dalitz-Plot Decomposition paper.
import qrules
from ampform_dpd.adapter.qrules import normalize_state_ids
reaction = qrules.generate_transitions(
initial_state="J/psi(1S)",
final_state=["p", "Sigma~-", "K~0"],
allowed_interaction_types="strong",
mass_conservation_factor=0,
)
reaction = normalize_state_ids(reaction)
Since AmpForm-DPD formulates amplitudes for three-body decays only, it defines its own topology class for encoding the decay chains. You can convert the QRules ReactionInfo object to such a ThreeBodyDecay as follows:
from ampform_dpd.adapter.qrules import to_three_body_decay
decay = to_three_body_decay(reaction.transitions)
The ampform_dpd.io module offers some handy functions for rendering the decay chains in different format.
We an now use the decay object to formulate amplitude models. Notice that we do not parametrize the dynamics for each of the chains yet, we only have a simple expression for the unpolarized differential cross-section (“intensity”) that contains the expected Wigner-\(d\) functions for the spin alignment.
Unfortunately, sympy.Indexed symbols with a superscript don’t render nicely. To fix this, AmpForm-DPD provides the following configuration function.
from ampform_dpd.io import simplify_latex_rendering
simplify_latex_rendering()
More interesting is to formulate a model with a dynamics parametrization. The most common choice is a relativistic Breit–Wigner function. It can optionally have form factors for the production vertex and/or for the decay vertex and have an energy-dependent width.
from ampform.dynamics.phasespace import PhaseSpaceFactor
from ampform_dpd.dynamics.builder import BreitWignerBuilder
dynamics_builder = BreitWignerBuilder(
energy_dependent_width=True,
production_form_factor=True,
decay_form_factor=True,
phsp_factor=PhaseSpaceFactor,
)
for chain in model_builder.decay.chains:
model_builder.dynamics_choices.register_builder(chain, dynamics_builder)
model = model_builder.formulate()
…or everything as a one-liner:
from ampform_dpd import DalitzPlotDecompositionBuilder
from ampform_dpd.dynamics.builder import BreitWignerBuilder
model_builder = DalitzPlotDecompositionBuilder(decay)
dynamics_builder = BreitWignerBuilder()
for chain in model_builder.decay.chains:
model_builder.dynamics_choices.register_builder(chain, dynamics_builder)
model = model_builder.formulate()
Here’s an example with a simple Breit–Wigner without form factors:
bw_builder = BreitWignerBuilder(
energy_dependent_width=False,
production_form_factor=False,
decay_form_factor=False,
)
for chain in model_builder.decay.chains:
model_builder.dynamics_choices.register_builder(chain, bw_builder)
simple_model = model_builder.formulate()
Custom dynamics can be defined by defining a dynamics builder function that has a signature that matches the DynamicsBuilder protocol.
import sympy as sp
from ampform_dpd import DefinedExpression
from ampform_dpd.decay import ThreeBodyDecayChain
from ampform_dpd.dynamics.builder import create_mass_symbol, get_mandelstam_s
def formulate_gaussian(decay_chain: ThreeBodyDecayChain) -> DefinedExpression:
resonance = decay_chain.resonance
s = get_mandelstam_s(decay_chain.decay_node)
x = sp.sqrt(s)
mu = create_mass_symbol(resonance)
sigma2 = sp.Symbol(Rf"\sigma_{{{resonance.latex}}}^2", nonnegative=True)
return DefinedExpression(
expression=sp.exp(-((x - mu) ** 2) / (2 * sigma2)),
parameters={
mu: resonance.mass,
sigma2: resonance.width**2,
},
)
for chain in model_builder.decay.chains:
model_builder.dynamics_choices.register_builder(chain, formulate_gaussian)
gaussian_model = model_builder.formulate()
The DalitzPlotDecompositionBuilder.formulate() method has other arguments that allow you to select a different reference subsystem and to use one complex-valued coefficient rather than two couplings per decay vertex. Also notice the tricks used here to hide the trivial summations over the helicities of spin-0 particles.
for chain in model_builder.decay.chains:
model_builder.dynamics_choices.register_builder(chain, bw_builder)
coefficient_model = model_builder.formulate(
cleanup_summations=True,
reference_subsystem=2,
use_coefficients=True,
)
coefficient_model.intensity.cleanup()