Analytical computation of process noise matrix in Kalman filter for fitting curved tracks in magnetic field within dense, thick scatterers

In the context of track fitting problems by a Kalman filter, the appropriate functional forms of the elements of the random process noise matrix are derived for tracking through thick layers of dense materials and magnetic field. This work complements the form of the process noise matrix obtained by Mankel[1].


Introduction
Kalman filter [2] is a versatile algorithm that has wide applications in various fields, like [3][4][5][6][7][8][9][10][11] etc. In 1987, Frühwirth [12] demonstrated its application to track fitting problems in high energy physics experiments for the first time. Since then, many experiments adopted this tool for track fitting purpose (for example, [13,14]) and various authors contributed to different aspects of the algorithm (for example, [15,16]). The problem is to estimate the charges, momenta, directions etc. of the observed particles from the measurements performed along their tracks.
These parameters are combined together to form a state vector. Usually, a Kalman filter based program (estimator) deduces the near-optimal values of the elements of the state vector iteratively, from the weighted averages of the predicted locations of the particle positions and the measured particle positions at the sensitive detector elements. In general, the prediction is done based on some analytical (or numerical) solution to the equation of motion of a charged particle passing through a dense material and magnetic field (see Ch. 3 of [13], or [16], for instance). However, the prediction represents the deterministic aspect of the particle motion. But the motion of the particle is also affected by the random processes like multiple Coulomb scattering [17] and energy loss fluctuations [18]. These are the stochastic perturbations to the deterministic motion of the particle, the latter being controlled by the magnetic field and the average energy loss. The estimator must take into account the random fluctuations appropriately, because precision of the filter estimation depends crucially on proper treatment of these random processes. Clearly, when the charged particle passes through thick layers of dense materials, the effects of such fluctuations are greater. This situation arises in case of track fitting in the Iron CALorimeter (ICAL) experiment, which is an upcoming neutrino oscillation experiment under the India-based Neutrino Observatory (INO) project [20]. It comprises a 50 kiloton magnetized iron calorimeter detector of dimension 48m × 16m × 14.4m, divided into three identical modules, as seen in Figure 1(a). The sensitive detector elements are made by 2m×2m resistive plate chamber detectors (RPCs), placed horizontally, which are sandwiched between 5.6 cm thick plates of iron. RPCs are planes of constant Z coordinates and any two RPCs are separated by vertical width of 9.6 cm. The iron plates are magnetized with current coils which generate up to 1.5T of magnetic field (Figure 1(b)). ICAL will try to resolve the neutrino mass hierarchy, by observing the earth matter effect on the neutrino oscillation. The experiment is most capable of the measurement of the properties of the muons, coming from the charged current interactions of the muon-neutrinos. These muons travel through the different layers of detector materials and leave electronic signals at the RPC planes. The position measurements done from these signals are used for track fitting. Since ICAL will observe atmospheric neutrinos of a wide energy range (E ν ∈ 1 − 15 GeV) coming from all directions, it is clear that a major fraction of muon tracks will be strongly affected by multiple scattering, while crossing the horizontal thick layers of iron at various angles. The thickness and the radiation lengths of the dense materials that the muons have to pass through within the ICAL detector are shown in the following The width of the scattering angle is related inversely to the particle momentum and the radiation length of the material it is passing through [21]. So, the muons will be subjected to significant amount of multiple scattering inside iron. The effect will clearly be more pronounced at lower energy. It is also important to note that the flux of atmospheric neutrinos is much higher at lower energy [22]. Thus, the ICAL track fitting program must account for the random effects in a proper fashion.
Let us consider the state vector x = (x, y, t x , t y , q/p) T which is used in many experiments like INO-ICAL [19,20,22,23], MINOS [24][25][26], LHCb [27][28][29][30] with forward detector geometry. Since the Kalman prediction is performed along an approximate particle trajectory, it introduces some deterministic uncertainties (dependent on magnetic field, average energy loss etc.) to the elements of the state vector. The random processes introduce additional uncertainties to these elements. These uncertainties are accounted for in a error covariance matrix C = (x−x)(x−x) T , wherex contains the true values of the elements of the state vector. The total error matrix propagated from a point l to the next l + dl along the track is given by: where F denotes the Kalman propagator matrix, encoding the deterministic factors between l and l + dl. F propagates the errors of the track parameters, represented by C matrix, deterministically, from l to l + dl. On the other hand, the matrix Q represents the error contributions from all the random processes to the total error C at l + dl. However, between two measurement sites, separated by some distance, the track fitting program should be sensitive to the possible variations of track parameters (momenta, direction etc.) and also to the possible variations of ambient parameters (materials, magnetic field components etc.). Then, one must apply Eq. (1) repeatedly, in small tracking steps, while approaching towards the next measurement site. Thus, the effective propagator matrix becomes F = Π N j=1 F j between two measurement sites [19]. Hence, the total propagated error at the next measurement site equals the sum of the (a) matrix representing deterministic error propagation (Π N j=1 F j ) C l 0 (Π N j=1 F j ) T and the (b) sum of the matrices of the deterministically propagated random uncertainties in all the tracking steps. It can be shown from Eq.(1) that this term becomes equal to (Eq. (3.16) of [13]): where F ms,k denotes the product of F j s between m s -th step and the final step. That is, to propagate the random uncertainties of a 'deeper' layer, a longer 'chain' of F j s is required.
The variances of the position, angle and the momentum elements of the state vector, arising from the multiple scattering and energy loss fluctuation in the thin layer of dense materials, have been investigated by various authors [17], [21,[31][32][33][34]. However, when passage of a particle through a thick layer of dense material is considered, one has to use effective variances and covariances, valid in the thick scatterer limit. These terms are obtained from a thorough study of Eq.(2) (see Appendix B of [1] written by Mankel). The author takes a simple form of the Kalman propagator matrix (F ) and obtains a set of 10 ordinary linear coupled differential equations. The solutions to these equations correspond to the elements of the random noise matrix in the thick scatterer limit.
However, the result of this work is not general in two respects: (1) the propagator matrix has been assumed to be constant and very simple in form (see section 2.3). This results in simple analytical form of elements of the random noise matrix Q (Eq. (12)). However, in many experiments, the Kalman propagator matrix may evolve significantly from iteration to iteration and may have a quite non-trivial form (for example, in ICAL track fitting program [19]). Naturally, in these cases, one needs to find the more appropriate form of the random noise matrix. (2) This work [1] concerns only the 4 × 4 block of the random noise matrix that corresponds to the position and the angular elements which directly suffer from multiple scattering. But the functional forms of the variance and covariance terms of q/p with other state vector elements which are affected by the fluctuations in energy loss, are not considered in this work.
The purpose of this paper is to derive the appropriate functional form of all the elements of the random process noise matrix for a curved track in magnetic field in the thick scatterer limit. We shall take a non-trivial and evolving propagator matrix for this purpose and ascertain what difference it makes to the track fitting performance. Even if the modification does not yield significant improvements in the track fitting performance, this exercise serves two purposes: (a) it completes the problem from a mathematical point of view and (b) it confirms that Mankel's approximate solutions are good enough. To the best of knowledge of the authors, no work has been done before which addresses these two issues.
The problem will be formulated mathematically in the next section 2. The desired elements of the random noise matrix will be seen to be solutions of a matrix differential equation. Then, we will describe two methods of obtaining its solutions in section 3. Among these methods, the first one (decoupling a set of linear coupled ODEs) is practical for implementation and will be used in the ICAL track fitting program in the presence of magnetic field. The details of implementation technique will be discussed in section 4. The relevant details on software supports will be covered in Appendix C. The reconstruction performance will be shown in section 5. We will conclude with a discussion of the merits and demerits of the approach in section 6.

Mathematical formalism
In case of the deterministic propagation of the random uncertainties, Kalman propagator matrix F transports these uncertainties at l to l + dl. The total random uncertainty matrix at l + dl has another term coming from the random uncertainties introduced to the direction and the momentum of the particle due to the multiple scattering and the energy loss fluctuations by the material between l and l + dl. We call this term δQ. The overall process noise matrix Q at l + dl is given by: where s = +1(−1) when the direction of propagation increases (decreases) the z coordinate while the tracking is carried out and I denotes the identity matrix. The elements of the F matrix (i.e. the dots within the parenthesis of the matrix in Eq.(4)) are the track length derivatives of the elements of the residual propagator matrix (F − I). These quantify the additional uncertainties introduced by the presence of the magnetic field etc. to the existing uncertainties at l [15, pp. [10][11][12]. Concrete examples of the elements can be found from [13], [15], [19] etc. We shall see that the nature of in Eq.(4) determines the functional forms of the elements of Q matrix.

Some comments on δQ
Since this uncertainty originates from a very small step of length dl, it may be assumed that the scattering took place in a plane of infinitesimal thickness. The elastic scattering with the Coulomb field of the nuclei of the dense detector material brings about a sudden change in the particle direction at the plane of the scattering. However, the particle position does not change laterally at that plane. Also, the magnitude of the momentum of the particle hardly changes as the energy imparted to these heavy nuclei is practically negligible [35, pp. 20]. If instead of q/p, q/p T is chosen to be a state element, where p T denotes the transverse momentum, it will change at that plane where the particle undergoes the scattering [1, pp. 9]. So, multiple scattering introduces uncertaintyonly in the particle direction and it is parametrized by two orthogonal angles θ 1 and θ 2 , defined with respect to the particle direction. On the other hand, the fluctuation in the energy loss happens due to uncertainty in the collision rate with the atomic electron when a high energy particle passes through a dense material. The physical mechanism of the ionization hardly changes the particle direction but surely changes the magnitude of the momentum. The fluctuation, therefore, is independent of multiple scattering angles, but dependent on particle momentum p. Now, the covariance between m th and n th elements of the state vector is given by: In Eq. (5), ξ i denotes any variable representing fluctuation due to the random processes (thus, ξ = θ 1 or θ 2 or p) and σ(ξ i ) is the width of that fluctuation. Since θ 1 , θ 2 and particle momentum p are independent parameters, one does not need to calculate the covariance terms between (ξ i , ξ j ) for i = j. Then, for the chosen state vector (x, y, t x , t y , q/p) T , the corresponding covariance elements may be calculated (for point scattering). All covariances with position coordinates (x or y) is zero according to our assumption that there is no horizontal shift of particle position in the infinitesimal plane of scattering. The covariances c(t x , q/p), c(t y , q/p) = 0, because: Now, change of direction due to multiple scattering does not change p and change of momentum due to energy loss fluctuation does not change direction t x or t y . As a result, ∂(q/p) Thus, over a tracking step length dl, the integrated random uncertainty matrix is given by: The nonzero variance and covariance elements of t x and t y are known in terms of the rms errors of the scattering angles [21,34]. The calculation of c(q/p, q/p) is available from [15].

Formulating the problem
In this section, we formulate the problem in the same way as indicated in Appendix B of [1]. However, we also take into account the effect of energy loss fluctuation on the q/p element of state vector. At a track length l, the random noise matrix Q(l) is given by: where in Eq. (8), the symmetric elements of the real symmetric matrix Q has been replaced by dots. This shows that there are exactly fifteen independent elements of the process noise matrix that need to be determined. If the propagator matrix F deviates from the identity matrix I by a matrix F s dl (see Eq. (4)), then we can say: From Eq.(9), one can deduce the differential equation of process noise Q: We note that dQ dl , δQ/dl and (F Q(l) + (F Q(l)) T in Eq.(10) are real symmetric matrices. This equation encodes a system of 15 coupled linear ODEs corresponding to the 15 independent elements of Q. The matrix (F Q(l) + (F Q(l)) T has been calculated with the help of Mathematica [36], assuming all the elements of F are nonzero. In fact, some elements of F were found to be rather high (of the order of one or more) depending upon the tracking directions and momenta. The functional form of every element of Q is the solution of the set of independent equations in Eq.(10).

Mankel's solution
In his work, Mankel [1] used a 4 × 4 block of random noise matrix whose elements were covariance terms of position and angular coordinates. The corresponding 4 × 4 block of the propagator matrix was given by: That is, except for F 13 = F 24 = 1, Mankel took all the other elements of the F matrix to be zero. In that case, the random noise matrix has 10 independent elements. Thus, 10 linear coupled ODEs are obtained. The simple form of the propagator (Eq. (11)) keeps the forms of the coupled equations simple. They can be easily solved and the resulting process noise matrix becomes: -where the symmetric counterparts are replaced by dots.

Solution of the problem
The matrix solution Eq. (12) is not valid in general, when all the elements of the propagator matrix are nonzero. From Eq.(10), if the matrix connecting the fifteen independent elements of Q (i.e. Q 11 to Q 55 ) to their derivatives is given as A 15×15 , then we can write: This matrix is real but not symmetric. From Appendix A, it is seen that 110 elements out of 225 elements of A 15×15 matrix are zero. Further simplifica-tions arise from the fact that only 4 elements of the 15 elements of δQ/dl vector are nonzero. Hence, Eq.(13) can be succinctly written as: where q is a column vector of the fifteen independent elements of the Q matrix (Q 11 , Q 12 , ...Q 55 ) and δq denotes the vector of the corresponding elements of δQ matrix (see Appendix A). Within the step of length dl, the elements of A remain unchanged, as they are obtained from the propagator matrix for that step. Hence, the problem is to solve non-homogeneous linear coupled system of differential equations with constant coefficients. Now, we shall investigate different approaches for solving this initial value problem and discuss their merits and demerits.

Solution by decoupling
The most elegant method to solve Eq. (14) is to decouple the equations by diagonalizing A. If A is diagonalizable (i.e. A = P DP −1 ) with an invertible P and a diagonal D, the system of equations can be decoupled through the substitution q = P u. In that case, Eq.(14) reduces to: Here P is the matrix of the eigenvectors of A; the corresponding eigenvalues are located at the diagonal position of the diagonal matrix D. As A is not necessarily real symmetric, the eigenvalues can be complex numbers as well and A may not be diagonalizable altogether in some cases. However, when it is diagonalizable, we can easily solve Eq.(15) for u from the fact that the j th component of the equation is just a first order linear ODE: where the set of {λ j } denotes the set of eigenvalues of A 15×15 . Eq.(16) can be solved by using the integrating factors and the solution to Eq.(14) becomes: We assume that P −1 δq varies very slowly over the small step of length l, so that it may be considered to remain constant while calculating the integral in Eq. (17). Thus, we get: In Eq. (18), there are 15 unknown coefficients u j (0) that must be deduced from the initial conditions. The initial condition is that at l = 0, all random noise errors are zero. We see that for l = 0, Eq.(18) reduces to: Eq. (19) is possible only if all u j (0)s are individually zero. Thus, we have: In the case when A is diagonalizable, the only difficulty of implementation is the occurrence of complex numbers in the result. In this case, we simply take the real parts of q i (l) to form the elements of the random noise matrix. The imaginary parts of q i (l) cannot be used, as the imaginary parts of q 1 (l), q 6 (l), q 10 (l), q 13 (l), q 15 (l) (that correspond to the diagonal elements of Q, i.e. the variance terms) are found to take negative values frequently. This inconsistency does not occur if real parts of q i (l) are used. As long as the matrix A is diagonalizable, P is invertible and λ j = 0, this method is observed to work. This typically happens inside the magnetized iron plates of ICAL detector. However, outside iron, the conditions are not satisfied (Det(A) → 0, one or more λ j are zero etc). As a result, Eq.(20) cannot be used there.

Reconciliation with the process noise matrix derived in [1]
The simple process noise matrix (Eq. (12)) derived in [1] is valid in a region with zero magnetic field where the simple form of the Kalman propagator matrix (Eq. (11)) is valid. On the other hand, Eq. (20) describes the form of every independent element of the process noise matrix in presence of magnetic field. It is not possible to directly reduce q i (l) of Eq. (20) to the corresponding elements of Eq. (12) in the absence of magnetic field to check whether the generalization has been consistent, since in that scenario Det(A) becomes zero (or very close to zero) which prohibits the computations of P matrix and the eigenvalues λ j . However, it is possible to reconcile Eq. (20) with Eq. (12) inside the magnetic field, by checking if the real parts of q i (l) are close to the corresponding elements in Eq. (12).
Although there is an exponential dependence in Eq. (20), the exponent can be replaced by its series in the limit of small step length l. As a result, each term in the summation becomes a power law in itself and can be represented as: So, one needs to check if the real parts of the coefficients of l, l 2 and l 3 in Eq.(21) are close to the corresponding coefficients in Eq. (12). This exercise has been performed and the results are shown in Appendix B.

Solution method without diagonalization
In this case, one first needs to solve the homogeneous equation dq dl = Aq (where the constant factor s is absorbed within the matrix A). The solutions for the vector q(l) are used to form a fundamental matrix solution M (l) [37], each column of which is independent and satisfies the homogeneous part of Eq. (14). Using the method of variation of parameters, the solution to the non-homogeneous initial value problem: can be given by [37]: When it is possible to find out all the possible eigenvalues and independent eigenvectors of A, construction of M (l) is straightforward [38,Ch.37]. However, matrices are not always diagonalizable. So, it is essential to have an alternative method of deriving M (l) when the calculation of all independent eigenvectors is not possible. This can be achieved by Putzer's algorithm [37]. The method is elegant in the sense that it does not require all the eigenvalues to be distinct or nonzero. However, in case of solving Eq. (14) it is seen that the calculation of M (l), a 15 × 15 matrix, becomes impractically lengthy, and therefore, the method has not been adopted. But if it is possible compute M (l) by this method, that may be used even outside magnetized iron plates.

Application to ICAL
In the track fitting program for ICAL [19], thick scatterer approximation has been used previously by implementing Mankel's form of random noise matrix [1]. Strictly speaking, this form of matrix is valid only if the track segment is linear, since the magnetic field dependent terms (that lead to curvature of the track) are assumed to be zero in the propagator matrix F (Eq.(11)). Therefore, it is a matter of interest to see how the performance of track fitting is affected, when a more appropriate solution (Eq. (20)) is applied to construct the random process noise matrix in the presence of the magnetic field.
This has been carried out through the use of a C++ based computational library it++ [39]. Details of the coding techniques etc. are given in Appendix C. It was seen that in all the cases where all the elements of A are non-trivial (which commonly happens within the magnetic field), the determinants of A assume large values (10 0 − 10 6 ) and the diagonalizations can be carried out quite easily. However, in the regions where the magnetic field is zero (outside the iron slabs in the ICAL detector) or its spatial derivatives are zero (inside iron), |Det(A)| = 0 (or |Det(A)| → 0) and Eq.(20) cannot be applied. This can be understood in the following way: outside iron, the propagator matrix reduces to Eq.(11), as all the magnetic field integrals vanish. Even inside iron, certain elements in the first two columns of F matrix (e.g. F 11 , F 12 etc. which depend on spatial derivatives of magnetic field components [19]) become zero occasionally. These zeros lead to additional zeros in the matrix A and the determinant of the latter becomes very small (close to zero) 1 . That the determinant is zero (or close to zero) suggests that one or more eigenvalues are zero (or close to zero). Hence, Eq.(20) cannot be evaluated properly and unphysical solutions are obtained if Eq.(20) is applied. Therefore, outside the iron plates (and occasionally inside the iron plates) where |Det(A)| is small (≤ 1), we used Mankel's form of the process noise matrix Eq. (12). In general, inside the magnetic field, where Det(A) is typically 1, we applied Eq. (20) to construct the elements of the process noise matrix by diagonalizing A through it++.
Since the solutions q i (l)s represent the terms of a covariance matrix, we expect that q 1 , q 6 , q 10 , q 13 , q 15 will be positive, because they correspond to the diagonal elements of the Q matrix (Q 11 , Q 22 , Q 33 , Q 44 , Q 55 respectively). However, the real parts of the solutions q i (l)s need not be positive. It is interesting to see that the computation automatically led to positive values of diagonal elements, as expected. No additional measure was needed to obtain these positive values. This shows that the analysis has been consistent.

Reconstruction performance
Since the method of computing the process noise matrix described in this paper is somewhat abstract, first we would like to show that the resulting Kalman filter works in a consistent fashion. Once that is done, we shall check if the Kalman filter, equipped with the random noise matrix developed in this paper, has better (or worse!) reconstruction performance compared to the one equipped with the random noise matrix derived by Mankel [1].
We used GEANT4 [40] to generate 5000 Monte Carlo muons (µ − ) inside the ICAL detector. The event vertices were smeared uniformly across a volume of (43.2m × 14.4m × 10.0m) around the center of the detector (see Figure 1(a)) in all φ directions (φ ∈ [0, 2π]). This ensures that the muon tracks from the inhomogeneous magnetic field region (Figure 1(b)) are also present in the total set of simulated tracks, in the same way it would happen in reality.
To show that the filter is working in the expected way, we shall present the 'goodness of fit' plots in the following. These are the pull distributions of the fitted variables and the reduced χ 2 distribution. The pull of a given variable ζ is defined as: where C ζζ denotes the error of the reconstructed ζ parameter, as estimated from the updated covariance matrix of the Kalman filter. In ICAL, we are mostly interested in the fitted parameters near the muon event vertex; hence, the pull is evaluated there only. For good fit, the pull distributions must have mean at zero and standard deviation equal to unity. In Figure 2, we show these plots for muons of momentum 5 GeV/c, with initial direction θ = cos −1 0.95 to the vertical. The reduced χ 2 of the model prediction of every event is obtained by dividing the total χ 2 p = k (r k−1 k ) T (R k−1 k ) −1 r k−1 k [12] by the no. of free parameters. Here r k−1 k denotes the residual of state prediction and R k−1 k is the corresponding error covariance matrix. The total no. of free parameters equals 2n − 5, found by subtracting 5 constraints (through initialization of the filter) from the total degrees of freedom (two times the no. of hits n along the track). The χ 2 p for prediction is equal to the χ 2 f [12] for the track fit. From Figure 2, it is observed that the pull distributions of the elements of the state vector have mean very close to zero and fitted width very close to unity. The reduced χ 2 plot (Figure 2(f)) peaks close to unity as well, as expected. Similar performance of track reconstruction is observed in a wide range of p and cos θ. At very low momenta (p < 2 GeV/c) and very large angles θ > 60 0 , gradual worsening of the reconstruction performance is observed. This degradation is intrinsic to the tracking problem, irrespective of whether or not the enhanced track scatterer treatment, described in this paper, is included. The gradual worsening is seen from the following momentum and direction resolution plots in Figure 3. Here, the momentum resolution has been defined as σ(p) p in , where σ(p) denotes the rms width of the reconstructed momentum distribution and p in denotes the input momentum. On the other hand, rms width of the reconstructed cos θ distribution (i.e. σ(cos θ)) has been chosen as the definition of the direction resolution.  At lower momenta and/or larger zenith angles, the muon tracks are affected by the multiple scattering to a greater extent. This affects the precision of muon momentum estimation, for which the resolution becomes poor (Figure 3(a)). On the other hand, at higher muon input momenta, the contribution to the momentum resolution from the spatial resolution component is higher [41,Eq. (3.5)] and that leads to gradual worsening of momentum resolution. The direction resolution steadily improves with increasing p in , but worsens as input θ is increased.
Let us now proceed to the comparative study of the Kalman filters equipped with two different process noise matrices: one derived in this paper and the other derived in [1]. To see if the former has any advantage or disadvantage over the latter, these two programs were used to fit the two copies of the simulated muon tracks of momentum 5 GeV/c and initial direction θ = cos −1 0.95. Thus, the two filters with the two different process noise matrices operated on identical sets of position measurements. It was observed that the quality and performance of reconstruction of these two programs are of the same order. For individual events, the correction in the reconstructed values of momentum or cos θ usually appeared at the second or third decimal places or beyond that. In fact, no significant improvement or deterioration was observed for any of the track parameters. Thus, no gross improvement was achieved by using the more appropriate form of random process noise matrix inside the iron plate equipped with magnetic field. This is shown in the following Figure 4(a) and 4(b).  This observation that there is hardly any difference in the reconstruction performance may raise some doubt about the validity of the process noise treatment. One may be interested to check how large the effect of the process noise treatment is in the first place. If the fitting is performed by switching off the process noise, we expect the fitting performance to deteriorate. This exercise was performed by setting all the elements of the process noise matrix to zero, but keeping the other of the program the same as before. The result is shown in Figure 5. It is seen that the fitting performance of the same set of events becomes less precise, due to inconsideration of the process noise treatment. In fact, many of the events are reconstructed with worse momenta values, as seen from the event count in Figure 5(a) (less compared to those in Figure 4(a)). The filter, however, converges to more or less accurate mean value, since it took into account the mean energy loss in correct manner. On the other hand, the direction estimation becomes very poor, as seen from the width of the distribution in Figure 5(b). This consistency check also confirms that the fitting performance improves significantly with respect to "no process noise treatment", when the process noise matrix is accounted for. Figure 4(a) shows that the formula of the process noise matrix developed in this paper, which was used for tracking inside magnetized iron plates, does not lead to gross improvement of track fitting performance. So, we conclude that Mankel's simple solution for random noise matrix is indeed a good approximation.

Summary
In this paper, a mathematical formalism has been developed for expressing the elements of the random noise matrix while performing track fitting with a Kalman filter through a thick scatterer and nonzero magnetic field. In this case, all the elements of the propagator are nonzero, unlike Mankel's approach [1] and we made use of the method of diagonalization (see 3.1) to construct the desired elements under such circumstances. Through this formalism, the elements cov(x, q/p) = cov(q/p, x) of the random noise matrix can also be calculated for a track deflected by a magnetic field in a thick scatterer. Evaluation of these elements was not included in Mankel's treatment [1]. Although no precaution was taken to render the real parts of q 1 , q 6 , q 10 , q 13 and q 15 positive (which correspond to the diagonal elements of the random noise matrix), they turned out to be positive in all the cases. However, this solution could not be used outside magnetic field region. Also, its use inside the magnetized iron plates did not improve the track fitting performance. The treatment by Mankel [1], derived under approximations, seems good enough for reconstruction of momentum, at least to the first or second decimal place. This is also clear from Table B.1 in Appendix B which shows that the corrections introduced to the elements of the process noise matrix are small. On the other hand, the mathematical form of the elements of the process noise matrix, derived in this paper, is quite general and can be used in the context of any state vector in other HEP experiments employing different state vectors (for example, those containing q/p T or curvature κ of the track as one of the elements).

Acknowledgements
K.B. expresses gratitude to the INO Cell, TIFR for providing short term visiting fellowship for supporting this research. The authors also acknowledge the comments of the anonymous referee for suggesting many improvements in the organization of the paper.

Appendix A.
Let us first formally define the vector q as: q = (Q 11 , Q 12 , Q 13 , Q 14 , Q 15 , Q 22 , Q 23 , Q 24 , Q 25 , Q 33 , Q 34 , Q 35 , Q 44 , Q 45 , Q 55 ) T Hence, q 1 , q 6 , q 10 , q 13 and q 15 represent the diagonal elements of Q matrix. Then, we formally define the vector δq as (see Eq. (7), Eq. (13), Eq. (14)): δq = (0, 0, 0, 0, 0, 0, 0, 0, 0, c(t x , t x ), c(t x , t y ), 0, c(t y , t y ), 0, c(q/p, q/p)) T Next, we construct the matrix A of Eq.(13). This is a 15 × 15 matrix with many non-trivial elements. Hence, we express it by dividing it into two blocks B 1 and B 2 such that A = (B 1 15×8 |B 2 15×7 ). By augmenting these two matrices, we can construct A. The matrix B 1 is given as: Similarly, the matrix B 2 is given as: A typical case of fitting an up-going muon track of momentum 2 GeV/c at initial direction θ = cos −1 0.75 is considered. In a small step inside the magnetized iron plate, we compare the coefficients of the powers of l, l 2 and l 3 to check the consistency between Eq. (12) and Eq. (20).  It is seen from Table B.1 that the corresponding coefficients are close enough. For example, the element Q 24 = 0.032906 l 2 when calculated from Eq. (12). In the modified approach, the value of Q 24 becomes: Q 24 = (10 −17 ) l + 0.032905 l 2 − 0.002956 l 3 where l ∼ 10 −3 m. So, the modification introduces very weak correction. In this case, the reconstructed momentum with the new process noise matrix became |p rec | = 2.259 GeV/c compared to |p rec | = 2.258 GeV/c, estimated by the standard approach. The table also shows the interesting point that the difference between Q 14 and Q 23 happens at the order of l 3 . In Eq. (12), these two independent elements were the same! Appendix C.
For computations in the high energy physics experiments, ROOT [42] is a widely accepted software. However, because of the inevitable occurrence of the complex numbers in this problem, it is rather difficult to implement the recipe of Eq.(20) using ROOT. The reason is following: the actual diagonal matrix becomes block-diagonal in the convention followed by ROOT, since it pushes the imaginary parts of the eigenvalues to the off-diagonal positions (see: 'Matrix Eigen Analysis' in Chapter 14 of [43]). The eigenvector matrix is also kept real in ROOT.
However, we wanted to proceed with the standard diagonalization method for which all the eigenvalues, real or complex, appear at the diagonal posi-tion. Therefore, we used a C++ based library it++ [39]. This library can be easily interfaced with existing code which is written in C++ by appending 'itpp-config --cflags' and 'itpp-config --libs' to LDFLAGS in the GNUMakefile. This library can be easily used to find eigenvalues and eigenvectors of A in the standard forms. The following it++ member function was used: itpp::eig(const mat &A, cvec &d, cmat &P ) to carry out the procedure. In this function, d denotes the complex vector of eigenvalues and P denotes the complex matrix obtained by augmenting the eigenvectors of A. This matrix is seen to have determinant nonzero and thus, is invertible. As required by Eq. (20), the inverse matrix P −1 is made to operate on δq and further computations are performed.
This package is based on external computational libraries, like BLAS [44] and LAPACK [45]. The level of accuracy of the computation is seen to be of the same order as of Mathematica [36]. For example, the eigenvalues and eigenvectors of a matrix computed by it++ and Mathematica are found to be consistent within ∼ 1%.