KinFit A Kinematic Fitting Package for Hadron Physics Experiments

A kinematic fitting package, KinFit , based on the Lagrange multiplier technique has been implemented for generic hadron physics experiments. It is particularly suitable for experiments where the interaction point is unknown, such as experiments with extended target volumes. The KinFit package includes vertex finding tools and fitting with kinematic constraints, such as mass hypothesis and four-momentum conservation, as well as combinations of these constraints. The new package is distributed as an open source software via GitHub. This paper presents a comprehensive description of the KinFit package and its features, as well as a benchmark study using Monte Carlo simulations of the pp → pK + Λ → pK + pπ − reaction. The results show that KinFit improves the parameter resolution and provides an excellent basis for event selection.


Introduction
Kinematic fitting is a powerful tool widely used in particle and nuclear physics analyses, recognised for its ability to improve the resolution of measured particle track parameters, suppress background, reconstruct undetected particles and to determine the position of vertices.Available information from measurements, such as momenta, angles and energy, combined with physics constraints, such as four-momentum conservation in a production or displaced decay vertex, or the mass of an undetected unstable particle through its decay products, are exploited.With this information, a mathematical minimization problem can be formulated and solved using e.g., the Lagrange multiplier technique.
Kinematic fitting techniques have been available in particle physics since the 1960's [1] and gained momentum during the bubble chamber days.However, existing packages like RAVE [2] from ILD or the full decay chain fitters from Belle II [3] or PANDA [4] are often embedded in the respective experiment software and specialized for the detector setup.In the measurement of complex decays of heavy particles, it is beneficial apply kinematic tree fits [3] to complex decay chains or jets [2,5].
In many hadron physics experiments like the HADES experiment [6] at GSI, the interaction point is not fixed to a point but can be located within a target volume that covers several centimeters of the path along the beamline.Hence, the interaction vertex position needs to be determined from the measured track parameters of the outgoing particles.In addition, weakly decaying particles such as hyperons produce daughter particles in a secondary vertex, located a measurable distance away from the interaction point.These are crucial to reconstruct since many hyperons, e.g.Λ and Σ 0 , are neutral and could otherwise escape detection.Figure 1 illustrates an example where a Λ hyperon is produced in the target volume and subsequently decays in a secondary vertex after travelling some distance.Tracks from secondary particles are more challenging to reconstruct and hence, the reconstructed track parameters have larger measurement uncertainties compared to those from primary particles.All this calls for kinematic fitting.A new package, KinFit, has been developed to facilitate the hyperon program at HADES [7] but is provided in an experiment-agnostic way and is therefore applicable also to other experiments.Hadron physics experiments at J-PARC [8], JLab [9] or COM-PASS [10] and AMBER at CERN, may profit from this package.
In addition to kinematic fitting, tools are provided to construct a particle candidate from its daughter tracks.A functionality for running a kinematic fit for event selection on a ROOT [11] file in an automated way is included as well.
This paper is outlined as follows: The methodology including the fitting procedure and available constraint equations are described in Section 2. In Section 3, the classes contained in the package are described and in Section 4, a benchmark study is provided for Λ hyperon reconstruction to demonstrate the performance of the fitter.In addition, a comprehensive user's guide is provided in Appendix A.

Methodology
The goal of a kinematic fitting algorithm is to find an improved set of track parameters as close as possible to the true values, such that they fulfill a set of kinematic and/or geometric constraints.This concept can be translated quantitatively into a χ 2 minimization expressed as follows: where ⃗ y is a vector of the N ∈ N measured track parameters provided by the tracking algorithm, V (⃗ y) is the corresponding covariance matrix and ⃗ η are the estimated track parameters.
In addition, the kinematic and geometric constraints are implemented in the procedure by constraint equations ⃗ f .These are in general represented by a set of K ∈ N continuously differentiable functions of the estimated parameters ⃗ η and a set of J ∈ N 0 unmeasured parameters combined in ⃗ ξ: The objective is to minimize the χ 2 from Equation (1) subject to these constraints.

