Crate laddu

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
Installation
Quick Start
- Writing a New Amplitude
- Calculating a Likelihood
Data Format
MPI Support
Future Plans
Alternatives

§Key Features

A simple interface focused on combining Amplitudes into models which can be evaluated over Datasets.
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

or

uv add 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::{
   AmplitudeID, Cache, DatasetMetadata, Expression, LadduResult, Mass,
   ParameterID, Parameter, Parameters, Resources,
};
use laddu::Event;
use laddu::amplitude::{IntoTags, Tags};
use laddu::resources::ScalarID;
use laddu::traits::*;
use laddu::math::{BarrierKind, QR_DEFAULT, Sheet, blatt_weisskopf_m, q_m};
use laddu::{Deserialize, Serialize, typetag};
use num::complex::Complex64;

#[derive(Clone, Serialize, Deserialize)]
pub struct MyBreitWigner {
    tags: Tags,
    mass: Parameter,
    width: Parameter,
    pid_mass: ParameterID,
    pid_width: ParameterID,
    l: usize,
    daughter_1_mass: Mass,
    daughter_2_mass: Mass,
    resonance_mass: Mass,
    daughter_1_mass_id: ScalarID,
    daughter_2_mass_id: ScalarID,
    resonance_mass_id: ScalarID,
}
impl MyBreitWigner {
    pub fn new(
        tags: impl IntoTags,
        mass: Parameter,
        width: Parameter,
        l: usize,
        daughter_1_mass: &Mass,
        daughter_2_mass: &Mass,
        resonance_mass: &Mass,
    ) -> LadduResult<Expression> {
        Self {
            tags: tags.into_tags(),
            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(),
            daughter_1_mass_id: ScalarID::default(),
            daughter_2_mass_id: ScalarID::default(),
            resonance_mass_id: ScalarID::default(),
        }
        .into_expression()
    }
}

#[typetag::serde]
impl Amplitude for MyBreitWigner {
    fn register(&mut self, resources: &mut Resources) -> LadduResult<AmplitudeID> {
        self.pid_mass = resources.register_parameter(&self.mass)?;
        self.pid_width = resources.register_parameter(&self.width)?;
        self.daughter_1_mass_id = resources.register_scalar(None);
        self.daughter_2_mass_id = resources.register_scalar(None);
        self.resonance_mass_id = resources.register_scalar(None);
        resources.register_amplitude(self.tags.clone())
    }

    fn bind(
        &mut self,
        metadata: &DatasetMetadata,
    ) -> LadduResult<()> {
        self.daughter_1_mass.bind(metadata)?;
        self.daughter_2_mass.bind(metadata)?;
        self.resonance_mass.bind(metadata)?;
        Ok(())
    }

    fn precompute(&self, event: &Event<'_>, cache: &mut Cache) {
        cache.store_scalar(self.daughter_1_mass_id, event.evaluate(&self.daughter_1_mass));
        cache.store_scalar(self.daughter_2_mass_id, event.evaluate(&self.daughter_2_mass));
        cache.store_scalar(self.resonance_mass_id, event.evaluate(&self.resonance_mass));
    }

    fn compute(&self, parameters: &Parameters, cache: &Cache) -> Complex64 {
        let m = cache.get_scalar(self.resonance_mass_id);
        let m0 = parameters.get(self.pid_mass).abs();
        let width0 = parameters.get(self.pid_width).abs();
        let m1 = cache.get_scalar(self.daughter_1_mass_id);
        let m2 = cache.get_scalar(self.daughter_2_mass_id);
        let q0 = q_m(m0, m1, m2, Sheet::Physical);
        let q = q_m(m, m1, m2, Sheet::Physical);
        let b0 = blatt_weisskopf_m(m0, m1, m2, self.l, QR_DEFAULT, Sheet::Physical, BarrierKind::Full);
        let b = blatt_weisskopf_m(m, m1, m2, self.l, QR_DEFAULT, Sheet::Physical, BarrierKind::Full);
        let width = width0 * (m0 / m) * (q / q0) * (b / b0).powi(2);
        1.0 / (Complex64::from(m0.powi(2) - m.powi(2)) - Complex64::I * m0 * width)
    }
}

