Expand description
laddu
(/ˈlʌduː/) is a library for analysis of particle physics data. It is intended to be a simple and efficient alternative to some of the other tools out there. laddu
is written in Rust with bindings to Python via PyO3
and maturin
and is the spiritual successor to rustitude
, one of my first Rust projects. The goal of this project is to allow users to perform complex amplitude analyses (like partial-wave analyses) without complex code or configuration files.
This crate is still in an early development phase, and the API is not stable. It can (and likely will) be subject to breaking changes before the 1.0.0 version release (and hopefully not many after that).
§Table of Contents
§Key Features
- A simple interface focused on combining
Amplitude
s into models which can be evaluated overDataset
s. - A single
Amplitude
trait which makes it easy to write new amplitudes and integrate them into the library. - Easy interfaces to precompute and cache values before the main calculation to speed up model evaluations.
- Efficient parallelism using
rayon
. - Python bindings to allow users to write quick, easy-to-read code that just works.
§Installation
laddu
can be added to a Rust project with cargo
:
cargo add laddu
The library’s Python bindings are located in a library by the same name, which can be installed simply with your favorite Python package manager:
pip install laddu
§Quick Start
§Rust
§Writing a New Amplitude
At the time of writing, Rust is not a common language used by particle physics, but this tutorial should hopefully convince the reader that they don’t have to know the intricacies of Rust to write performant amplitudes. As an example, here is how one might write a Breit-Wigner, parameterized as follows:
I_{\ell}(m; m_0, \Gamma_0, m_1, m_2) = \frac{1}{\pi}\frac{m_0 \Gamma_0 B_{\ell}(m, m_1, m_2)}{(m_0^2 - m^2) - \imath m_0 \Gamma}
where
\Gamma = \Gamma_0 \frac{m_0}{m} \frac{q(m, m_1, m_2)}{q(m_0, m_1, m_2)} \left(\frac{B_{\ell}(m, m_1, m_2)}{B_{\ell}(m_0, m_1, m_2)}\right)^2
is the relativistic width correction, $q(m_a, m_b, m_c)
$ is the breakup momentum of a particle with mass $m_a
$ decaying into two particles with masses $m_b
$ and $m_c
$, $B_{\ell}(m_a, m_b, m_c)
$ is the Blatt-Weisskopf barrier factor for the same decay assuming particle $a
$ has angular momentum $\ell
$, $m_0
$ is the mass of the resonance, $\Gamma_0
$ is the nominal width of the resonance, $m_1
$ and $m_2
$ are the masses of the decay products, and $m
$ is the “input” mass.
Although this particular amplitude is already included in laddu
, let’s assume it isn’t and imagine how we would write it from scratch:
use laddu::{
ParameterLike, Event, Cache, Resources, Mass,
ParameterID, Parameters, Float, LadduError, PI, AmplitudeID, Complex,
};
use laddu::traits::*;
use laddu::utils::functions::{blatt_weisskopf, breakup_momentum};
use laddu::{Deserialize, Serialize, typetag};
#[derive(Clone, Serialize, Deserialize)]
pub struct MyBreitWigner {
name: String,
mass: ParameterLike,
width: ParameterLike,
pid_mass: ParameterID,
pid_width: ParameterID,
l: usize,
daughter_1_mass: Mass,
daughter_2_mass: Mass,
resonance_mass: Mass,
}
impl MyBreitWigner {
pub fn new(
name: &str,
mass: ParameterLike,
width: ParameterLike,
l: usize,
daughter_1_mass: &Mass,
daughter_2_mass: &Mass,
resonance_mass: &Mass,
) -> Box<Self> {
Self {
name: name.to_string(),
mass,
width,
pid_mass: ParameterID::default(),
pid_width: ParameterID::default(),
l,
daughter_1_mass: daughter_1_mass.clone(),
daughter_2_mass: daughter_2_mass.clone(),
resonance_mass: resonance_mass.clone(),
}
.into()
}
}
#[typetag::serde]
impl Amplitude for MyBreitWigner {
fn register(&mut self, resources: &mut Resources) -> Result<AmplitudeID, LadduError> {
self.pid_mass = resources.register_parameter(&self.mass);
self.pid_width = resources.register_parameter(&self.width);
resources.register_amplitude(&self.name)
}
fn compute(&self, parameters: &Parameters, event: &Event, _cache: &Cache) -> Complex<Float> {
let mass = self.resonance_mass.value(event);
let mass0 = parameters.get(self.pid_mass);
let width0 = parameters.get(self.pid_width);
let mass1 = self.daughter_1_mass.value(event);
let mass2 = self.daughter_2_mass.value(event);
let q0 = breakup_momentum(mass0, mass1, mass2);
let q = breakup_momentum(mass, mass1, mass2);
let f0 = blatt_weisskopf(mass0, mass1, mass2, self.l);
let f = blatt_weisskopf(mass, mass1, mass2, self.l);
let width = width0 * (mass0 / mass) * (q / q0) * (f / f0).powi(2);
let n = Float::sqrt(mass0 * width0 / PI);
let d = Complex::new(mass0.powi(2) - mass.powi(2), -(mass0 * width));
Complex::from(f * n) / d
}
}
While it isn’t shown here, we can often be more efficient when implementing
Amplitude
s by precomputing values which do not depend on the
free parameters. See the Amplitude::precompute
method for more details.
§Calculating a Likelihood
We could then write some code to use this amplitude. For demonstration purposes, let’s just calculate an extended unbinned negative log-likelihood, assuming we have some data and Monte Carlo in the proper parquet format:
use laddu::{Scalar, Mass, Manager, NLL, parameter, open};
let ds_data = open("test_data/data.parquet").unwrap();
let ds_mc = open("test_data/mc.parquet").unwrap();
let resonance_mass = Mass::new([2, 3]);
let p1_mass = Mass::new([2]);
let p2_mass = Mass::new([3]);
let mut manager = Manager::default();
let bw = manager.register(MyBreitWigner::new(
"bw",
parameter("mass"),
parameter("width"),
2,
&p1_mass,
&p2_mass,
&resonance_mass,
)).unwrap();
let mag = manager.register(Scalar::new("mag", parameter("magnitude"))).unwrap();
let expr = (mag * bw).norm_sqr();
let model = manager.model(&expr);
let nll = NLL::new(&model, &ds_data, &ds_mc);
println!("Parameters names and order: {:?}", nll.parameters());
let result = nll.evaluate(&[1.27, 0.120, 100.0]);
println!("The extended negative log-likelihood is {}", result);
In practice, amplitudes can also be added together, their real and imaginary parts can be taken, and evaluators should mostly take the real part of whatever complex value comes out of the model.
§Data Format
The data format for laddu
is a bit different from some of the alternatives like AmpTools
. Since ROOT doesn’t yet have bindings to Rust and projects to read ROOT files are still largely works in progress (although I hope to use oxyroot
in the future when I can figure out a few bugs), the primary interface for data in laddu
is Parquet files. These are easily accessible from almost any other language and they don’t take up much more space than ROOT files. In the interest of future compatibility with any number of experimental setups, the data format consists of an arbitrary number of columns containing the four-momenta of each particle, the polarization vector of each particle (optional) and a single column for the weight. These columns all have standardized names. For example, the following columns would describe a dataset with four particles, the first of which is a polarized photon beam, as in the GlueX experiment:
Column name | Data Type | Interpretation |
---|---|---|
p4_0_E | Float32 | Beam Energy |
p4_0_Px | Float32 | Beam Momentum (x-component) |
p4_0_Py | Float32 | Beam Momentum (y-component) |
p4_0_Pz | Float32 | Beam Momentum (z-component) |
aux_0_x | Float32 | Beam Polarization (x-component) |
aux_0_y | Float32 | Beam Polarization (y-component) |
aux_0_z | Float32 | Beam Polarization (z-component) |
p4_1_E | Float32 | Recoil Proton Energy |
p4_1_Px | Float32 | Recoil Proton Momentum (x-component) |
p4_1_Py | Float32 | Recoil Proton Momentum (y-component) |
p4_1_Pz | Float32 | Recoil Proton Momentum (z-component) |
p4_2_E | Float32 | Decay Product 1 Energy |
p4_2_Px | Float32 | Decay Product 1 Momentum (x-component) |
p4_2_Py | Float32 | Decay Product 1 Momentum (y-component) |
p4_2_Pz | Float32 | Decay Product 1 Momentum (z-component) |
p4_3_E | Float32 | Decay Product 2 Energy |
p4_3_Px | Float32 | Decay Product 2 Momentum (x-component) |
p4_3_Py | Float32 | Decay Product 2 Momentum (y-component) |
p4_3_Pz | Float32 | Decay Product 2 Momentum (z-component) |
weight | Float32 | Event Weight |
To make it easier to get started, we can directly convert from the AmpTools
format using the provided [amptools-to-laddu
] script (see the bin
directory of this repository). This is not bundled with the Python library (yet) but may be in the future.
§MPI Support
The latest version of laddu
supports the Message Passing Interface (MPI) protocol for distributed computing. MPI-compatible versions of the core laddu
methods have been written behind the mpi
feature gate. To build laddu
with MPI compatibility, it can be added with the mpi
feature via cargo add laddu --features mpi
. Note that this requires a working MPI installation, and OpenMPI or MPICH are recommended, as well as LLVM/Clang. The installation of these packages differs by system, but are generally available via each system’s package manager.
To use MPI in Rust, one must simply surround their main analysis code with a call to laddu::mpi::use_mpi(true)
and laddu::mpi::finalize_mpi()
. The first method has a boolean flag which allows for runtime switching of MPI use (for example, disabling MPI with an environment variable).
§Future Plans
- GPU integration (this is incredibly difficult to do right now, but it’s something I’m looking into).
- As always, more tests and documentation.
§Alternatives
While this is likely the first Rust project (aside from my previous attempt, rustitude
), there are several other amplitude analysis programs out there at time of writing. This library is a rewrite of rustitude
which was written when I was just learning Rust and didn’t have a firm grasp of a lot of the core concepts that are required to make the analysis pipeline memory- and CPU-efficient. In particular, rustitude
worked well, but ate up a ton of memory and did not handle precalculation as nicely.
§AmpTools
The main inspiration for this project is the library most of my collaboration uses, AmpTools
. AmpTools
has several advantages over laddu
: it’s probably faster for almost every use case, but this is mainly because it is fully integrated with MPI and GPU support. I’m not actually sure if there’s a fair benchmark between the two libraries, but I’d wager AmpTools
would still win. AmpTools
is a much older, more developed project, dating back to 2010. However, it does have its disadvantages. First and foremost, the primary interaction with the library is through configuration files which are not really code and sort of represent a domain specific language. As such, there isn’t really a way to check if a particular config will work before running it. Users could technically code up their analyses in C++ as well, but I think this would generally be more work for very little benefit. AmpTools primarily interacts with Minuit, so there aren’t simple ways to perform alternative optimization algorithms, and the outputs are a file which must also be parsed by code written by the user. This usually means some boilerplate setup for each analysis, a slew of input and output files, and, since it doesn’t ship with any amplitudes, integration with other libraries. The data format is also very rigid, to the point where including beam polarization information feels hacked on (see the Zlm implementation here which requires the event-by-event polarization to be stored in the beam’s four-momentum). While there isn’t an official Python interface, Lawrence Ng has made some progress porting the code here.
§PyPWA
PyPWA
is a library written in pure Python. While this might seem like an issue for performance (and it sort of is), the library has several features which encourage the use of JIT compilers. The upside is that analyses can be quickly prototyped and run with very few dependencies, it can even run on GPUs and use multiprocessing. The downside is that recent development has been slow and the actual implementation of common amplitudes is, in my opinion, messy. I don’t think that’s a reason to not use it, but it does make it difficult for new users to get started.
§ComPWA
ComPWA
is a newcomer to the field. It’s also a pure Python implementation and is comprised of three separate libraries. QRules
can be used to validate and generate particle reaction topologies using conservation rules. AmpForm
uses SymPy
to transform these topologies into mathematical expressions, and it can also simplify the mathematical forms through the built-in CAS of SymPy
. Finally, TensorWaves
connects AmpForm
to various fitting methods. In general, these libraries have tons of neat features, are well-documented, and are really quite nice to use. I would like to eventually see laddu
as a companion to ComPWA
(rather than direct competition), but I don’t really know enough about the libraries to say much more than that.
§Others
It could be the case that I am leaving out software with which I am not familiar. If so, I’d love to include it here for reference. I don’t think that laddu
will ever be the end-all-be-all of amplitude analysis, just an alternative that might improve on existing systems. It is important for physicists to be aware of these alternatives. For example, if you really don’t want to learn Rust but need to implement an amplitude which isn’t already included here, laddu
isn’t for you, and one of these alternatives might be best.
Re-exports§
pub use typetag;
Modules§
- amplitudes
Amplitude
s and methods for making and evaluating them.- breit_
wigner - The Breit-Wigner amplitude.
- common
- Common amplitudes (like a scalar value which just contains a single free parameter).
- data
- Methods for loading and manipulating
Event
-based data. - experimental
- extensions
- Module for likelihood-related structures and methods
- ganesh_
ext - A module containing the
laddu
interface with theganesh
library- kmatrix
- Amplitudes related to the K-Matrix formalism.
- likelihoods
- Extended maximum likelihood cost functions with support for additive terms
- phase_
space - A phase space factor for
$a+b\to c+d$
with$c\to 1+2$
.- piecewise
- Piecewise functions as amplitudes.
- resources
- Structures for manipulating the cache and free parameters.
- traits
- Useful traits for all crate structs
- utils
- Utility functions, enums, and traits
- ylm
- A spherical harmonic amplitude.
- zlm
- A polarized spherical harmonic amplitude.
Structs§
- AmplitudeID
- A tag which refers to a registered
Amplitude
. This is the base object which can be used to buildExpression
s and should be obtained from theManager::register
method. - Angles
- A struct for obtaining both spherical angles at the same time.
- Binned
Dataset - A list of
Dataset
s formed by binningEvent
s by someVariable
. - Breit
Wigner - A relativistic Breit-Wigner
Amplitude
, parameterized as follows: - Cache
- A single cache entry corresponding to precomputed data for a particular
Event
in aDataset
. - Complex
- A complex number in Cartesian form.
- Complex
Scalar - A complex-valued
Amplitude
which just contains two parameters representing its real and imaginary parts. - CosTheta
- A struct for obtaining the $
\cos\theta
$ (cosine of the polar angle) of a decay product in a given reference frame of its parent resonance. - Dataset
- A collection of
Event
s. - Ensemble
- A collection of
Walker
s - Evaluator
- A structure which can be used to evaluate the stored
Expression
built on registeredAmplitude
s. This contains aResources
struct which already contains cached values for precomputedAmplitude
s and any relevant free parameters and constants. - Event
- A single event in a
Dataset
containing all the relevant particle information. - Likelihood
Evaluator - A structure to evaluate and minimize combinations of
LikelihoodTerm
s. - LikelihoodID
- An identifier that can be used like an
AmplitudeID
to combine registeredLikelihoodTerm
s. - Likelihood
Manager - A
Manager
but forLikelihoodTerm
s. - Likelihood
Scalar - A
LikelihoodTerm
which represents a single scaling parameter. - MCMC
Options - A set of options that are used when Markov Chain Monte Carlo samplings are performed.
- Manager
- A manager which can be used to register
Amplitude
s withResources
. This structure is essential to any analysis and should be constructed using theManager::default()
method. - Mandelstam
- A struct used to calculate Mandelstam variables ($
s
$, $t
$, or $u
$). - Mass
- A struct for obtaining the mass of a particle by indexing the four-momenta of an event, adding together multiple four-momenta if more than one index is given.
- Minimizer
Options - A set of options that are used when minimizations are performed.
- Model
- A struct which contains a set of registerd
Amplitude
s (inside aManager
) and anExpression
. - NLL
- An extended, unbinned negative log-likelihood evaluator.
- Parameters
- This struct holds references to the constants and free parameters used in the fit so that they
may be obtained from their corresponding
ParameterID
. - Phase
Space Factor - An
Amplitude
describing the phase space factor given in Equation A4 here - Phi
- A struct for obtaining the $
\phi
$ angle (azimuthal angle) of a decay product in a given reference frame of its parent resonance. - Piecewise
Complex Scalar - A piecewise complex-valued
Amplitude
which just contains two parameters representing its real and imaginary parts. - Piecewise
Polar Complex Scalar - A piecewise complex-valued
Amplitude
which just contains two parameters representing its magnitude and phase. - Piecewise
Scalar - A piecewise scalar-valued
Amplitude
which just contains a single parameter for each bin as its value. - PolAngle
- A struct defining the polarization angle for a beam relative to the production plane.
- PolMagnitude
- A struct defining the polarization magnitude for a beam relative to the production plane.
- Polar
Complex Scalar - A complex-valued
Amplitude
which just contains two parameters representing its magnitude and phase. - Polarization
- A struct for obtaining both the polarization angle and magnitude at the same time.
- Resources
- The main resource manager for cached values, amplitudes, parameters, and constants.
- Scalar
- A scalar-valued
Amplitude
which just contains a single parameter as its value. - Status
- A status message struct containing all information about a minimization result.
- Vec3
- A vector with three components
- Vec4
- A vector with four components (a Lorentz vector)
- Ylm
- An
Amplitude
for the spherical harmonic function $Y_\ell^m(\theta, \phi)
$. - Zlm
- An
Amplitude
representing an extension of theYlm
Amplitude
assuming a linearly polarized beam as described in Equation (D13) here
Enums§
- Expression
- An expression tree which contains
AmplitudeID
s and operators over them. - Laddu
Error - The error type used by all
laddu
internal methods - Likelihood
Expression - A combination of
LikelihoodTerm
s as well as sums and products of them. - ParameterID
- An object which acts as a tag to refer to either a free parameter or a constant value.
- Parameter
Like - An enum containing either a named free parameter or a constant value.
Constants§
- PI
- The mathematical constant $
\pi
$.
Traits§
- Deserialize
- A data structure that can be deserialized from any data format supported by Serde.
- Serialize
- A data structure that can be serialized into any data format supported by Serde.
Functions§
- constant
- Shorthand for generating a constant value (which acts like a fixed parameter).
- open
- Open a Parquet file and read the data into a
Dataset
. - parameter
- Shorthand for generating a named free parameter.
Type Aliases§
- DVector
- A dynamically sized column vector.
- Float
- A floating-point number type (defaults to
f64
, seef32
feature).
Derive Macros§