The Iterative Lagrange Multiplier Method
The method of Lagrange multipliers [12] is wellsuited for addressing such problems.This technique enables the transformation of the constrained minimization into the minimization of a single Lagrange function, L: where the K additional variables summarized in ⃗ λ are introduced, referred to as Lagrange multipliers.
Minimizing the Lagrange function (Equation ( 2)) involves finding the derivatives of L with respect to all unknowns ⃗ η, ⃗ ξ, ⃗ λ.By setting these derivatives to zero and subsequently solving for ⃗ η and ⃗ ξ, the minimization can be achieved.As this problem is generally non-linear, an iterative procedure is employed, with each iteration yielding improved approximations for ⃗ η and ⃗ ξ.Assuming the values of all quantities of the iteration ν have already been extracted, we want to express the quantities of the subsequent iteration ν + 1 in terms of the values of iteration ν.We then proceed as follows [12]: 1. First, the following notations are introduced where the derivative of the constraint equations with respect to ⃗ η. 2. The updated, unmeasured variables, ⃗ ξ ν+1 , are obtained by where F ξ is a K × J Jacobian matrix (F ξ = D ξ ⃗ f ), i.e. the derivative of the constraint equations with respect to ⃗ ξ. 3. The updated Lagrange multipliers, ⃗ λ ν+1 , are obtained from 4. The updated measured parameters, ⃗ η ν+1 , are calculated as 5. Finally, the new L is calculated and the results compared with the previous iteration.
To decide when the solution is sufficiently close to the minimum, convergence criteria are defined, see Section 2.4.
In the end, the new covariance matrix, V ν+1 , is calculated This ansatz assumes a quadratic minimum of the χ 2 function.The function will look different if it deviates too much from its minimum.Currently there is no guard against this but the user should be mindful of potential deviations and assess the validity of the results accordingly.

Track Parametrization
The input objects for KinFit include track parameters such as momentum and angles, defined in polar coordinates, as well as the point of closest approach to the beam axis.Functions for converting from Cartesian coordinates are provided to ensure compatibility with other experiments (see the User's Guide in Appendix A).In addition, the covariance matrix is required as input.The particles are assumed to move in a region free of magnetic fields, meaning they are propagated as straight tracks.The parameters employed here to uniquely describe a track include: • Polar angle θ [radians] defined from 0 to π.
• Azimuthal angle ϕ [radians] defined from −π to π relative to the beam (z) axis.• R [mm], the distance between the beam axis and the point of closest approach of the track to the beam axis.• Z [mm], the z coordinate of the point of closest approach of the track to the beam axis.

and
• Covariance matrix, where the diagonal entries correspond to the uncertainties in the parameters 1/p, θ, ϕ, R and Z.
For the purely kinematic fits, only the track parameters 1/p, θ and ϕ are required.

Constraints
KinFit offers a variety of constraints that can be selected for kinematic fitting.

1C: Vertex Constraint
This purely geometrical vertex constraint minimizes the distance between two tracks, ensuring they originate from the same vertex.A straight line in 3D is uniquely defined by a base vector that points from the origin of a chosen coordinate system to a coordinate on the line and the direction of the line.The components of these vectors are: and (5) In the base vector, π/2 is added to ϕ to ensure the base vector is constructed in the HADES coordinate system.A suitable constraint equation can then be formulated as follows:

) This expression is proportional to the minimum distance between the straight lines parameterized by the respective base and direction vectors
The vertex fit can be performed either separately or as part of a series of consecutive fitting procedures.In addition, there is the possibility to perform a mass and a geometrical vertex fit simultaneously.

4C: Four-momentum conservation in the beam-target interaction
The 4C constraint demands that the sum of the four-momenta of the final state particles equals that of the initial beam-target system (in the following denoted with the suffix ini).Since all parameters are treated as measured parameters with uncertainties, there are four over-constraints.
For N particles, the constraint equations are

