Fadi Zeini Paul Avery CBX 95–34 July 23, 1995

UFTRAK: New Energy Correction for Charged Tracks

1 Introduction

When a particle traverses a detector, its path deviates from a true helix because of interactions with material, most notably energy loss and multiple scattering. Track fitting programs which neglect these effects do progressively worse at low momenta where the effects are most pronounced. Even DUET, which accounts for multiple scattering at the VD–DR interface and includes a post-hoc correction to account for other sources of multiple scattering, performs poorly at low momentum, partly because of its neglect of energy loss.

Multiple scattering and energy loss are functions of velocity, not momentum, so a correct accounting of material interactions necessitates fitting the track 5 times, once each for the hypotheses e, μ, π, K, and p. An efficient algorithm is the Kalman filter[1], which uses an iterative procedure to trace the track from the outermost hit back to the point of production. The algorithm is surprisingly simple, and consists of the following repeated steps: (1) transport the track parameters and covariance matrix to a hit or material layer; (2) if material, add the effects of energy loss to the track parameters and multiple scattering to the error matrix; (3) if a hit, average the track parameters with the hit to get an updated track and covariance matrix.

Several experiments have adopted some form of Kalman filter algorithm, and CLEO is expected to have one for the next recompress, thanks to Rob Kutschke’s work. In the meantime, we must correct the DUET track parameters for energy loss after the fact, the subject of this note. The routines described here improve the algorithm[2] used since 1991, since they correct, in addition to momentum, the impact parameter (DBCD) and angle (FICD) as a function of transverse momentum. The new subroutines are in libuftrak.a.

2 The previous dE/dx correction

A full description of the previous dE/dx correction can be found in CBX 92–40[2], while the routines are in the ftlib library. The graphs in the CBX note showing the K_{s} and Λ masses before and after the correction demonstrate quite clearly its importance. For the K_{s}, the mass shifts by 5 MeV! Comparable mass shifts were seen for the D^{0} and D^{+}. The plots also show that the masses appeared to be relatively constant with momentum after correction.

Although it was reasonably successful, the energy correction ignored effects on phi shift and impact parameter. These routines called the low level function DKSDEL (written elsewhere) which actually carried out the correction. This routine solved the dE/dx differential equation in a slightly different way. Unfortunately DKSDEL contains two errors. First there is a discontinuity at p/m = 1.8, where different algorithms are used for momenta on either side of that point. The discontinuity occurs at p = 0.25, 0.89 and 1.69 GeV, for π, K and p, respectively. Second, the correction incorrectly saturates for p/m > 10 because the table entries run out for these momenta. Thus, for pions with p > 1.4 GeV the routine underestimates the true correction.

3 Energy Losses

3.1 Effect and overall strategy

The objective of LIBUFTRAK is to provide the user with a collection of subroutines that calculate energy loss corrections to the DUET track parameters. The interface between the end user and the library subroutines is carried out through the subroutine UFTCOR which, given a track CD number, returns the corrected parameters, namely, DBCD, FICD, CUCD, Z0CD,

and CTCD.

As a particle traverses a material thickness, its energy losses lead to a decrease in its velocity. Consequently, the change in momentum results in an observed curvature that deviates from the true curvature of the track. As Figure 1 shows, a track with a certain measured curvature also has a measured impact parameter with respect to the beam spot, DBCD, that differs from its true value[3]. This difference constitutes the shift in DBCD. Other measured parameters, not shown in the plot, are also shifted from their true values.

Figure 1: effects of energy loss on the track curvature and impact parameter. The momentum-dependent shift also depends on the sign of the track; a track with a positive impact parameter is shifted towards the beam spot while a negative-DBCD track shifts away from the beam spot. These effects are more pronounced at low momentum.

For our analysis, we categorized tracks based on their “miss-category”; the fraction missing of all DR hits used in reconstructing the track. All tracks were divided into four groups,

