Implementation of GENFIT2 as an experiment independent track-fitting framework
Tadeas Bilka, Nils Braun, Thomas Hauth, Thomas Kuhr, Lia Lavezzi, Felix Metzner, Stephan Paul, Elisabetta Prencipe, Markus Prim, Johannes Rauch, James Ritman, Tobias Schlueter, Stefano Spataro
IImplementation of
Genfit2 as anexperiment independent track-fittingframework
Tadeas Bilka [1] , Nils Braun [2] , Thomas Hauth [2] , Thomas Kuhr [3] , LiaLavezzi [4 , , Felix Metzner [2] , Stephan Paul [6] , Elisabetta Prencipe [7 , a ] ,Markus Prim [2] , Johannes Rauch [6] , James Ritman [7 , , Tobias Schlüter [9] ,Stefano Spataro [5] [a] = corresponding author[1] Charles University of Prague (CZ);[2] Karlsruher Institut für Technologie (DE); [3] Ludwig-Maximilians-Univesität München (DE); [4] Institute of High Energy Physics (CN); [5] Universitádi Torino and INFN (IT); [6] Technische Universität München (DE); [7] ForschungszentrumJülich (DE); [8] Ruhr-Universität Bochum (DE); [9] LP-Research Inc. (JP) Abstract
The
Genfit toolkit, initially developed at the Technische Universität München,has been extended and modified to be more general and user-friendly. The new
Genfit , called
Genfit2 , provides track representation, track-fitting algorithms andgraphic visualization of tracks and detectors, and it can be used for any experimentthat determines parameters of charged particle trajectories from spacial coordinatemeasurements. Based on general Kalman filter routines, it can perform extrapo-lations of track parameters and covariance matrices. It also provides interfaces toMillepede II for alignment purposes, and RAVE for the vertex finder. Results ofan implementation of
Genfit2 in basf2 and PandaRoot software frameworks arepresented here.
February 14, 2019 a r X i v : . [ phy s i c s . d a t a - a n ] F e b ONTENTS Contents
Genfit2 standalone toolkit . . . . . . . . . . . . . . . . . . . 188.2 Visualization of tracks in basf2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 188.3 Performance of
Genfit2 in basf2 . . . . . . . . . . . . . . . . . . . . . . . . . . . 248.4 Performance of
Genfit2 in PandaRoot . . . . . . . . . . . . . . . . . . . . . . . 27
INTRODUCTION Track reconstruction is an essential part of most nuclear and particle physics experiments.
Genfit2 [1] is an experiment-independent modular framework for track-fitting and otherrelated tasks. The project is open-source and available on GitHub .All particle physics experiments need to identify and classify processes based on detectorsignals. Combining these signals to recover particle trajectories is the task called “track recon-struction”: suitable collections of measurements must be combined into track candidates. Tracksmust be fitted, and points of common origin or exit, vertices, must be found and fitted. Trackfinding and track fitting are not independent from each other. Steps of different levels of trackfinding and refining alternate with track-fitting steps. Fitted tracks from different tracking sub-detectors can be matched and combined into larger tracks; and single measurements can beappended to existing tracks, which might come from another subdetector. With provided suit-able collections of measurements and detector geometry, all these tasks can be performed by Genfit2 . Genfit2 extends and improves the work previously performed by
Genfit [2], the externaltracking tool available in the PandaRoot releases since 2009. PandaRoot [3, 4] is the officialsoftware framework of the PANDA Collaboration [5], which is building a future ¯ pp experimentinvestigating reactions where an antiproton beam with momentum between 1.5 GeV/ c and 15GeV/ c impinges on a fixed target. PANDA, which will be located at the international FAIRfacility in Darmstadt (Germany), will consist of two main sections: the target spectrometer,which will have a solenoidal magnetic field ( B = 2 T); and a forward spectrometer, which willinclude a dipole field with 2 T · m maximum bending power.The new Genfit2 toolkit was already successfully used during a test of the Belle II high-level trigger architecture and the combined Belle II vertex detector readout architecture [6]. Itis used for recorded Phase 2 data and huge MC productions. Belle II [7] is a major upgrade ofthe Belle detector [8], which operates in a solenoid magnetic field of 1.5 T at an e + e − collider inTsukuba (Japan), running at a center of mass energy equal to that of the Υ(4 S ) . During theaforementioned test it served for online track reconstruction in the data reduction stage of the https://github.com/GenFit/GenFit. The Belle II experiment concluded the Phase 2 data taking period on July 17 th , 2018. During thePhase 3, planned on spring 2019, the experiment will run at the energy in the center of mass of the Υ(4 S ) and Υ(3 S ) . INTRODUCTION high-level trigger, as well as for offline analysis. It was also used to supply input to the MillepedeII alignment software [9].This paper documents the implemented classes, how to use them, and also presents perfor-mance results with the PandaRoot and the basf2 software frameworks. The latter is the officialAnalysis and Software Framework [10, 11] of the Belle II Collaboration. This document is orga-nized as follows: after this general introduction to Genfit2 , the motivation to build this toolkitis presented in the second Section. In Section 3 and 4, generalities about track parameteriza-tion, extrapolation and linearization are provided; the documentation of the
Genfit2 classes isgiven in Section 5, with the description of the general structure of this modular code. Section 6describes the different track fitters available in
Genfit2 ; Section 7 provided a short descriptionof the interface between
Genfit2 and RAVE [12]; Section 8 shows the performace in basf2 andPandaRoot; finally, the work is summarized in the Conclusion.
MOTIVATION FOR AN EXPERIMENT-INDEPENDENT TRACK-FITTING TOOLKIT In order to perform an analysis it is necessary to measure the position of a particle with itsuncertainty, and have a knowledge of the particle momentum.A measurement is a processed electronic signal recorded by a detector, which serves as inputto the reconstruction algorithms. Physics analyses generally work with particle tracks . Tracks can e.g. be parametrized by an Euclidean position (cid:126)x and momentum (cid:126)p with the correspondingcovariance matrix. To form a track from a set of measurements, two tasks have to be performed: track-finding and track-fitting . Track-finding groups measurements which are likely to originatefrom the same particle into track candidates. This is a task on its own, and will not be detailedhere.
Track-fitting takes the track candidates as input, and computes track parameters at anygiven point along the track, together with the covariance matrix of the track parameters, takingin consideration effects due to particle interactions with materials and the magnetic fields.Several high-energy-physics experiments implement their own track fitters; however, theyuse similar algorithms, such as the well known
Kalman filter [13]. The idea of
Genfit2 isto provide an open-source, modular and extensible framework which is capable of performingtrack-fitting and other related tasks, and can easily be adapted to various experimental setups.Smaller experiments, which do not have the manpower to develop their own track fitter, or newexperiments which need a working tool to do research and development, are encouraged to use
Genfit2 . In the absence of interaction, the trajectory of a particle with charge q in a magnetic field B isdetermined by the equations of motion, given its momentum vector p , in a single point alongthe track. In order to describe a track, its trajectory has to be parameterized. In a generalgeometry, and inside a magnetic field, it is most useful to use the distance along the track s asfree parameter, and then give the values: ( q/p, u, v, u (cid:48) , v (cid:48) ) (1)where u, v are the rectangular coordinates of a plane which the track intersects at distance s ,and u (cid:48) , v (cid:48) are the direction cosines, i.e. the projections of the direction of momentum on the TRACK PARAMETERIZATION AND EXTRAPOLATION coordinate axes. If the plane is that of a planar detector, it is most useful to have u (and v )coincide with the coordinate(s) measured by the detector. Genfit2 extrapolates tracks using the standard Runge-Kutta-Nyström method [14] as mod-ified by Bugge and Myrheim to carry along the Jacobian matrix [15, 16]. The code shares thesame pedigree as the STEP extrapolation code used in the ATLAS experiment [17, 18]. Thecorresponding class in
Genfit2 is RKTrackRep .During the Runge-Kutta extrapolation, global track parameters are used, namely q/p , thespatial coordinates x = ( x, y, z ) T and the unit tangent vector pointing along the track, i.e. thederivative of x with respect to the track length s .In order to reduce bias of the final fit result due to inaccurate seed values, which can leadto failed linearization of the Kalman filter, usually several fitting steps are done until the fitconverges. For example, we can assume here that the track is e.g. processed three times inforward and backward directions, giving an accurate estimation of the track parameters at thestarting position of the track. The last backward update of one iteration is usually taken as astarting value for the next iteration. In order to minimize the bias of the starting values, theircovariance has to be very large. However, if it is too large, numerical problems can arise. It ispractical to multiply the covariance of the starting value by a factor of 500 up to 1000 between2 iterations. This factor can of course be adjusted by the user, depending on the specific case.The fitter has a minimum and maximum number of iterations, which are assigned defaultvalues of 2 and 4, respectively, or the user can assign other values. As soon as the minimumnumber of iterations has been reached, the tool checks if the p -value has changed by less than acertain amount with respect to the previous iteration. The default value to check if the fit shouldperform a next iteration is set to 10 − . However, tracks with a p -value close to zero are oftenconsidered as converged with this criterion, while the χ , albeit big, is still changing significantly,indicating that the fit is still improving. This situation occurs often for tracks which are givenbad seed values. Thus, a non-convergence criterion has been introduced: if the relative changein χ from one iteration to the next is more than 20%, the fit will continue. The value of thisnon-convergence criterion can be defined by the user.After the extrapolation has reached the target surface, global coordinates are converted tolocal coordinates (see Eq. (1)), and the Jacobian for the extrapolation is calculated from theJacobian matrices of the extrapolation steps J i together with those of the coordinate changes LINEARIZATION local → global on the starting plane A , and global → local on the target plane B as: F A → B = J B global → local J n J n − · · · J J A local → global . (2)More details about the global → local coordinate system transformation are given in the nextSec. 4 and Appendix A. All implemented track fitting algorithms use a discrete model, where the particle travels fromdetector plane A to detector plane B such that the parameters on plane B are given as a functionof the parameters on plane A , i.e. p B = h A → B ( p A ) , and the Jacobian of this transport is thederivative of this function. Following the notation in Ref. [27], we use the symbol: F A → B ≡ ∂h A → B ( p ) ∂ p (cid:12)(cid:12)(cid:12) p = p A .The covariance of the track parameters on plane B can then expressed as: C B = F A → B C A F TA → B + N A → B , (3)where N A → B is the noise matrix evaluated by the track extrapolation algorithm. The evaluationof ( p B , C B ) is referred to as prediction.Each measurement m i gives rise to a residual r i , the difference between the measurementexpected based on the track parameters p i and the measured value, which in general is given as: r i = m i − f i ( p i ) , (4)where f is the measurement function which calculates the expected measurement given a set oftrack parameters. The measurement itself has a covariance matrix V i . Since our track param-eterization uses planar coordinates, this can be brought into a simple linear form when the u (and v for two-dimensional detectors) direction and zero point of the plane - in which the trackparameters are given - coincide with that of the measured coordinates, r i = m i − H i p i . (5)In this way, no further linearization is needed in order to apply the Kalman filter to this mea-surement model, and H i is a constant projection matrix for each detector. In general, H i has tobe replaced with the derivative of f i . STRUCTURE OF THE CODE The evaluation of F A → B and N A → B on the other hand cannot be eliminated by such a judi-cious choice of coordinates, and its value depends on the unknown track parameters. Especiallythe first few hits pose a problem: the track parameters only become determined successivelywhile the transport matrix already needs to be evaluated. Genfit2 has two approaches to this problem: • use the current result of the track fit, which may lead to imprecise or false results (it isthe SimpleKalmanFitter described in Sec. 6); • use information from the track finders to estimate an initial set of transport and noisematrices, then after each iteration of the fit prepare a new set of matrices which is used inthe next iteration (it is the ReferenceKalmanFitter described in Sec. 6).
Genfit2 is a C++ object-oriented modular tool, structured as shown in Fig. 1. It providesgeneral track-fitting tools, which are based on three pillars:
Figure 1: Structure of the
Genfit2 package [1].
STRUCTURE OF THE CODE
Genfit2 class
RKTrackRep . • measurements; • track representations; • fitting algorithms. M easurements serve as objects containing measured coordinates from a detector. Theyprovide functions to construct a (virtual) detector plane and to provide measurement coordinatesand covariance in that plane. The abstract base class
AbsMeasurement defines the interface. Ameasurement can be one-dimensional ( e.g. a line segment measured by a semiconductor stripdetector), two-dimensional ( e.g. a position on a semiconductor pixel detector or a drift isochroneof a wire detector), or three-dimensional ( e.g. a point in space measured by a Time ProjectorChamber (TPC)).
Genfit2 comes with predefined measurement classes for various detectortypes, which are called
PlanarMeasurement , WireMeasurement and
SpacePointMeasurement .In
Genfit2 , state vectors p and measurements m are given in local coordinate systems (seeformulae in Appendix A), which are defined by detector planes and implemented in the DetPlane class. A detector plane has an origin (cid:126)o , and is spanned by two orthogonal unit vectors (cid:126)u and (cid:126)v .The normal vector (cid:126)n is then equal to (cid:126)u × (cid:126)v (see Fig. 2). Measurements provide a matrix thatcan project the state p and covariance C into the measurement coordinate system of m and itscovariance V . For the Kalman filter mathematics, which are detailed in Appendix A, the state p and covariance C have to be projected into the measurement coordinate system of m and V .These projections are described by matrices, as in Eqs. 6-8.If we consider a silicon strip detector with strips parallel to (cid:126)u , the coordinate v in the local STRUCTURE OF THE CODE plane system reads as: m = ( v ) . (6)The H matrix now simply has to project out the v component of p : H = (cid:0) (cid:1) . (7)Now, Hp is the projected state, and HCH T is the projected covariance.For a detector which measures u and v , the following applies: m = (cid:0) u v (cid:1) (8) H = (cid:18) (cid:19) . (9)For planar detectors, the detector plane is given by the detector geometry, whereas for wire and spacepoint measurements, so-called virtual detector planes are constructed by extrapolating thetrack to the point of closest approach (POCA ) to the spacepoint or wire. Additional detectorseffects can be used here, allowing one to compensate for detector deformations like plane bending,wire sagitta, and misalignment during this procedure.Drift-time corrections and cluster fits (track-dependent clustering) are possible. All per-track data is kept in the
T rack object. It holds a set of
TrackPoint objects, which cancontain
M easurements , F itterInf o objects, which hold all fitter-specific information, and thinscatterers (currently only used in the GBL fitter, see Section 6.3).A
T rack contains one or more track representations, representing the particle hypotheses thatshould be used for the fit. One of them must be selected as the cardinal track representation.This can be done by the user or by
Genfit2 , which selects the track representation that best fitsthe measurements ( e.g. it has the lowest χ ). The T rack also contains a
F itStatus object, whichstores general information (number of iterations, convergence, etc.) and fit properties ( χ , NDF, p -value, track length, etc.). The track candidate ( TrackCand ) serves as a helper class, basicallystoring indices of raw detector hits in
T rackCandHit objects, which can also be overloaded bythe user to store additional information.
W ireT rackCandHit objects can store the left-rightambiguity. POCA stands for Point of Closest Approach.
STRUCTURE OF THE CODE A possible use case for this possibility to incorporate within these objects is the patternrecognition: algorithms can work on raw objects that have a lighter weight content than a
M easurement , or implement other features needed by the pattern recognition, but not neces-sarily by
Genfit2 . The
Genfit2 class
MeasurementFactory can build a
T rack from a
TrackCand . This
T rack canthen be processed by the various fitting algorithms. After fitting,
T rack objects contain a lotof data: track representations,
TrackPoint objects with measurements,
F itterInf o objects, etc.Usually, not all of this information needs to be stored on disk. For example, for physics analysisonly parameters related to the vertex and information on the track quality are necessary. Theuser can decide which data to keep, e.g. only the fitted state of the first
TrackPoint , and onlyfor the cardinal track representation.
Wire detectors measure a drift isochrone (a cylinder, in the simplest case). Intersected withthe detector plane, this gives two (one dimensional) measurements, since the particle could havepassed on either side of the wire.
Genfit2 provides the possibility to store several measurementsof the same type in one
TrackPoint . These tracks can also be fitted with the Kalman filter,therefore Genfit2 provides several options how to handle multiple measurements: • average : the average of the individual measurements is calculated. This option is primarilyused for the DAF (see Sections 6.1 and 6.2). • closest to prediction : the measurement which is closest to the state prediction is selected. • closest to ref erence : the measurement which is closest to the reference track is selected.This can only be used with the Kalman filter with reference track. • closest to prediction/ref erence f or wire measurements : if the TrackPoint has one
W ireM easurement , select the side which is closest to the prediction/reference. Otherwise,it makes use of the average . This feature allows the DAF to assign weights to the
TrackPoint . More information related to theDAF are in Sec. 6.2.
STRUCTURE OF THE CODE For wire measurements, it turned out that selecting the side which is closest to the statepredictions is often the best option. The average weigthed measurement m and its covariance V are calculated from the individual measurements i as follows: V = (cid:32)(cid:88) i (cid:0) w i V − i (cid:1)(cid:33) − m = V · (cid:88) i (cid:0) w i V − i m i (cid:1) (10) The Runge-Kutta track representation (
RKTrackRep ) is based on a Runge-Kutta extrapolatoradapted from Geant3 [19]. An abstract interface class interacts with the detector geometry.Implementations using the Root class
TGeoManager [20] and Geant4 class
G4Navigator [21] areavailable.The state p is parametrized with 5 coordinates in a local plane coordinate system, as alreadyshown in Fig. 2. The cartesian position (cid:126)x and direction (cid:126)a translate into plane coordinatesaccording to the following equations: p = (cid:0) q/p, u (cid:48) , v (cid:48) , u, v (cid:1) T u (cid:48) = (cid:126)a · (cid:126)u(cid:126)a · (cid:126)nv (cid:48) = (cid:126)a · (cid:126)v(cid:126)a · (cid:126)nu = ( (cid:126)x − (cid:126)o ) · (cid:126)uv = ( (cid:126)x − (cid:126)o ) · (cid:126)v. (11)where q is the charge of the track, and p its momentum. During fitting, material properties areused to calculate the following effects: energy loss and energy-loss straggling for charged particlesaccording to the Bethe Bloch formula (the code is based on the code ported from Geant3);multiple scattering (according to Ref. [22] or using the Highland-Lynch-Dahl formula [23]), wherethe full noise matrix is calculated; and soft Bremsstrahlung energy loss and energy-loss stragglingfor e − and e + (the code is based on the code ported from Geant3). The step sizes used for theRunge-Kutta extrapolation should be as large as possible to avoid unnecessary computation,while still being small enough to keep errors reasonably small. An adaptive step-size calculationis done in the RKTrackRep , taking field inhomogeneities and curvature into account. To calculatematerial effects correctly, extrapolation stops at material boundaries, and the step size is limited
FITTER IMPLEMENTATION so that a maximal relative momentum loss in the material is not exceeded. RKTrackRep providesdifferent methods to find the POCA of the track to non-planar measurements. These are usedto construct virtual detector planes: • extrapolate to point finds the POCA of the track to a given spacepoint. The virtualdetector plane contains the spacepoint and the POCA, and it is perpendicular to the trackdirection. A weight matrix can be used as a metric, defining the space in which the POCAwill be calculated. By default, the inverted 3D covariance of a spacepoint measurementis used as a metric, which gives correct fitting results also for spacepoints with arbitrarycovariance shapes. Again, the virtual detector plane contains the spacepoint and thePOCA, and it is perpendicular to the track direction in the space defined by the metric; • extrapolate to line finds the POCA of the track to a given line or wire. The virtual detectorplane contains the line and the POCA. This routine is used for fitting wire measurements.The intersection of the virtual detector plane and the measurement covariance gives thecovariance in the plane. TrackRep combines track parameterization and track extrapolation code. The state can beextrapolated through material and the magnetic field, to: • detector planes; • the POCA to points, lines, cylinders, and cones; • or by a certain distance.The state and covariance can be converted to/from a state defined in the global cartesian coor-dinate system, ( (cid:126) x (cid:126) p ) , as detailed in Appendix A. The abstract base class AbsTrackRep definesthe interface, and
Genfit2 ships with the
RKTrackRep as a concrete implementation.
Track representations combine track-parameterization and track-extrapolation code. Fittingalgorithms use the measurements and track representations to calculate fit results, which arestored in corresponding objects in the
T rack objects.Four different track-fitting algorithms are currently implemented in
Genfit2 : FITTER IMPLEMENTATION • two smoothing Kalman filters, one which linearizes the transport around the state predic-tions, and one which linearizes around a reference track; • a deterministic annealing filter ( DAF ); • a general broken lines (GBL) fitter. Genfit2 provides the possibility to store several measurements of the same type in one
TrackPoint , mainly for using the DAF to assign weights to them. Wire measurements alsoproduce two
MeasurementOnPlane objects, representing the passage of the particle on either sideof the wire. These tracks can also be fitted with the Kalman filters.
The Kalman filter is a progressive algorithm that successively adds the information from eachdetector hit to the estimated track parameters using a linear model of measurement and trackpropagation. It uses a series of measurements observed over time, containing “noise” (randomvariables) and other inaccuracies, and produced estimates of unknown variables; then iterates aleast square method in 3 steps: prediction, update, and smoothing.Smoothing [26] is a standard technique: the track is fitted in forward and backward direc-tion, where predictions and updates are saved in the
F itterInf os at each
TrackPoint . Withboth of these, smoothed track states can be calculated. The weighted average between forwardand backward prediction gives the unbiased state, whereas the average between prediction inone direction and update of the opposite direction results in the biased smoothed state. Thissmoothed track gives more accurate states than either forward or backward updates alone.The Kalman filter method introduced in
Genfit2 is widely documented. We refer specificallyto Ref. [27], and here we only document the specifics of our implementation which deviatesin two ways. First, all matrix operations are performed in a way which explicitly maintainsthe symmetry of the covariance matrices at each step in the calculation, thus minimizing thenumber of operations and avoiding inconsistencies due to round-off errors. Second, an alternativeimplementation is available which implements a square-root filter which additionally guaranteesthat the covariance matrices remain positive-definite while improving the robustness againstround-off errors or vastly different matrix entries which can especially happen at the begginingof the track fit, where individual track parameters only become known successively.
FITTER IMPLEMENTATION The Kalman filter performs a succession of fits (called “updates”), moving from the innermosthit along the series of hits identified by the track finders. Then the result of this first iterationof fits in the outwards direction is used as the starting value for the next iteration in which thefit proceeds in the backwards direction. In order to avoid biasing the inwards fit, the covariancematrix of the result of the outwards fit is multiplied by a large factor and the off-diagonal entriesare discarded. At each
TrackPoint the results from inward and outward fit are combined byaveraging the updated state of the outward fit and the prediction of the inward fit. This averagingis called smoothing, but it should be kept in mind that also for the smoothed track, the trackparameters change discretely at each measurement. Averaging is thus an important operationin
Genfit2 , and some care has gone into a numerically stable, efficient implementation. Thisimplementation is detailed in Appendix B.Given the track parameters p i − | i − , C i − | i − at TrackPoint i − after the preceding updatestep, we evaluate the prediction at plane i : p i | i − = h ( i − → i ( p i − | i − ) , (12) C i | i − = F ( i − → i C i − | i − F T ( i − → i + N ( i − → i , (13)and the residual of the prediction: r i | i − = m i − H i p i | i − . (14)Taking care to maintain the symmetry properties of the matrices, we can then write the updatedstate and covariance as: p i | i = p i | i − + C i | i − H ti ( V i + H i C i | i − H Ti ) − r i | i − , (15) C i | i = C i | i − − C i | i − H Ti ( V i + H i C i | i − H Ti ) − H i C i | i − . (16)This summarizes the formalism of the Kalman filter. The only difference between the so-called SimpleKalmanFitter and the
ReferenceKalmanFitter lies in the evaluation of the transport andnoise matrices, as explained above. Nevertheless, we have implemented an alternative equivalentformalism, using square-roots of covariance matrices to further improve numerical accuracy andresilience. This implementation is described in Appendix C. Employing it allows significantincreases in the blow up factor applied to the covariance matrix on fitter returns.During the first iteration, the reference track is constructed by extrapolating the seed statefrom measurement to measurement. In the following iterations, the smoothed states of the
FITTER IMPLEMENTATION previous iterations are taken as a new reference states. A χ value is calculated according to: r k = p k | n − p k,r χ k = (cid:88) i ( r k | n,i ) V k | n,i,i . (17)If the χ is below a certain value (default assigned value is 1), the corresponding reference statewill not be updated to save computing time.The same convergence criteria as for the standard Kalman filter are used. Additionally, thefit converges if the reference states do not change anymore, according to Eq. [17].The choice of the linearization point can significantly affect the performance of the track fit.In particular, if the first few hits lead to a large misestimate of the curvature, state predictionsmay stay very far from the actual trajectories. Consequently, linearizing around them is notoptimal; material and magnetic field lookups are also not done at the proper places. It istherefore common to linearize around reference states [24], which are calculated by extrapolatingthe start parameters to all TrackPoint objects. At later iterations, the smoothed states fromthe previous iterations are used as linearization points. However, if the change would be verysmall, reference states are not updated, saving computing time. It is also possible to let thefitter sort the measurements along the reference track, which can improve the fitting accuracy.In addition to the convergence criteria detailed above, the fit is regarded as converged if none ofthe reference states has been updated since the previous iteration.
The deterministic annealing fitter (DAF) [27] is a powerful tool for the rejection of outlyingmeasurements. It is a Kalman filter that uses a weighting procedure between iterations basedon the measurement residuals to determine the proper weights. The user can select which of
Genfit2 two Kalman filter implementations should be used, and specify the annealing scheme.The DAF is also perfectly suited to resolving the left-right ambiguities of wire measurements.However, a problem can occur: the weights of the
MeasurementOnPlane objects must be initial-ized. The basic solution is to initialize both left and right measurements with a weight of 0.5.Effectively, the wire positions are taken as measurements in the first iteration, and their covari-ance is twice the mean of the individual covariances. So, all the wire positions have the samecovariance, no matter how far from the actual trajectory they are. This systematically false
VERTEX RECONSTRUCTION estimate of the covariances biases the fit. Genfit2 implements a novel technique to initializethe weights that improves the fitting efficiency: measurements with larger drift radii are assignedsmaller weights, leading to larger covariances since the wire position is expected to be fartheraway from the trajectory. The mathematical espression used in the code is: w = 12 (cid:18) − r drift r drift, max (cid:19) , (18)where w is the given weight, and r drift is the distance of the measurement from trajectory.In contrast, measurements with smaller drift radii, which are closer to the trajectory, getlarger weights. When the annealing scheme has been passed, convergence is checked: if the ab-solute change of all weights is less than − (user configurable), the fit is regarded as converged.Otherwise another iteration is done, until the fit is converged or a maximum number of iterationsis reached. The generalized broken lines (GBL) method of track fitting [28, 29] was implemented especiallyfor the purpose of alignment with Millipede II [30]. It is mathematically equivalent to the Kalmanfilter (with thin scatterers instead of continuous materials), and works with a large number ofparameters, but fits tracks in their entirety in one step. For alignment purposes,
Genfit2 alsoprovides a set of interfaces for alignment parameters and derivatives which can be implementedby the detector classes.
GFRave , an interface to the vertex-fitting framework RAVE [12], has been implemented in
Gen-fit2 . RAVE is a detector-independent toolkit for vertex reconstruction, originally developed forthe CMS experiment [31].
GFRave takes full advantage of the
Genfit2 material model and thesophisticated algorithms of RAVE, allowing for precise and fast vertex reconstruction. GFRaveis part of the basf2 framework, where it is actively used. However, there is still no correspondingimplementation in PandaRoot, which is using for the time being other vertex and kinematicfitters.
PERFORMANCE Genfit2 standalone toolkit
Examples of what can be the result of a standalone
Genfit2 simulation are given in Figs. 3-4, which display the screenshots of an event reconstructed with
Genfit2 , performed by analgorithm based on
Teve : the reader can see that many features are available in the left-menu.Figures 5, 6, 7 show tests performed with a demonstration macro provided together with the
Genfit2 toolkit ( makeGeom.C ). By default, this test macro creates a geometry which comprisesa 1 m cube filled with air, and 4 layers of cylindrical silicon material with a radius of 1.0 cm,a length of 10 cm and a wall-thickness of 0.05 cm. The creation of measurements in thetest macros is independent on the geometry, but the geometry is basically provided to evaluatematerial effects. One million events with a single pion in a solenoidal field of 2 T have beensimulated. Different pion momenta have been simulated: 50 MeV/c, 100 MeV/c, 400 MeV/c, 1GeV/c. The information illustrated in the standalone Genfit2 event visualization is related tothe detector, hits, track with errors and track markers, the reference track, the result of the fitin forward and backward direction.
Genfit2 features a sophisticated 3D event display, which allows to visualize fitted tracks. De-tector geometry, measurements, detector planes, reference tracks, forward and backward fits(predictions and updates), smoothed tracks, and covariances of measurements and tracks canbe drawn. Tracks can be refitted with different algorithms and settings, and fit results can beviewed instantly. Examples have been shown in Figs. 3-4.
PERFORMANCE
Genfit2 tool. Here apion of 400 MeV/c is generated in a solenoidal magnetic field B = 2 T. PERFORMANCE (a) Measurements with covariance(yellow), planar detectors and driftisochrones (cyan), respectively, andreference track (blue). (b) Detecor planes (grey). For thespacepoint- and wire-hits, virtual detec-tor planes have been constructed.(c) Forward (cyan) and backward (ma-genta) fit with covariances of the stateupdates. (d) Smoothed track with covariance(blue). Figure 4:
Genfit2 event display screenshots [32]. The fit of a set of hits with the Kalmanfilter with reference track is shown. For demonstration purposes, the different hit typessupported by
Genfit2 are used, starting from the origin: planar pixel hit, spacepointhit, prolate spacepoint hit, two perpendicular planar strip hits, double sided planar striphit, wire hit, wire hit with second coordinate measurement.
PERFORMANCE Entries 999996 / ndf χ ± ± − Sigma 0.001 ± − − − Entries 999996 / ndf χ ± ± − Sigma 0.001 ± q/p pull Entries 999996 / ndf χ ± ± − χ ± ± − p value Entries 999996 / ndf χ ± ± ± − − − Entries 999996 / ndf χ ± ± ± u’ pull Entries 999996 / ndf χ ± ± − Sigma 0.001 ± − − − Entries 999996 / ndf χ ± ± − Sigma 0.001 ± v’ pull Entries 999996 / ndf χ ± ± − Sigma 0.001 ± − − − Entries 999996 / ndf χ ± ± − Sigma 0.001 ± u pull Entries 999996 / ndf χ ± ± ± − − − Entries 999996 / ndf χ ± ± ± v pull Figure 5: Pull of the track parameters ( u, v, u ’ , v ’ and q/p ), with the pion momentum =1 GeV/ c . Gaussian fits to these pull distributions are consistent with the expectation ofa normal distribution. PERFORMANCE Entries 999835 / ndf χ ± ± − Sigma 0.001 ± − − − Entries 999835 / ndf χ ± ± − Sigma 0.001 ± q/p pull Entries 999835p0 2.028e+01 ± ± − Entries 999835p0 2.028e+01 ± ± − p value Entries 999835 / ndf χ ± ± ± − − − Entries 999835 / ndf χ ± ± ± u’ pull Entries 999835 / ndf χ ± ± ± − − − Entries 999835 / ndf χ ± ± ± v’ pull Entries 999835 / ndf χ ± ± − Sigma 0.001 ± − − − Entries 999835 / ndf χ ± ± − Sigma 0.001 ± u pull Entries 999835 / ndf χ
157 / 159Constant 2.858e+01 ± ± ± − − − Entries 999835 / ndf χ
157 / 159Constant 2.858e+01 ± ± ± v pull Figure 6: Pull of the track parameters ( u, v, u ’ , v ’ and q/p ), with the pion momentum =100 MeV/ c . Gaussian fits to these pull distributions are consistent with the expectationof a normal distribution. PERFORMANCE Entries 998018 / ndf χ ± ± − Sigma 0.0007 ± − − − χ ± ± − Sigma 0.0007 ± q/p pull Entries 998018 / ndf χ ± ± − Entries 998018 / ndf χ ± ± − p value Entries 998018 / ndf χ ± ± ± − − − Entries 998018 / ndf χ ± ± ± u’ pull Entries 998018 / ndf χ ± ± ± − − − Entries 998018 / ndf χ ± ± ± v’ pull Entries 998018 / ndf χ ± ± − Sigma 0.001 ± − − − Entries 998018 / ndf χ ± ± − Sigma 0.001 ± u pull Entries 998018 / ndf χ ± ± ± − − − Entries 998018 / ndf χ ± ± ± v pull Figure 7: Pull of the track parameters ( u, v, u ’ , v ’ and q/p ), with the pion momentum =50 MeV/ c . Gaussian fits to these pull distributions are consistent with the expectation ofa normal distribution. PERFORMANCE Genfit2 in basf2
Figure 8: Scheme of the track reconstruction design in basf2 [33].
Genfit2 performance depends on many parameters, such as computer hardware and com-piler, number and types of measurements, momentum, geometry, magnetic field, number ofiterations, convergence settings, and so on. Nevertheless, a study has been conducted in both,basf2 and PandaRoot, to give the reader an idea of the performance expected in high energyphysics using
Genfit2 . Figure 8 shows the track reconstruction design in Belle II. The BelleII detector is structured with an internal vertex detector (VXD), which is made by 2 inner lay-ers of silicon Pixel Detector (PXD) and 4 double-side strips of Silicon Detector (SVD); then aCentral Drift Chamber (CDC) is given for tracking purposes and particle identification. A time-of-propagation (TOP) detector is also furnished, surrounded by an electro-magnetic calorimeter(EM), while the outer detector is for µ and K L detections. A Combinatorial Kalman Filter(CKF) was implemented at software level for track finding issues. Genfit2 takes care of track
PERFORMANCE fitting in basf2: it acts in a way to bind together track segments and information of PXD (planarhits), SVD (position along strips), and CDC (wire and drift time). The execution time of thebasf2 GenFitter module has been measured on a current 3.4 GHz office PC in single threadedoperation, an Intel machine with 4 cores. All code has been compiled with -O3 optimizationsettings.The module performs the fitting and a few more other features, like producing the
Gen-fit2 track from a
TrackCand , storing the fitted track in an output array, creating relations etc.,resulting in an overhead of a little less than 1 ms. Tracks are generated with a particle gun,with θ = 100 o and a momentum of 0.9 GeV/ c in a constant magnetic field equal to 1.5 T. Theresulting tracks have ± TrackPoints . Genfit2 is configured to do from 3 to 10 iterationswith default convergence criteria. From the results shown in Tab. 1, one can see that the
KalmanFitter with reference track needs less iterations to converge. Material lookup takes around 2.2ms per iteration for the Kalman filter. The
Reference Kalman fitter is also faster here, sincethe reference states are not recalculated if they are close to the smoothed states of the previ-ous iteration. This is also the reason why the
Reference Kalman fitter can often finish after2 iterations. The reference track is already so close to the smoothed track, that it would notchange in subsequent iterations. Table 1 shows that the calculation of the material effects takesa large proportion of the computing time (of course, it strongly depends on the complexity ofthe detector geometry under exam), and the Kalman fitter looks reasonably fast.
Table 1: Execution time of the GenFitter module in basf2.Fitter no material effects with material effects n. iterationsKalman 3.4 ms 10 ms 3Reference Kalman 4.0 ms 8.2 ms 2DAF 7.0 ms 11 ms 4
Examples of the good performance of
Genfit2 in basf2 are shown in Figs. 9-10: in theformer a study with cosmics is performed, showing the good shape of the track reconstructionsoftware code; in the latter a study showing the efficiency of tracks, reconstructed with basf2, isshown. The pull distributions of those parameters, defined as the ratio between the differenceof reconstructed and generated values, and the error of such parameter, is parameterized by adouble Gaussian function, and displayed in Fig. 9.
PERFORMANCE (a) (b)(c) (d) Figure 9: Study of track reconstruction in Belle II, using cosmics collected duringFebruary-March 2018: z (a) and d (b) resolution vs momentum. z and d are de-fined as the projections of the z axis of the POCA, and that of the x-y plane of thePOCA, respectively. Fits to real data (cosmics, red plot) and to the MC simulated singletrack events (blue plot) are displayed. Figures (c) and (d) show respectively the fits tothe pull distributions of those parameters [34]-[35]. PERFORMANCE
Genfit2 in basf2 [34], withand without machine background simulation.
Genfit2 in PandaRoot
An interface has been built between
Genfit2 and PandaRoot:
GenfitTools . The basic structureof
GenfitTools is the same used to interface the previous
Genfit package, but with a fewmodifications due to the new classes introduced in
Genfit2 (see Fig. 11). The new
GenfitTools structure has been created to allow both,
Genfit and
Genfit2 , to be used as fitting tool, andthen compare their results.
GenfitTools includes several objects: adapters (2); recohits (2); recotasks (2); trackrep .The objects adapters , recohits , recotasks refer to the Genfit interface to the PandaRootclasses, and the objects adapters recohits recotasks Genfit2 interface to thePandaRoot classes.
PndTrack is the standard object for PandaRoot tracks, and it is convertedinto a
Genfit2 or Genfit
Track by means of the adapters (2), which are able to deal withchanges in the
Genfit2 code without changing additional packages.In order to test the performance of the
Genfit2 and the new
GenfitTools package, a so-called particle-gun generator has been used. Five particle hypoteses have been tested: electrons,muons, pions, kaons and protons. A MC particle-gun generator is used to produce single par-
PERFORMANCE
GenfitTools package, in the PandaRoot trunk 28747 [36]. ticle beam for each different particle hypotheses, with different p T values, running at the sameantiproton beam momentum . For the example mentioned in this section a beam momentumof 15 GeV/c was used. We tested particles emitted at a fixed polar angle ( θ = 60 ). We alsochecked the Genfit2 performance in PandaRoot by changing the beam momentum, to evaluatethe resolution and pull of the tracking parameters. The magnetic field map used in these tests isa realistic one, as shown in Fig. 12(a). Since the PANDA detector is composed by two spectrom-eters, one with a solenoidal field (on the left part of Fig. 12(a)) and another with dipole field (onthe right of Fig. 12(a)), the intermediate region is characterised by strong inhomogeneities whichhas to be dealt by the extrapolation and fitting code. For comparison, the Belle II magnetic fieldscheme is shown in Fig. 12(b). More details related to this study are published in Ref. [36].The track representation used in the following tests is not
RKTrackRep , because PandaRootcurrently uses GEANE [38]; therefore an estimate on how long it takes to reconstruct an eventwith
Genfit2 in PandaRoot framework is biased, and not comparable with the numbers mea-sured for basf2.Figure 13 shows the performance of tracking tools in 2 different PandaRoot trunk revisions:the old tool (red squares) in comparison with the new
Genfit2 tool (blue squares). The an- A change in the antiproton momentum beam changes the flux density of the dipole magnetic field inthe PANDA forward spectrometer.
CONCLUSION tiproton beam momentum is set to 15 GeV/c. No significant differences are observed in the2 different PandaRoot trunk revisions under investigation, for high p T ; however, for p T < c the Genfit2 performance is now better. In the PandaRoot trunk revision 28747 a switchis provided to use both,
Genfit2 and
Genfit . Genfit2 has been introduced in the officialPandaRoot releases since 2017.The result of our simulation shows a clear improvement for reconstruction of low momentumtracks (see Fig. 13). A minimum transverse momentum threshold of 50 MeV/c is imposed, sincethat is required for at least one hit in the Straw-Tube-Tracker (STT).
Genfit2 is built in a wayto track low momentum tracks down 0 MeV/c.We summarize in Fig. 14 a preliminary achievement of this work in PandaRoot: it shows thepull distribution of the spatial tracking parameters d and z , running the particle-gun beam at p T = 1 GeV/ c . As expected, the width of a Gaussian fit to the distribution is consistent with1.0, and the mean with zero. Small variations are observed if we decrease the p T value, but theyare not significant, and the results confirm the stability of the fitter [40]. In summary, we have demonstrated that
Genfit2 can successfully provide track-fitting toolswhen the magnetic field is nearly homeogeneous (such as in Belle II) and when it is non-homogeneous (such as in PANDA).
Genfit2 has shown good track-fitting performance andhas been validated with MC samples. It has successfully been integrated into PandaRoot andbasf2. It easily adapts to different detector setups and geometries. With the new fitter imple-mentation classes, compared to the former
Genfit2 , the performance of low momentum trackreconstrucion is improved. Tests have been performed changing particle hypotheses: electron,muon, pion, kaon, and proton. The Runge-Kutta track representation works well in basf2, and isavailable also in PandaRoot now. The newly available interfaces to Millepede II and RAVE arealso useful and easy to utilize, as it is the new event display feature included in the toolkit. Pre-liminary results are promising, as reported in Ref. [35] (cosmics, and first data collected duringthe early time of Phase 2 of Belle II for a total integrated luminosity equal to 5 pb − ). Gen-fit2 is now used in several projects, including Belle II, COMET, PANDA, FOPI, SHIP, FOOTand BGO-OD. We encourage scientists to make use of our generic tracking tool and participatein our ongoing work to further improve and extend
Genfit2 . CONCLUSION (a)(b) Figure 12: (a) Magnetic field scheme in PANDA [36]: in the central part of the detectorit is constant and equal to 2T; in the dipole area the maximum bending power is 2 T · m.(b) Magnetic field scheme in Belle II [37]: it is a solenoidal field of 1.5 T at the center ofthe detector. CONCLUSION p T ) of charged pions, emitted atthe polar angle value θ = 60 , in PandaRoot. A comparison between the old tool (redsquares) and the new Genfit2 tool (blue points) is provided. The track fitting efficiencyis evaluated in a 3 σ window around the reconstructed momentum [36]. (a) (b) Figure 14: Generated values are subtracted from the reconstructed values of the variables d and z , and then divided by their errors. These plots have been obtained for p T =400 MeV/ c and they show the pull distributions [36] using Genfit2 in PandaRoot. Thewidth of the Gaussian fit to both distributions is consistent with being 1.0, as expectedfor a correct pull distribution, and the mean is consistent with being zero.
10 BibliographyReferences [1] J. Rauch and T. Schlüter, arXiV:1410.3698 [physics.ins-det].[2] C. Höppner et al. , NIM A 620, 518-525 (2010).[3] M. Al-Turany et al. , J. Phys.: Conf. Ser. 119, 032011 (2008).[4] S. Spataro, J. Phys.: Conf. Ser. 396, 022048 (2012).[5] W. Erni et al. , arXiv:0903.3905 [hep-ex].[6] T. Bilka et al. et al. , Progr. Theor. Exp. Phys. 04D001 (2012).[9] V. Blobel, NIM A 566, 5 (2006).[10] T. Abe et al. , KEK Report 2010, arXiV:1011.0352 [physics.ins-det].[11] T. Kuhr et al. , arXiv:1809.04299 [physics.comp-ph][12] W. Waltenberger. IEEE Trans. Nucl. Sci., 58, 434-444 (2011).[13] R. E. Kalman, Transcription of ASME - Jou. Basic Eng. 82, D35 (1960).[14] M. Abramowitz et al. , Handbook of mathematical functions with formulas, graphs, and math-ematical tables , Dover, New York (1964).[15] J. Myrheim et al. , NIM 160, 43-48 (1979).[16] L. Buggie et al. , NIM 179, 365-381 (1981).[17] E. Lund et al. , JINST 4, P04001 (2009).[18] E. Lund et al. , JINST 4, P04016 (2009).[19] GEANT - CERN Program Library Long Write-up, W5013 (1993).
EFERENCES [20] R. Brun and F. Rademakers, AIHENP96 Workshop (Lausanne) 389, 81-86 (1996).[21] GEANT4 Coll., GEANT 4 NIM A 506, 250-303 (203).[22] A. Fontana et al. , JPCS 117, 032018 (2008).[23] V.L.Highland NIM 129, 479 (1975); G.R. Lynch and O.I. Dahl NIM B58, 6 (1991).[24] R. Frühwirth et al. , Data Analysis T echniques f or High − Energy P hysics . CambridgeUniversity Press (2000).[25] B. D.O. Anderson and J. B. Moore,
Optimal F iltering . Dover Publications (1979).[26] R. Frühwirth, NIM A 262,444-450 (1987).[27] R. Frühwirth and A. Strandlie,
T rack F itting with Ambiguities and N oise : A Studyof Elastic T racking and N onlinear F ilters . Computer Physics Communications, 197-214(1999).[28] V. Blobel and C. Kleinwort, NIM A 673, 107-110 (2012).[29] V. Bloblel, NIM A 566 (2006) 14-17.[30] V. Bloblel and C. Kleinwort, Proc. Conf. ASTPP (2002).[31] CMS Collaboration, JINST, 0803:S08004 (2008).[32] J. Rauch and T. Schlüter, JPCS 608, 012042 (2015).[33] Scheme presented on behalf of the Belle II Coll. by S. Spataro at the CHEP conference,Sophia, 2018.[34]
The Belle II Physics Book , arXiV:1808.10567 [hep-ex], pag.60-61.[35] E. Prencipe, arXiV: 1808.04721 [hep-ex] (submitted to EPJ).[36] E. Prencipe, JPCS 898, 042037 (2017).[37]
The Belle II Physics Book , arXiV:1808.10567 [hep-ex], pag.56.[38] V. Innocente et al. , GEANE: Average Tracking and Error Propagation Package, CERN ,Program Library W5013-E (1991)
EFERENCES [39] D. N. Brown et al. , The BaBar Track Fitting Algorithm , Proc. CHEP Conf. (2000).[40] E. Prencipe, talk given at the
Common Track Reconstruction Software Forum , Vienna(2016).[41] G. H. Golub, Num. Math. , 206 (1965).[42] P. Businger et al. , Num. Math. , 269-276 (1965).[43] B.D.O. Anderson and J.B. Moore, Optimal filtering , Dover, 2005, pag. 5.
COORDINATE SYSTEM A Coordinate system
GENFIT2 uses three different sets of variables to characterize trajectories. Their definitions andthe conversions between them is discussed below.We start by introducing the conventions for labelling planes. A plane in three-space is definedby three vectors: we call its origin (cid:126)o (0,0,0), and define the two orthogonal unit vectors u , v ,giving the directions of the coordinate axes; the normal of the plane points in the direction of (cid:126)n = (cid:126)u × (cid:126)v . (cid:126)v , (cid:126)u and (cid:126)n identify a left-handed system, respectively, as shown in Fig. 2, Sec. 2.The three sets of variables are: • . This is the common set used for user interfacing. It contains the cartesiancoordinates for the position and the momentum vector, labeled as:( x, p ).These, together with an independently provided charge q , specify a helix uniquely. • local coordinates . These are the five coordinates used in track fitting. They are givenrelative to a plane with local orthogonal coordinates u and v . The coordinates are labeledas: ( q/p, u (cid:48) , v (cid:48) , u, v ) Here q is the charge of the particle; p is its momentum; u (cid:48) , v (cid:48) are the direction cosinesrelative to the plane; u , v are the local coordinates on the plane.These variables are supplemented by a variable σ = ± , indicating whether the verticalcomponent of the track enters the supporting plane parallel or antiparallel to its normalvector. Conversion to the 6D-coordinates is achieved as follows: q/pu (cid:48) v (cid:48) uv → (cid:32) O + uU + vV σ N + u (cid:48) U + v (cid:48) V √ u (cid:48) + v (cid:48) × p (cid:33) COORDINATE SYSTEM • global coordinate . These are the seven coordinates used in track extrapolation: the position x , the direction unit vector T and charge over momentum of the track, written as:( x , T , q/p)Conversion to the 6D-coordinates is performed by the mapping:( x , T , q/p) → ( x , p T ).These coordinates contain the complete helix information. Additionally, we can representa track by its helix parameters, which are labeled as:( d , φ, ω, z , tanλ )and which are defined following the conventions used at BaBar [39].Jacobians are used for the various possible coordinate transformations. For the mapping betweentwo sets of vectors, x → y , we define the Jacobian matrix as: J ( x → y ) = ( ∂y∂x ) , (19)where the y labels the rows and the x labels the columns. This, together with matrix multipli-cation, gives the usual chain rule for the mapping x → y → z , J ( x → z ) = J ( y → z ) J ( x → y ) (20)and we obtain the usual error propagation formula: C y = J ( x → y ) · C x · J ( x → y ) T (21)which expresses the covariance matrix of the y in terms of the error matrix of the x togetherwith the Jacobian of the mapping between them. J ( local → global ) = U x V x U y V y U z V z − u (cid:48) σN x − u (cid:48) v (cid:48) V x +(1+ v (cid:48) ) U x (1+ u (cid:48) + v (cid:48) ) / − v (cid:48) σN x − u (cid:48) v (cid:48) U x +(1+ u (cid:48) ) V x (1+ u (cid:48) + v (cid:48) ) / − u (cid:48) σN y − u (cid:48) v (cid:48) V y +(1+ v (cid:48) ) U y (1+ u (cid:48) + v (cid:48) ) / − v (cid:48) σN y − u (cid:48) v (cid:48) U y +(1+ u (cid:48) ) V y (1+ u (cid:48) + v (cid:48) ) / − u (cid:48) σN z − u (cid:48) v (cid:48) V z +(1+ v (cid:48) ) U z (1+ u (cid:48) + v (cid:48) ) / − v (cid:48) σN z − u (cid:48) v (cid:48) U z +(1+ u (cid:48) ) V z (1+ u (cid:48) + v (cid:48) ) / (22) COORDINATE SYSTEM J (6 D → global ) = p y + p z p / − p x p y p / − p x p z p / − p x p y p / p x + p z p / − p y p z p / − p x p z p / − p y p z p / p x + p y p / − qp x p / − qp y p / − qp z p / (23) J ( global → D ) = p − p qT x p − p qT y p − p qT z (24)The important special cases, where the plane is orthogonal to the track, i.e. u (cid:48) = v (cid:48) = 0, delivers: J ( local → D ) = J ( global → D ) · J ( local → global ) = U x V x U y V y U z V z − p qσN x U x V x − p qσN y U y V y − p qσN z U z V z (25)In order to express efficiently the Jacobian for conversion from 6D to local coordinates, we expressthe momentum vector in its components ( p U , p V , p N ), where, e.g. p U ≡ p · U , and σ = sgn( p N ).Then: J (6 D → local ) = − q ( p N N x + p U U x + p V V x ) p − q ( p N N y + p U U y + p V V y p − q ( p N N z + p U U z + p V V z p σ ( p N U x − p U N x ) p N σ ( p N U y − p U N y p N p N U z − p U N z p N σ ( p N V x − p v N x p N σ ( p N V y − p V N y p N p N V z − p V N z p N U x U y U z V x V y V z (26) J ( global → local ) = ( T · N ) U x − ( T · N ) N x ( T · N ) σ ( T · N ) U y − ( T · N ) N y ( T · N ) σ ( T · N ) U z − ( T · N ) N z ( T · N ) σ
00 0 0 ( T · N ) V x − ( T · N ) N x ( T · N ) σ ( T · N ) V y − ( T · N ) N y ( T · N ) σ ( T · N ) V z − ( T · N ) N z ( T · N ) σ U x U y U z V x V y V z (27)Again, the special case where momentum is orthogonal to the plane in question ( T · N = p p · N AVERAGING = σ = ± J (6 D → local ) = − qp N x − qp N y − qp N z σ U x p σ U y p σ U z p σ V x p σ V y p σ V z p U x U y U z V x V y V z (28)and: J ( global → local ) = σU x σU y σU z
00 0 0 σV x σV y σV z U x U y U z V x V y V z (29)As consistency check, one can try to calculate the transformation: local → global → local :the Jacobian of this transformation is indeed equal to 0; the converse transformation, e.g.global → local → global , has singular Jacobian. Its null space is spanned by: ( , , , T x , T y , T z ),( N x , N y , N z , , , , ), which comes about because T is normalized (this brings about the firstdirection of the null space) and because a point in local coordinates is constrained to a plane(the second direction). B Averaging
Given two measurements of the same n -dimensional quantity and their respective covariances, ( x , C ) and ( x , C ) , the two measurements can be combined to a least-likelihood estimate ˜ x byminimizing the χ of the combination, i.e.χ = (˜ x − x ) T C − (˜ x − x ) + (˜ x − x ) T C − (˜ x − x ) ≡ min ˜ x . (30)The solution is the well-known formula: ˜ x = ( C − + C − ) − ( C − x + C − x ) , (31)and the covariance of ˜ x is by linear error propagation: ˜ C = ( C − + C − ) − . (32)We can simplify the evaluation of these textbook formulas significantly in the following manner,applying a technique from Ref. [41]: covariance matrices are positive definite matrices. As such, AVERAGING they can be written C = S T S where S is some matrix which can be chosen to be upper triangular.In the case where C is diagonal, the entries of C correspond to the squared errors σ and theentries of S correspond to the errors σ . It is in this sense that S is referred to as a squareroot of C . It is commonly evaluated by Cholesky decomposition. Since S is a upper triangularmatrix its inverse, and thus the inverse of C , can be evaluated directly by means of forwardsubstitution. Likewise products of the form S − x can be evaluated by forward substitutionwithout first evaluating the inverse. Our goal is to replace the matrices in Eq. (31) and Eq.(32)by their square roots and use those directly, thus operating with the errors instead of the squarederrors.Starting with the inner part of (32), we can write: C − + C − = S − T S − + S − T S − (33) = ( S − T , S − T ) (cid:18) S − S − (cid:19) , (34)where on the second line we use block matrix notation, i.e. the matrix on the right is n × n .As such, there is an n × n -dimensional orthogonal matrix Q , QQ T = 1 , which converts it toan upper triagonal matrix, written explicitly: (cid:18) S − S − (cid:19) = Q (cid:18) R (cid:19) , (35)where R is n × n . This operation is the QR decomposition of the left-hand side, and we canwrite: C − + C − = ( R T , QQ T (cid:18) R (cid:19) (36) = R T R. (37)Going back to Eq. (32), we have: ˜ C = R − R − T . (38)In other words, we can obtain the square root of ˜ C directly from the square roots of C , C bystacking them and performing a QR decomposition. This increases the numerical range in thepresence of vastly different errors, as we no longer operate on the squared errors but on the errorsthemselves. Additionally, any matrix of the form in Eq. (38) is always invertible and positivedefinite which prevents catastrophic conditions in subsequent numerical operations. PROOF OF SQUARE-ROOT KALMAN FILTER FORMULA It is perhaps revealing to consider the similarity to the Pythagorean theorem which gives thelength of a vector ( a, b ) T as a + b = c . The corresponding orthogonal transformations in thiscase are the rotations which take the vector ( a, b ) T to ( ±√ a + b , T . The computational gaincomes from the fact that the orthogonal operation disappears because the relevant quantities aresquares which always lead to factors QQ T or Q T Q (for the inverse matrices).In order to evaluate ˜ x , we insert those into Eq. (31), and obtain: ˜ x = R − R − T ( R T , Q (cid:18) S − x S − x (cid:19) (39) = R − ¯ Q (cid:18) S − x S − x (cid:19) , (40)where ¯ Q refers to the upper n rows of Q . It never explicitly appears in the implementation,instead the right-hand multiplication is evaluated simultaneously with the QR decomposition inEq. (35), using the algorithm of Ref. [42].In practice this algorithm is not only preferable because of its numerical properties, we alsomeasured it to perform slightly faster than a naive implementation of Eqs. (31) and (32). C Proof of Square-root Kalman filter formula