3C: Four-momentum conservation in a displaced vertex
The 3C constraint utilizes four-momentum conservation at a given decay vertex, where a mother particle (M) decays weakly into N particles.This procedure relies on measured information about the primary vertex and the decay vertex, and on a mass hypothesis of the mother particle as e.g.Λ or K s .The angles of the mother particle, θ M and ϕ M , are determined by the direction of the vector pointing from the primary to the decay vertex and are therefore treated as measured parameters, whereas the momentum p M is unmeasured and hence obtained from the fit.This results in three over-constraints (3C).The initial value is estimated from energy conservation of the decay products using In this expression, the subscript M represents the mother particle, while i refers to the decay products.The user must provide the mass of the mother particle as a mass hypothesis.The uncertainties in the vertex positions must estimated in the previous step and are propagated to the uncertainties in the parameters.Since the momentum of the mother particle p M is an unmeasured quantity, its uncertainties are not known a priori but are, as the momentum itself, obtained from the fit.The constraint equations ensure four-momentum conservation in the decay of the mother:

1C: Four-momentum conservation with a missing particle
Four-momentum conservation at the interaction point enables the reconstruction of one undetected or missing (indicated by the suffix miss) particle with a mass hypothesis m miss .I all other inital and final state particles are known or measured, the four-momentum conservation results in the following constraint equations: Here one has measured particles with indices i, a missing particle with mass m miss and the initial beam-target four vector p µ ini .The mass, m miss needs to be provided by the user as a hypothesis as well as p µ ini .The four-momentum of the missing particle can be extracted after the fit.Since there are three unmeasured variables, i.e. the three-momentum of the missing particle, only one over-constraint (1C) remains.

1C: Missing Mass Constraint
The missing mass constraint is suitable in the case when there is one undetected final state particle whose momentum is not of interest.Instead, a missing mass constraint can be formulated from the mass hypothesis of the missing particle:

1C: Invariant Mass Constraint
The mass constraint can be used to constrain a set of particles to originate from a common mother particle, whose mass is known.The invariant mass of this set of particles is then constrained to the mass of the hypothetical mother particle:

Convergence
The iterative fitting procedure is terminated when convergence is achieved.There are three different convergence criteria: The difference in χ 2 of the Lagrange function between consecutive iterations, the Euclidean norm of the sum of all constraint equations and the Euclidean norm of the difference of all track parameters, normalized to the initial measurement, between two consecutive iterations.The default value for the convergence criteria is 10 −4 but can be customized.If convergence is not reached before a specified number of iterations (default 20), the fitting procedure is exited.

Goodness of fit
The performance of the fit is assessed by the final χ 2 f inal of the fit and the so-called pull distributions for all fitted variables.The value of χ 2 f inal should be small, on the order of the number of overconstraints, N C .For an ensemble of events, N C defines the shape of the χ 2 f inal distribution.There is a one-to-one relation between χ 2 f inal for a given N C and the corresponding probability, as defined by the probability density function [13].Ideally, the probability distribution should be uniform if the correct particle hypotheses have been used in the fit and if the covariance matrices accurately describe the measurement precision.However, if the particle hypotheses are incorrect, then a peak towards zero probability should be discernible.If elements of the covariance matrix are over-or underestimated, the distribution will decrease or increase towards larger values of probability.
The pull distributions are defined as: for each variable, where η/y are the fitted/measured variables, respectively, and σ(η)/σ(y) are the corresponding uncertainties.The pull distribution should follow a normal distribution with a mean value of 0 and a standard deviation of 1.

Class Descriptions
The KinFit package is written in C++ and based on ROOT [11] (version 6) and uses CMake [14] (version 3.0 or newer) for the installation.It is available online at: https : // github .com / KinFit / KinFit .git Documentation about the usage of the provided tools can be found in the README of the git repository and in the User's Guide in Appendix A. The properties of the particle candidates needed as input are the particle's track parameters, mass, and covariance matrix.Functions are provided to transform Cartesian track parameters to the R, Z track parameters.

KFitParticle
The track parameters of particle candidates are organized in KFitParticle objects that provide all information needed by the fitter class.It inherits from the ROOT class TLorentzVector.For each candidate, the values for the track parameters described in Section 2.2 have to be set.Optionally, an arbitrary particle ID and track ID can be chosen by the user for later reference.

KFitDecayCandFinder
KFitDecayCandFinder calculates properties of an unmeasured candidate that decays in a displaced vertex to be used later in the fit.The angles, θ and ϕ are calculated from the line segment connecting the primary to the displaced decay vertex.The uncertainties for these angles are propagated from the uncertainties in the vertex positions in X, Y and Z using the matrix formalism for error propagation.The user needs to provide the uncertainties for the vertex positions.