(a) 0 to 25% missing hits in the DR, (b) 25% to 50%, (c) 50% to 75% and, (d) tracks with a miss category greater than 75%. Figure 2 is a general plot of the fraction of hits missing for tracks of all momenta. The fraction of hits missing was calculated as the ratio of the number of missing DR hits to the total number of layers in the Drift Chamber. To determine the number of missing hits, a maximum track radius was found based on the track’s curvature, CUCD, or, for tracks exiting the endplates, the cotangent of the polar angle, CTCD. The layers with no hits (using KBITCD) are then counted up to, and including, the layer closest to the maximum radius (to within 2mm). Figure 3 shows the probability that a track falls into one of the four categories discussed above, in bins of log transverse momentum. This was determined by averaging the number of tracks in each ’miss- category’ as a function of logp_{⊥}. The plot clearly shows that the low momentum region is dominated by tracks reconstructed from a high ’miss’ ratio (more than 75% missing hits), whereas, as expected, high-momentum tracks register more hits on average. The higher probability of category 2 tracks (25% to 50% missing hits) over category 1 tracks (0% to 25% missing hits) in regions of high momentum has to do with hardware efficiency (drift distances, pulse times, ... etc).

GEANT calculates energy losses based on the particle’s momentum. To account for the deviation from a constant ionization energy, we compared the observed momentum of a track to its Monte-Carlo generated value, using the difference as an indication of the decrease in the particle’s momentum as it crosses the beam material. Figure 4 shows the difference between the average observed and generated momenta for a sample of 4s3 continuum MC tracks. Consequently, a minimum energy loss that is consistent with the momentum decrease is found, as explained in the next section. This (dE/dx)_{min} is a function of momentum and is grouped into one of the four miss-categories discussed above.

As explained in section 3.3, the Bethe-Bloch formula provides a relation between Δp, the change in momentum due to energy losses, and (dE/dx)_{min}.

The track is projected onto the x-y plane and the Bethe-Bloch formula is integrated to obtain a value for the track’s transverse momentum. Once Δp_{⊥} was calculated, the corresponding shifts in the track parameters can be determined. For a helical trajectory, as shown in ref. 4, the momentum components can be written as

p_{0x} = p_{⊥} cos ϕ_{0} = \frac{a}{2c} cos ϕ_{0}
p_{0y} = p_{⊥} sin ϕ_{0} = \frac{a}{2c} sin ϕ_{0}
p_{0z} = p_{⊥} cot θ = \frac{a}{2c} cot θ

where p_{0x}, p_{0y}, and p_{0z} are the components at the point of closest approach to the origin, c is the curvature and a is a constant equal to -0.2998BQ where B is the magnetic field and Q is the charge of the track. ϕ_{0} and θ are the azimuthal and polar angles, respectively.

From the above equations, the curvature can be related to the transverse momentum,

p_{⊥} = \frac{a}{2c}
Δp_{⊥} = −\frac{a}{2c^{2}} Δc
(Δp_{⊥}/p) = −(Δc/c)

Thus, the change in a track’s curvature, Δc, is found from the change in transverse momentum. Based on the change in curvature, changes in the azimuthal angle, FICD and the impact parameter, DBCD, are calculated. The changes are then added to the old values to give the corrected parameters returned by UFTCOR.

3.2 Analysis procedure

A sample of approximately 10^{6} 4s3 continuum Monte-Carlo events was used. Track selection was based on “good” tracks, with KINCD=0 and a cut on the Z-coordinate at point of closest approach of |Z0CD − Z0B| ≤ 0.02 − 0.04, where Z0B is the beam’s z-position, bmpos(3,1).

The momentum dependence of the shift in the impact parameter with respect to the beam spot was determined by plotting the mean value of DBCD as a function of momentum. The plots were fitted to a single Gaussian plus a Chebychev polynomial. The results of the fits were then plotted in bins of logp_{⊥}. Figure 5 compares the shifts in the DBCD mean value, as a function of logp_{⊥}, between 4s3 continuum MC and 4s3 data as well as 4s5 data. The three values are close enough to justify the use of Monte-Carlo generated tracks in determining corrections to energy losses. A similar procedure was carried out to determine the shift in the azimuthal angle, Δϕ.

