Laura++ : a Dalitz plot fitter

The Dalitz plot analysis technique has become an increasingly important method in heavy flavour physics. The Laura++ fitter has been developed as a flexible tool that can be used for Dalitz plot analyses in different experimental environments. Explicitly designed for three-body decays of heavy-flavoured mesons to spinless final state particles, it is optimised in order to describe all possible resonant or nonresonant contributions, and to accommodate possible CP violation effects.


Introduction
Decays of unstable heavy particles to multibody final states can in general occur through several different intermediate resonances. Each decay channel can be represented quantummechanically 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) [1,2] was introduced originally to describe the phase space of K 0 L → πππ decays, but is relevant for the decay of any spin-zero particle to three spin-zero particles, P → d 1 d 2 d 3 . In such a case, energy and momentum conservation give where m(d i d 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 i d 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 → 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. 1 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 where A is the amplitude for the three-body decay. Thus, any deviation from a uniform distribution is due to the dynamical structure of the amplitude. Examples of the kinematic boundaries of a DP, and of resonant structures that may appear in this kind of decay, are shown in Fig. 1. 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 π 0 π 0 π 0 system produced in pp annihilation [3]. Recently, it has also been used to study three-body η c decays [4,5]. 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 [6,7] then contributions beyond the Standard Model must be present. Until around the year 2000, most DP analyses of charm decays were focussed on understanding hadronic structures at low ππ or Kπ mass. In particular, pioneering analyses of D → Kππ decays were carried out by experiments such as MARK-II, MARK-III, E687, ARGUS, E691 and CLEO [8][9][10][11][12][13]. These analyses revealed the existence of a broad structure in the Kπ 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π) [14]. Subsequent uses of this approach include further studies of the Kπ S-wave [5,[15][16][17] as well as the K + K − [18] and π + π − [19] S-waves, in various processes. Similarly, DP analyses of decays such as D + → π + π + π − [20][21][22][23] indicated the existence of a broad low-mass ππ S-wave known as the σ pole [24].
With the advent of the e + e − B-factory experiments, BaBar [25,26] and Belle [27], DP analyses of B meson decays became feasible. The method was used to obtain insights into charm resonances through analyses of B + → D − π + π + [28,29] and B 0 → D 0 π + π − [30,31] decays. Studies of B meson decays to final states without any charm or charmonium particles also became possible [32][33][34]. 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 B + → ρ(770) 0 K + decay [35,36]. 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 [37,38] and B decays [39][40][41][42][43][44][45][46].
With the availability of increasingly large data samples at these experiments and, more recently, at the Large Hadron Collider experiments (in particular, LHCb [47]), 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 α, β and γ of the CKM Unitarity Triangle with low theoretical uncertainty from, respectively B 0 → π + π − π 0 [48], B 0 → Dπ + π − [49,50] and B 0 → DK + π − decays [51,52], 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 C ++ 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 [53]. In particular, Root is used to handle data file input/output, histogrammed quantities, and the minimisation of negative log-likelihood functions with Minuit [54]. Further documentation and code releases (distributed under the Apache License, Version 2.0 [55]) are available from http://laura.hepforge.org/ .
The description of the software given in this paper corresponds to that released in Laura ++ version v3r4.
In Sec. 2, a brief summary of the Dalitz-plot analysis formalism is given, and the conventions used in Laura ++ are set out. Section 3 describes effects that must also be taken into account when performing an experimental analysis. Sections 4, 5 and 6 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. 7, where the overall work flow in Laura ++ is described. The performance of the software is discussed in Sec. 8, ongoing and planned future developments are briefly mentioned in Sec. 9, and a summary is given in Sec. 10.

Dalitz-plot analysis formalism
Given two variables that describe the Dalitz plot of the P → d 1 d 2 d 3 decay, all other kinematic quantities can be uniquely determined for fixed initial-and final-state (subsequently referred to as parent and daughter) particle masses. The convention adopted in Laura ++ is that the DP is described in terms of m 2 13 ≡ m 2 (d 1 d 3 ) and m 2 23 ≡ m 2 (d 2 d 3 ). Hence, these two variables are required to be present in any input data provided to Laura ++ .
The description of the complex amplitude is based on the isobar model [56][57][58], which describes the total amplitude as a coherent sum of N amplitudes from resonant or nonresonant intermediate processes. 2 This means that the total amplitude is given by where c j are complex coefficients, discussed further in Sec. 4.5, giving the relative contribution of decay channel j. There are several different choices used in the literature to express the resonance dynamics contained within the F j (m 2 13 , m 2 23 ) terms. Here, one common approach, which is the default in Laura ++ , is outlined; other possibilities are discussed in Appendices A and B. For a resonance in m 13 , the dynamics can be written as F m 2 13 , m 2 23 = N × R (m 13 ) × T ( p, q) × X(| p| r P BW ) × X(| q| r R BW ) , where the functions R and T describe the invariant mass and angular dependence of the amplitude, the X functions are form factors, and N is a normalisation constant. In Eq. (4) only the kinematic -i.e. DP -dependence has been specified; the functions may also depend on properties of the resonance such as mass, width and spin. The arguments q and p are the momentum of one of the resonance decay products (d 3 in this case, see Sec. 2.4 for further information) and that of the so-called "bachelor" particle (i.e. the particle not associated with the decay of the resonance; d 2 in this case), both evaluated in the rest frame of the resonance. The parameters r P BW and r R BW are characteristic meson radii described below. The resonance dynamics are normalised in Laura ++ such that the integral over the DP of the squared magnitude of each term is unity DP F j m 2 13 , m 2 23 2 dm 2 13 dm 2 23 = 1 .
Although not strictly necessary, as only the total probability density function (PDF) needs to be normalised, this allows a meaningful comparison of the values of the c j coefficients.

Resonance lineshapes
In Eq. (4), the function R (m 13 ) is the resonance mass term. The detailed forms for all available shapes in Laura ++ are given in App. A. Here the most commonly used relativistic Breit-Wigner (RBW) lineshape is given as an example where m 0 is the nominal mass of the resonance and the dependence of the decay width of the resonance on m is given by where Γ 0 is the nominal width of the resonance and q 0 denotes the value of q when m = m 0 . In Eq. (7), L is the orbital angular momentum between the resonance daughters. Note that since all the initial-and final-state particles have zero spin, this quantity is the same as the spin of the resonance and is also the same as the orbital angular momentum between the resonance and the bachelor. It is relevant to note that Eq. (6) can be written where tan φ = m 0 Γ(m) m 2 0 −m 2 . This shows the characteristic phase rotation of a resonance as m 2 increases from far below to far above m 2 0 .

Angular distributions and Blatt-Weisskopf form factors
Using the Zemach tensor formalism [59,60], the angular probability distribution terms T ( p, q) are given by L = 4 : T ( p, q) = 16 35 35( p · q ) 4 − 30( p · q ) 2 (p q) 2 + 3(p q) 4 , L = 5 : T ( p, q) = − 32 63 63( p · q ) 5 − 70( p · q ) 3 (p q) 2 + 15( p · q )(p q) 4 , (14) where q ≡ | q | and p ≡ | p |. These have the form of the Legendre polynomials P L (cos θ), where θ is the "helicity" angle between p and q, multiplied by the appropriate power of −2 p q. These factors act to suppress the amplitude at low values of the break-up momentum in either the decay of the parent or the resonance -the so-called "angular momentum barrier". However, these factors on their own would cause the amplitude to continue to grow with increasing break-up momentum even once the barrier was exceeded.
The terms X(z), where z = p r P BW or q r R BW , are Blatt-Weisskopf form factors [61], which act to cancel this behaviour once above the barrier. They are given by where z 0 represents the value of z when m = m 0 . The radius of the barrier, denoted r P BW or r R BW where the superscript indicates that the parameter is associated with either the parent or resonance in the decay chain, is usually taken to be 4.0 GeV −1 ≈ 0.8 fm [34]. Alternative descriptions of the angular distributions and Blatt-Weisskopf form factors are given in App. B.

Fit fractions
In the absence of any reconstruction effects, the DP PDF would be In a real experiment, the variation of the efficiency across the DP and the contamination from background processes must be taken into account; these details are discussed in Sections 3, 5 and 6. Typically, the primary results -i.e. the values obtained directly in the fit to data -of a DP analysis include the complex amplitude coefficients, given by c j in Eq. (3), that describe the relative contributions of each intermediate process. These results are dependent on a number of factors, including the amplitude formalism, choice of normalisation and phase convention used in each DP analysis. This makes it difficult to make useful comparisons between complex coefficients obtained from different analyses using different software. Fit fractions provide a convention-independent method to make meaningful comparisons of results from different fits. The fit fraction is defined as the integral of a single decay amplitude squared divided by that of the coherent matrix element squared for the complete DP, The sum of these fit fractions is not necessarily unity due to the potential presence of net constructive or destructive interference. Such effects can be quantified by defining interference fit fractions (for i < j only) as The interference fit fractions describe the net interference between the amplitudes of two intermediate processes.
Interference effects between different partial waves in a given two-body combination cancel when integrated over the helicity angle. Therefore, non-zero interference fit fractions should arise only between contributions in the same partial wave of one two-body combination, or between contributions in different two-body combinations. Large interference fit fractions, or equivalently a sum of fit fractions very different from unity, can often be an indication of inadequate modelling of the Dalitz plot.

Helicity angle convention
In the formalism just described, there is a choice as to which of the two resonance daughters the momentum q should be attributed (and hence to attribute the momentum − q to the other). This choice, although arbitrary, will affect the values of the measured phases and hence it is important that it is documented to allow comparisons between results. The convention used in Laura ++ is as follows: • θ 12 is defined as the angle between d 1 and d 3 in the rest frame of the d 1 and d 2 system, i.e. q is the momentum of d 1 ; • θ 23 is defined as the angle between d 3 and d 1 in the rest frame of the d 2 and d 3 system, i.e. q is the momentum of d 3 ; • θ 13 is defined as the angle between d 3 and d 2 in the rest frame of the d 1 and d 3 system, i.e. q is the momentum of d 3 .
This convention is illustrated in Fig. 2. One important point to note is that it is not a cyclic permutation. Rather it is designed such that for decays where d 1 and d 2 are identical particles the formalism is already symmetric under their exchange, as required. For decays of neutral particles to flavour-conjugate final states containing two charged daughters, e.g.

( )
B → π + π − K 0 S , there is a further complication that must be considered. In the example given, if one chooses for the B decay that d 1 would be π + , d 2 would be π − and d 3 would be K 0 S , one should then define the B decay using the conjugate particles, i.e. d 1 would be π − , d 2 would be π + and d 3 would be K 0 S . In practice, however, one often has an untagged data sample that contains both B and B decays that are not distinguished and so a single definition of the DP must be used. (Flavour-tagged analyses are discussed in Sec. 9.) Consequently, the amplitude model must account for the incorrect particle assignments for one of the flavours. Due to the choice of convention in Laura ++ this can be handled in a straightforward manner as long as the self-conjugate particle (the K 0 S in the example given) is assigned to be d 3 . Under these circumstances, the relation F (m 2 13 , m 2 23 ) = F (m 2 23 , m 2 13 ) can be restored simply by multiplying by −1 the cosine of the helicity angle θ 12 in the amplitude calculations for either the particle or antiparticle decay -we choose to do this for the particle decay (the B decay in the example given). So, in the considered example, a contribution from B 0 → ρ(770) 0 K 0 S would have its helicity angle definition inverted with respect to that of B 0 → ρ(770) 0 K 0 S . Figure 2: Values of the cosine of the helicity angles (left) θ 13 , (middle) θ 23 , and (right) θ 12 as functions of DP position. The kinematic boundary of this DP corresponds to that for the B 0 → π + π − π 0 decay.

Experimental effects
In order to extract physics results from reconstructed P → d 1 d 2 d 3 candidates using real experimental data, several effects need to be taken into account. One major concern will be that any backgrounds that fake the signature of signal decays need to be removed, which is usually done by imposing various selection criteria that exploit differences in the kinematics and topologies between signal and background events. The effect of applying selection criteria invariably means that the probability of reconstructing signal decays will not be 100% and furthermore could vary as a function of DP position. Along the boundaries, at least one of the daughter particles has low momentum, which typically reduces the reconstruction efficiency compared to decays at or near the DP centre. To account for this effect properly, the signal PDF, originally defined in Eq. (21), needs to be modified to where the signal efficiency (m 2 13 , m 2 23 ) is defined as the fraction of signal decays at the given DP position that are retained after all selection criteria have been applied.
For certain modes, there can be decay channels that can mimic the properties of the signal mode under study. For example, there may be significant backgrounds to the charmless decay Such backgrounds can be removed, or at least suppressed, by applying a "veto", which means excluding candidates that lie within, typically, three to five widths on either side of the mass peak. Alternatively, they can be accounted for within the signal model.
Another issue that needs to be considered is the effect of finite experimental resolution in the determination of the momentum of the parent P and its daughter particles. This leads to imperfect measurements of the invariant mass-squared combinations of the daughters, and also causes uncertainty on the invariant mass of the reconstructed P candidate. To avoid creating a DP with a fuzzy boundary, the mass of P can be fixed to its expected value, with adjustments made to the four-momenta of its daughter particles to ensure momentumenergy conservation. This will in general improve the resolution of the measurement of the DP co-ordinates, however there may still be significant effects related to events migrating from one region of the DP to another, especially near the corners of the kinematic boundary. This effect is usually ignored if the size of the migration/resolution is smaller than the width of the narrowest resonance under consideration, or if the largest migration probability is below 10% or so (although in the latter case, it is likely that systematic uncertainties on the physics results would need to be evaluated). If particularly narrow resonances contribute to the decay or if the final state particles under study suffer from significant misreconstruction effects (as is often the case for decays involving neutral pions), these effects can be taken into account rather generically by adding a "self cross-feed" component to the signal PDF. In this component, the true PDF is smeared by the resolution function w scf (s reco , s true ), given in terms of the reconstructed and true DP positions s reco ≡ (m 2 13 , m 2 23 ) and s true . The total signal PDF is then where for the first component the resolution is negligible (equivalent to w scf being a delta function), and the level of the second is determined by the self cross-feed fraction f scf (s true ), i.e. the fraction of reconstructed events with true DP position s true that are misreconstructed. The integral is over all true DP positions, although in practice only those with non-zero values of w scf for the given s reco need to be included. The fraction f scf can vary between 0 and 1, which correspond to the cases that resolution is negligible or that it must be considered for all signal events. The map w scf is sufficiently flexible to account for the fact that resolution may be more important to consider in some regions of the DP than others.
Despite imposing selection requirements in order to select signal candidates, there can remain significant fractions of various backgrounds in the DP analysis sample. This means that an extended likelihood function L needs to be employed in order to include these additional contributions: where N is equal to k N k , N k is the yield for the event category k (signal or background), N c is the total number of candidates in the data sample, and P j k is the PDF for the category k for event j, which consists of the product of the DP PDF and any other (uncorrelated) PDFs that are used to discriminate between signal and background. The function −2 ln L is minimised in an unbinned fit to the data in order to extract all of the parameters.
Amplitude analyses often feature multiple solutions, which are local minima of the negative log-likelihood function. For example, in the case that two broad overlapping resonances appear in the same partial wave, it may be possible to have solutions with either constructive or destructive interference with similar log-likelihood values. In order to find the true global minimum, the fit should be repeated many times with randomised initial values of the free parameters. The best solution can then be found by taking that with the minimum negative log-likelihood. In addition to providing confidence that the result obtained corresponds to the global minimum, this procedure is helpful to understand the sensitivity of the data to rejecting the secondary solutions.

Implementation of the signal component
In this section we begin to describe the code structure of the Laura ++ package by first outlining the classes and methods used to build up the total DP amplitude of the signal, given in Eq. (3). Furthermore, we describe how this is normalised in order to form the signal PDF defined in Eq. (21).

Particle definitions and kinematics
The most fundamental parts of the code define the properties of the parent particle P and its three daughters d 1 , d 2 and d 3 and their associated kinematics. Allowed types for P are The information on the decay that is to be modelled is encapsulated within the LauDaughters class, which is constructed by providing the names or PDG codes [62] of the parent and daughters. The particle properties are retrieved using the LauDatabasePDG singleton class, which extracts and supplements information from the Root TDatabasePDG particle property class.
The LauDaughters object then instantiates a LauKinematics instance, supplying to it the masses of the parent and its daughters. Instances of LauKinematics are used throughout the Laura ++ code to calculate and store all of the required kinematic variables for a given position in the DP (usually supplied as m 2 13 and m 2 23 ). These kinematic variables include the two-body invariant masses and helicity angles, the momenta of the daughters in the parent rest frame and in each of the two-body rest frames. In addition, there is the option to calculate the co-ordinates of the so-called "square Dalitz plot".

Square Dalitz plot
Since, particularly in B decays, signal events tend to populate regions close to the kinematic boundaries of the DP, it can be convenient to use a co-ordinate transformation into the so-called square Dalitz plot (SDP) [33]. The SDP is defined by variables m and θ that have validity ranges between 0 and 1 and are given by where m max 12 = m P − m d 3 and m min 12 = m d 1 + m d 2 are the kinematic limits of m 12 allowed in the P → d 1 d 2 d 3 decay, while θ 12 is the helicity angle between d 1 and d 3 in the d 1 d 2 rest frame, as explained in Sec. 2.4. Similar to how a choice of DP variables must be made, the SDP can be defined in several different ways. The expressions of Eq. (27) correspond to the choice used in Laura ++ , which must be employed consistently whenever a SDP is used.
To transform between DP and SDP representations, it is necessary to ensure correct normalisation. This is achieved by including the determinant of the Jacobian of the transformation, which is given by where p and q are evaluated in the d 1 d 2 rest frame and the partial derivatives evaluate to The SDP coordinate system has several advantages that apply whenever there is a need to bin the phase space, which are illustrated in Fig. 3. Firstly, the regions near the kinematic boundary are spread out, which means that these regions where the signal is often concentrated and where also there can be rapid variation in efficiency and background distributions can be treated with a much finer resolution, even while maintaining a uniform binning. Secondly, the kinematic boundary is perfectly aligned with the bin edges.  Figure 3: Illustration of the transformation between conventional and square Dalitz plot representations, for resonances in the B 0 s → D 0 π + K − decay (here the final-state particles are ordered following the d 1 d 2 d 3 convention of Laura ++ ). Compared to the same DP shown in Fig. 1, a fake D 0 π + resonance, with parameters corresponding to those of the D * 2 (2460) + state (blue points) has been added in order to better visualise the transformation in relevant DP regions.

Isobar dynamics and resonances
Once the kinematics of a particular decay mode have been established, the structure of the signal DP model can be by defined by creating a LauIsobarDynamics object. Components of the model are specified using the addResonance member function, which requires: • the name of the resonance, • an integer that specifies which of the daughters is the bachelor particle, and hence in which invariant mass spectrum this resonance will appear (1 for m 23 , 2 for m 13 , 3 for m 12 or 0 for some nonresonant models), • an enumeration to select the form of the dynamical amplitude.
Appendix C contains lists of the names of the allowed resonances along with their nominal mass, width, spin, charge and Blatt-Weisskopf barrier radius. This information is all automatically retrieved from LauResonanceInfo records that are stored in the LauResonanceMaker class. Appendix C also provides information on how to account for a state that is not already included in Laura ++ , and how to change the nominal values of the properties of any resonance. Appendix A gives details of all the dynamical amplitude forms that are currently implemented in the package and in Table A1 supplies the corresponding enumeration types. Examples of usage are given in Sec. 7.1. The signal model may also include contributions that do not interfere with the other resonances in the DP. These may arise due to decays that proceed via intermediate long-lived, i.e. negligible natural width, states; an example is the contribution from Experimentally, such components can be considered either as signal or background, and selection requirements (e.g., based on the consistency of the three daughter tracks originating from the same vertex position) may be used to suppress them, but in certain cases some contribution will remain. Within Laura ++ , the user can choose how to treat such contributions. When considered as part of the signal model, non-interfering components can be added with the addIncoherentResonance member function, which has the same number and type of arguments as the addResonance function. In this case the form of the dynamical amplitude should be specified as GaussIncoh, and the width should be changed to correspond to the experimental resolution. Note that the resolution for incoherent contributions is handled in a different way to the approach described in Sec. 3 and Sec. 5.2. As part of the signal model, a non-interfering component will contribute to the denominator of the fit fractions, but it is simple for the user to subtract it from the results since there is no interference with other components. When building the model, Laura ++ performs a simple check of charge conservation, while angular momentum is conserved by construction. However, the onus is on the user to make sure that the strong decays of resonances included in the model respect conservation of parity and flavour quantum numbers, since these are not checked by the code. In order to help with this, a summary of all resonances used in the model is printed out during the initialisation.
The complex dynamical amplitudes R(m) of the various resonance forms are defined using classes that inherit from the abstract base class LauAbsResonance. For example, the relativistic Breit-Wigner lineshape is defined within the LauRelBreitWignerRes class. All such classes implement the resAmp member function that returns a LauComplex class that represents the amplitude at the given value of the relevant two-body invariant mass. The LauAbsResonance base class implements the calculation of the angular distribution factor. Those amplitude forms that require the calculation of the Blatt-Weisskopf factors make use of the LauBlattWeisskopfFactor helper class.

Symmetry
If the decay of P contains two identical daughters, such as in B + → π + π + π − , then the DP will be symmetric. As mentioned in Sec. 2.4, the identical particles should be positioned as d 1 and d 2 . This situation is automatically detected by LauDaughters and the information propagated to the amplitude model. In this case, it is required only to define the resonances for the pair d 1 d 3 ; the amplitude is automatically symmetrised by LauIsobarDynamics by flipping the invariant-mass squared variables m 2 13 ↔ m 2 23 , recalculating the amplitude and summing.
When all of the daughters are identical, for example in

Normalisation of signal model
Various integrals of the dynamical amplitude across the DP need to be calculated in order to normalise the signal PDF given by Eq. (21), as well as for calculating the fit fractions for individual resonances defined in Eqs. (22) and (23). Since the c j coefficients are constant across the DP, only the amplitude terms F j need to be integrated. In general, these integrals cannot be found analytically, so Gauss-Legendre quadrature is used to evaluate them numerically. This is achieved by dividing the DP into an unequally spaced grid whose points correspond to the abscissa co-ordinates from the quadrature procedure. The granularity of the grid is chosen to ensure sufficiently precise integration, as discussed in more detail below. The F j terms are then multiplied by the quadrature weights and summed over all grid points that lie within the DP boundary. To remove the quadratic variable dependence, and to improve numerical precision, the DP area element dm 2 13 dm 2 23 is replaced by dm 13 dm 23 multiplied with the Jacobian factor 4 m 13 m 23 . This means that the normalisation of the total amplitude A is given by where w a (w b ) are the weights for the Gauss-Legendre quadrature abscissa values h a (h b ), which correspond to the grid points along the m 13 (m 23 ) axis, and the amplitude A is evaluated for all abscissas inside the DP kinematic boundary. An equivalent expression is used for normalisation of the experimental signal PDF of Eq. (24). The number of points N a (N b ) is set by dividing the m 13 (m 23 ) mass range by a default "bin width" δm of 5 MeV/c 2 , which can be changed using the setIntegralBinWidths function in LauIsobarDynamics, giving ∼ 1000 integration bins along each mass axis for B decays. It is important to realise that δm is not, in general, equal to the separation between neighbouring abscissa points. The LauIntegrals class handles the calculation of the general weights and abscissas for the integration range (−1, 1), while the LauDPPartialIntegralInfo class scales these using the half-ranges and mean values of m 13 and m 23 , following the numerical recipe given in Ref. [63]. This information is then used within the LauIsobarDynamics class to find the normalisation integral for the total amplitude, as well as the integrals for the fit fractions.
If the DP contains narrow resonances with widths below a threshold value, which defaults to 20 MeV/c 2 but can be changed using the setNarrowResonanceThreshold function in LauIsobarDynamics, then the quadrature grid is split up into smaller regions to ensure that the narrow lineshapes are integrated correctly. The range, in the invariant mass of the resonance daughters, for these sub-regions is taken to be m 0 ± 5Γ 0 , where m 0 and Γ 0 are the nominal mass and width of the resonance. The number of quadrature points (∼ 1000) along each sub-grid axis is set by dividing the resonance mass range (ensuring the limits stay within the DP boundary) by a δm value which defaults to 1% of Γ 0 and is configurable using the setIntegralBinningFactor function in LauIsobarDynamics. Grid regions that are outside the narrow resonance bands use the default bin width. When there are narrow resonances along the diagonal axis m 12 , the integration scheme switches to use the SDP defined in Section 4.1.1. The number of points on the integration grid defaults to 1000 for each of the m and θ axes. This can be tuned using the setIntegralBinWidths function.
The fact that resonance parameters can float in the fit means that the integrals will need to be recalculated if and when those values change. In order to minimise the amount of information that needs to be recalculated at each fit iteration, a caching and bookkeeping system is employed that stores the amplitudes of each component of the signal model as well as the Gauss-Legendre weights and the efficiency for every point on the integration grid. At each fit iteration it then checks which, if any, resonance parameter values have changed and with which amplitudes those parameters are associated. Only those affected amplitudes are recalculated (for each integration grid point and for each event in the data sample). While greatly improving the speed of fits, this comes at some cost in terms of memory usage, in particular if the integration grid is very fine. If the number of grid points is extremely large a warning message is printed, which recommends that the integration scheme be tuned using the setNarrowResonanceThreshold and setIntegralBinningFactor functions of LauIsobarDynamics to reduce the number of points.

Signal model and amplitude coefficients
Having defined the signal PDF for P → d 1 d 2 d 3 , as in Eq. (24), it is necessary to define the parameterisation of the complex coefficients c j defined in Eq. (3). Several different parametrisations have been used in the literature and are available in Laura ++ -a complete list is given in Table 1. These can be separated into two categories: cases in which it is assumed that there is no difference between the decay of P and its CP -conjugate P and cases in which such CP -violating differences are accommodated.
In the former case, the signal model is constructed by passing the corresponding LauIsobarDynamics instance to a LauSimpleFitModel object, which inherits from the abstract base class LauAbsFitModel. The fit model classes implement the functions needed to generate events according to the DP model as well as to perform fits to data. In the latter case, the decays of each of P and P should be represented by their own LauDaughters instance. These are used to construct two instances of LauIsobarDynamics, which in turn are used to construct an instance of the LauCPFitModel class that, like LauSimpleFitModel, also inherits from LauAbsFitModel.
The F j terms in Eq. (21) are calculated using the amplitudes that make up the LauIsobarDynamics model. The complex coefficients c j are each represented by an object inheriting from the LauAbsCoeffSet base class, which provides an abstract interface for combining a set of real parameters to form the complex c j . Each coefficient is constructed by providing the name of the resonance j and a series of parameter values that will be used to form the complex number c j . They are then applied to the model using the setAmpCoeffSet function of the LauSimpleFitModel or LauCPFitModel class. Checks are made to ensure that the coefficient name matches that of one of the components of the isobar model. The coefficient is then assigned to that component. As such the various coefficients are reordered to match the ordering in the LauIsobarDynamics model(s). Table 1: List of coefficient sets to represent c j in Eq. (3), separated into cases where CP conservation is assumed and those where CP violation is accommodated in the model. Where parameters are preceded by ± signs in the expressions for c j , the + (−) sign corresponds to the usage for P (P ) decays. The corresponding class for each set is LauLabel CoeffSet, where Label is the set label given below.

Set label
Parameters

Implementation of efficiency and resolution effects
As discussed in Sec. 3, in an experimental analysis it is usually necessary to modify the signal PDF in order to account for effects such as the variation of the reconstruction and selection efficiency over the DP and detector resolution or misreconstruction. In this section we describe the classes and methods in the Laura ++ package that are used to implement these modifications to the pure physics PDF described previously. All experimental effects are handled within Laura ++ through histograms. Optional spline interpolation is available to smooth effects due to the finite size of samples used to obtain the histograms. In some analyses it may be possible to obtain reliable parametric descriptions of, for example, the efficiency variation over the DP. Where this is a possible and desirable approach the user can simply generate a histogram from the parametric shape and use that as input. The granularity of the histogram can be as fine as necessary to describe local variations; often the limiting factor on the binning will be the size of the sample used to obtain the shape.

Efficiency
The variation of the signal efficiency over the DP, (m 2 13 , m 2 23 ), is implemented by the LauEffModel class. Its constructor requires a LauDaughters object, which defines the kinematic boundary, as well as a LauVetoes object, which is used to specify any region in the DP that has been excluded from the analysis (perhaps to remove particular sources of background or to exclude a region of phase space where the efficiency variation is poorly understood). The signal efficiency is set to zero inside a vetoed region; the resulting discontinuity at the boundary motivates different treatment of vetoes to other sources of inefficiency that vary smoothly across the DP. Vetoes can be added using the addMassVeto or addMassSqVeto functions of the LauVetoes class, which require the bachelor daughter index as well as the lower and upper invariant mass (or mass-squared) values for each exclusion region. Since version v3r2 of Laura ++ , vetoes are automatically symmetrised as appropriate for DPs containing two or three identical particles in the final state.
All other information on the efficiency variation over the DP needs to be supplied in the form of a uniformly binned two-dimensional Root histogram. Alternatively, a set of histograms can be provided; in this case the total efficiency is obtained by multiplying the efficiencies at the appropriate position in phase space from each of the components. This provides a convenient way to assess the impact of systematic uncertainties from different contributions to the total efficiency. Each of these component efficiency histograms can have different binning.
Efficiency histograms will usually be constructed by applying all selection requirements to a simulated sample of signal decays that has been passed through a full detector simulation. A ratio is then formed of all decays that survive the reconstruction and selection to all those that were originally generated. Since the effect of explicit vetoes in the phase space is separately accounted for, it is advised that these are not applied when constructing the numerator of the efficiency histogram.
As previously mentioned in Sec. 4.1.1, signal events often occupy regions close to the kinematic boundaries. The SDP transformation defined in Eq. (27) can be used to spread out these regions so that the efficiency variation can be modelled more accurately. As such, the LauEffModel class will accept histograms that have been created in either m 2 13 -m 2 23 or m -θ space. The histogram can then be supplied to the LauEffModel class via the setEffHisto or setEffSpline function as appropriate. Where the total efficiency is to be obtained from the product of several components, further histograms can be included with the addEffHisto and addEffSpline functions. In each case, a boolean argument squareDP is used to indicate the space in which the histogram has been defined. For symmetric DPs, there is also the option to specify that the histogram provided has already been folded and hence only occupies the upper half of the full DP (or the corresponding lower half of the SDP). Internally, each histogram is stored as a Lau2DHistDP or Lau2DSplineDP object, which implement (optional) bilinear or cubic spline interpolation methods, respectively.
Functionality is also available to help estimate systematic uncertainties due to imperfect knowledge of the efficiency variation by creating Gaussian fluctuations in the bin entries. The fluctuations are based on the uncertainties provided by the user, which, depending on the function used to provide the histogram, can be asymmetric. It is up to the user as to whether the provided uncertainties are simply due to the limited size of the simulated sample or whether they also account for effects such as possible disagreements between data and simulation. The fluctuations are activated by providing optional boolean arguments to the function where the histogram is provided.

Resolution
The self cross-feed component of the signal likelihood, modelled as described in Eq. (25), is implemented using information from the LauScfMap class, which stores all possible values of w scf via its setHistos function. As for the description of the efficiency, the histograms used to describe w scf and f scf can be constructed in either m 2 13 -m 2 23 or m -θ space. However, to simplify the implementation in Laura ++ it is currently required that all histograms related to the description of resolution have the same binning.
The implementation proceeds as follows. Consider a SDP describing the true position, s true , divided into a uniformly binned two-dimensional histogram. Each bin will have associated with it another two-dimensional histogram, with identical binning, whose entries give the migration probability w scf (s reco , s true ) of the true position (given by the original bin centre) being reconstructed in the bin that contains s reco . These can be constructed as follows: • Each histogram contains only the events that were generated in a given true bin.
• The events are plotted at their reconstructed co-ordinates.
• The histogram is then normalised.
Some histograms may be empty if there were no events generated in that bin, although this is of course dependent on the size of the samples used and the size of the bins. The order of the histograms in the vector supplied to the setHistos function should be in terms of the Root "global bin number".
The splitSignalComponent method of LauSimpleFitModel and LauCPFitModel takes a LauScfMap object as an argument, as well as a two-dimensional Root histogram whose entries give f scf (s true ). The fit models then use this information to evaluate the self cross-feed contribution to the likelihood. In the ideal case of Eq. (25) this is described by an integral; in practice this becomes a summation, where the hat (ˆ) notation is used to indicate quantities evaluated in the SDP (as is the case in this example), and the bracket ( ) notes quantities obtained from histograms. The pure signal PDF P sig is as in Eq. (24), evaluated at the position corresponding to the SDP point at the centre of theŝ true i bin. The phase space factors ∆Ω are equal to ∆m ∆θ |J| where the SDP bin size is given by ∆m ∆θ and the Jacobian of the SDP transformation is that of Eq. (28); since equal binning is required, the ratio of phase space factors reduces to a ratio of Jacobians.
6 Implementation of background components 6

.1 Dalitz-plot distributions
Backgrounds in the data sample can be taken into account by including them in the total likelihood defined in Eq. (26). This means that the DP distributions of all background categories need to be provided. In an analogous way to the implementation of the signal efficiency described in the previous section, the DP distribution of each background category is represented with a uniformly binned two-dimensional Root histogram. Background contributions are handled in Laura ++ as follows. Firstly, the names of all background categories need to be provided to the fit model via the setBkgndClassNames function of the appropriate LauAbsFitModel class. Then, each named category needs to have its DP distribution defined and added to the fit model. This is achieved by supplying one (or more if the model subdivides the data to account for effects such as CP violation) LauBkgndDPModel object. The constructor of the LauBkgndDPModel class requires pointers to the usual LauDaughters and LauVetoes objects. The histogram is supplied to the LauBkgndDPModel object via its setBkgndHisto (setBkgndSpline) function, and is then internally stored as a Lau2DHistDPPdf (Lau2DSplineDPPdf) object that implements bilinear (cubic spline) interpolation. The PDF value is then calculated as the interpolated number of background events B(m 2 13 , m 2 23 ) divided by the total integrated area of the histogram: Like signal events, backgrounds also tend to populate regions close to the kinematic boundaries of the DP. Therefore, the use of histograms in the SDP space can improve the modelling of backgrounds. This is achieved by providing a histogram in m -θ space and setting the squareDP boolean flag to true in the setBkgndHisto or setBkgndSpline functions of LauBkgndDPModel. The normalisation of these PDFs then automatically includes the Jacobian for transforming from normal to square Dalitz-plot space.
Some special treatment is necessary for backgrounds modelled from sources that contain contributions that are vetoed in the DP fit. Following the example of Sec. 3, the combinatorial background to B − → K − π + π − decays may be modelled from a sideband in the B candidate mass distribution that also contains some genuine D 0 → K − π + decays. The histogram binning will introduce some smearing of such contributions, so that once the veto is applied later some residual background may remain. In principle this can be avoided with sufficiently fine histogram bins, but this will often be impractical due to finite sample sizes. Instead, and in contrast to the procedure for efficiency histograms described in Sec. 5.1, the veto should be applied when the background histogram is made. This will, however, lead to an underestimation of the background that is being modelled (in the example above, of the random combinations of three tracks) within bins that lie partially inside the vetoed regions. To correct for this effect, each background histogram needs to be divided by another histogram (with the same binning) whose entries contain the fraction of events, generated from a high statistics sample that is uniform in phase space, that lie outside any veto region; bins that are completely outside (inside) a veto have a weight of unity (zero), and division by zero is interpreted as zero weight. The histogram supplied to LauBkgndDPModel should already have had this correction applied to it.

Other discriminating variables
When fitting a data sample that contains (significant) backgrounds, additional discriminating variables can be included in the total likelihood function defined in Eq. (26), in order to provide improved separation between signal and background categories. Examples of such discriminating variables include the mass of the parent particle P candidate and the output of a multivariate discriminant to separate signal from background. Assuming that all of the variables x = (x 1 , x 2 , ..., x n ) are uncorrelated, the PDF P j k of the signal or background category k, for event j, is given by the product of the individual variable PDFs (including that of the DP distribution): . Additional PDFs are represented by classes that inherit from LauAbsPdf, whose constructor requires the variable name, a vector of the PDF parameters, as well as minimum and maximum abscissas that specify the variable range. Each class implements the member function for evaluating the PDF function at a given abscissa value and that for evaluating the maximum value of the PDF in the fitted range. A list of classes that can be used for additional PDFs is given in App. D. Each PDF needs to be added to the fit model; for LauSimpleFitModel (LauCPFitModel), signal PDFs are added using the setSignalPdf (setSignalPdfs) function, self cross-feed PDFs are included via setSCFPdf (setSCFPdfs) function, and background PDFs are added with the setBkgndPdf (setBkgndPdfs) function.
An alternative approach is to perform first a fit to some or all of the other discriminating variables to determine the yields of signal and background. These values can then be fixed in the fit to the DP to extract the amplitude parameters. Requirements can also be imposed on those variables, before or after such a fit, in order to enhance the signal purity among the events entering the DP fit. This approach will have slightly reduced sensitivity but avoids the need to model correlations between the DP and the other variables (although classes are provided that allow the modelling of such correlations, see App. D.13). It is possible to fit only other discriminating variables, i.e. excluding the DP variables, in Laura ++ through the useDP function of LauAbsFitModel, which takes a boolean argument.

Work flows
The Laura ++ package is designed to perform three main work flows: • Monte Carlo generation of toy datasets from a defined PDF, • fitting a defined PDF to an input data set, • calculating per-event weights for an input data sample (usually full detector simulation) from a defined DP model.
In this section we describe the setup stages that are common to all tasks, then each of the three work flows in turn, and finally some variations that are available. Code taken from the various examples included with the package is used to illustrate many of these points. These code snippets are in C ++ but it is also possible to use the python bindings that are generated by rootcling when the Laura ++ library is compiled and to write the application code in python. An example that demonstrates how to do this, GenFit3pi.py, is included in the package since version v3r3.

Common setup
The first setup task is to define the DP with which we are working. As mentioned in Section 4.1 this is achieved through the declaration of a LauDaughters object: bool squareDP = true ; LauDaughters * daughters = new LauDaughters ( " B + " , " pi + " , " pi + " , " pi -" , squareDP ) ; specifying the types of the parent particle and the three daughters. The final argument specifies whether or not we will be using the variables of the SDP (defined in Section 5) and as such whether these should be calculated -if they are not required then it is more efficient to switch off this calculation. Next to be defined is the efficiency model, including any explicit vetoes in the Dalitz plane: LauVetoes * vetoes = new LauVetoes () ; LauEffModel * effModel = new LauEffModel ( daughters , vetoes ) ; Without further specification, this will give a uniform efficiency (of unity) over the DP. This can be modified by specifying vetoes and/or supplying the efficiency variation in the form of a histogram, e.g.
vetoes -> addMassVeto (2 , 1.7 , 1.9) ; // D0 veto , covering 1.7 < m13 < 1.9 GeV / c^2 TH2 * effHist = ...; effModel -> setEffHisto ( effHist ) ; The next step is to define the signal amplitude model. Before defining any resonances the first thing to be done is to specify the Blatt-Weisskopf barrier radii to be used for the different resonances. This is done by accessing the singleton LauResonanceMaker factory object: reson -> s e t R e s o n a n c e P a r a m e t e r ( " g1 " ,0.2) ; reson -> s e t R e s o n a n c e P a r a m e t e r ( " g2 " ,1.0) ; reson = sigModel -> addResonance ( " f_2 (1270) " , 1 , La u Ab sR es o na nc e :: RelBW ) ; reson = sigModel -> addResonance ( " rho0 (1450) " , 1 , La uA bs R es on an c e :: RelBW ) ; reson = sigModel -> addResonance ( " NonReson " , 0 , L au Ab s Re so na n ce :: FlatNR ) ; Note how the returned pointer to the added resonance can be used to modify parameters of the resonance. It is also possible to specify that resonance parameters should be floated in the fit using the LauAbsResonance::floatResonanceParameter member function. The fit model can now be created based on the isobar model that has just been defined. And the complex coefficients for each resonance defined using one of the various parametrisations defined in Table 1, in this case using a simple magnitude and phase form: The boolean arguments to the constructor indicate whether the value should be fixed in the fit. Since Laura ++ version v3r1, the ordering of the components is defined by the isobar model (first all coherent resonances in order of addition to the isobar model, then all incoherent resonances in order of addition), and the coefficients are automatically rearranged to match that order. Thus, the user does not need to worry about the order of the components given here.
To either generate toy datasets or to fit data, it is necessary to specify the number of events in each of the signal and (if appropriate) background categories. The Laura ++ work flow requires that these values are specified also in the case of weighting data, although they are not used. In the case of fitting data the given values may be floated, as is often appropriate when additional discriminating variables are included in the fit (see Sec. 6.2), or fixed, as is usually necessary when only the DP is being fitted. In the former case, if the doEMLFit option is specified the fitted values of the yields can be expected to have the correct Poisson uncertainties. The user should ensure that the specified number of events in the signal and background categories correspond to the number of events in the data file that will be read in to that fit. The number of signal events is set as follows: LauParameter * nSigEvents = new LauParameter ( " nSigEvents " ,500.0 , -1000.0 ,1000.0 , false ) ; fitModel -> setNSigEvents ( nSigEvents ) ; The number of experiments, as well as which experiment to start with, must also be specified: const int nExpt (1) ; const int firstExpt (0) ; fitModel -> setNExpts ( nExpt , firstExpt ) ; The above example values are appropriate when fitting a single data sample. If generating and/or fitting toy pseudoexperiments then a larger number of experiments can be specified.
Optionally, models for one or more background categories can also be specified along with the corresponding event yields.
Finally there are various control and configuration options that can be set, which are listed in Table 2. Once all the configuration is completed, the required operation (toy generation, fitting, or event weighting) can be run: where all the arguments are strings, with the following purposes • command is either "gen", "fit" or "weight"; • dataFile and treeName specify the name of the Root file and tree from which the data should be read or to which the generated data should be written (depending on the mode of operation); • rootFileName specifies the name of the Root file to which the results of the fit (fit convergence status, likelihood, parameter values, uncertainties, correlations, etc.) should be saved; • tableFileName specifies the name of the L A T E X file to which the results of the fit or generation should (optionally) be written.
The run function performs the initialisation of all the PDFs that have been defined and checks the internal consistency of the model, e.g. that all event categories have a PDF defined for each variable. All the parameters of the model are gathered into a single list and the bookkeeping of those parameters that affect the DP normalisation is set up.
Any constraints on the model parameters, both on individual and on particular specified combinations of parameters, are also recorded. Following these initialisation steps, the particular routines for toy generation, fitting, or weighting of events are called, depending on the specified command. Each of these routines are described in the following sections.

Toy generation
The generation of ensembles of simplified pseudoexperiments, or toys, is of great importance in Laura ++ . Not only is toy generation used to test the performance of a fitter, to test its stability and study potential biases and/or correlations between fitted parameters, it can also be used to determine uncertainties on parameters (see Sec. 7.3.1). Large samples of toy experiments, generated according to the results of a fit, can also be used to display the results of the fit projected onto any variable in any region of the phase space. This is often a convenient way of drawing the result of the fit, including contributions from separate components. The first action specific to the generation of toy pseudoexperiments is to ensure that the value of each fit parameter is set to that to be used for generation, i.e. the genValue of the LauParameter. Then an ntuple (specifically, a Root TTree) is created in which the generated events will be stored. The branches to be stored are determined based on the PDFs that have been configured in the initialisation step, and will include: • iExpt and iEvtWithinExpt: these indicate to which pseudoexperiment this event belongs and the position of the event within that experiment; • evtWeight: a weight for the event that is normally unity but can take other values 4 ; • a set of MC truth branches that indicate whether a given event was generated as signal or one of the background categories, e.g. genSig; • efficiency: the value of the signal efficiency at the generated DP position (only stored for signal events); • the DP position in terms of the invariant masses and helicity angle cosines: m12, m13, m23, m12Sq, m13Sq, m23Sq, cosHel12, cosHel13, cosHel23; • optionally, the position in the SDP coordinates: mPrime and thPrime; • the values of all other variables used by additional PDFs (see Section 6.2).
Each experiment is then generated in turn, with the following procedure being followed for each: 1. Determine the number of events to generate for each category. This will either be the value of the yield parameter for the category or, if doPoissonSmearing is specified (see Table 2), will be sampled from a Poisson distribution whose mean is equal to the yield parameter value. Calculate total intensity and multiply by efficiency: 2. For each category, generate the number of requested events. Truth information is stored as to which category the event belongs in branches named "genCAT", where "CAT" is replaced by the name of the category. For a given event, the variable corresponding to the category being generated will have value 1, while all others will have value 0.
3. For signal events, the generation of the DP position is performed by the LauIsobarDynamics::generate function, the operation of which is shown in Fig. 4. In essence, it is an "accept/reject" method in which the value of the signal model intensity at a given point (I tot in Fig. 4) is compared with a "ceiling value" (I max in Fig. 4). Note that the signal model intensity includes effects due to efficiency, vetoes and resolution. Since, in general, it is not computationally efficient to evaluate analytically the maximum value of the intensity, an approximate value is provided by the user via the LauIsobarDynamics::setASqMaxValue function. During the generation process it is checked that a larger value is not encountered. If it is, which would indicate that the generation is biased, the ceiling value is increased and the generation processes starts again from the beginning -all previously generated experiments are discarded. Similarly, too large a ceiling value can make the generation extremely inefficient, so it is also attempted to detect this condition and correct it. If a self cross-feed component has been included in the signal description, once a DP point has been generated it may then be smeared according to the procedure outlined in Section 5.2. Finally, the values of any other discriminating variables can be generated from the corresponding PDFs.
4. For background events, the DP position is generated by the LauBkgndDPModel::generate function.
This performs a similar accept/reject routine to that used to generate signal events but here the (optionally interpolated) histogram height is compared with the maximum value of the histogram. There is therefore no need to verify that the encountered value is less than the maximum since this is true by construction. Generated values within a veto region are automatically rejected. Finally, the values of any other discriminating variables are generated from the corresponding PDFs.
Once all experiments have been generated successfully, all events are written to the ntuple.

Embedding events from a file
One possible variation of the scheme just outlined is that events for one or more components of the model are loaded from a Root file rather than being generated from the model PDFs. These events could, for example, be taken from samples of full detector simulation or from data control samples. This procedure can be useful to check for biases that can result due to certain experimental effects not being taken into account in the model (for example, experimental resolution or correlations between fit variables). The events are sampled randomly from the provided file. Depending on the number of events within the samples it may be necessary to sample with replacement, i.e. to allow events to be used more than once.
The exact procedure to enable the embedding of events from a Root file varies slightly depending on the fit model class being used. For LauSimpleFitModel there are two functions embedSignal and embedBkgnd, which can be used to embed signal and background events. For LauCPFitModel there are four functions, to allow separate sources of events in the case that the parent candidate is the particle or antiparticle. The arguments taken by the various functions are, however, essentially the same: • bkgndClass: in the case of the functions for embedding background events, the name of the background class must be specified as the first argument; • fileName: the name of the Root file from which the events should be loaded; • treeName: the name of the TTree object within the file in which the events reside; • the arguments reuseEventsWithinEnsemble and reuseEventsWithinExperiment control whether or not the sampling of events should allow replacement, firstly within the context of the entire ensemble of experiments being generated (i.e. replacement occurs once each experiment has been generated) and secondly within each experiment (i.e. replacement occurs immediately); • useReweighting: in the case of embedding signal events, this argument controls whether an accept/reject routine should be used (based on the configured amplitude model) when sampling the events. As discussed in Sec. 7.4, reweighting allows an existing sample to be reused with an alternative model, which can avoid high computational costs associated with obtaining a new sample. This option requires that the both the generated and reconstructed DP coordinates are available in the provided sample since the generated point is used to find the amplitude value to be used in the accept/reject routine but the reconstructed point is the one saved to the output ntuple. This ensures that the effects of experimental resolution and misreconstruction are preserved.
For example, to embed events for the signal component, sampling with immediate replacement, and performing the accept/reject routine: fitModel -> embedSignal ( signalSampleFileName , signalSampleTreeName , kTRUE , kTRUE , kTRUE ) ; As for other control and configuration functions, such as those in Table 2, these should be called prior to the invocation of the LauAbsFitModel::run function.

Fitting
While the primary goal of the Laura ++ package is to facilitate fits to the DP distributions of various decays of interest in data, it is also essential to be able to fit simulated samples such as ensembles of toy experiments generated as discussed above. The same work flow is used for both types of fits. Unbinned maximum likelihood fits are the default in Laura ++ , with extended maximum likelihood fits also an option as noted in Table 2.
The first action specific to the routine for fitting is to create the ntuple in which to store the results of the fit. The branches to be stored are determined based on the PDFs that have been configured in the initialisation step, and will include: • iExpt: an integer that indicates to which pseudoexperiment the results in this entry in the ntuple belong (zero for a fit to data); • fitStatus: an integer that indicates whether or not the fit has converged and the degree of accuracy of the resulting covariance matrix (0: fit has failed to converge; 1: fit has converged but covariance matrix is only approximate; 2: the covariance matrix is a full matrix but is forced to be positive definite; 3: the covariance matrix is full and accurate); • EDM: the estimated distance of the negative log likelihood to the true minimum of the function; • NLL: the minimised value of the negative log likelihood; • for each fit parameter: its fitted value, its value at initialisation (for toy data, this will be the true value) and, if it was floated in the fit, its uncertainty (including asymmetric uncertainties, if calculated), pull and global correlation coefficient; • for each pair of floated fit parameters: their correlation coefficient; • any extra parameters defined by the fit model, for example, the fit fractions of each component and the interference fit fractions.
If information on the per-event likelihood values is to be stored in order to allow the calculation of sPlot weights [64], this ntuple is also created. The file containing the data to be fitted is opened and the ntuple retrieved. An error will be returned if the file cannot be opened, if the ntuple cannot be found or if it does not have a flat format, i.e. if it contains stored arrays. Each experiment is then fitted in turn, with the following procedure being followed for each: 1. The data for the given experiment is read into memory. The fit model then passes the data to each of the PDFs such that they can store the values of any variables that they require and also calculate and cache as much information as they can towards (and possibly including) the likelihood value for each event.
2. Optionally (see according to the specification of useRandomInitFitPars; see Table 2), the initial values of the isobar coefficient parameters are randomised.
3. The fitter is initialised by providing it with the list of fit parameters and two boolean options to control whether a two-stage fit (see Sec. 7.3.2) is to be performed and whether asymmetric uncertainties on the fit parameters should be determined.
4. The minimisation of the negative log likelihood is then performed and the uncertainties on the parameters and their correlations are determined. The fit status information and the covariance matrix is stored, and the final values and uncertainties are written back to the fit parameter objects.
5. Some final manipulation of the fit parameter values is performed, e.g. such that phases lie within a particular range, and their pull values are calculated. Quantities derived from the fitted parameters, such as fit fractions, are also calculated. All this information is written to the fit results ntuple.
6. Optionally (see writeSPlotData in Table 2), per-event likelihood information is stored in a separate ntuple.
7. Optionally (see compareFitData in Table 2), and only if the fit was successful, toy MC pseudoexperiments are generated based on the fitted values of the parameters.
Once all experiments have been fitted, information is printed on the success rate of the fits and the fit results ntuple is written to file. Optionally (see writeSPlotData in Table 2), per-event weights are calculated from the likelihood information, using the sPlot method [64], and written to an ntuple. Due to the often complicated dependence of the likelihood on many parameters in a typical DP analysis, generally there will appear local minima in the likelihood as well as the global minimum. Depending on the starting values of the fit parameters it is possible that the fit will converge to one of the local minima. As such it is highly advisable to perform multiple fits to each data sample using randomised starting values for the parameters of the isobar coefficients. As mentioned in Table 2, the control function useRandomInitFitPars can be used to toggle this behaviour. It is then required only to re-run the executable to obtain a fit based on a new randomised starting point. The applications in the examples directory of the package all include a compulsory command-line option to provide an integer to label the particular fit (or set of fits, if fitting multiple experiments), which is then incorporated into the name of the file containing the ntuple of results for that fit or set of fits. Once all fits have been performed, the fit that returns the best likelihood value can then be selected as the likely global minimum for each experiment. In the examples directory of the package a utility class ResultsExtractor (with accompanying application code) is provided to simplify the process of checking for multiple solutions and extracting the results of the fit that is the candidate global minimum for each experiment. This application writes out a Root file that contains an ntuple that is filled with the candidate global minimum results for each experiment, and a histogram for each experiment that shows the frequency with which each negative log likelihood value is obtained.

Determination of uncertainties
In addition to central values, the fit will return estimates of the uncertainties of each fitted parameter. However, it is typical to want to obtain results not only for the fitted parameters (e.g. complex coefficients for each component of the model) but also for derived parameters (e.g. fit fractions and interference fit fractions). The central values of such derived parameters can be obtained algebraically, however non-linear effects typically render unreliable the determination of uncertainties by standard error propagation. Instead, both central values and uncertainties can be obtained from ensembles of pseudoexperiments generated according to the result of the fit to data. In this approach the uncertainty is obtained from the spread of the distribution of the obtained values in the pseudoexperiments of the parameter of interest. Since this procedure guarantees correct coverage of the uncertainties, and also allows study of effects such as biases and non-linear correlations in the results, it is often preferable to use it also for directly fitted parameters. A further alternative to determining uncertainties on derived parameters is by bootstrapping [65,66].

Two-stage fit
Fits with large numbers of free parameters can be slow to converge, and may be unstable. These problems can be ameliorated with a two-stage fit, in which certain parameters are initially fixed to specified values while all other parameters are floated. In the second stage, all parameters are floated ensuring that correlations between parameters are correctly accounted for. Note that parameters which are floated only in the second stage cannot have randomised initial values.
One situation in which this fit procedure has been found to be particularly useful is when resonance parameters such as masses and widths are to be floated. Another is when CP -violation parameters are to be determined. It may be noted that in the CartesianCP coefficient set, one can simultaneously swap x ↔ δ x , y ↔ δ y with an effect that is equivalent to rotating all amplitudes for P decays by π. Such a transformed set of parameters may be hard for the fit to distinguish from the untransformed case. As such a two-stage fit, in which CP -violation parameters are initially fixed to zero before being floated in the second stage, can help to resolve ambiguities and reject unrealistic multiple solutions.

Relations between fit parameters
There are a number of situations where a given fit parameter can appear in more than one PDF within the likelihood function. Such cases must be flagged in the construction of the fit model in order to ensure that such a parameter is provided only once to the minimiser. Consider the following example that concerns a fit to the invariant mass distribution of candidate B + → K + π + π − decays. The signal is modelled as a Gaussian function, which has two parameters: the mean µ (corresponding to the mass of the B + ) and width σ (corresponding to the experimental resolution). A possible background arises from the decay B + → η K + ; η → π + π − γ, which is modelled using an ARGUS threshold function (see App. D.1). The threshold parameter of this function is the B + mass and as such it needs to be encoded in the likelihood function that values of the mean of the signal distribution and the threshold of this background are due to the same parameter. This can be achieved in Laura ++ using the LauParameter::createClone function, as follows: const Double_t mbMin (5.10) ; const Double_t mbMax (5.60) ; LauParameter * sig_mb_mean = new LauParameter ( " sig_mb_mean " , 5.28 , 5.26 , 5.30) ; LauParameter * sig_mb_sigma = new LauParameter ( " sig_mb_sigma " , 0.20 , 0.10 , 0.30) ; std :: vector < LauAbsRValue * > mbPars ; mbPars . push_back ( sig_mb_mean ) ; mbPars . push_back ( sig_mb_sigma ) ; LauAbsPdf * sig_mb_pdf = new LauGaussPdf ( " mB " , mbPars , mbMin , mbMax ) ; LauParameter * bkg1_mb_m0 = sig_mb_mean -> createClone () ; LauParameter * bkg1_mb_xi = new LauParameter ( " bkg1_mb_xi " , 20.0 , 0.0 , 50.0) ; mbPars . clear () ; mbPars . push_back ( bkg1_mb_m0 ) ; mbPars . push_back ( bkg1_mb_xi ) ; LauAbsPdf * bkg1_mb_pdf = new LauArgusPdf ( " mB " , mbPars , mbMin , mbMax ) ; For the isobar coefficients, it is possible to clone the parameters related to a particular coefficient using the LauAbsCoeffSet::createClone function. One can specify precisely which parameters should be cloned; for example, to allow only CP -violating parameters to be cloned while the CP -conserving parameters are still free to vary independently.
Extending the previous example to consider a second possible background from B + → K + π − π + π 0 decays, where the π 0 is not reconstructed, we encounter a case where the threshold parameter is now the difference between the B + and π 0 masses. This scenario can be accommodated in Laura ++ using a LauFormulaPar object, the constructor of which takes as arguments two strings, which are the name of the compound parameter and the formula to be used to combine the values of the input parameters, and a std::vector of LauParameter objects, which are the input parameters. For the example in question this would be: Note the syntax of the formula expression (from the Root TFormula class) where the input parameters are referred to by their index in the std::vector. This syntax allows the usual arithmetic operators to be used as well as functions such as those contained in the TMath namespace.

Blinding of fit parameters
In analyses of data when one is searching for a new signal or searching for CP violation it is quite a common practice to "blind" the value of the observables of interest until the event selection and analysis procedures have been finalised [67,68]. This is commonly achieved by applying an offset, whose value is unknown to the analysts but is determined uniquely from a so-called "blinding string", when the parameter value is passed to the minimiser or is printed out or otherwise saved at the end of the fit but the true, "unblind", value is used when calculating the value of the likelihood function. This same technique is used in Laura ++ in the LauBlind class. The hash of the blinding string is used as the seed for a random number generator, using which a value is sampled from a normal distribution. The offset is then calculated by multiplying this random number by a scaling factor. The blinding string and scaling factor are supplied by the user as follows: LauParameter * nSig = new LauParameter ( " signalEvents " , 1000.0 , -1000.0 , 2000.0 , kFALSE ) ; nSig -> blindPa rameter ( " dalitzplot " ,1000.0) ; Some care must be taken when defining the scaling factor and the allowed range for the parameter to try and avoid situations where the parameter hits one of the limits. It is also advisable to use different blinding strings for each observable that is to be blinded.

Fitting with external constraints
When the value of a parameter is known from other measurements to within some precision it can be useful to incorporate that information into the likelihood function. This is most commonly achieved by multiplying the likelihood by a Gaussian term, where the mean and width of the Gaussian are the central value and uncertainty from the external source and the abscissa is the value of the parameter in the fit. The addition of such terms to the likelihood expression can be carried out in Laura ++ using the function LauParameter::addGaussianConstraint, where the arguments are the central value and uncertainty from the external measurement.
In the event that the constraint that needs to be applied is not just on the value of a single parameter but on some combination of fit parameters, this can be achieved via the addConstraint function that is available in all fit model classes. In addition to the mean and width of the Gaussian function, this function takes the following arguments: • a formula string that specifies how to combine the parameters, which should use the syntax specified by LauFormulaPar (see Sec. 7.3.3); • a std::vector of strings that are the names of the fit parameters (i.e. the LauParameter objects in the list of fit parameters) whose values are to be used in the formula to calculate the abscissa of the Gaussian function. Note that the order of the names within the vector must match the ordering specified in the formula.
This feature can also be used to perform a scan of the likelihood as a function of any combination of fit parameters. This can be achieved by setting the mean of the Gaussian constraint to one of a range of fixed values, in turn, and setting the width to a sufficiently small value that the quantity is effectively fixed. The negative log-likelihood values obtained from independent fits at each of the scan points can be converted into whatever format is most convenient for the user.

Fitting with background subtraction via sPlot weights
It is possible to perform fits with background subtraction via sPlot weights as proposed in Ref. [69]. In this approach, the weights are obtained from a fit to a discriminating variable in which the signal and any background categories can be distinguished [64]; for the purposes of this discussion this will be considered to be the mass of a B meson candidate for a particular decay to a three-body final state. The advantage of this approach is that it becomes unnecessary to have models for the DP distributions of the background components, and therefore the discussion of Sec. 6 becomes irrelevant. However, the formalism relies on the discriminating variable being uncorrelated with the DP variables. This is likely to be a good approximation for the signal, and also reasonable for combinatorial background. However, additional sources of background such as those involving misidentified decays are likely to have significant correlations between the reconstructed mass and the DP position. Therefore this approach is only valid in the case that such backgrounds are negligible across the whole B candidate mass range -not only the signal region (in this approach to background subtraction, there is no concept of signal region).
Within Laura ++ this approach to fitting can be implemented by calling the doSFit member function of LauAbsFitModel, which takes as argument the name of the branch that contains the relevant weights, as shown in Table 2. It is also possible to pass an optional second argument, which is a scaling factor that is needed in order to obtain correct uncertainties from the likelihood; this is unnecessary if the uncertainties will be evaluated from pseudoexperiments.

Simultaneous fitting of independent data samples
Simultaneous fits to independent data samples have a variety of applications. For example, the samples may correspond to different experimental conditions and so each may require a different model of the efficiency and background distributions. A simultaneous fit to the various sub-samples allows the statistical power of the full sample to be used to determine the parameters of interest, while accounting as accurately as possible for the changes in the experimental environment. Furthermore, it simplifies the evaluation of systematic uncertainties, in particular the treatment of sources that are correlated among the samples. The technique can also be employed in order to extract information that would otherwise not be feasible, for example to perform a coupled-channel analysis or to exploit flavour symmetries or other relations between decay modes, in order to extract information on CP violation observables, as recently carried out in Ref. [70].
The implementation of simultaneous fitting in Laura ++ is based on the Jfit framework, an overview of which is given here and explained in more detail in Ref. [71]. The framework is based around a master-worker architecture in which the master drives the minimiser, combining the likelihood values for each category that have been calculated by a number of workers. The master and each of the workers run as separate processes that communicate via network sockets. This means that the calculation of the likelihood for each category is performed in parallel, increasing the speed of the computations if sufficient CPU cores are available. Note that the framework even allows the master and various worker processes to run on separate hosts with a modest performance penalty due to the increased latency in the communication between the processes, which depends on the network speed between the hosts, as shown in Table 5.
In order to modify an existing fit code to act as a worker process it is sufficient to modify only the final line: where the additional arguments to the runSlave function are a string and an unsigned integer that specify the hostname (e.g. "localhost" when all processes are running on the same host) and port number on which the master process is listening for connections. The code to start a master process is extremely simple: L au Si mF i tM as te r master ( nSlaves , port ) ; master . runSimFit ( ntupleName , nExpt , firstExpt , useAsymmErrors , twoStageFit ) ; The arguments to the constructor specify the number of worker processes that are expected to connect and the port on which it should listen for connections (a value of 0 indicates that it should use the first available port). The arguments to the runSimFit function specify the name of the file to which the fit results ntuple should be written, the number of experiments to be run and the ID of the first experiment, whether or not asymmetric uncertainties should be determined, and whether or not the fit consists of two stages. A shell script can then be used to launch a master process and the appropriate number of worker processes. An example script, runMasterSlave.sh, and source code for master and worker executables, Master.cc and Slave.cc, are included in the package examples.
During the set-up phase, each of the worker processes provides to the master the list of fit parameters needed for that worker to calculate the value of its likelihood function. The master stores this information and configures the minimiser. As part of this procedure, it must decide which parameters are common to some or all workers, such that each of these are provided only once to the minimiser. This decision is based purely on the parameters' names, so it is vital that these are correctly set when building the fit model in each of the workers. In order to aid the user, a summary is printed by the master process at the end of the set-up of the fit parameter lists.

Weighting events
Typically when generating samples of full detector simulation for three-body decays they are generated uniformly either in phase space or in the SDP, such that they can be used to describe the variation of the efficiency as a function of phase-space position. While it would often be useful to also have samples that are generated according to an amplitude model or indeed a selection of models, it can be infeasibly expensive in terms of CPU time and disk space, particularly at hadron collider experiments. It is therefore convenient to be able to apply per-event weights to such samples so that the weighted distributions reproduce those that would be obtained if the samples had been generated according to a particular model. The workflow to achieve this goal is described here.
The common parts of the workflow are the same as if preparing to generate toy MC samples or to fit a data sample. The only exception is that the argument to the LauDaughters constructor that determines whether or not the SDP coordinates are calculated takes on an additional significance -it indicates to the weighting procedure whether the provided sample was originally generated uniform in the conventional or square Dalitz plot. In case the original sample was not generated as a uniform distribution in either DP or SDP variables, this must be accounted for by applying a user-calculated correction to the weighting factor.
The weights are calculated using the MC-truth coordinates of the events (so as to preserve any resolution/migration effects from the full simulation). Hence, it is required that these coordinates are available in the input file with the following variable names: "m13Sq MC" and "m23Sq MC". The weight is calculated as the total amplitude-squared divided by the maximum value of the amplitude squared, which should be set by the user (as described in Sec. 7.2). If the input sample was generated uniformly in the SDP then the weight is multiplied by the Jacobian of Eq. (28). The weights are written to an ntuple in a new file, which has the same name as the input data file with " DPweights" appended. The variables iExpt and iEvtWithinExpt are also written to the ntuple and an event index is built from their values. This allows the weights ntuple to be used in conjunction with the input data ntuple via the Root "friend tree" mechanism.

Performance
A DP fitting package such as Laura ++ must be highly performing in several different ways. First, it must allow the user to obtain a good fit to data, and to evaluate the goodness of the fit. Several different methods to quantise the goodness of a multidimensional fit have been proposed in the literature and those which are available in Laura ++ are described below. A selection of examples of results obtained using Laura ++ is then given. Another important performance metric is the speed of execution; this is discussed for a number of example uses of the package.

Goodness of fit
Evaluating the level of agreement between an amplitude fit and the data can be difficult. Three methods to perform this task are discussed further; a two-dimensional binned χ 2 test, a mixed-sample test and a point-to-point dissimilarity test. These methods are described in detail in Ref. [72]. It should be noted that all methods test whether the data and the model are consistent to within the statistical uncertainty of the data sample; in some cases it may be necessary to consider in addition whether systematic effects could lead to differences between the data and the model.

Binned χ 2 method
The SDP distribution of the data is divided into bins with approximately equal bin content using an adaptive binning technique. The same binning distribution is applied to a sample of toy events that are generated from the amplitude fit model. A standard χ 2 test is then performed to compare the data and toy distributions within the chosen binning scheme. The relevant test statistic is where d i and t i are the number of events in the i th bin from data and toy, respectively, and n bins is the number of bins. Note that the generated toy sample can be much larger than the data, with the t i values obtained by scaling appropriately to correspond to the expectation for the data in each bin from the result of the fit. A drawback of this method is that the minimum number of events in each bin should not be too small, in order to have a reliable test statistic. But it should also not be too large, as this will cause the bins sizes to increase, leading to a loss of sensitivity to the variation of the amplitude over small scales. Typically a minimum number of events per bin of around 20 can be used, but the user should verify for themselves if this is appropriate in their case.

Mixed-sample method
The mixed-sample method tests how likely it is that the data and toy samples, produced from the fit model, come from the same parent distribution by evaluating where n data and n toy are the number of data and toy events, respectively. The number of nearest-neighbours to each data or toy data point considered by the test is given by n k . The term I (i, k) is equal to 1 if the i th event and its k th neighbour belong to the same sample and is 0 otherwise. Reference [72] suggests that n k = 10 and n toy = 10n data are sensible values. The statistic T H can be calculated many times, by using subsamples of the data and toy events, to build up a distribution of values. The quantity used to evaluate goodness-of-fit is (T H − µ T ) /σ T , where µ T and σ T are the mean and standard deviation of T H , respectively. Thus, by definition, the distribution of (T H − µ T ) /σ T has a mean of 0 and a width of 1 in the case that the data and toy samples are identical.

Point-to-point dissimilarity method
The consistency of a data sample and a toy sample generated from a model obtained by fitting the data can also be assessed using the following test statistic, Here, x denotes DP position, and ψ(| x i − x j |) is a weighting function. It can be shown that, in the limit of infinite statistics, with the choice ψ(| x i − x j |) = δ(| x i − x j |), the expression of Eq. (35) is equivalent to a χ 2 statistic [72]. In realistic scenarios, a form for the weighting function must be chosen, and the most appropriate choice may depend on the specific use case. The choice has been shown to work well in DP analysis [72]. The term is the value of the model at position x and dx is the area of the DP (which is included so that the mean value of the denominator becomes 1). The optimal value of the nuisance parameterσ is expected to be around the square of the typical width of the resonances in the DP in question, and thus usually ∼ 0.01 GeV 2 /c 4 . Unlike the mixed-sample test, the number of toy events should be large (n toy n data ) to avoid statistical fluctuations.
To compute a p-value using this test statistic, first the test statistic T h is calculated using the full available statistics. Then, a permutation test is performed as follows. The data and toy samples are pooled together and a new sample of size n data is randomly selected from the pooled sample. The new sample is then treated as the data sample, while the remaining events become the toy sample, and a new value of T h = T perm is calculated. This is repeated many times, and the p-value of the test is obtained from the fraction of times that T perm > T h . This can then be repeated with additional toy samples to build up a distribution of p-values.

Examples
The Laura ++ package has been used for numerous publications by several experimental collaborations and groups of phenomenologists. Below, several examples that demonstrate the features of the package are discussed. In addition, Laura ++ has also been used for various other studies of three-body charmless B meson decays by the BaBar collaboration [73][74][75][76][77], studies of charm decays by the LHCb collaboration [78], unpublished studies by several collaborations (for example Refs. [31,79]), and investigations into the phenomenology of different three-body decays [50,80,81].
Studies of charmless three-body B meson decays provide interesting opportunities to investigate the dynamics of hadronic B decays including potential CP violation effects. The B + → π + π + π − [33,82] and B + → K + π + π − [34,36] decays have been investigated by the BaBar collaboration using the Laura ++ package. In the most recent amplitude analysis of B + → π + π + π − decays [82], the amplitude model includes contributions from the ρ(770) 0 , ρ(1450) 0 , f 2 (1270), f 0 (1370) resonances and a nonresonant component. In the most recent amplitude analysis of B + → K + π + π − decays [36], the amplitude model includes the K * (892) 0 , K * 2 (1430) 0 , ρ(770) 0 , ω(782), f 0 (980), f 2 (1270), f X (1300) and χ c0 resonances together with Kπ S-wave and nonresonant components. In both cases, CP violation is allowed in the amplitudes. Projections of the fit results around the ρ(770) 0 resonance are shown in Fig. 5. The B + → π + π + π − data are consistent with CP conservation while there is evidence for CP violation in B + → ρ(770) 0 K + decays, which becomes more evident when inspecting the data in different regions of the π + π − helicity angle, θ π + π − . As model-independent analyses of larger data samples of these decays by the LHCb collaboration [83][84][85] have revealed large CP violation effects that vary significantly across the DP, there is strong motivation for updated amplitude analyses.  Figure 5: Projections of the data and fit results onto the π + π − invariant mass in the ρ(770 0 ) region, for (left) B + → π + π + π − [82] and (right) B + → K + π + π − [36] candidates observed by the BaBar collaboration. In both cases the top row shows all candidates, the middle row shows those with cos θ π + π − > 0 and the bottom row shows those with cos θ π + π − < 0, while in each row the left (right) plot is for B − (B + ) candidates. The data are the points with error bars, the red/dark filled histogram shows the continuum background component, the green/light filled histogram shows the background from other B meson decays, and the blue unfilled histogram shows the total fit result.
Understanding the origin of these CP violation effects requires related modes to also be studied. The Laura ++ package has also been used by the BaBar collaboration for a time-dependent DP analysis of B 0 → K 0 S π + π − decays [44], as well as for an amplitude analysis of B + → K 0 S π + π 0 decays [86]. In the latter, the modelling of the large background contribution, as well as of the smearing of the DP position due to the limited resolution of the neutral pion momentum (self cross-feed), is particularly important. In addition, correlations between the DP position and the variables that are used to discriminate signal decays from background contributions are taken into account as described in App. D.13. The amplitude model includes components from the K * (892) resonance and Kπ S-wave (both appearing in both charged and neutral channels) as well as the ρ(770) + resonance. Projections of the fit results are shown in Fig. 6. The analysis reveals evidence for a CP asymmetry in B + → K * (892) + π 0 decays.  Figure 6: Projections of the data and fit results onto (top) K 0 S π ∓ , (middle) K 0 S π 0 and (bottom) π ∓ π 0 invariant mass distributions for B + → K 0 S π + π 0 candidates observed by the BaBar collaboration [86]. Background from D 0 → K 0 S π 0 has been vetoed. In each row the left (right) plot is for B − (B + ) candidates. The data are the points with error bars, the (black) dash-dotted curves represent the signal contribution, the dotted (red) curves to the continuum background component, the dashed (green) curves to the total background contribution and the solid (blue) curves the total fit result.
The LHCb collaboration has used the Laura ++ package for several studies of multibody B meson decays to charmed final states, with important results for charm spectroscopy and CP violation measurements. For example, the B 0 s → D 0 K − π + decay was found to have a DP structure that contains effects from overlapping spin-1 and spin-3 resonances with masses around m(D 0 K − ) ∼ 2.86 GeV/c 2 [87,88]. The neutral charm meson is reconstructed through its D 0 → K + π − decay. The model contains contributions from the K * (892) 0 , K * (1410) 0 , K * 2 (1430) 0 , K * (1680) 0 resonances as well as a K − π + S-wave component, and  A similar DP analysis with the Laura ++ package has been performed by the LHCb collaboration for the B 0 → D 0 K + π − mode [89]. The model obtained is an essential input into a subsequent analysis of the same decay with the neutral charm meson reconstructed through D decays to the CP eigenstates K + K − and π + π − [70]. In the latter case, contributions from both B 0 → D 0 K + π − and D 0 K + π − amplitudes can interfere, giving sensitivity to the angle γ of the CKM unitarity triangle. A DP analysis allowing for CP violation effects provides a powerful way to determine γ without ambiguities [51,52]. In this analysis, a simultaneous fit to the final states with different D decays is implemented using the Jfit method described in Sec. 7.3.7 and Ref. [71]. Projections of the fit result onto the data (weighted by the signal purity) are shown in Fig. 8. No significant CP violation effect is observed, and the resulting limits on γ are not strongly constraining. The method is, however, expected to give competitive constraints on γ as larger data samples become available and as additional D meson decay modes are included in the analysis. Moreover, the analysis also gives results for hadronic parameters that must be known in order to interpret results from quasi-two-body analyses of B 0 → DK * (892) 0 decays in terms of γ. As such, the results have an important impact in combinations of results to obtain the best knowledge of γ [90,91].  [70]. A legend describing the various contributions is also given.
Other decays of the type B → D ( * ) Kπ have sensitivity to γ, and are also important to study as possible background contributions to the two-body B → D ( * ) K type decays that are more conventionally used for this purpose. The LHCb collaboration has also published results on the B + → D − K + π + [92] and D + K + π − [93] decays, which were obtained from analyses using the Laura ++ package. The higher-yield B → D ( * ) ππ channels provide some of the most interesting possibilities to explore charm spectroscopy. An amplitude analysis of B + → D − π + π + [94] has been performed by the LHCb collaboration, using the Laura ++ package, in which the model contains contributions from the D * states. In the absence of sufficiently detailed models for the D − π + S-wave, a quasi-modelindependent description based on spline interpolation is used. Projections of the results of the fit onto the data are shown in Fig. 9.   Figure 9: Projections of the data and fit results onto m(D − π + ) min for B + → D − π + π + candidates observed by the LHCb collaboration [94] on (top left) linear and (bottom right) logarithmic y-axis scales. Here, m(D − π + ) min is the smaller of the two values of m(D − π + ) for each B + → D − π + π + candidate. A legend describing the various contributions is also given. The (bottom right) Argand diagram of the D − π + S-wave amplitude shows the expected phase motion corresponding to the D * 0 (2400) 0 resonance. The numbered points correspond to the spline knots.

Speed
It is essential for the Laura ++ amplitude analysis package to run quickly, since otherwise the large data samples available in modern experiments can lead to unmanageably long execution time. In this subsection some performance benchmarks are provided. More specifically, a selection of the examples that are provided with the package (several of which are based on analyses presented in the previous subsection) are run out of the box on the same machine (an Intel Core i5-3570 3.4 GHz quad-core CPU with 8 Gbytes of RAM). In each case, timings for both generation of 50 toy datasets and for fitting those same 50 datasets are provided in Table 3. The fitting times are averaged over 20 fits with randomised starting parameters. The scenario demonstrated in each example is as follows: • GenFit3pi.cc Example analysis of the symmetric final state B + → π + π + π − , using LauSimpleFitModel  Table 4.
• GenFit3KS.cc Example analysis of the fully-symmetric final state B 0 → K 0 S K 0 S K 0 S , using LauSimpleFitModel (i.e. not including effects of CP violation). By default there are 10000 signal events per experiment and the signal isobar model contains four components: f 0 (980), f 0 (1710), f 2 (2010) and χ c0 . All of the resonance parameters are fixed in the fit. No background contributions are included.

• GenFitDs2KKpi.cc
Example analysis of the decay D + s → π + K + K − , using LauSimpleFitModel (i.e. not including effects of CP violation). By default there are 10000 signal events and the signal isobar model contains three components: φ(1020), K * (892) 0 and a nonresonant component. All of the resonance parameters are fixed in the fit. No background contributions are included.

• GenFitEFKLLM.cc
Example analysis of the decay B 0 → D 0 K + π − , using LauSimpleFitModel (i.e. not including effects of CP violation). By default there are 5000 signal events per experiment and the signal isobar model contains two components, which are the two parts of the EFKLLM model for the K + π − S-wave (see Eq. (52) in App. A). This example mainly serves to demonstrate how to use this particular model.
• GenFitBelleCPKpipi.cc Example analysis of the decay B + → K + π + π − , using LauCPFitModel, which includes effects from CP violation. By default there are 5000 signal events per experiment and the signal isobar model contains seven components: K * (892) 0 , K * 0 (1430) 0 , ρ(770) 0 , f 0 (980), χ c0 , and two nonresonant components. The slope of the π + π − exponential nonresonant model and the effective range and scattering length of the K + π − LASS nonresonant component are floating parameters by default. The isobar coefficients use the LauBelleCPCoeffSet form to parameterise the CP violation (see Table 1) and a two-stage fit is employed. No background contributions are included.

• GenFitKpipi.cc
Example analysis of the decay B + → K + π + π − , using LauCPFitModel, which includes effects from CP violation. This example is based closely on the BaBar analysis from Ref. [36]. By default there are 4585 signal events per experiment and the signal isobar model contains seven components: K * (892) 0 , K * 0 (1430) 0 , K * 2 (1430) 0 , ρ(770) 0 , ω(782), f 0 (980), f 2 (1270), f 0 (1300), χ c0 , and a nonresonant component. All of the resonance parameters are fixed in the fit. The isobar coefficients use the LauCartesianCPCoeffSet form to parameterise the CP violation (see Table 1) and a two-stage fit is employed. Background contributions for combinatorial candidates and those from other B decays are included.

• KMatrixExample.cc and KMatrixDto3pi.cc
These examples demonstrate how to use the K-matrix description of the S-wave. The first scenario is for the B 0 → π + π − K 0 S DP but the only terms included in the signal model are ρ(770) 0 , K * (892) + and a K-matrix component for the π + π − S-wave. By default there are 5000 signal events per experiment. The second scenario is for the symmetric D + → π + π + π − DP and includes in the signal model the ρ(770) 0 , f 2 (1270), and the K-matrix. By default there are 20000 signal events per experiment. In both scenarios there are no background contributions and all resonance parameters are fixed in the fit.

• GenFitNoDP.cc and GenFitNoDPMultiDim.cc
These examples demonstrate how to use the package to perform fits to variables other than the DP. The first case fits only a single variable (the invariant mass of the B + candidate), while the second performs a 2D fit. In both cases there are 5000 signal events. In the first case there are two background components included, which have yields of 7000 and 3000 events. In the second case there is a single background component that has a yield of 7000 events. All yield parameters are floating, along with some of the shape parameters of the PDFs. Asymmetric uncertainties are evaluated in the first case.
• runMasterSlave.sh, Master.cc and Slave.cc This example demonstrates how to perform a simultaneous fit to two categories of data in the decay channel B 0 → π 0 π 0 K 0 S . The data are split based on the reconstruction category of the K 0 S candidate. In one category there are 500 signal and 1200 background events, while in the second there are 750 signal and 2500 background events. The common signal DP model contains contributions from the f 0 (980), f 2 (1270), K * (892) 0 , and K * 0 (1430) 0 resonances. The mass of the K * (892) 0 resonance is a floating parameter of the fit and a two-stage fit is employed. The background component is distributed uniformly in the DP by default. To demonstrate further the impact on the performance of simultaneous fitting, this example is run with: the Master and the two Slave processes all running on the same host, the Master and the two Slave processes running on three separate hosts, and the timings for each of these scenarios are provided in Table 5. To run this particular example, the hosts have dual Intel Xeon E5-2620v2 2.1 GHz 6-core CPUs and 64 Gbytes of RAM and are connected via 10 Gbits/s ethernet.

Future developments
There are several features that would help to improve and extend the functionality of Laura ++ in the future. Some of these potential future developments are described below. In addition, it is anticipated that the Laura ++ code will be continually updated to take advantage of features in the latest C ++ standards.

Plotting the amplitude
Currently the user can easily make plots of the DP distribution or its projections from the result of the fit (see Sec. 7.2). It can also be of interest to draw amplitude-level quantities, for example to show the phase variation with two-body invariant mass of a particular partial wave. While this is currently possible in Laura ++ (see Fig. 9 for an example), it would be desirable to provide an interface to simplify matters for the user.

Decay-time-dependent fits
As mentioned in Sec. 1, the evolution of DP structure with decay time of a B or D meson can be of interest to study CP violation effects. For example, studies of B 0 → π + π − π 0 and Dπ + π − decays are of interest to measure the angles α and β of the CKM Unitarity Triangle with low theoretical uncertainty. Studies of B 0 → K 0 S π + π − and K 0 S K + K − decays enable determinations of β that are not as clean, but are potentially sensitive to effects of physics beyond the Standard Model. Also, the D 0 → K 0 S π + π − channel appears the most sensitive to possible CP violation effects associated with charm mixing. Many other channels are potentially of interest.
The Laura ++ package has been used for decay-time-dependent DP analysis, for example, in Refs. [44,50]. However, this implementation was experiment-specific and therefore unsuitable for use in the more general case. Further development is necessary to establish a model that can deal generically with issues such as flavour tagging, decay time resolution and acceptance as well as production and detection asymmetries.

Alternative handling of resolution effects
The treatment of resolution discussed in Sec. 3 and Sec. 5.2 is completely general, but is likely to be inefficient in certain cases. For example, in the B + → K + K − K + decay, there are narrow contributions from the φ(1020) and χ c0 resonances in specific regions of the DP, for which resolution effects may be important. The rest of the DP is populated with broad or nonresonant states so that resolution can be safely neglected. An approach in which Gaussian (or non-Gaussian) smearing of Dalitz plot position can be used in only selected regions of the phase space may be useful. In such a case it will be necessary to take care to avoid issues due to edge effects, including the possibility of an event being smeared to positions outside the kinematic boundary of the DP.

Non-zero spins
There are many interesting three-body decays that include particles with non-zero spin in the initial and final states, which cannot be fitted using the current version of Laura ++ . Large samples of b-baryons are available in the LHCb data samples, which have baryons in the initial and final states, for example Λ 0 b → D 0 pπ − and Ξ − b → pK − K − decays. Decays of B → J/ψ hh , where h and h are charged pions or kaons, are also interesting and contain the spin-1 J/ψ particle. Another group of decay modes, B → D * hh , which include the vector D * meson, are interesting for spectroscopy of D * * and D * * s states. To enable the Laura ++ package to cope with the decays above, several things would require updating or changing. Firstly, the phase space of the problem is expanded from two to five dimensions where additional degrees of freedom would be some angular variables. The helicity of particles of non-zero spin, like the J/ψ , must also be considered, requiring a sum over the helicity states. For the initial and final state this sum must be incoherent and for the intermediate states a coherent sum is required. The Zemach spin terms currently implemented must also be changed for the non-zero spin particles. One potential way to calculate the spin factors would be to interface Laura ++ with the qft++ package [95], which allows the spin terms to be calculated for any process.

Genetic algorithms
To ensure that the global minimum in the NLL has been found, Laura ++ allows many fits to be performed with randomised starting values for the floated fit parameters. This method works well, but can become time consuming if the global minimum is not found in a high proportion of fit attempts. Genetic algorithms could provide a method to find sensible starting values for the floated parameters such that the global minimum is always found. This could be achieved by interfacing to existing software packages with implementations of genetic algorithms.

Interface to EvtGen
The EvtGen package [96] is designed to predominantly simulate the decays of b-and c-hadrons. It is important in experimental particle physics to produce simulated data samples that are realistic and match the true data samples. Typically in three-body decays the simulated events are produced flat in the DP, without resonant contributions. It would be beneficial if Laura ++ could be used directly to provide realistic DP distributions for simulated samples of three-body decays.

Summary
The Laura ++ package provides a flexible and optimised framework for Dalitz-plot analysis. While it can be used for the decay of any stable spin-zero particle to any final state containing three stable spin-zero particles, it has until now been most widely used for decays of B or D mesons to three pseudoscalars. Features included in Laura ++ allow the physics of such decays to be probed in detail, including studies of the resonances appearing in the contributing partial waves, and investigations of CP -violating effects. Use of the Laura ++ software has resulted in numerous publications to date, with many more expected in future with the increasingly large data samples available at LHCb, Belle II and other experiments.

Acknowledgements
The Laura ++ package has been under development for many years, with support principally from the Science and Technology Facilities Council (United Kingdom) and by the European Research Council under FP7. Individual authors have received support from Marie Sk lodowska-Curie Actions and from the University of Warwick. We thank the Laura ++ user community for feedback, bug reports, testing and other contributions to improve the package, in particular Louis Henry, Adlene Hicheur, Patricia Magalhães, Jussara Miranda, Sian Morgan and Charlotte Wallace. We also acknowledge productive discussions regarding aspects of Dalitz plot analysis with many members of the BaBar, Belle and LHCb collaborations, notably Eli Ben-Haim, Alex Bondar, Jeremy Dalseno, Bill Dunwoodie, Brian Meadows and Anton Poluektov. Similarly, instructive communications with members of the theory community have been of great benefit; in particular we thank Vladimir Anisovich, David Bugg, Leonard Lesniak, Benoit Loiseau, Mike Pennington, Alessandro Pilloni, Andrey Sarantsev and Adam Szczepaniak. Furthermore, we thank Bertram Kopf and Matthias Steinke for providing access to the PAWIAN software (PANDA collaboration) for cross-checking the K-matrix implementation. Finally, we acknowledge useful input on technical features of the package from René Brun and Bertrand Echenard.

Appendices A Formulae for available lineshapes
This section presents the complete formulae for all resonance shapes implemented in Laura ++ . Table A1 gives the list of shapes, together with the corresponding LauResonanceModel enumeration integer that is required to specify the resonance type for the LauIsobarDynamics addResonance function, as well as the equation number(s) that provide the expression for the resonance mass term R(m) used in Eq. (4). The K-matrix shape is particularly complicated and is therefore described in a dedicated subsection. The simple Breit-Wigner lineshape is given by where m 0 is the pole mass and Γ 0 is the resonance width. The more commonly used relativistic Breit-Wigner lineshape is described in Sec. 2.1. 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 eff 0 within the kinematically allowed range. This is accomplished with the ad-hoc formula where m max and m min are the upper and lower limits of the kinematically allowed mass range. For virtual contributions, only the tail of the RBW function enters the Dalitz plot. The Gounaris-Sakurai form of the Breit-Wigner lineshape [97] is usually used as an alternative model for the ρ resonance. It is given by where q is the magnitude of the momentum of one of the daughter particles in the resonance rest-frame, and The normalisation condition at R(0) fixes the parameter D = f (0)/(Γ 0 m 0 ), and is found to be [97] The Flatté [98], or coupled two-channel Breit-Wigner, lineshape is usually used to model f 0 (980), K * 0 (1430) and a 0 (980) states. It was originally introduced in the form The decay widths in the two systems are usually represented by products of couplings and dimensionless phase-space factors: Here the fractional coefficients come from isospin conservation, m i,j 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 A2. The Clebsch-Gordan coefficients in Eqs. 45 and 46 are not guaranteed to be correct for every possible resonance that could be modelled with the Flatté lineshape, but are appropriate for every case considered in Table A2. The expressions for the widths are continued analytically (Γ → i|Γ|) 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 [99] (otherwise f A is set to unity).
Variants of the Flatté lineshape have been used in the literature. In some cases, e.g. Refs. [100,101], the constant m 0 that multiplies the widths in the denominator of Eq. (44) is absorbed into the couplings. As a consequence the couplings have dimensions of mass-squared, and are sometimes denoted as g i [100] and sometimes as g 2 i [101]. In Table A2 all values have been converted to be consistent with Eqs. (44)- (46). In Laura ++ it is only possible to use the Flatté lineshape for the systems specified in Table A2. At construction time the resonance name is checked and the corresponding parameter values are set; these can be modified by the user if desired. Table A2: 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é lineshapes. Units of GeV for g 1,2 (or GeV 2 for m 0 g 1,2 when taken from Refs. [100,101]) and GeV 2 /c 4 for s A are implied.

Resonance Channel 1
Channel 2 g 1 or m 0 g 1 g 2 or m 0 g 2 s A Reference f 0 (980) π 0 ,π 0 ,π ± ,π ± K ± ,K ± ,K 0 ,K 0 0.165 4.21g 1 - The σ or f 0 (500) → ππ and κ or K * 0 (800) → Kπ low-mass scalar resonances can be described using the form where M is the mass where the phase shift goes through 90 • for real s ≡ m 2 , and the width where the square-root term is the phase space factor, which requires the invariant masses of the daughter particles m 1 and m 2 , s A is the Adler-zero constant, while b 1 , b 2 and A are additional constants [99]. Table A3 gives the default values of the parameters. Table A3: Default values of the parameters for the σ and κ lineshapes based on BES data [99].
The Dπ S-wave can be parameterised using the form provided by Bugg [102], who labels the pole state as "dabba": where ρ is the Lorentz invariant phase space factor 1 − s 0 /m 2 , s 0 is the square of the sum of the invariant masses of the D (m D ) and π (m π ) daughters, s A is the Adler-zero term m 2 D − 0.5m 2 π that comes from chiral symmetry breaking [103], while b = 24.49, α = 0.1 and β = 0.1.
The RBW function is a very good approximation for narrow resonances well separated from any other resonant or nonresonant contribution in the same partial wave. This approximation is known to be invalid in the Kπ S-wave, since the K * 0 (1430) resonance interferes strongly with a slowly varying nonresonant term [104]. The so-called LASS lineshape [105] has been developed to combine these amplitudes, with cot δ B = 1 aq where m 0 and Γ 0 are now the pole mass and width of the K * 0 (1430), and a and r are parameters that describe the shape. Most implementations of the LASS shape in amplitude analyses of B meson decays [34,106] apply a cut-off to the slowly varying part close to the charm hadron mass (∼ 1.7 GeV/c 2 ).
An alternative representation of the Kπ S-wave amplitude can be made using the EFKLLM model described in Ref. [107] (the acronym comes from the surnames of the authors of that paper), which uses a tabulated form-factor f Kπ 0 (m 2 ) that is interpolated using two splines (one each for the magnitude and phase parts), multiplied by a scaling power-law mass-dependence m , leading to where suggested values for the exponent are zero for κ (this is also the default value) and −2 for K * 0 (1430). Because of the large phase-space available in three-body B meson decays, it is possible to have nonresonant amplitudes (i.e. contributions that are not associated with any known resonance, including virtual states) that are not constant across the Dalitz plot. One possible parameterisation, based on theoretical considerations of final-state interactions in B ± → K ± π + π − decays [108], uses the form where with the constant parameters d 0 = 1.3232 × 10 −3 GeV −8 , a 1 = 0.65 GeV −2 , b 1 = 18 GeV 2 , a 2 = 0.55 GeV −2 and b 2 = 15 GeV 2 in natural units. There are several empirical methods that can be used to model the nonresonant contributions. One is to use an exponential form factor [32] R(m) = e −αm 2 , while another form is simply a power-law distribution where in both cases α is a parameter that must be determined from the data. For symmetric DPs, the exponential form is modified to while a Taylor expansion up to first order can also be used: where m P is the invariant mass of the parent P . Another possible description for nonsymmetric DPs is based on the polynomial function [109] where m k is the invariant mass of daughter particle k and n is the integer order equal to 0, 1 or 2; a quadratic dependence in m can be constructed by using up to three polynomial R(m) terms, one for each order along with their individual c j amplitude coefficients. We next come to the model that implements the ρ−ω mass mixing amplitude described in Ref. [110] where A ρ is the ρ lineshape, A ω is the ω lineshape, |B| and φ B are the relative magnitude and phase of the production amplitudes of ρ and ω, and ∆ ≡ δ(m ρ + m ω ), where δ governs the electromagnetic mixing of ρ and ω (with pole masses m ρ and m ω ). Here, the amplitude A ω is always given by the RBW form of Eq. (6), while the amplitude A ρ can either be represented using the Gounaris-Sakurai formula given in Eq. (39) or the RBW form; the required shape is selected using either the RhoOmegaMix GS or RhoOmegaMix RBW enumeration integer labels given in Table A1. When ignoring the small ∆ 2 term in the denominator of Eq. (60), this is equivalent to the parameterisation described in Ref. [111]; this option can be chosen using either the RhoOmegaMix GS 1 or RhoOmegaMix RBW 1 enumeration labels, depending on what lineshape is needed for the ρ resonance. From SU(3) symmetry, the ρ and ω are expected to be produced coherently, which gives the prediction |B|e iφ B = 1. To avoid introducing any theoretical assumptions, however, it is advisable that these parameters are left floating in the fit. In general δ is complex, although the imaginary part is small so this is neglected. The theory prediction for δ is around 2 MeV [112], and previous analyses have found |δ| to be 2.15 ± 0.35 MeV [110] and 1.57 ± 0.16 MeV, and arg δ to be 0.22 ± 0.06 [111]. These parameters can be also be floated in the fit. If the dynamical structure of the DP cannot be described by any of the forms given above, then the LauModIndPartWave class can be used to define a model-independent partial wave component, using splines to produce an amplitude. It requires a series of mass points called "'knots", in ascending order, which sets the magnitude r(m) and phase φ(m) at each point m that can be floated when fitted to data: The amplitude for points between knots is found using cubic spline interpolation, and is fixed to zero at the kinematic boundary. There are two implementations for representing the amplitudes: one uses magnitudes and phases (MIPW MagPhase), while the other uses real and imaginary parts (MIPW RealImag). Finally, the incoherent Gaussian lineshape form is given by where m 0 is the mass peak and G 0 is the width. This can be used to parameterise the amplitude for a very narrow resonance, where the measurement of the width is dominated by experimental resolution effects, producing a lineshape that is indistinguishable from a Gaussian distribution. The narrow width ensures that the resonance will not interfere significantly with other resonances in the DP, i.e. it will be incoherent. This form could also be used to parameterise narrow background resonance contributions that would otherwise be excluded with mass vetoes, such as the D 0 meson decay to K − π + in the charmless mode B − → K − π + π − , when used with the addIncoherentResonance function of LauIsobarDynamics.

A.1 K-matrix
The isobar model, described earlier in Sec. 2, can be used to describe the dynamics of three-body decays when the quasi two-body resonances are relatively narrow and isolated. However, this model does not satisfy scattering (S-matrix) unitarity, thereby violating the conservation of quantum mechanical probability current, when there are broad, overlapping resonances (with the same spin-parity), such as the intermediate S-wave states σ for ππ and κ for Kπ channels. Assuming that the dynamics is dominated by two-body processes, meaning that the S-wave does not interact with the rest of the products in the final state, then unitarity is naturally conserved within the K-matrix approach [113], which was originally developed for two-body scattering [114] and the study of resonances in nuclear reactions [115,116], but was extended to describe resonance production in a more general way [117]. The scattering matrix S, describing the general transformation of an initial state to a final state, can be defined as where I is the identity matrix, representing the case when the initial and final states do not interact at all, and T is the physical (observable) transition matrix. Conservation of scattering probability means that the n × n S matrix, where n is the number of channels, is unitary (SS † = S † S = I). The factor 2i is introduced so that the transition amplitude for a single resonance channel corresponds to a circle of unit diameter centred at (0, i/2) in the complex plane; physically allowed values of T will be along the boundary (elastic scattering) or inside (inelastic scattering) this unitarity circle. Using Eq. (63), it can be shown that the n × n K matrix operator, defined as is Hermitian (K = K † ) [113]. Furthermore, K is real and symmetric, owing to the time-reversal invariance of the S and T matrices. Rearranging Eq. (64) gives the following expression for the scattering transition operator in terms of the K matrix: The normalisation of the two-body wave functions requires the inclusion of phase-space factors in both the initial and final states [118]. This then leads to the following definition of the Lorentz-invariant transition amplitudeT : where u and v indicate the channel indices (from 1 to n) and ρ is the normalised diagonal n × n phase space matrix, whose elements are equal to 2q/m, where q is the magnitude of the momentum of either daughter in the rest-frame of the two-body state that has invariant mass m = √ s. In general, the phase space element of channel u is given by where m 1u and m 2u are the rest masses of the two daughters [62], and is continued analytically by setting ρ u to be i|ρ u | when the channel is below its mass threshold, provided it does not cross into another channel. The Lorentz-invariant form of the K matrix, which will also be real and symmetric, can be written asK which then implies that the Lorentz-invariant transition amplitude is given bŷ We can now use this expression to give the general amplitude of the production of overlapping resonance states. This model or ansatz [117] describes the amplitude of channel u in terms of the initialP -vector preparation of channel states v, that has the same form asK, transforming ("scattering") into the final state u via the propagator term (I − iKρ) −1 : The scatteringK matrix can be parameterised as a combination of the summation of N poles with real bare masses m α , together with nonresonant slowly-varying parts (SVPs), so-called since they essentially have a 1/s dependence, with real (and symmetric) coupling constants f scatt uv [119]: where g (α) u denotes the real coupling constant of the pole m α to the channel u, the factor is the Adler zero term that suppresses the false kinematic singularity when s goes below the ππ production threshold [103], while m 2 0 , s scatt 0 , s A and s A0 are real constants of order unity; typical values are m 2 0 = 1 GeV 2 /c 4 , s scatt 0 = −5 GeV 2 /c 4 , s A = 1, and s A0 = 0 GeV 2 /c 4 . Note that the real poles m α are the masses of the so-called bare states of the system, which do not correspond to the masses and widths of resonances (mixtures of bare states) from the complex poles in the physical T matrix. The production vectorP is parameterised in an analogous form to theK matrix: where β α and f prod v (which both depend on the final state channel u) are complex production constants for the poles and nonresonant SVPs, respectively, and s prod 0 is of order unity and is usually taken to be approximately equal to s scatt 0 . It is important that the production and scattering processes use the same poles m α , otherwise the transition amplitude would vanish (diverge) at theK-matrix (P -vector) poles; the singularities need to cancel out for the total amplitude. Also note that the Adler zero suppression factor given in Eq. (72) is generally not needed forP , since its inclusion does not improve the description of S-wave amplitudes found in experimental data [16,22,120,121].
In order to clarify what amplitudes are used, we can separate out the production pole and SVP terms shown above. The amplitude of each production pole m α for the final state u is given by where we need to sum the propagator contributions over the channels v, while the SVP production amplitudes are separated out for each individual v → u channel as We can then sum all of these contributions to give the total S-wave amplitude The elements ρ u of the diagonal phase space matrix depend on the channel type u. For ππ systems, the five available channels are ππ, KK, 4π, ηη and ηη multimeson states [119]; note that η η is above the open charm threshold and is not considered. The phase space factor for ππ (u = 1), KK (u = 2) and ηη (u = 4) is given by Eq. (67), with m 1u and m 2u equal to the rest masses of the two pseudoscalars forming channel u (m 1u = m 2u ). For ηη (u = 5), the second multiplicative term involving m η − m η is ignored (set to unity), since below threshold this crosses channels and we cannot continue this analytically in the usual way. As given in Ref. [119], the phase space term for the 4π multimeson state (u = 3) is .
Here, s 1 and s 2 refer to the invariant mass-squared of the two di-pion states (which are simply considered as integration variables), M 0 is the pole mass of the ρ resonance (775 MeV) and Γ(s) = Γ 0 [1 − (4m 2 π /s)] 3/2 is the energy-dependent width, where Γ 0 is taken to be 0.3 GeV, which is approximately 75% of the total width of the f 0 (1370) → 4π state [62]. The constant factor ρ 0 ensures continuity at s = 1 GeV 2 /c 4 , while the limits of integration are 4m 2 π to ( √ s − 2m π ) 2 for s 1 and 4m 2 π to ( √ s − √ s 1 ) 2 for s 2 in order to satisfy kinematic constraints. The ρ 31 term needs to be evaluated numerically and is approximated very well by a 6 th order polynomial in s: ρ 31 (s) = 1.0789s 6 + 0.1366s 5 − 0.2974s 4 − 0.2084s 3 + 0.1385s 2 − 0.0193s + 0.0005. (78) For Kπ systems, we have the three channels Kπ, Kη and Kπππ multimeson states [122]. The phase space factors for the first two channels (u = 1,2) are again given by Eq. (67), while the multimeson phase space element is given by where r 0 is a constant of continuity at s = 1.44 GeV 2 /c 4 . The K-matrix formalism is a way to describe the dynamics of a set of broad, overlapping resonances with the same isospin I s , spin L and parity P . Final states with different I s L P values would require the appropriate number of K-matrices. To avoid overcomplicating the Dalitz plot analysis, the usual procedure is to parameterise only the S-wave (L P = 0 + ) components with the K-matrix approach, and then combine the other (narrow) resonances with the isobar model. This means that the total amplitude would be given by where F u,Is (s) is the K-matrix amplitude defined in Eq. (70) for the channel u and isospin state I s . The recommended procedure would then be to first use scattering data to completely define the K-matrix elements in Eq. (71), such as using the values quoted in Ref. [121] which are obtained from a global analysis of ππ scattering data [119]. Subsequently in the DP analysis the user can fit for the coefficients β α and f prod v of thê P -vector used in Eqs. (73), (74) and (75).

A.1.1 Implementation details for K-matrix
Special commands are required in order to use the K-matrix amplitude defined in Eqs. (70) and (76), which is combined automatically with the other isobar resonances to produce the total dynamical amplitude given by Eq. (80).
First, the (I − iKρ) −1 propagator term is formed using the defineKMatrixPropagator function in LauIsobarDynamics, which requires a descriptive name, a text file containing a keyword-defined list of the scattering and Adler zero coefficients, as well as an integer to specify which daughter is the bachelor particle. This function also requires the total number n of K-matrix scattering channels (sum over v = 1 to n), the number of bare poles N (sum of m α terms), as well as the final channel index u. Note that the completê K matrix in Eq. (71), which is real and symmetric, is found for all possible values of u and v in order to find the propagator; the specific u index is only needed for finding the final F u amplitude. For ππ S-wave, we normally have five channels (ππ, KK, 4π, ηη and ηη multimeson states), five poles and the index u is equal to 1 (ππ). The production vectorP defined in Eq. (73) is then formed using the addKMatrixProdPole and addKMatrixProdSVP functions of LauIsobarDynamics that create the β α pole term given by Eq. (74) (which internally sums the propagator function over the initial channels v owing to the g v coupling dependence) and the slowly-varying f prod v term given by Eq. (75), respectively. They each require a descriptive name, the name of the propagator term defined earlier and the pole or channel integer number (starting from 1). These functions also accept a boolean useProdAdler to specify if the Adler zero suppression factor given in Eq. (72) is also used for the production vectorP ; by default useProdAdler is set to false. Additional K-matrix amplitudes (e.g. for different isospin settings) can be included by simply defining additional propagators with unique names together with their required production terms.
Internally, the K-matrix propagator is defined by the LauKMatrixPropagator class, in which each unique propagator is created using an instance of the LauKMatrixPropFactory factory method, while the LauKMatrixProdPole and LauKMatrixProdSVP classes represent the production pole and slowly-varying terms, respectively. Since Root does not implement complex matrices, the K-matrix propagator is expanded into real and imaginary parts using the following method. If A, B, C and D are real matrices (TMatrix objects), then the propagator can be expressed as where A is equal to I +KIm(ρ) and B is −KRe(ρ). Both A and B are completely determined if we know the real, symmetricK matrix and the diagonal phase space matrix ρ (which can have imaginary terms if the invariant mass is below threshold). The real and imaginary propagator terms are then given by

A.1.2 Pedagogical K-matrix plots
In order to better understand the properties of the K-matrix description we will now show a series of instructional plots. The first of these is Fig. 10 which shows the transition amplitude of the ππ → ππ S-wave, corresponding to the first element T 11 of the T matrix defined in Eq. (66), using the parameters given in Table A4 and where we are not considering the effect of the production vectorP . Figure 10a) shows the phase motion of the amplitude, which lies within a circle of unit diameter centred on (0, i/2), while Fig. 10b) is the equivalent intensity or amplitude squared. First, we can see that the amplitude follows the unit circle anticlockwise, corresponding to the very broad σ or f 0 (500) resonance structure, until we reach an invariant mass near to the threshold of the f 0 (980) resonance, where its interference with the σ produces a striking dip in the intensity; the amplitude is purely elastic until we reach the f 0 (980). As we follow the phase motion counterclockwise, new channels such as KK open up at higher energies, producing other resonance structures that give extra interference terms and so the scattering process remains inelastic. A more detailed discussion of these features is given in Ref. [123], which has a slightly different amplitude intensity distribution at high invariant mass owing to different scattering data being considered. If we now imagine a vector starting from the centre of the unitarity circle (0, i/2) and ending on the position of the complex amplitude T 11 , then the phase shift δ is defined as half of the angle that subtends with the imaginary Table A4:K-matrix parameters taken from Ref. [121], which are obtained from a global analysis of ππ scattering data by Anisovich and Sarantsev [119]. Only f 1v parameters are listed here (ππ S-wave). Masses m α and couplings g (α) u are given in GeV/c 2 , while units of GeV 2 /c 4 for s-related quantities are implied; s prod 0 is taken from Ref. [22].
while the inelasticity η is defined as twice the length of Figure 10c) shows the evolution of δ with invariant ππ mass, while Fig. 10d) shows the variation of the inelasticity η, where a purely elastic (inelastic) process has η = 1 (η = 0); rapid changes in δ and η are observed at the thresholds of various resonance structures. We now move onto theK matrix itself, which is the main ingredient of the scattering propagator. Figure 11 shows the first row of theK matrix (K 1j ), where we have split up the various components that make up each matrix element. The red lines show only the bare pole m α contributions, given by the first summation term on the right hand side of Eq. (71), where u = 1, v = 1 − 5 and we sum over all five poles (α = 1 − 5), and the Adler zero suppression factor f A0 (s) is set to unity. All of the plots show the strong effect of the bare pole singularities. The blue lines show the rather small SVP contributions, corresponding to the second term on the right hand side of Eq. (71) (with f A0 (s) = 1). The summation of these various pole and SVP contributions is given by the magenta lines, while the black lines show the inclusion of the Adler zero suppression factor of Eq. (72), which only significantly effects the overall shapes at low invariant mass (m 0.5 GeV/c 2 ). The absence of singularities forK 13 below 1.2 GeV/c 2 is due to the fact that the coupling of the first two bare poles in the 4π channel is zero. Similar distributions are obtained for the other rows of theK matrix.  Next, let us look at the s-dependence of the first row of the propagator matrix [I − iKρ] −1 1v , since this will effectively modulate the individual production pole and SVP shapes that are combined to form the total ππ S-wave amplitude F 1 (s) using Eq. (70). Figure 12 shows the real and imaginary components, as well as the magnitude, of the propagator elements. The overall impression we get is that the propagator amplitudes have non-trivial variations as a function of √ s, owing to the matrix inversion process mixing and transposing the superposition of the bare pole states m α . These poles essentially produce the various cusps and peaks in the propagator amplitude, where the channel couplings g (α) u,v and phase space ρ u,v (mass threshold) weighting factors shift and distort these features away from the original m α values. The corresponding Argand diagrams show similar behaviour as the transition amplitude T 11 shown in Fig. 10a), although in general they exhibit distortions due to the channel-dependent couplings and mass thresholds. In particular, the ππ → ππ propagator (u = 1, v = 1) exactly matches the phase motion of T 11 if we first rotate T 11 by 90 degrees anticlockwise (δ = 45 degrees) around the centre of the unitarity circle at (0, i 2 ) and then shift it by the translation ( 1 2 , − i 2 ). Discussing the features in Fig. 12 in detail, the first pole (m 1 = 0.651 GeV/c 2 ) produces the first cusps around 0.65-0.8 GeV/c 2 in all of the channels, except for 4π (v = 3) which has a coupling of zero. The second pole (m 2 = 1.2036 GeV/c 2 ) produces cusps at 1, 1.1 and 1.5 GeV/c 2 for the KK (v = 2), ηη (v = 4) and ηη (v = 5) channels, respectively, and also generates a very broad dip centred near 1.2 GeV/c 2 for the ππ (v = 1) channel. The third pole (m 3 = 1.55817 GeV/c 2 ) generates broad peaks near 1.55 GeV/c 2 , where the low mass end terminates in a cusp at the threshold of the given channel, except for the ηη mode where a narrow cusp at 1.5 GeV/c 2 is generated against a smoothly varying amplitude. The fourth pole (m 4 = 1.21 GeV/c 2 ) creates structures very similar to the second pole owing to their almost degenerate mass values, with an additional broad peak around 1.1-1.3 GeV/c 2 present for the 4π channel, while the fifth pole (m 5 = 1.82206 GeV/c 2 ) generates the broad peaks near 1.8 GeV/c 2 for all channels, with additional cusps at 1 (KK), 1.1 (ηη) and 1.5 GeV/c 2 (ηη ) that are very similar to those found for the second and fourth poles.   Figure 13 shows the pole production amplitudes A α,u=1 (s) defined in Eq. (74), which are formed by modulating the pole singularity term 1/(m 2 α − s) with a weighted sum of the s-dependent propagator distributions for all channels v = 1 to 5 shown in Fig. 12, along with β α ≡ 1. Concentrating on the magnitudes, we can see that the first pole (m 1 = 0.651 GeV/c 2 ) has an amplitude that begins rather flat from the ππ threshold until it starts to peak near 1 GeV/c 2 before rapidly falling at higher invariant mass. Even though the ππ propagator (ππ → ππ) is slowly decreasing from unity at the ππ threshold down to zero near 0.8 GeV/c 2 , the presence of both the rapid rise of the amplitude from the pole singularity at 0.65 GeV/c 2 as well as the increasingly influential KK propagator (KK → ππ), owing to its larger coupling constant, ensures that the amplitude remains fairly constant in the region below 1 GeV/c 2 . The propagators for all of the channels (except 4π) have a peak near 1 GeV/c 2 , and these combine to give the same local peak for the first production pole. As we increase the invariant mass, the falling shape of the pole singularity starts to dominate the amplitude modulation, and so we get a rapid reduction in the magnitude no matter what shapes the propagators have. The amplitude for the second production pole (m 2 = 1.2036 GeV/c 2 ) strongly depends upon the ππ and KK propagator shapes, where the former has a coupling constant almost double that of the latter. As we decrease the invariant mass from 1.2 GeV/c 2 , the pole amplitude would be very small at the strong KK dip at 1 GeV/c 2 if not for the compensating sharp peak in the ππ propagator. Likewise, the zero ππ propagator amplitude at 0.8 GeV/c 2 is nullified by the non-zero KK contribution. These two effects conspire to shift the location of the sharp dip in the production pole amplitude by 50 MeV/c 2 from 1 to 0.95 GeV/c 2 . As we decrease the invariant mass, the pole amplitude becomes more influenced by the rising ππ propagator until it starts to fall as we move further away from the pole mass. Above 1.2 GeV/c 2 the production amplitude tends to follow the undulations of the ππ and KK propagators, producing broad local peaks centred on 1.3 and 1.9 GeV/c 2 as well as a more narrow one at 1.5 GeV/c 2 . The modulation of the third pole (m 3 = 1.55817 GeV/c 2 ) amplitude is dominated by the 4π channel, although the other channels give significant contributions, varying from 33% (ηη) up to 66% (ππ). The 4π propagator has a very broad maximum centred very close to the pole mass, and so we would expect the production pole peak to remain very close to 1.56 GeV/c 2 (with an asymmetric width that is slightly narrower on the low side). However, as we decrease the invariant mass, the rising contributions from the ππ and KK propagators effectively shift the production peak by 60 MeV/c 2 down to 1.5 GeV/c 2 . Below 1.2 GeV/c 2 , the modulations from the ππ and KK propagators become washed out since they are too far away from the pole position. Above 1.5 GeV/c 2 , the width of the production amplitude remains very wide owing to the dominant 4π propagator. The fourth pole located at 1.21 GeV/c 2 is almost degenerate with the second pole (1.2036 GeV/c 2 ) and so we would naively expect them to have essentially identical production shapes. However, the coupling coefficients are completely different, where now the 4π channel propagator dominates, with a factor of two or more reduction in the other contributions, leading to the production of the two broad peaks centred around 1.4 and 1.7 GeV/c 2 . For invariant masses at 1 GeV/c 2 and below, the amplitude does indeed closely follow the shape of the second production pole since the 4π propagator becomes negligible and the ππ and KK contributions dominate. The fifth and last pole (m 2 = 1.82206 GeV/c 2 ) has an amplitude that is strongly influenced by the 4π channel, with much smaller contributions from the others. The large, broad 4π propagator maximum near 1.6 GeV/c 2 shifts the production peak by around 70 MeV/c 2 down to 1.75 GeV/c 2 , while the other smaller production peaks near 1.4 and 1.1 GeV/c 2 match those seen in the 4π shape. The magenta lines in Fig. 13 show the effect of multiplying the Adler zero suppression factor f A0 (s) to the production pole amplitudes, where we can see that it only reduces the magnitudes for invariant masses below 1 GeV/c 2 .
The mass distributions of the SVP production amplitudes A SVP,u=1,v (s) defined in Eq. (75) follows very closely the shape modulations of the propagator terms shown in  The final set of pedagogical plots are given in Fig. 14, which show the normalised mass projections of the individual production pole and SVP amplitudes, as well as the ρ(770) resonance for comparison, using events generated uniformly across the Dalitz plot. In general, we see that the peaking structures observed in Fig. 13 are replicated here, with a reduction in the intensity as the invariant mass approaches the ππ threshold (fewer events are generated at the kinematic boundary), with an almost flat intensity seen for invariant masses above 2 GeV/c 2 , which is the effective cut-off for the K-matrix parameterisation owing to the fact that there are no bare poles defined in this region. The first production pole generates the peak corresponding to the f 0 (980) along with a broad shoulder on the low mass side, which can be referred to as the f 0 (500) or σ resonance. The third pole generates the f 0 (1500) peak, the fifth pole is the main contributor for the f 0 (1710), while a combination of the second and fourth poles (which have almost degenerate bare masses) produces peaks in the f 0 (1370) region. All of these peaks are generated dynamically by the K-matrix amplitude. The first two SVPs have rather oscillatory shapes in the region around 1 to 2 GeV/c 2 ; when they are combined we can approximately obtain a very broad bump between the f 0 (980) and f 0 (1370) peak locations. The third SVP generates a large peak very near 1.6 GeV/c 2 , in between the f 0 (1500) and f 0 (1710) regions. Lastly, the fourth and fifth SVPs are essentially degenerate, peaking at the same position of the f 0 (980) from the first production pole. This means that we can ignore these two contributions, or at least remove the fifth SVP, since nothing is gained by their inclusion in the total amplitude description.

