FUMILI-based minimization with constraints using method of elimination of differentials

Kinematic fitting is one of the popular particle physics problems where constraint minimization is used. The constraints setting additional relations between parameters p can be given in form of equations φ(p1, . . . , pn) = 0. Often these equations are non-linear and complicated, and thus it is impossible or impractical to eliminate redundant parameters directly. The article covers employing of the minimization approach called a method of elimination of differentials, that is being developed at JINR as an extension to the FUMILI minimizer, and is intended for kinematic fitting in particle physics experiments. 1 Kinematic fitting A track reconstruction problem in general can be understood as a problem of multidimensional linear regression. Here, for every track j we are looking for a such initial momenta pj that the reconstructed coordinates of detector hits ĉ j i (p j), where i is a hit number, would be as close as possible to the registered hit coordinates c j i , taking into account the known detector uncertainties δc j i . In the other words, the task to be solved here is a minimization of the functional: F(pj, c j i ) = ∑ i c j i − ĉ j i (pj) δc j i 2 . (1) Further we would omit the indices i, j for simplicity, so that c and p would be vectors of all hits and momentum components correspondingly. Often one could get an additional information from reaction kinematics in a form of the energy-momentum conservation laws ∑ Pinitial = ∑ Pfinal, (2) or, in case of one of the final particles is not registered, ∣∣∣∣∑Pinitial −∑Pfinal∣∣∣∣2 = M2 X , (3) where P is a four-momentum of a particle, and MX is a reaction missing mass. Using (2), (3) allows us to narrow down the multidimensional domain of functional (1) and improves the precision of the result. A procedure that takes the kinematics constraints ∗e-mail: tokareva.vict@gmail.com © The Authors, published by EDP Sciences. This is an open access article distributed under the terms of the Creative Commons Attribution License 4.0 (http://creativecommons.org/licenses/by/4.0/). EPJ Web of Conferences 201, 07001 (2019) https://doi.org/10.1051/epjconf/201920107001 AYSS-2018 into account is called kinematic fitting. This method is widely employed in high energy physics, e.g. in experiments like BES-III [1], CLAS [2], CMS [3]. In this paper we are going to cover a FUMILI-based framework for kinematic fitting we have developed, which is currently used for data analysis in the ANKE experiment (Jülich, Germany) [4]. 2 Constrained minimization The conservation laws (2–3) could be expressed as φk(p) = 0, 1 ≤ k ≤ m, (4) where m is a number of constraints. Then (1), (4) is a constraint minimization problem written in the augmented form. The cases of interest are those where system (1), (4) is non-linear and could not be practically solved in a straightforward way, with taking in account practical time and memory limitations. The first ideas on overcoming this issue were suggested in several works starting from early sixties. Solmitz and Berge proposed to use the method of Lagrange multipliers[5], which is still widely used in problems of this type [6]. Approximately at the same time the penalty method for minimization was proposed at JINR by Moroz [7]. Later, one of the coauthors of this work suggested the method of elimination of differentials for kinematic fitting with constraints [8]. 2.1 Method of elimination of differentials In a neighborhood of a fixed point p∗ functional (1) can be represented by a power series F(p∗ + ∆p) = F(p∗) + n ∑ i=1 ∂F(p∗) ∂pi ∆pi + 1 2 n ∑


Kinematic fitting
A track reconstruction problem in general can be understood as a problem of multidimensional linear regression. Here, for every track j we are looking for a such initial momenta p j that the reconstructed coordinates of detector hitsĉ j i ( p j ), where i is a hit number, would be as close as possible to the registered hit coordinates c j i , taking into account the known detector uncertainties δc j i . In the other words, the task to be solved here is a minimization of the functional: (1) Further we would omit the indices i, j for simplicity, so that c and p would be vectors of all hits and momentum components correspondingly. Often one could get an additional information from reaction kinematics in a form of the energy-momentum conservation laws or, in case of one of the final particles is not registered, where P is a four-momentum of a particle, and M X is a reaction missing mass. Using (2), (3) allows us to narrow down the multidimensional domain of functional (1) and improves the precision of the result. A procedure that takes the kinematics constraints into account is called kinematic fitting. This method is widely employed in high energy physics, e.g. in experiments like BES-III [1], CLAS [2], CMS [3].
In this paper we are going to cover a FUMILI-based framework for kinematic fitting we have developed, which is currently used for data analysis in the ANKE experiment (Jülich, Germany) [4].

Constrained minimization
The conservation laws (2-3) could be expressed as where m is a number of constraints. Then (1), (4) is a constraint minimization problem written in the augmented form.
The cases of interest are those where system (1), (4) is non-linear and could not be practically solved in a straightforward way, with taking in account practical time and memory limitations.
The first ideas on overcoming this issue were suggested in several works starting from early sixties. Solmitz and Berge proposed to use the method of Lagrange multipliers [5], which is still widely used in problems of this type [6]. Approximately at the same time the penalty method for minimization was proposed at JINR by Moroz [7]. Later, one of the coauthors of this work suggested the method of elimination of differentials for kinematic fitting with constraints [8].

Method of elimination of differentials
In a neighborhood of a fixed point p * functional (1) can be represented by a power series and constraints (4) by Equations (5-6) can be rewritten in the matrix form as: where G( p * ) is the gradient, and Z( p * ) is the Hesse matrix. Here the rectangular matrix D( p * ) has m rows and n columns, since we have n parameters and m constraints. The vector ∆ p could be split into ∆p c that has m components, and ∆p f that has n − m components; the same could be done with the matrix D. Then (8) could be rewritten as Since (9) is a linear equation, it could be solved in relation to ∆p c , resulting in Returning to initial functional (7), we could employ (10) to eliminate the sub-vector ∆p c , and obtain a similar functional, in contrast depending on only n − m increments ∆p f :

FUMILI-based realization of the method
FUMILI, a package for minimizing χ 2 -like functionals, is one of the first minimizers included into ROOT release. It has been showing it's reliability, stability and high convergence rate while it had been being used by scientific community for decades. The greedy minimization algorithm which is employed in FUMILI was first proposed at JINR by S.N. Sokolov [9] and developed by I.N. Silin and V.S. Kurbatov [10]. It linearizes the second derivatives, so the error matrix is always positively defined, and thus each step leads to a minimum. Since FUMILI is designed specifically to minimize χ 2 -like functionals, it finds a minimum faster and requires significantly fewer calls to the computationally intensive objective function than non-specialized minimization algorithms.
An implementation for the method of elimination of differentials has been designed as an extension of FUMILI. A corresponding class inherits from the ROOT class TFumili and adds the methods for setting constraints and their derivatives.

Example: Monte-Carlo simulation
The algorithm implementation was tested on a simulated set of events. For this purpose 500 000 (a, b) events were Monte-Carlo generated according to the function with the parameters {x 1 = 0.5, 1}. Then the same equation (12) was fitted to them using an event-by-event log. likelihood method with constraints The results of the fit are presented in table 1. They show that constraints could lead to significant improvement in fit precision, that is particularly pronounced in the present case for the parameters x 2 and x 3 .

Kinematic fitting at ANKE
The developed FUMILI extension was used by the authors to perform the kinematic fitting of the data obtained at the ANKE spectrometer (Jülich, Germany) [4] for the pp → pp and pp → ppπ 0 reactions.
The early realization of the method, supporting only a single constraint equation, was applied during earlier ANKE studies where one constraint could be sufficient [11]. Then the penalty method [7] was used at ANKE for a long time. The recent developments of the realization described in this paper have allowed to employ it for measuring the pp → ppπ 0 Figure 1. Scheme of the pp → pp reaction being registered in the ANKE forward detector. Other ANKE detectors are not shown. forward differential cross section; the obtained results have been recently presented at [12] and now being prepared for publication.
Let us consider the results of the simulation that was performed for [12] for the pp → pp reaction registered in the ANKE forward detector ( fig. 1). In this case the detector acceptance allows to register only one of the final protons, so we use equation (3) for missing mass as a constraint. Tracking is being done in the forward detector by three multiwire gas chambers, each consisting of several multiwire planes. A set of signals produced by adjacent wires when a particle passes nearby is called a cluster. So, an event would contain a set of cluster coordinates {c 1 , . . . , c K } with their errors {δc 1 , . . . , δc K }, and the minimized functional is defined by equation (1). For the study under consideration, the reconstructed coordinatesĉ k are calculated on the basis of the known field map of the analyzing magnet using the Runge-Kutta method, assuming that these particles originated from a point-like source at the center of the target-beam interaction volume.
In order to check the improvement of the momentum precision, about 700 000 pp → pp events were simulated and then reconstructed either with or without kinematic constraint (3).
To check the precision of the reconstruction, we have calculated the difference between the original Monte-Carlo generated momentum components and the ones reconstructed by the detector simulation program. The reconstruction precision for momenta in polar coordinates in the center of mass system is shown in fig. 2. Since pp → pp is a binary reaction, condition (3) fixes the absolute value of the momentum |p cm p | to the delta function. Kinematic fitting improves the precision of the polar angle θ cm p nearly twice. For the azimuthal angle φ cm p the improvement not so significant and can be noticed only at the edges of the acceptance.

Conclusion and outlook
An approach has been proposed for performing minimization of χ 2 -like functionals with an arbitrary number of constraints called a method of elimination of differentials. A realization of the method has been implemented as an extension to the FUMILI minimizer and employed for studying various processes at the ANKE spectrometer. Simulation results are presented for the method showing significant improvement in fit precision. During the work on the code for the ANKE project, the idea arose of generalizing this code and publishing it in open access for use by the entire community. Our immediate plans include refactoring and modifying the existing code, setting up continuous integration and development for the project. In the long term, it seems interesting to create a universal framework that would use a single interface to perform kinematic fitting not only with the method  Figure 2. Errors of reconstructed proton momentum in polar coordinates (|p cm p |, θ cm p , φ cm p ) in the center of mass system for the pp → pp reaction at ANKE, simulated for T beam = 700 MeV either with (a) or without (b) kinematic fitting. The error bars show the average difference between the momentum components used as the simulation input and the ones obtained by the detector simulation program during the reconstruction. of elimination of differentials, but also using the penalty function or the Lagrange multiplier method.