While it isn’t shown here, we can often be more efficient when implementing Amplitudes 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, Dataset, DatasetReadOptions, Mass, NLL, parameter};
let p4_names = ["beam", "proton", "kshort1", "kshort2"];
let aux_names = ["pol_magnitude", "pol_angle"];
let options = DatasetReadOptions::default().p4_names(p4_names).aux_names(aux_names);
let ds_data = laddu::io::read_parquet("test_data/data.parquet", &options).unwrap();
let ds_mc = laddu::io::read_parquet("test_data/mc.parquet", &options).unwrap();

let resonance_mass = Mass::new(["kshort1", "kshort2"]);
let p1_mass = Mass::new(["kshort1"]);
let p2_mass = Mass::new(["kshort2"]);
let bw = MyBreitWigner::new(
    "bw",
    parameter!("mass"),
    parameter!("width"),
    2,
    &p1_mass,
    &p2_mass,
    &resonance_mass,
).unwrap();
let mag = Scalar::new("mag", parameter!("magnitude")).unwrap();
let expr = (mag * bw).norm_sqr();

let nll = NLL::new(&expr, &ds_data, &ds_mc, None).unwrap();
println!("Parameters names and order: {:?}", nll.parameters());
let result = nll.evaluate(&[1.27, 0.120, 100.0]).unwrap();
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, optional auxiliary scalars (for example, polarization magnitudes/angles), and a single column for the weight. Four-momenta are described by choosing a unique particle identifier and appending the suffixes _px, _py, _pz, and _e. Auxiliary values are listed explicitly by name. All numeric columns may be stored as Float32 or Float64; they are promoted to f64 on read. For example, the following columns 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
`beam_px`	`Float32` or `Float64`	Beam momentum (x-component)
`beam_py`	`Float32` or `Float64`	Beam momentum (y-component)
`beam_pz`	`Float32` or `Float64`	Beam momentum (z-component)
`beam_e`	`Float32` or `Float64`	Beam energy
`pol_magnitude`	`Float32` or `Float64`	Beam polarization magnitude
`pol_angle`	`Float32` or `Float64`	Beam polarization angle
`proton_px`	`Float32` or `Float64`	Recoil proton momentum (x-component)
`proton_py`	`Float32` or `Float64`	Recoil proton momentum (y-component)
`proton_pz`	`Float32` or `Float64`	Recoil proton momentum (z-component)
`proton_e`	`Float32` or `Float64`	Recoil proton energy
`kaon1_px`	`Float32` or `Float64`	Decay product 1 momentum (x)
`kaon1_py`	`Float32` or `Float64`	Decay product 1 momentum (y)
`kaon1_pz`	`Float32` or `Float64`	Decay product 1 momentum (z)
`kaon1_e`	`Float32` or `Float64`	Decay product 1 energy
`kaon2_px`	`Float32` or `Float64`	Decay product 2 momentum (x)
`kaon2_py`	`Float32` or `Float64`	Decay product 2 momentum (y)
`kaon2_pz`	`Float32` or `Float64`	Decay product 2 momentum (z)
`kaon2_e`	`Float32` or `Float64`	Decay product 2 energy
`weight`	`Float32` or `Float64`	Event weight

AmpTools-format ROOT tuples can also be loaded through the Python bindings of laddu by calling laddu.io.read_amptools(...), which performs the conversion automatically. The Rust API currently supports Parquet and standard ROOT TTrees.

§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).

The current ROOT backend always materializes the entire TTree on rank 0 before partitioning, so distributed runs do not yet see the same memory savings they do with Parquet inputs.

§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§