B Formulae for available angular distributions
The angular distributions and Blatt-Weisskopf form factors set out in Sec. 2.2 are the default settings in Laura ++ . However, other formalisms to describe the angular distributions are also implemented in the package and it is straightforward to switch between them. This appendix details these alternative formalisms and illustrates the few additional lines of code required to use them.
The four spin-factor formalisms are defined in the enumeration LauAbsResonance::LauSpinType, which can take the values Zemach P (the default setting), Zemach Pstar, Covariant, and Legendre. The simplest description of the spin factors is that of the Legendre formalism, where the spin factors are simply the Legendre polynomials (with some additional numerical constants in order to maintain consistency of the phase conventions among the various formalisms) The spin factors for Zemach P are those given in Eqs. (9)- (14), which differ from the expressions of Eqs. (85)-(90) by factors of (p q) L . Similarly, those for Zemach Pstar are the same as those for Zemach P but with the bachelor momentum evaluated in the rest frame of the parent particle (p * ), rather than that of the resonance (p). The angular distributions have been implemented in Laura ++ up to L = 5, which is two units larger than the maximum spin of any resonance observed to be produced in any Dalitz plot to date [87,88,124]. The angular distributions discussed above are based on a non-relativistic assumption. For certain channels, this may not be sufficiently precise, and therefore the Covariant formalism is also made available. This is given by The first three of these expressions are derived in Ref. [125] and, based on that work, the last two were derived in Ref. [124]. As can be seen from the expressions, the differences between formalisms are more significant for higher spin resonances, and particularly affect tails of the distributions. To give an idea of the effect, the lineshapes for the f 2 (1270) and ρ 3 (1690) 0 resonances decaying to π + π − are shown in Fig. 15.  Figure 15: Lineshapes for the (left) f 2 (1270) and (right) ρ 3 (1690) 0 resonances decaying to π + π − (in the B + → K + π + π − Dalitz plot) with the (blue) Legendre, (red) Zemach P, (green) Zemach Pstar and (magenta) Covariant spin formalisms. In all cases the relativistic Breit-Wigner description is used, with mass and width parameters as given in App. C and the two Blatt-Weisskopf factors set to unity.
It is possible to switch between these different formalisms via a function of the LauResonanceMaker factory object. For example, to use the Covariant formalism one would do: L a u R e s o n a n c e M a k e r & resMaker = L a u R e s o n a n c e M a k e r :: get () ; resMaker . s e t S p i n F o r m a l i s m ( L au A bs Re so n an ce :: Covariant ) ; It is important to note that any such operation must be performed prior to constructing any resonances, i.e. before calling LauIsobarDynamics::addResonance or LauIsobarDynamics::addIncoherentResonance for the first time.
As the angular and Blatt-Weisskopf factors are strongly coupled, it is also possible to straightforwardly modify the form of the Blatt-Weisskopf factors. In particular, the momentum value used for the factor that is related to the decay of the parent particle into the resonance and the bachelor can be selected from the following options (defined in the LauBlattWeisskopfFactor::RestFrame enumeration): • LauBlattWeisskopfFactor::ResonanceFrame, the momentum of the bachelor in the rest frame of the resonance, p (the default setting), • LauBlattWeisskopfFactor::ParentFrame, the momentum of the bachelor in the rest frame of the parent, p * , • LauBlattWeisskopfFactor::Covariant, the product of the momentum of the bachelor in the rest frame of the parent, p * , and a function of the ratio of the energy and mass of the resonance in the rest frame of the parent, 1 + p 2 /m 2 P . More precisely, this function is the expression in the middle term in Eqs. (91) to (95)  where in this example the momentum of the bachelor in the rest frame of the parent (p * ) is to be used. Again, this operation must be performed before constructing any resonances.
In addition, it is possible to change the form of the Blatt-Weisskopf factors, with the different types being defined by the LauBlattWeisskopfFactor::BarrierType enumeration.
The default setting, corresponding to Eqs. (15)- (20), is given by LauBlattWeisskopfFactor::BWPrimeBarrier and is recommended when the angular terms contain momentum factors.
An exponential form for these factors, LauBlattWeisskopfFactor::ExpBarrier, which has been used in some analyses for virtual contributions has also been implemented, To change the form of the barrier factors for all resonances, the following lines are required L a u R e s o n a n c e M a k e r & resMaker = L a u R e s o n a n c e M a k e r :: get () ; resMaker . setBWType ( L a u B l a t t W e i s s k o p f F a c t o r :: BWBarrier ) ; where in this example the forms in Eqs. (96)-(101) are to be used. Again, this operation should be performed before constructing any resonances. As for the T ( p, q) terms, the differences between Blatt-Weisskopf form factor formalisms are more significant for higher spin resonances, and far from the peak of the resonance. An illustrative comparison of the shapes is given in Fig. 16.  Figure 16: Lineshapes for the (left) f 2 (1270) and (right) ρ 3 (1690) 0 resonances decaying to π + π − (in the B + → K + π + π − Dalitz plot) with (blue) no Blatt-Weisskopf factors, and with the (red) ResonanceFrame, (green) ParentFrame and (magenta) Covariant settings for evaluating the momentum that enters the Blatt-Weisskopf factor associated with the decay of the parent particle. In all cases the relativistic Breit-Wigner description is used, with mass and width parameters as given in App. C and the Zemach P formalism for the spin factors.
It is possible to make all of the changes discussed in this Appendix at the level of individual resonances, using the functions LauAbsResonance::setSpinType and LauAbsResonance::setBarrierRadii, but this requires much care to be taken and is not generally recommended.