3.3 dE/dx Correction

3.3.1 dE/dx thickness

The formula used by GEANT to calculate energy losses gives a result different from the constant (dE/dx)_{min} of 1.92 MeV. In order to account for that difference, (dE/dx)_{min} is modified to a slowly varying function of momentum. This is achieved by using the Bethe-Bloch formula

dE/dx = (A/β^{2})(log(\frac{β^{2}γ^{2}m_{e}c^{2}}{I_{0}} ) − β^{2})

where β and γ are the usual relativistic parameters, I_{0} is the ionization potential parametrized by I_{0} = 16Z^{0.9} where the atomic number Z = 9 is used as an approximate average for the materials used, m_{e} is the electron mass and A is a constant that depends on (dE/dx)_{min}, the quantity listed in the Particle Data Book. The minimum of this function occurs at βγ = 3.4, which yields an A value of A = (1/11.528)(dE/dx)_{min}. By using (dp/dx) = (1/β)dE/dx we get the following relation between (dE/dx)_{min} and the change in momentum

(dp/dx) = (1/β^{3})(dE/dx)_{min} \frac{1}{11.528}(2logβγ + 9.0872 − β^{2})

which we can use to solve for (dE/dx)_{min} by integrating (dp/dx) using Simpson’s rule. The integration is performed over the range between the observed momentum and the Monte-Carlo generated value. By plotting the (dE/dx)_{min} solution and fitting the plots with a third order Chebychev polynomial we parametrize the momentum-dependent variation. Figure 6 (a,b,c and d) is a plot of the (dE/dx)_{min} fits in bins of logp_{⊥}, as a function of the four miss categories discussed earlier.

The effective dE/dx thickness is thus calculated based on

ΔE_{min} = 1.92MeV/ sin θ ∗ F (p_{⊥})

where

F (p_{⊥}) = a + b × logp_{⊥} + c × (2(logp_{⊥})^{2} − 1) + d × logp_{⊥} × (4(logp_{⊥})^{2} − 3)

and a, b, c and d are coefficients determined from the fit (see figure 6). The function P is a third-order Chebychev polynomial. The 1/ sin(θ) term corrects for the increase in the path length due to non-normal tracks.

3.3.2 Effective Radius, R_{eff}

The DKSDEL routine used previously applied dE/dx track corrections at the origin with no change in the track’s position and direction. The new LIBUFTRAK subroutines, instead, apply the corrections at R_{eff}, the effective point at which the particle is assumed to undergo energy losses. The effective radius was calculated using

Δϕ = −2RΔC

where Δϕ is the shift in the azimuthal angle of the track and ΔC is the change in curvature caused by energy losses (see section 3.1). As mentioned earlier, The average change in the azimuthal angle, Δϕ, was determined by fitting plots of FICD (in bins of logp_{⊥}) and plotting the results of the fits divided into the four miss-categories, as shown in Figure 7.

As a check, the effective radius was calculated using

ΔD = R^{2}ΔC

where ΔD is the change in the impact parameter with respect to the beam spot, DBCD. Both calculations agreed allowing us to find the effective radius using the shifts in both the impact parameter and the azimuthal ϕ,

R_{eff} = \frac{−2δD}{δϕ}

Figure 8 is a plot of the effective radius as a function of the log of the transverse momentum, logp_{⊥}, for all tracks divided into the four miss categories. As shown in the Figure, R_{eff} is momentum-dependent and not a fixed quantity which justifies our deviation from previous corrections. This modification compensates for the change in a track’s position and direction depending on it’s momentum.

3.3.3 Correcting the momentum

The Bethe-Bloch formula is integrated to find ΔP, the change in momentum due to energy loss, as illustrated in Appendix A. The dE/dx thickness calculated in section 3.3.1 is used for the ΔE_{min} mentioned in the appendix. A corrected momentum is then calculated using the formula

P_{corr} = P_{old} + ΔP

This is the actual value of the track’s momentum before energy losses. The B field correction is already applied in GETCD and is included in the P_{old} term.