amplitude: Amplitudes and methods for making and evaluating them.
amplitudes: Concrete amplitude constructors and physics-model building blocks.
angular: Angular, barrier, and density-matrix factor amplitudes. Angular, barrier, and density-matrix factor amplitudes.
data: Methods for loading and manipulating event datasets.
distributions
execution: Execution-policy and thread-pool coordination helpers.
experimental
expression: Expression trees, compiled diagnostics, and evaluator interfaces.
extensions: Module for likelihood-related structures and methods
generation: Monte Carlo event generation tools.
io: Format-specific IO helpers for Datasets.
kinematics: Kinematic frame helpers and angle containers.
kmatrix: Amplitudes related to the K-Matrix formalism. K-matrix amplitudes.
likelihood: Extended maximum likelihood cost functions with support for additive terms. Extended likelihood terms and composite likelihood expressions.
lookup: Lookup-table amplitudes.
math: Special functions and numerical helpers.
optimize: A module containing optimization interfaces and algorithm integrations. Optimization helpers and algorithm integrations for likelihood terms.
parameters: Parameter handles, identifiers, and assembled parameter storage.
quantum: Quantum-number helpers and discrete analysis enums.
random: Randomized helpers used by stochastic likelihood terms. Randomized helpers used by stochastic likelihood terms.
reaction: Reaction topology, particles, and decay-node helpers.
resonance: Resonance line shapes and related factors. Resonance line shapes and related factors.
resources: Structures for manipulating the cache and free parameters.
scalar: Scalar-valued amplitude components. Scalar-valued amplitude components.
topology
traits: Useful traits for all crate structs
variables: Event variables derived from reactions and particle selections.
vectors: Three- and four-vector types used throughout the library.

Macros§

parameter: Convenience macro for creating parameters. Usage: parameter!("name") for a free parameter, or parameter!("name", 1.0) for a fixed one.

Structs§

AmplitudeID: A tag set which refers to a registered Amplitude.
AmplitudeSemanticField: A single named field in an AmplitudeSemanticKey.
AmplitudeSemanticKey: A semantic identity key used to opt into deduplicating equivalent amplitude computations.
Angles: A struct for obtaining both spherical angles at the same time.
BinnedDataset: A list of Datasets formed by binning EventData by some Variable.
BlattWeisskopf: A Blatt-Weisskopf barrier-factor amplitude for a two-body decay.
BreitWigner: A relativistic Breit-Wigner Amplitude.
BreitWignerNonRelativistic: A non-relativistic Breit-Wigner Amplitude.
Cache: A single cache entry corresponding to precomputed data for a particular EventData in a Dataset.
ClebschGordan: A Clebsch-Gordan coefficient expression.
ComplexScalar: A complex-valued Amplitude which just contains two parameters representing its real and imaginary parts.
CompositeGenerator: Generator settings for a generated composite particle.
CosTheta: A struct for obtaining the cosine of the polar angle of a decay product in a given frame of its parent resonance.
Dataset: A collection of events with optional metadata for name-based lookups.
DatasetMetadata: A collection of EventData.
DatasetReadOptions: Options for reading a Dataset from a file.
DatasetWriteOptions: Options for writing a Dataset to disk.
Decay: An isobar decay view used to derive variables and decay-local amplitudes.
Evaluator: Evaluator for Expression that mirrors the existing evaluator behavior.
Event: Borrowed, metadata-aware event access for variable and amplitude evaluation.
EventData: Raw event data in a Dataset containing all particle and auxiliary information.
EventGenerator: Event generator for generated reactions.
Expression: A holder struct that owns both an expression tree and the registered amplitudes.
ExpressionIntensity: Evaluates a real-valued Expression as a generated-batch intensity.
Flatte: A Flatte Amplitude.
GeneratedBatch: A generated dataset batch plus the metadata needed to interpret it.
GeneratedEventLayout: Metadata describing the columns and generated particles in a generated event batch.
GeneratedParticleLayout: Metadata for one generated particle in a generated event layout.
GeneratedReaction: A generated reaction layout.
GeneratedTwoToTwoReaction: A generated two-to-two reaction preserving p1 + p2 -> p3 + p4 role semantics.
GeneratedVertexLayout: Metadata for one generated vertex in a generated event layout.
HistogramSampler: Sampler for drawing values from a weighted histogram.
InitialGenerator: Generator settings for an initial-state particle.
LikelihoodExpression: A combination of LikelihoodTerms as well as sums and products of them.
LikelihoodScalar: A LikelihoodTerm which represents a single scaling parameter.
LikelihoodTermObserver: An observer which calls LikelihoodTerm::update on each step of the algorithm.
LookupAxis: A lookup-table axis defined by monotonically increasing coordinates.
LookupTable: A lookup table over one or more event Variables.
Mandelstam: A struct used to calculate Mandelstam variables (s, t, or u).
Mass: A struct for obtaining the invariant mass of a selected or reaction-defined particle.
NLL: An extended, unbinned negative log-likelihood evaluator.
OwnedEvent: Owned metadata-aware event data.
P4Selection: A reusable selection that may span one or more four-momentum names.
Parameter: An enum containing either a named free parameter or a constant value.
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.
Particle: A kinematic particle or composite system used to define a reaction.
ParticleGraph: Owned particle definitions and decay edges for a reaction.
PhaseSpaceFactor: A phase-space Amplitude for $a+b\\to c+d$ with $c\\to 1+2$ .
Phi: A struct for obtaining the azimuthal angle of a decay product in a given frame of its parent resonance.
PhotonHelicity: A real-photon helicity using the physical values +1 and -1.
PhotonSDME: A photon spin-density matrix element.
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.
PolPhase: An Amplitude representing the polarization phase $P_\gamma e^{2 i \Phi}$.
PolarComplexScalar: 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.
Reaction: A reaction with owned particle definitions and topology roles.
RejectionSampleIter: Iterator over accepted generated batches.
RejectionSampler: Rejection sampler over generated batches.
RejectionSamplingDiagnostics: Rejection-sampling diagnostics accumulated while sampling.
RejectionSamplingOptions: Options for rejection sampling generated events.
ResolvedTwoToTwo: Resolved event-level momenta for a two-to-two reaction.
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.
StableGenerator: Generator settings for a stable generated particle.
TwoToTwoReaction: A direct two-to-two reaction preserving p1 + p2 -> p3 + p4 semantics.
VariableScalar: A real-valued Amplitude which evaluates an event Variable.
Vec3: A vector with three components.
Vec4: A four-vector (Lorentz vector) whose last component stores the energy.
Voigt: A Voigt line-shape Amplitude.
Wigner3j: A Wigner-3j symbol expression.
WignerD: An amplitude evaluating a Wigner-D matrix element from decay angles.
Ylm: An Amplitude for the spherical harmonic function $Y_\ell^m(\theta, \phi)$.
Zlm: An Amplitude representing an extension of the Ylm Amplitude for a linearly polarized beam.