C Standard resonances
This section provides the complete set of available resonances, indicating the name, mass m 0 , width Γ 0 , spin, charge and Blatt-Weisskopf barrier radius r R BW . Table C1 contains  information for light meson resonances, Table C2 for charm, charmonium, strange-charm, beauty and strange-beauty resonances, Table C3 for K * resonances and Table C4 for nonresonant terms. Most data are taken from Ref. [62]. The tables list the information contained in the information records for both neutral and positively-charged resonances. Negatively-charged resonance records are implemented as charge-conjugates of the positively charged ones; the plus sign in the name is replaced with a minus sign.
In case a user wishes to modify the values of the parameters from those given in the tables, the LauAbsResonance::changeResonance function, which takes the mass, width and spin as arguments, can be used. The Blatt-Weisskopf barrier radius can be changed with the LauAbsResonance::changeBWBarrierRadii function, and other parameters specific to particular lineshapes can be changed with the LauAbsResonance::setResonanceParameter function. The same approach can be used to include a resonance that is not available in these tables, by using any of the existing states of appropriate charge and redefining its properties.

D PDF classes
This section details the formulae within classes that can be used to parameterise additional PDFs P(x; p 1 , p 2 , ..., p n ) for the likelihood function. Here, x denotes the dependent variable, while p 1 , p 2 , ..., p n is the list of parameters in the form of a vector of LauParameter objects, each containing a descriptive name (which must contain the case-sensitive word shown in quotes), the value of the parameter and optionally its validity range, uncertainty and constantness. All PDFs used in Laura ++ are normalised to unity, although the normalisation factors are omitted in many equations in this section for brevity.