KFitVertexFinder
KFitVertexFinder finds the vertex by calculating the point of closest approach between at least two tracks by a matrix formalism.For more than two tracks, the vertex is taken as the center of gravity, i.e. the point that is simultaneously closest to all tracks..This is calculated from a least square method.The vertex finding code is based on a procedure in HYDRA1 , the HADES software.

KinFitter
KinFitter is the main class containing the fitting functions.It can perform a vertex, 4C, 3C, missing particle, missing mass, mass and a combined vertex+mass fit (see Section 2.3).The constraint equations and Jacobi matrices are implemented here.This is where the iterative fitting procedure is carried out.The maximum number of iterations or convergence criteria can be changed to custom values.

KFitAnalyzer
KFitAnalyzer is a user interface class.It contains user settings and performs the event loop with a fit of choice.The input particles need to be stored as KFitParticles in a TClonesArray .For each event, particles are selected based on their particle ID (PID), which is also provided to the KFitAnalyzer.A KFitDecayBuilder object is created which takes the selected particles and the desired constraint as input.In the end the fitted particles are retrieved from the KFitDecayBuilder and stored in an output file together with the fit probability.The KFitAnalyzer currently performs all fits except the 3C fit.

KFitDecayBuilder
This class is responsible for building all possible combinations of the particles within one event.Each combination is passed to the KinFitter which performs the selected fit.The KFitDecayBuilder selects the combination with the best fit probability.

Performance Study
The performance of KinFit is benchmarked using Monte Carlo (MC) simulations.These represent an ideal scenario where a finite detector resolution is added by smearing, and thus known, whereas in real data, it might be difficult to estimate the covariance matrix exactly.This approach allows for a focused investigation of the quality of the KinFit tools.The reaction investigated is pp → pK + Λ, Λ → pπ − , as illustrated in Figure 1.This reaction provides an opportunity to examine all tools contained in the package.All fits performed in this section use the default settings of KinFitter, described in Section 2.4.
The event generator Pluto [15] was used to generate 100 000 events of the pp → pK + Λ Λ → pπ − reaction at a beam kinetic energy of T=4.5 GeV.In our simulations, the produced particles and their decay products are distributed isotropically across the available phase space.The primary vertex is generated at the origin, i.e. (0,0,0).A full 4π acceptance is assumed and no material effects are included.The Λ hyperons decay at a displaced decay vertex according to the mean Λ hyperon life time c τ (Λ) = 7.98 cm [13].The track parameters of the final state particles are smeared according to a Gaussian distribution, simulating uncertainties from e.g., detector resolution in a controlled way.The uncertainty for each track parameter is listed in Table 1.Since the results presented here do not depend on a specific experiment, the uncertainties are set to be constant, except the momentum uncertainty, which is momentum dependent.

Mass Fit using KFitAnalyzer
In each event, there are four measured particles: a proton (p 1 ) and a kaon from the interaction point and a pion and a proton (p 2 ) from the Λ hyperon decay.This means there are two possibilities to combine a pion with a proton, either p 1 π − and p 2 π − , of which the latter is the correct Λ decay candidate.The KFitAnalyzer is used to find the pion-proton combination that comes from the Λ hyperon decay through the application of a mass fit.The pion-proton pair giving the largest mass fit probability (see Fig. 2) are selected as Λ daughters.In 99 % of the events, the correct proton p 2 is chosen.The probability distribution is uniform, as expected.

Reconstruction of Λ Hyperon from Vertex Positions and a 3C Fit
The Λ hyperon reconstruction procedure includes three steps: 1. Vertex finding As a first step the position of the primary vertex, where the Λ hyperon was produced is determined from the point of closest approach between the other tracks coming from the production vertex.In this case, these are the proton (p 1 ) and kaon tracks.This task is performed by KFitVertexFinder.In a similar way, the Λ hyperon decay vertex is obtained from the point of closest approach between the proton (p 2 ) and pion tracks originating   from the Λ hyperon decay.The vertex positions are shown in Fig. 3.All primary particles are produced at the origin which means that the distribution in the top panel reflects the finite resolution of the vertex estimation.The dominating effect in the Λ decay vertex distribution is the distance travelled by Λ before decaying.

