Decays of unstable heavy particles to multibody final states can in general occur through several different intermediate resonances.
Each decay channel can be represented quantum-mechanically by an amplitude, and the total density of decays across the phase space is represented by the square of the coherent sum of all contributing amplitudes.
Interference effects can lead to excesses or deficits of decays in regions of phase space where different resonances overlap.
Investigations of such dynamical effects in multibody decays are of great interest to test the Standard Model of particle physics and to investigate resonant structures.
The Dalitz plot (DP)~\cite{Dalitz:1953cp,Fabri:1954zz} was introduced originally to describe the phase space of $\KL \to \pi\pi\pi$ decays, but is relevant for the decay of any spin-zero particle to three spin-zero particles, $P \to d_{1} d_{2} d_{3}$.
In such a case, energy and momentum conservation give
where $m(d_id_j)$ is the invariant mass obtained from the two-body combination of the $d_i$ and $d_j$ four momenta.
Consequently, assuming that the masses of $P$, $d_1$, $d_2$ and $d_3$ are all known, any two of the $m^2(d_id_j)$ values -- subsequently referred to as Dalitz-plot variables -- are sufficient to describe fully the kinematics of the decay in the $P$ rest frame.
This can also be shown by considering that the 12 degrees of freedom corresponding to the four-momenta of the three final-state particles are accounted for by two DP variables, the three $d_i$ masses, four constraints due to energy-momentum conservation in the $P \to d_{1} d_{2} d_{3}$ decay, and three co-ordinates describing a direction in space which carries no physical information about the decay since all particles involved have zero spin.
A Dalitz plot is then the visualisation of the phase space of a particular three-body decay in terms of the two DP variables.\footnote{
The phrase ``Dalitz plot'' is often used more broadly in the literature.
In particular, it can be used to describe the projection onto two of the two-body invariant mass combinations of a three-body decay even when one or more of the particles involved has non-zero spin.
}
Analysis of the distribution of decays across a DP can reveal information about the underlying dynamics of the particular three-body decay, since the differential rate is
(Top left) kinematic boundaries of the three-body phase space for the decay $\Bs\to\Dzb\Km\pip$.
The insets indicate the configuration of the final-state particle momenta in the parent rest frame at various different DP positions.
(Top right) examples of the resonances which may appear in the Dalitz
plot for this decay: (red) $\D^*_2(2573)^-$, (orange) $\kaon^*(892)^0$,
(green) $\kaon\pion$ S-wave.
(Bottom) projections of this DP onto the squares of the invariant masses (from left to right): $m^2_{\Dzb\Km}$, $m^2_{\Dzb\pip}$, $m^2_{\Km\pip}$.
}
\label{fig:dpexample}
\end{figure}
The Dalitz-plot analysis technique, usually implemented with model-dependent descriptions of the amplitudes involved, has been used to understand hadronic effects in, for example, the $\pi^0\pi^0\pi^0$ system produced in $p\bar{p}$ annihilation~\cite{Amsler:1995bf}.
Recently, it has also been used to study three-body $\eta_c$ decays~\cite{Lees:2014iua,Lees:2015zzr}.
However, DP analyses have become particularly popular to study multibody decays of the heavy-flavoured $D$ and $B$ mesons.
Not only do the relatively large masses of these particles provide a broad kinematic range in which resonant structures can be studied but, since the decays are mediated by the weak interaction, there may be \CP-violating differences between the DP distributions for particle and antiparticle.
Studying these differences can test the Standard Model mechanism for \CP violation: if the asymmetries are not consistent with originating from the single complex phase in the Cabibbo-Kobayashi-Maskawa (CKM) quark mixing matrix~\cite{Cabibbo:1963yz,Kobayashi:1973fv} then contributions beyond the Standard Model must be present.
% as predicted for certain decay modes in various theories.
Until around the year 2000, most DP analyses of charm decays were focussed on understanding hadronic structures at low $\pi\pi$ or $K\pi$ mass.
In particular, pioneering analyses of $D \to K\pi\pi$ decays were carried out by experiments such as MARK-II, MARK-III, E687, ARGUS, E691 and CLEO~\cite{Schindler:1980ws,Adler:1987sd,Anjos:1992kb,Albrecht:1993jn,Frabetti:1994di,Kopp:2000gv}.
These analyses revealed the existence of a broad structure in the $K\pi$ S-wave that could not be well described with a Breit--Wigner lineshape.
In later analyses, it was shown that this contribution could be modelled in a quasi-model-independent way, in which the partial wave is fitted using splines to describe the magnitude and phase as a function of $m(K\pi)$~\cite{Aitala:2005yh}.
Subsequent uses of this approach include further studies of the $K\pi$ S-wave~\cite{Bonvicini:2008jw,Link:2009ng,delAmoSanchez:2010fd,Lees:2015zzr} as well as the $\Kp\Km$~\cite{Aubert:2007dc} and $\pip\pim$~\cite{Aubert:2008ao} S-wave{s}, in various processes.
Similarly, DP analyses of decays such as $\Dp \to \pip\pip\pim$~\cite{Aitala:2000xu,Muramatsu:2002jp,Link:2003gb,Bonvicini:2007tc} indicated the existence of a broad low-mass $\pi\pi$ S-wave known as the $\sigma$ pole~\cite{Pelaez:2015qba}.
With the advent of the $e^+e^-$ \B-factory experiments, \babar~\cite{Aubert:2001tu,TheBABAR:2013jta} and Belle~\cite{Abashian:2000cg}, DP analyses of \B meson decays became feasible.
The method was used to obtain insights into charm resonances through analyses of $\Bp\to \Dm\pip\pip$~\cite{Abe:2003zm,Aubert:2009wg} and $\Bz \to \Dzb\pip\pim$~\cite{Kuzmin:2006mw,delAmoSanchez:2010ad} decays.
Studies of \B meson decays to final states without any charm or charmonium particles also became possible~\cite{Garmash:2004wa,Aubert:2005sk,Aubert:2005ce}.
-Once baseline DP models were established, it was then possible to search for \CP violation effects, with results including the first evidence for \CP violation in the $\Bp \to \rho^0 \Kp$ decay~\cite{Garmash:2005rv,Aubert:2008bj}.
+Once baseline DP models were established, it was then possible to search for \CP violation effects, with results including the first evidence for \CP violation in the $\Bp \to \rho(770)^0 \Kp$ decay~\cite{Garmash:2005rv,Aubert:2008bj}.
Moreover, analyses that accounted for possible dependence of the \CP violation effect with decay time as well as with DP position were carried out for both $D$~\cite{Asner:2005sz,Abe:2007rd} and $B$ decays~\cite{Aubert:2007jn,Kusaka:2007dv,Kusaka:2007mj,Aubert:2007sd,:2008wwa,Aubert:2009me,Nakahama:2010nj,Lees:2013nwa}.
With the availability of increasingly large data samples at these experiments and, more recently, at the Large Hadron Collider experiments (in particular, LHCb~\cite{Alves:2008zz}), more detailed studies of these and similar decays become possible.
In addition, many ideas for DP analyses have been proposed, since they provide interesting possibilities to provide insight into hadronic structures, to measure \CP violation effects and to test the Standard Model.
These include methods to determine the angles $\alpha$, $\beta$ and $\gamma$ of the CKM Unitarity Triangle with low theoretical uncertainty from, respectively $\Bz \to \pip\pim\piz$~\cite{Snyder:1993mx}, $\Bz \to D\pip\pim$~\cite{Charles:1998vf,Latham:2008zs} and $\Bz \to D\Kp\pim$ decays~\cite{Gershon:2008pe,Gershon:2009qc}, among many other potential analyses.
Thus, it has become increasingly important to have a publicly available Dalitz-plot analysis package that is flexible enough both to be used in a range of experimental environments and to describe many possible different decays and types of analyses.
Such a package should be well validated and have excellent performance characteristics, in particular in terms of speed since complicated amplitude fits can otherwise have unacceptable CPU requirements.
This motivated the creation, and ongoing development, of the \laura\ package, which is described in the remainder of the paper.
\laura\ is written in the \cpp\ programming language and is intended to be as close as possible to being a standalone package, with a sole external dependency on the \root\ package~\cite{Brun:1997pa}.
In particular, \root\ is used to handle histogrammed quantities, and to implement minimisation of negative log-likelihood functions with \minuit~\cite{James:1975dr}.
Further documentation and code releases, distributed under the Boost Software License~\cite{boost}, are available from
The description of the software given in this paper corresponds to that
released in \laura\ version {\tt v3r2}.
In Sec.~\ref{sec:DalitzGeneralities}, a brief summary of the Dalitz-plot analysis formalism is given, and the conventions used in \laura\ are set out.
Section~\ref{sec:expt-effects} describes effects that must also be taken into account when performing an experimental analysis.
Sections~\ref{sec:signal},~\ref{sec:eff-resol} and~\ref{sec:backgrounds} then contain discussions of, respectively, the implementation of the signal model, efficiency and resolution effects, and the background components in \laura, including explicit classes and methods with high-level details given in Appendices.
These elements are then put together in Sec.~\ref{sec:workflow}, where the overall work flow in \laura\ is described.
The performance of the software is discussed in Sec.~\ref{sec:performance}, ongoing and planned future developments are briefly mentioned in Sec.~\ref{sec:future-devel}, and a summary is given in Sec.~\ref{sec:summary}.
where $m_0$ is the pole mass and $\Gamma_0$ is the resonance width.
The more commonly used relativistic Breit--Wigner lineshape is described in Sec.~\ref{sec:lineshapes}.
We note here that the relativistic Breit--Wigner lineshape can also describe so-called virtual contributions, from resonances with masses outside the kinematically accessible region of the Dalitz plot, with one modification:
in the calculation of the momenta, the mass $m_0$ is set to a value $m_0^{\rm{eff}}$ within the kinematically allowed range.
This is accomplished with the {\it ad-hoc} formula
where the fractional coefficients come from isospin conservation,
$m_{ij}$ denotes the invariant mass of the daughter particle $j$ (1-4) in channel $i$ (1-2),
and $g_1$ and $g_2$ are coupling constants whose values are assumed to be those
provided in Table~\ref{tab:flattepars}.
The Clebsch-Gordan coefficients in Eqs.~\ref{eq:Gamma1} and~\ref{eq:Gamma2} are not guaranteed to be correct for every possible resonance that could be modelled with the Flatt\'e lineshape, but are appropriate for every case considered in Table~\ref{tab:flattepars}.
The expressions of the widths are continued analytically ($\Gamma \ra i |\Gamma|$) when $m$ is below any of the specific channel thresholds,
contributing to the real part of the amplitude, while the Adler-zero
term $f_A = (m^2 - s_A)/(m^2_0 - s_A)$ can be used to suppress false kinematic
singularities when $m$ goes below threshold~\cite{Bugg:2003kj}
(otherwise $f_A$ is set to unity).
%%
\begin{table}[!htb]
\caption{
The four daughter particles used for each channel term $m_{ij}$, as well as the coupling ($g_1$, $g_2$) and Adler-zero ($s_A$) constants for the Flatt\'e lineshapes.
Units of \gev\ for $g_{1,2}$ and \gevgevcccc\ for $s_A$ are implied.