3.3.4 Correcting for shifts in CUCD, DBCD and FICD

Once ΔP is found, the change in curvature is calculated from the relation

Δp/p = −Δc/c

where C and P are the measured curvature and momentum of the track, respectively (post dE/dx). From the change in curvature, the azimuthal angle correction is then computed using

Δϕ = −2R_{eff} × Δc

where R_{eff} is the effective radius discussed in section 3.3.2 . Using the same curvature correction and effective radius, the impact parameter correction is then found from

DDA = Δ(DBCD) = Δc × R_{eff}^{2}

Figure 9 (a, b, c and d) is a comparison of the uncorrected and corrected impact parameter, DBCD, in bins of logp_{⊥} (divided into four miss-category plots), for 4s3 continuum MC. The procedure for obtaining the uncorrected DBCD mean value was discussed earlier (see sect. 3.2). The corrected values were determined by adding the corrections obtained (as described above) to the old values and then fitting the plots as was done before.

Figure 10 is a general plot comparing 4s3-data uncorrected and corrected DBCD mean values as a function of logp_{⊥} for all miss-categories combined.

4 Subroutines

* indicates output from the subroutine (write)

SUBROUTINE UFTCOR(ICD, MASS, PCR, CUC, FIC, DBC, Z0C, CTC)

Apply dE/dx corrections to a CD track and return the corrected momentum, curvature, phi, impact parameter, point of closest appraoach in Z-axis, and the cotangent of the polar angle, theta.

ICD integer CD track number
MASS real  Assumed mass of particle
PCR* real  Corrected momentum
CUC* real  Corrected curvature
FIC* real  Corrected phi angle
DBC* real  Corrected impact parameter
Z0C* real  Corrected point of closest appraoach in Z-axis
CTC* real  Corrected cotheta

SUBROUTINE DBSHFT(ICD, MASS, RADIUS, P, DCU, DPHI, DDA)

Apply dE/dx corrections to a CD track and compute the updated momentum, and the changes in curvature, phi and impact parameter.

ICD  integer CD track number
MASS real assumed mass of the particle
RADIUS real  reserved for future use - radius at which to apply dE/dx corrections
P*  real  corrected momentum
DCU*  real  change in curvature
DPHI* real  change in azimuthal angle, ϕ
DDA*  real  change in the impact parameter

SUBROUTINE DBSHF4(ICD, MASS, RADIUS, P4, XPOS)

Apply dE/dx corrections to a CD track and put the answer in a 4-vector and a position.

ICD  integer  CD track number
MASS  real  assumed mass of the particle
RADIUS real  reserved for future use - radius at which to apply dE/dx corrections
P4*  real array corrected 4-momentum
XPOS* real array Output 3-position relative to BMPOS(1-3)

SUBROUTINE DBSMOM(PIN, MASS, CZ, POUT, RUNN)

Applies dE/dx corrections to the momentum of a particle, returning the new momentum.

PIN real Input momentum
MASS real  assumed mass of the particle
CZ real Z direction cosine used to compute path length through material
POUT* real  Output momentum after energy loss
RUNN integer Run number (material may be time dependent)

REAL FUNCTION DBDPUP(P, XMAS, DEDXMN, FRAC)

Returns momentum before encounter with dE/dx material (i.e., the correct momentum before any energy losses). The momentum is calculated using the Runge-Kutte method, discussed in Appendix B.

P  real Momentum after energy loss
XMAS  real Assumed mass of particle
DEDXMN real (dE/dx)_{min} × thickness of material
FRAC  real If fractional momentum change is < FRAC then just calculate momentum

assuming constant dP/dx, i.e., avoid full Runge-Kutte calculation.

REAL FUNCTION DBDPDN(P, XMAS, DEDXMN, FRAC)

Returns momentum after encounter with dE/dx material. The momentum is calculated using the Runge-Kutte method, discussed in Appendix B.

P  real Momentum before energy loss
XMAS  real Asssumed mass of particle
DEDXMN real (dE/dx)_{min} × thickness of material
FRAC  real If fractional momentum change is < FRAC then just calculate momentum