Reconstruction of the neutral candidate
The direction of the Λ candidate is given by the vector pointing from the primary to the decay vertex.The magnitude of the Λ hyperon momentum is estimated by the of the sum of the three-momenta of its daughter particles, see Equation ( 8).The KFitDecayCandFinder carries out this task and calculates the uncertainties of the track parameters of the Λ candidate.The latter requires information on the vertex uncertainties which were estimated by Gaussian fits to the vertex resolution in x, y and z for each vertex.

3C fit
A kinematic fit ensuring four-momentum conservation at the Λ hyperon decay vertex is the final step and is executed by KinFitter.This improves the resolution of the Λ hyperon track parameters significantly (Figure 4) and can be used to reject false combinations of protons in the vertices.
Only the correct combination of particles was used in this example, i.e. correct assignment of the protons to the vertices.The pull distributions for the pion track are shown as an example in Fig. 5.These nearly follow a normal distribution, as expected.Fig. 6 shows the probability distribution of the 3C fit.The probability deviates slightly from a uniform distribution and there are some events for which the fit converged, but gives a low probability.This effect and the slightly larger deviation of the pull distributions from a normal distribution are only seen in the 3C fit and originate from the estimation of the vertex resolutions, which are not exactly described by a Gaussian.

4C Fit of the Reaction
pp → pK + Λ The 4C fit is used to constrain the four-momenta of all final state particles to that of the initial beam-target system.In this example, the final state particles are a proton and a kaon from the primary vertex and a proton and a pion from the Λ hyperon decay vertex.Fig. 7 shows the momentum resolution for the kaon and the pion before and after the fit.The improvement is more substantial for the kaon momentum resolution compared to the other particles.This is due to the larger uncertainty of the kaon associated with its larger total momentum.The probability distribution is uniform and the pull distributions follow a normal distribution, as shown in Figure 8.In this scenario, it is assumed that the kaon is not detected.The missing particle fit is employed to constrain the four-momenta of the detected particles along with the undetected particle to match the beam-target system, as well as to estimate the momentum of the missing particle.Figure 9 presents the probability distribution of the fit and compares the kaon momentum resolution from the initial guess and after the fit.The initial guess for the kaon momentum is calculated from threemomentum conservation in the interaction point.The momentum resolution before the fit is worse than in the example in Section 4.3, where the momentum was measured directly.A noticeable improvement in the momentum resolution can be observed after the fit.However, the resolution is still worse than in the case of the 4C fit, as the missing particle fit has fewer over-constraints.The resolutions of the other track parameters and the pull distributions appear quite similar to those of the 4C fit.

Vertex Fit
The vertex fit aims to constrain the track parameters of final state particles that originate from the same single point in space.To perform the vertex fit, at least two outgoing particles from the same vertex must be measured.In this example, these particles are the kaon and proton (p 1 ) produced at the beam-target interaction point.Fig. 10 illustrates the resolution of the estimated R-parameter, defined in Section 2.2, for the proton.The probability distribution for this vertex fit, along with an example pull distribution of the R-parameter of p 1 , is depicted in Fig. 11.As anticipated, the probability is uniformly distributed across its range, while the pull follows a normal distribution.The KinFitAnalyzer was used to choose the proton track that has the largest probability to come from the same vertex as the kaon track.The correct proton (p 1 ) is chosen in 85 % of the events, which is lower than the fraction of correctly identified combinations by the mass fit presented in Section 4.1.This shows that the mass constraint is more powerful than the vertex fit.

Runtime Performance
The runtime performance of various components of KinFit was evaluated on a laptop with the following hardware specifications: • Processor: Intel Core i7-1185G7, 3.00 GHz • Memory: 32 GB RAM The results are displayed in Table 2.The average runtime per 1000 events was calculated for those in which the fit converged.A low overhead of approximately a microsecond was observed.

Conclusions
The KinFit package is a versatile kinematic fitting package that has been developed to provide tools for reconstruction of vertices, momenta and masses of particles in a generic hadron physics experiment.In particular, it is capable of fitting