Enums§

Distribution
GeneratedParticle: A generated particle with generation and reconstruction metadata.
GeneratedReactionTopology: A generated reaction topology.
GeneratedStorage: Selects which generated particle four-momenta are written into generated datasets.
GeneratedVertexKind: The semantic kind of a generated vertex.
LadduError: The error type used by all laddu internal methods
LookupBoundaryMode: Out-of-range behavior used by LookupTable.
LookupInterpolation: Interpolation scheme used by LookupTable.
MandelstamTDistribution
ParameterID: An object which acts as a tag to refer to either a free parameter or a constant value.
ParticleSource: Source of a particle four-momentum.
ParticleSpecies: Experiment-neutral metadata describing a generated particle species.
PhotonPolarization: Photon polarization state used by PhotonSDME.
ReactionTopology: A reaction topology.
Reconstruction: Reconstruction interpretation for a generated particle.
RejectionEnvelope: Envelope strategy used by rejection sampling.
SimpleDistribution

Constants§

PI: The mathematical constant $\pi$.

Traits§

BatchIntensity: Evaluates unnormalized intensities for generated batches.
DecayAmplitudeExt: Extension methods that build decay-local angular expressions from Decay objects.
Deserialize: A data structure that can be deserialized from any data format supported by Serde.
IntoP4Selection: Helper trait to convert common particle specifications into P4Selection instances.
LadduGenRngExt
ProductionAmplitudeExt: Extension methods that build production-local angular expressions from Production objects.
RngSubsetExtension: An extension to Rng which allows for sampling from a subset of the integers [0..n) without replacement.
Serialize: A data structure that can be serialized into any data format supported by Serde.
VariableExpressionExt: Extension methods for building expressions from event Variables.

Type Aliases§

LadduResult: A Result type alias for LadduErrors.

Derive Macros§

Deserialize
Serialize

Crate laddu

Crate laddu Copy item path

§Table of Contents

§Key Features

§Installation

§Quick Start

§Rust

§Writing a New Amplitude

§Calculating a Likelihood

§Data Format

§MPI Support

§Future Plans

§Alternatives

§AmpTools

§PyPWA

§ComPWA

§Others

Re-exports§

Modules§

Macros§

Structs§

Enums§

Constants§

Traits§

Type Aliases§

Derive Macros§

Crate laddu