assuming constant dP/dx, i.e., avoid full Runge-Kutte calculation.

REAL FUNCTION DBPCOR(P, MASS, DEDXMN)

Applies dE/dx correction to the momentum of a particle, returning the new momentum.

P  real Input momentum
MASS  real Assumed mass of the particle
DEDXMN real (dE/dx)_{min} × thickness of material

REAL FUNCTION DBRCOR(PT)

Compute effective radius vs PT

PT real Input transverse momentum

In addition to the above subroutines, libuftrak also includes the following four utility subroutines written by Jorge Rodriguez:

REAL FUNCTION DBCDCT(ITRAK)

Returns the number of sigmas away from the mean, the value of ITRAK’s DBCD is. The value is normalized by the measured spread in DBCD.

ITRAK  integer Track number
DBCDCT* real  Number of sigmas calculated

REAL FUNCTION DBSIG(ITRK)

Returns the value of the rms spread in DBCD as a function of momentum at FICD = 0. The parameterization that gives the momentum dependence is done separately for each data set (4s2-4s8). If the run number of ITRAK is not found the parameterization for the 4s2 data is used. The function also checks the Monte Carlo flag LMONTE. If LMONTE = TRUE the parameterization used is that found by the fit to 4s2-4s4 Monte Carlo data.

ITRK  integer Track number
DBSIG* real  Value of the spread for PQCD(ITRK) at FICD=0.

REAL FUNCTION DBSIGFI(ITRAK)

Returns the spread of DBCD as a function of momentum and track angle FICD.

ITRAK  integer Track number
DBSIGFI* real  Value of the spread for PQCD(ITRK) as a function of FICD.

REAL FUNCTION DBMEAN(ITRAK)

Returns the mean value of DBCD for track ITRAK as a function of PQCD.

ITRAK  integer Track number
DBMEAN* real  Mean value of DBCD

A Elementary Theory of Charged Particle Energy Loss

The corrected momentum p_{cor} is calculated from the initial momentum p_{old} by the formula

p_{cor} = Rp_{old} + Δp  (1)

where R is a B field correction (it is applied in GETCD) and Δp is the dE/dx correction. Δp is calculated by integrating the Bethe-Bloch formula

Δp = ∫_{0}^{t} \frac{dp}{dx}dx = ∫_{0}^{t} \frac{1}{β} \frac{dE}{dx} dx = ∫_{0}^{t} \frac{A}{β^{3}} [ln (\frac{2m_{e}c^{2}β^{2}γ^{2}}{I_{0}} ) − β^{2}] dx  (2)

where t is the thickness of the material to be traversed (including the increased path length due to non-normal tracks), A is a constant that depends on the material, Z is the atomic number, m_{e} is the electron mass, I_{0} is the ionization potential and β and γ are the usual relativistic parameters. I use the approximation I_{0} = Z^{0.9}eV, and take Z = 9 for representative media. A is determined from our knowledge of (dE/dx)_{min}, the quantity listed in the Particle Data Book. The minimum occurs at βγ = 3.40, which gives A = (dE/dx)_{min}/11.528. It is easier to write the differential equation in terms of τ ≡ βγ = p/m and the normalized thickness x′ = x/t:

\frac{dτ}{dx′} = \frac{ΔE_{min}}{m} F (τ)
F (τ) = \frac{(1 + τ^{2})^{3/2}}{11.528τ^{3}} [9.0872 + 2 ln(τ) − \frac{τ^{2}}{1 + τ^{2}} ]  (3)
τ = \frac{p}{m}

where ΔE_{min} = (dE/dx)_{min}t is the average energy loss suffered by a minimum ionizing particle passing through a distance t of material (including path length corrections from non-normal tracks). The change in momentum can then be calculated from the integral

Δp = mΔτ = ΔE_{min} ∫_{0}^{1} F (τ)dx′  (4)