Stage
Fitting Initialization Vertex Finding Mother Candidate Finding Time / event [µs] 10 1 1 1    an unknown production vertex from an extended beam-target interaction volume, and of combining kinematic and geometric constraints.The tools provided by KinFit have been tested with toy MC simulations, demonstrating outstanding performance in terms of improved track parameter resolutions and the ability to select correct particle hypotheses and -combinations.
While the focus has been on providing suitable tools for hadronic interactions and specifically hyperon decays, the package can be applied for other types of reactions as well.Though KinFit uses the same track parametrization as HADES, but it is generally considered experimentindependent.
The minimal overhead introduced by adding KinFit to an analysis suggests that the package is suitable for running on a local machine.

Outlook
Although the majority of particle physics analyses are based on C++, the demand for similar Python tools is currently growing due to the rapid expansion of the Python ecosystem.Therefore, the KinFit package is planned to be integrated into the SciKit-HEP project [16], which is a community-driven Python ecosystem for data analysis.candidates from each.

Handling the KinFitter
To use fitting functions, the KFitParticle objects that will be used need to be placed in a vector of the type std::vector from the C++ STD library.This vector needs to be passed to the constructor of KinFitter.How this is done is illustrated in the examples below.For most fit options, an arbitrary number of particles can be added.However, for the vertex fit, exactly two particles need to be added.Following this, the user must choose one of the constraints, either; the 3C, the 4C, mass, missing mass, the vertex or the missing particle.Adjusting the number of iterations as well as the convergence criteria is possible but optional; if not set, the default criteria are applied.Then the fit() function is called.This performs the actual fitting and returns true if the fit converged.Some information can be obtained independently of the individual fit function that is used, e.g.getChi2() returns the χ 2 of the fit and getProb() returns the corresponding probability.getPull(int val) returns the pull of variable v for a daughter particle p, val = 5 • d + v.The function isConverged() returns a boolean; true if the fit has converged and false otherwise, this can be used for event or sample selections.The initialization of KinFitter settings that can be used before calling fit() are summarized below.

2 IPFig. 1
Fig.1Illustration of a Λ hyperon produced by a proton beam hitting a proton target.The position of the interaction point (IP, marked in blue) and the Λ decay vertex (marked in green) are determined by the point of closest approach of the p 1 and K and p 2 and π − tracks, respectively.

Fig. 3
Fig. 3 Reconstructed production vertex from proton and K + (top) and location of reconstructed decay vertices of the Λ hyperon from the secondary proton and π − (bottom) as reconstructed by the KFitVertexFinder.The positions are projected onto the x-z plane.

Fig. 4
Fig. 4 Resolution of the Λ hyperon track parameters before and after the 3C fit.

Fig. 5
Fig. 5 Pull distributions for the π − track parameters after the 3C fit with respective mean and standard deviation.

4. 4
Missing Particle Fit of K + in the Reaction pp → pK + Λ.

Fig. 7
Fig. 7 Momentum resolution for the K + (top) and π − (bottom) before and after the 4C fit.

Fig. 8
Fig. 8 4C fit.Probability distribution (top) and example pull distribution of the π − momentum with mean and standard deviation (bottom).

Fig. 9
Fig. 9 Missing particle fit.Probability distribution (top) and resolution of the estimated and fitted K + momentum (bottom).

Fig. 10
Fig. 10 Vertex fit.Resolution of the R-parameter of the primary proton (p 1 ) before and after the fit.

Fig. 11
Fig. 11 Vertex fit.Probability distribution (top) and example pull distribution of the R-parameter for the primary proton (p1) (bottom).

Fig. 12
Fig.12Pull distributions for the track parameters of all particles after the 4C fit with respective mean and standard deviation.

Fig. 14
Fig.14Pull distributions for the track parameters of all particles after the 3C fit with respective mean and standard deviation.

Fig. 15
Fig.15Pull distributions for the track parameters of the primary proton tracks after the vertex fit with respective mean and standard deviation.

Table 1
Gaussian standard deviations of the track parameters in the MC sample after smearing.

Table 2
Runtime per event for different parts of KinFit.