D.1 LauArgusPdf
The ARGUS threshold function [126] can be used to parameterise the shape of combinatorial or partially-reconstructed backgrounds of the invariant mass m of parent candidates: where x = m/m 0 , m 0 is the end-point of the curve ("m0"), ξ is the shape parameter ("xi"), while θ(x ≤ m 0 ) = 1 and θ(x > m 0 ) = 0.

D.5 LauCrystalBallPdf
The Crystal Ball function [127] is a PDF that contains a Gaussian core with a continuous power-law tail on one side: where µ and σ are the Gaussian "mean" and width ("sigma"), respectively, α ("alpha") is the positive or negative distance from the mean in which the Gaussian and the tail parts match up, n ("order") is the power exponent for the tail, while t is equal to (x − µ)/σ, which changes sign if α is negative.

D.6 LauExponentialPdf
This exponential function is simply given by where λ is the "slope" parameter.

D.8 LauLinearPdf
This linear function only needs the gradient ("slope") λ: where c is the intercept and is evaluated using the range of the abscissa x.

D.10 LauParametricStepFuncPdf
This parametric step function is a binned distribution whose parameters are the contents of each bin (except one), essentially representing a histogram with variable bin content. The content of the remaining bin is determined from that of the others and the requirement of normalisation. The constructor requires two vectors of LauParameter objects; the first stores the bin contents or weights (which are parameters that can be fitted), while the second stores the lower edge abscissa limits of each bin (in ascending order) as well as the upper edge limit of the last bin. It can also be specified whether the normalisation bin is the first or last; it is advisable to use the bin with the larger content.