For moderate velocities (β > 0.7), F (τ) is approximately constant, giving Δp ≃ ΔE_{min}F (τ). However, when β, τ ≪ 1 and F (τ) changes rapidly, it is usually necessary to integrate the dτ/dx′ differential equation numerically using a method such as the Runge-Kutte technique. The integration is aided by transforming to a new variable, ξ = τ^{4}, which removes the 1/τ^{3} term from the differential equation and makes the numerical solution more accurate. Using this variable the dE/dx differential equation becomes

\frac{dξ}{dx′} = \frac{ΔE_{min}}{m} F′(ξ)
F′(ξ) = \frac{(1 + \sqrt{ξ})^{3/2}}{4 · 11.528} [9.0872 + \frac{1}{2} ln(ξ) − \frac{\sqrt{ξ}}{1 + \sqrt{ξ}} ]  (5)
ξ = (\frac{p}{m})^{4}

B Calculation of energy loss using Runge-Kutta method

The Fourth-order Runge-Kutta method requires four evaluations of the function per step. The function is evaluated once at the beginning of the interval, x_{n}, twice at some trial midpoints, x_{n} + \frac{1}{2}h, and once at the endpoint of the interval, x_{n} + h. A weighted mean of the four values is then computed and a value of the function is found based on the weighted mean. As an illustration, assume we have some function representing y as a function of x. For convenience, let K_{1} = F (x_{n}, y_{n}), the value of the function at the initial point. Next, we choose

K_{2} = F (x_{n} + \frac{1}{2}h, y_{n} + \frac{1}{2}hK_{1})

where h is the interval size. We define a third value K_{3}, such that

K_{3} = F (x_{n} + \frac{1}{2}h, y_{n} + \frac{1}{2}hK_{2})

The function is evaluated a fourth time at the endpoint of the interval. Define a new K_{4}

K_{4} = F (x_{n} + h, y_{n} + hK_{3})

to be the value of the function at the fourth point. Thus, K_{1}, K_{2}, K_{3} and K_{4} represent the values of the function at four different points, The weighted mean of the four values is

K = \frac{1}{6}(K_{1} + 2K_{2} + 2K_{3} + K_{4})

Note that the error is of fifth order in the step size, h, which makes the fourth-order Runge-Kutta a very good approximation to use (a method is called n^{th} order if its error term is O(h^{n+1})). In our application, the method is used to compute ΔP in order to get P + ΔP. The change in momentum as a result of a change in energy (energy loss) is calculated. Thus, using the above illustration, y is P and dy is dP. The value (1/m)(dE/dx)_{min} is used as a step (h in the illustration above), where m is the assumed mass of the particle and (dE/dx)_{min} is the average energy loss per unit distance suffered by a minimum ionizing particle, as quoted in the Particle Data Book (see Appendix A).

References

[1] P. Avery, CBX 92–39, August, 1992.

[2] P. Avery, CBX 92–40, August 1992.

[3] J. L. Rodriguez, CBX 94–70a, February 1995.

[4] P. Avery, CBX 92–45, August 1992.

[5] R. Kutschke, CSN 94–334, December 1994.

Figure 2: Tracks vs. fraction of DR hits missing when reconstructing a track. A peak around 0.23 indicates that most of the tracks were reconstructed with 23% of the hits missing in the drift chamber.

Figure 3: The fraction of tracks in each of the four miss-categories as a function of momentum 13 Figure 4: The difference between MC generated momentum, PGEN, and a track’s momentum, PQCD, as a function of logp_{⊥}. This difference is used as a guide to find the minimum energy loss for a particle traversing material Figure 5: A comparison of DBCD shifts between 4s3 continuum MC, 4s3 data and 4s5 data

15 Figure 6: the variation of (dE/dx)_{min} with logp_{⊥} divided into the 4 miss-categories

16 Figure 7: change in phi as a function of logp_{⊥}, categorized by miss-category

17 Figure 8: Effective radius as a function of logp_{⊥} for tracks with a miss category of 25% to 50% Figure 9: 4s3-continuum-MC DBCD, before and after dE/dx correction, in bins of logp_{⊥} divided into the four miss categories Figure 10: A comparison of uncorrected and corrected DBCD for 4s3 data in bins of logp_{⊥}

20