D.11 LauSigmoidPdf
The sigmoid function follows an "S" shape and is defined as with parameters "a" and "b" defining the steepness of the slope and the shift of the distribution, respectively. A negative value of a will flip the distribution around the ordinate axis.

D.12 LauSumPdf
This class implements the summation of two PDFs, P 1 and P 2 , with a relative fraction ("frac") f :

D.13 Dalitz-plot-dependent PDFs
What follows are descriptions of PDFs with parameters that can vary across the Dalitz plot via power-law scaling with Laurent polynomials containing either positive or negative exponents. Here, the parameterisation of a given PDF parameter g is given by where D is the invariant mass-squared variable representing the DP position and c i is the coefficient of expansion for the power term i. The variable D can be either m 2 13 , m 2 23 , or m 2 12 , the minimum or maximum values of m 2 13 or m 2 23 , or the distance from the DP centre. The constructors of these PDFs require vectors of the coefficients c i for each function parameter g, as well as the pointer to the LauDaughters object in order to find D from the kinematics. The LauDPDepBifurGaussPdf, LauDPDepCruijffPdf and LauDPDepGaussPdf classes represent bifurcated Gaussian, Cruijff and normal Gaussian PDFs, given in Eqs. (104), (107) and (110) respectively, with parameters that can vary according to Eq. (115).
The LauDPDepMapPdf class can be used to define a PDF that requires different functions depending on the DP region. It uses a Root histogram that divides the DP, or its projection onto one axis, into ascending, numbered regions. The region number (starting from zero) for a given value of D is then used to choose the corresponding PDF from the ordered list of functions provided as a vector in the constructor.
The LauDPDepSumPdf class implements the sum of two PDFs, defined by Eq. (114), in which the fraction f ("frac") depends on the variable D, using either the contents from a two-dimensional histogram directly or a vector of coefficients c i for the positive-power Laurent polynomial shown in Eq. (115).