Tracking in Cylindrically Symmetric Magnetic Fields
By Jim Thomas
Introduction
Tracking charged particles in a cylindrically symmetric magnetic field relies upon our knowledge of the integrals and evaluated along the paths of the particle trajectories. I will show that the angle between the particles final direction of motion and its initial direction of motion is proportional to while the angle between the particles final direction of motion and a line drawn between the origin and the end point of the track (e.g. a detector hit) is proportional to . Clearly, is relevant when a detector configuration includes particle tracking before and after the magnet while applies to detector configurations where the tracking detectors are placed after the magnet.
Tracking before and after a magnetic field region is the most powerful method of spectroscopy. In this case, the displacement of the particle’s endpoint from the initial trajectory endpoint (i.e. the infinite momentum ray) is proportional to and this is an excellent experimental observable. (R is the value of the radius coordinate at the endpoint of the track.)
These observations are important in designing tracking algorithms for modern high energy and nuclear physics detectors. They are also important in designing magnet mapping devices for these detectors because it makes it clear what to examine when asking the question "how well must the field be mapped to achieve the required spectrometer resolution?".
Feynman Finger Physics in a Cylindrically Symmetric Field
A cylindrically symmetric field is defined by the quantities (r,f,z) and B(r,f,z). The only symmetry I assume is in f. It is not necessary to postulate a field uniform in z.
Consider Fig 1. Detectors are placed inside and outside of a cylindrically symmetric magnetic field volume. In this example, the strongest field is near the center of the volume ... but this is not required in the following discussion. The field is uniform in f but not necessarily uniform in z. The inner detectors are at » 0 radius. The outer detectors are placed at radius R.
Figure 1: A typical magnetic spectrometer
A particle is produced at the center of the detector and travels outwards. It passes through the outer detectors at an angle a; which is the angle between the final momentum vector and a line drawn from the origin to the point where the particle hit the detector.
Similarly, g is the angle between the final momentum vector and the initial trajectory of the particle. The displacement of the particle’s track, measured on the outer detectors, is R f or approximately R (g-a). (The latter equality is strictly true only for the case where the magnetic field is uniform in z. See below.)
Now consider a differential element of the magnetic field volume as shown in Fig 2.
Fig. 2: Particle trajectories through a differential element of the field volume
The particle receives an infinitesimal pt kick from the field element located at radius l. We can calculate the kick using the equations of motions of a particle in a uniform field:
where B^ is the component of the field that is perpendicular to the particles instantaneous direction of motion. Thus, assuming the field is constant inside the differential volume element, it is clear that dg = r / dl . Strictly speaking, l is the path length along the particle trajectory and will differ from the radius, r, especially for low momentum particles or very strong magnetic fields. dl is the differential path length which connects the incoming track with the outgoing track at the boundaries of the differential field volume. It is a segment of a circle. So,
Using the law of sine’s, it is also easy to show that
and the total displacement of the track, from the infinite momentum ray, at radius R is R(g -a).
These equations are correct in the limit of small angular deflections of the track and for high momentum particles. They are only approximations in the case of large angular deflections and lower momentum particles. But the equations show that there are well-defined relationships between a, g and the integrals of the field. So a knowledge of a and/or g is equivalent to knowing the transverse component of the particles momentum in all cases.
These equations are also good approximations for high pt particles in the case where the B field is changing slowly in Z. This is because both the Z component of the field and the R component of the field give a kick in the f direction. So a particle follows a helical trajectory where the radius of the helix changes slowly with Z. This simple helical motion is only disrupted when the f component of the velocity builds up to an appreciable level and the Vf x B terms compete significantly with the VR x B and VZ x B terms. And since Vf = 0 at the origin, by definition, it takes a while for this velocity component to build up to an appreciable level.
The equations are correct in the general case of a cylindrically symmetric magnetic field with BR and BZ components, however, we must be careful to integrate over the component of the B field that is perpendicular to the particles instantaneous direction of motion. Thus B^ is used in the equations above as a reminder of these general considerations.
Finally, I should note that a and g are the angles between vectors and I have not projected them into the ( r, z ) plane. Some authors do this.
A More Rigorous Approach
There is a simple way to prove that Beta is proportional to Int(B.dl) using a blackboard and finger physics. You can also show that Alpha is proportional to Int(B.l.dl) with this technique.
But I thought it would be fun to try to prove the case using analytic methods. Starting from the Lorentz force equation, it is simple to show the following two relations.
The position vector is defined to be:
R(s) = R0(s) + s x t0 + Q/pc x INT{ds' x (s-s') x t(s') X B(s')}
and the tangent vector t is given by:
t(s) = t0 + Q/pc x INT{ds' t(s') X B(s') }
where cap X means cross product. The "tangent vector" is shorthand for p-bar/p or v-bar/v where p-bar stands for the vector P and p is its magnitude.
Now take the cross product between the initial direction tangent vector and the final direction tangent vector to define the angular shift between the initial and final directions:
t0 X t(s) = Q/pc x INT{ds' t0 X [t(s') X B(s')] }
and note that t(s) X B(s) is what we normally call B-perp, the component of the field that is perpendicular to the direction of motion of the particle.
= Q/pc x INT{ds' t0 X B-perp}
Due to the phi symmetry of our field, t0 (the initial direction unit-vector) is perpendicular to B-perp to second order and so the magnitude of the cross product is merely B-perp. You will note that the magnitude of t0 X t(s) defines beta (I stick to this notation, with apologies to Shiva and Alexi) and so if:
sin(beta) = Mag{ t0 X t(s) }
then
beta = Q/pc x Int{ds' B-perp}
QED. This is otherwise known as Int(B.dl). This conclusion is true for any phi symmetric field because both the Z component of the field and the R component of the field give kicks in the phi-hat direction. And it only breaks down when the velocity of the particle becomes large in the phi-hat direction. In this case the phi-hat component of the velocity interacts with the Z field (and R field) yielding some strange focussing at low momentum in the CM. In the MM, these effects tend to pull the muons off a cone of constant theta. But these effects are small at energies above 400 MeV or so and so I am comfortable neglecting them.
Next, take the cross product between t(s) and R(s). This defines the angular displacement between the final tangent vector and a line drawn between the vertex and the point the particle hit the detector. Be careful about the order of the cross products and many terms will cancel:
t(s) X R(s) = Q/pc x Int{ds' x s' x [ t0 X [t(s') X B(s')] ]}
+ Q**2/(pc)**2 Int{ds" Int{ds' (s-s') [B-perp(s") X B-perp(s')]}}
I got lazy there at the end and identified t(s) X B(s) as B-perp. The interesting thing about the second term is that it goes to zero in the case of phi symmetry because B-perp(s") and B-perp(s') point in the same direction until the phi-hat velocity builds up to an appreciable amount. This term also goes to zero at large momentum. Thus,
t(s) X R(s) = Q/pc x Int{ds' x s' x [ t0 X B-perp ] }
and if you use all the same arguments, and limitations, noted above you will immediately recognize that
sin(alpha) = Mag{ t(s) X R(s) }
and
alpha = Q/(Rpc) x Int{B.l.dl}
Where R is the magnitude of R(s-final), the total distance along the path, or, to a good approximation the straight line distance from the vertex to the hit on the detector. QED II.
This establishes a direct relationship between Beta (the angle between the initial direction and the final direction) and Int(B.dl). It also shows that Alpha (the angle between the final direction and a line drawn between the vertex and the hit on the detector) is proportional to Int(B.l.dl). Both angle relationships are valid to second order. Note that alpha depends on R and that R is constant on a cone of constant theta, where theta is the emission angle from the vertex with respect to the beamline. To a good approximation, our particles are confined to a cone of constant theta and only the phi-hat velocity components work to pull trajectories off that cone.
The utility of these angles is that they are amenable to Feynman finger physics. As such, I like them. It also turns out that Alpha is what we measure in the CM. We don't know the initial direction vector and so we can only measure the angle between the track, as it goes through the detector, and a line between the vertex and the hit point. Fortunately, this angle is proportional to Int{B.l.dl} and we have a good spectrometer after all.
I also like to quote alpha and beta in the Muon spectrometer because they tell a story about the magnet. One is proportional to B.dl and the other to B.l.dl. You would be surprised how different these are due to the large changes in the B field through the MM. One emphasizes the front of the spectrometer and the other emphasizes the back of the spectrometer. In this capacity, they are independent and useful information and a knowledge of both improves our ability to measure the momentum of particles inside the MM.
Another discussion
Let a particle traverse a magnetic field region. Starting from the Lorentz force equation:
F = dP/dt = Ze/c ( V X B )
Integration of this equation along a path yields:
t(L) = t(0) + Ze/pc INT{ dl' t(l') X B(l') }
Cap X means cross product. "t" is the tangent direction vector and it is shorthand for v-bar/v; where v-bar stands for the velocity vector V, and v is its magnitude. The integrals are evaluated from 0 to L, where L is equal to INT{ dl' } or equivalently L is the time-of-flight divided by the invariant velocity.
Since [ t(l) X B(l) ] is what we usually call B-perp (the component of the field that is perpendicular to the direction of motion of the particle), we can write this in the more familiar form:
1> t(L) = t(0) + Ze/pc INT{B dl'}
It is interesting to note that the angle between the initial direction vector and the final direction vector is "gamma" and is precisely defined by:
( 1 - cos(gamma) ) * 2 = Magnitude[ Ze/pc INT{ B dl' } ]**2
or in the small angle approximation (magnitude of the vector implied below):
2> gamma = Ze/pc INT{ B dl' }
Integrating equation 1> again, and interchanging the order of integration yields the position vector R(L):
R(L) = R(0) + L t(0) + Ze/pc INT{ dl' (L-l') t(l') X B(l') }
or more colloquially:
3> R(L) = R(0) + L t(0) + Ze/pc [ L INT{ B dl'} - INT{ B l' dl'} ]
Using equations 1> and 3>, it is easy to show that the angle between the final direction vector and the vector defined by the initial and final points ( R(0) and R(L) ) is "alpha":
( R(L)-R(0) )**2 + L**2 - 2 L ( R(L)-R(0) ) cos(alpha) =
Magnitude[ Ze/pc INT{ B l' dl' } ]**2
or in the small angle approximation
4> alpha L = Ze/pc INT{ B l' dl' }
Similarly, the displacement of the track from its initial direction vector endpoint can be calculated using the angle "phi":
( R(L)-R(0) )**2 + L**2 - 2 L ( R(L)-R(0) ) cos(phi) =
Magnitude[ L Ze/pc INT{ B dl' } - Ze/pc INT{ B l' dl' } ]**2
so in the small angle approximation the displacement is:
5> phi L = L Ze/pc INT{ B dl' } - Ze/pc INT{ B l' dl' }
Equations 1-5 form the basis for all tracking and momentum reconstruction algorithms; either directly, or as sums and differences of these equations, or projected onto detector planes. Note that the magnitude of the momentum vector is conserved in a static magnetic field and so p has been placed outside the integrals .... but this is true only for a particle in a vacuum! For particles traversing real detectors there will be multiple coulomb scattering and the momentum is progressively decreasing with each time step. Fortunately, MCS is not difficult to calculate with modern computational methods and p (the magnitude of P) can be placed inside the integrals. Then these equations are correct even with MCS.
Examples: Consider measurements made in an apparatus where there is only one detector and it is outside the magnetic field region. Only the final position and direction of the particle is measured at the external detector; the initial direction vector is not known. By substituting equation 1> into 3> we can derive an equation that includes only experimental data and no unknowns.
6> R(L) = R(0) + L t(L) - Ze/c INT{ 1/p B l' dl'}
Note that alpha, in equation 4>, is the angle between ( R(L) - R(0) ) and L t(L). This is the angle between the final direction vector and a line drawn between the vertex and the "hit" on the external detector. In real applications, t(L) is usually projected into the (x,y) plane and the (y,z) planes yielding two parameter fits to the momenta of tracks.
A second example: Equation 3> is the defining equation for the "floating-wire" method where a wire under tension is floated inside a magnet. The deflection of the wire is proportional to the tension and current in the wire according to equation 3>. However, the effects of MCS are not measured by this technique.
As a final application of these equations, consider the case where there are three detectors in a magnetic field. One before the field region, one in the middle, and one at the end. Further, suppose that these detectors do not measure the angles of the particles as they pass through the detectors but only report the position where the track passed through the detector planes. In this case neither equations 1> or 3> are helpful. However, equation 3> can be applied again to a track traversing from detector 1 to detector 2 and then applied to the case of a track traversing from detector 3 to detector 2 (in the reverse direction). This allows us to eliminate the direction vector at detector 2 and yields:
R(3) - (1+L2/L1) R(2) + L2/L1 R(1) =
Ze/pc INT(3,2){B l' dl'} + L2/L1 Ze/pc INT(1,2){B l' dl'}
This is the expression for two times sagitta; the displacement of a track, at detector 2, from the infinite momentum ray. Note that this elegant definition is the distance between the "hit" on detector 2 and the mid point of a straight line drawn between the "hit" on detector 1 and the "hit" on detector three. The difference between the integrals of B.l.dl and B.dl is twice the sagitta if L1 is equal to L2 and the bending power of the magnets is the same between detectors 1 and 2 and between 2 and 3. In all other cases, the intuitive definition of sagitta (i.e. projected onto detector 2) is momentum dependent. As a further complication, we cannot evaluate the integrals "in the reverse direction" from which the particle flies because the effects of MCS are not time reversal invariant. We can, however, rearrange terms and evaluate the integrals in the forward direction to find:
7> R(3) - (1+L2/L1) R(2) + L2/L1 R(1) =
L2 Ze/c INT(2,3){1/p B dl'} - Ze/c INT(2,3){1/p B l' dl'}
+ L2/L1 Ze/c INT(1,2){1/p B l' dl'}
Therefore, sagitta depends on knowing the integrals of B.dl and B.l.dl independently.
Remember that B inside the integrals really means B-perp and the integrals are vector equations. This makes it tedious to calculate the integrals directly in a tracking code because B-perp needs to be evaluated at each time step.
Instead of calculating INT{B dl} and INT{B l dl}, here is a better idea. Use a tracking code with a good B field map to simulate tracks. Then use the detector "hits" to calculate the field integrals via equations 1> or 6> or some variation on this theme.
For example, the PHENIX Central Magnet spectrometer tracks particles, which are outside the magnetic field region. There is no tracking before, or inside, the magnetic field region. In particular, we don't know t(0), the initial direction vector. The experimental data is R(L), R(0), and t(L). Equation 6> tells us that
8> (p(0)c/Ze) (R(L) - R(0) - L t(L)) = (-1) INT{ k(l) B l' dl'}
where k(l) is p(0)/p(l) (the normalized inverse momentum as a function of l which is changing due to the energy loss of the particles in air and in the detectors).
Using simulated tracks, it is therefore a simple matter to generate extensive tables of INT{ k(l) B l' dl'} for all possible paths through the magnet.
Equation 8 is also a very good momentum reconstruction algorithm. In practice, you can solve 8> assuming a high momentum track in the evaluation of INT{B l' dl'} to get a crude estimate of the momentum and then iteratively search the table of integrals in a small window around this estimated momentum to get a better estimate of p(0).
Note that 8> is a vector equations and we really have to match magnitudes and directions. In practice, we will have real data on the left hand side of these equations and "mapped" data on the right; meaning that the equality won't be exact and we have to find the best value of p(0) to minimize the vector residual.
The solution to this minimization problem is left to the reader. However, the answer is given as follows: suppose that you have two vectors A, and B, and a scalar p. p times A is approximately B. The best value for p is p = [ A DOT B ] / [ A DOT A ]. Where A DOT A is the dot product.
Note, also, that the left-hand side of equation 8> assumes we know R(0), the vertex location. If this information is unknown, then we must assume the vertex location because the problem is under-constrained. There are six unknown parameters at the vertex and only four parameters are measured by a detector outside the field region. This comment makes it clear that we can't find the vertex for a single track using a single tracking station. We need several tracking stations inside the magnetic field, multiple tracks to show that pt is conserved, or multiple tracks in combination with high quality time-of-flight information to reconstruct a pair mass. This is important for Lambda decays and particles with significant c-tau decay lengths.
OK, now what the bloody heck does all this have to do with magnet mapping?
Consider the PHENIX Central Magnet. You need to run a good tracking code with a high-resolution field map. (e.g. PISA coupled with Orrin's 1 cm by 1 cm map of the CM, for example. The map is available on the PHENIX WWW homepage under PHENIX/magnet/maps.) This will give you "perfect" data to save. Next, you can bugger the magnetic field (with "bfield_shift.f" available at the same location) and recalculate the data assuming the field map was accidentally shifted a few millimeters by some sleepy graduate student. Do this for the same vector momentum particle in both cases.
Next, put the perfect data on the left-hand side of equation 8>, the buggered map data on the right, and calculate p(0) to see how much it has changed from the "perfect" case. I have finessed the fact that you should re-arrange 8> to solve for p(0) and then take the dot product between "buggered" and "perfect" data to get the best estimate of p(0).
For the Central Magnet the master equation becomes:
Equation CM>
p-new [R(L) - R(0) - L" t(L)]dat DOT [R(L) - R(0) - L" t(L)]map
------ = ---------------------------------------------------------
p-map [R(L) - R(0) - L" t(L)]dat DOT [R(L) - R(0) - L" t(L)]dat
In magnetic field mapper simulations, [ ]map means to use values derived from the buggered map studies. [ ]dat means to use "perfect" data from a simulation without artificial errors. If you want to use this master equation as a momentum reconstruction algorithm then [ ]map means to use values derived from the best available field map and [ ]dat means to use real experimental data (i.e. detector "hits").
However, note that L" is L-map in all cases. Most spectrometers don't measure L (except through time of flight studies) and so you have to use the value estimated from the simulations as input to these equations. This is a subtle, not very big effect, but the correct imitation of the real world.
For the PHENIX Muon Magnet mapper simulation you could use either equations 1> or 3> to derive
Equation MM1>
p-new [ t(L) - t(0) ]dat DOT [ t(L) - t(0) ]map
------ = -----------------------------------------
p-map [ t(L) - t(0) ]dat DOT [ t(L) - t(0) ]dat
or MM2>
p-new [R(L) - R(0) - L" t(0)]dat DOT [R(L) - R(0) - L" t(0)]map
------ = ---------------------------------------------------------
p-map [R(L) - R(0) - L" t(0)]dat DOT [R(L) - R(0) - L" t(0)]dat
Remember that t(0) is V-vector(0)/v-scalar(0) and so it represents the initial direction vector and R(0) is the location of the starting point on tracking station number one in the MM.
These equations give you a quantitative way to understand the effect of systematic errors on the tracking resolution. You can explore the errors introduced by grid sizes in the maps, rotations of the measuring devices, or displacements of the measuring devices. All of this is easier to do than to describe. You should be able to implement these methods with existing codes fairly easily. The hard part is to visualize the data after you get it!
Finally, one last example. For momentum reconstruction in the Muon Magnet, you might want to choose equation 7>, instead, because it uses the hits in all three detector stations simultaneously and eliminates the lower quality direction vector information needed in MM1 and MM2.
Re-arranging 7 and solving for the best estimate of p(0) yields:
Equation MM3>
p-est [ R(3) - (1+L2/L1) R(2) + (L2/L1) R(1) ]dat DOT [...]map
----- = ----------------------------------------------------------
p-map [...]dat DOT [...]dat
Equation MM3> could be used for mapper studies, too, but it is based on an unusual combination of integral data rather than B.dl or B.l.dl, alone. There is nothing wrong with this except that you have to know what you are doing in relation to what you actually measure in the magnet map.
Equation MM3> is a very good algorithm for momentum reconstruction, however, because it over-constrains the solution for the momentum. R(1), R(2), R(3) represent 6 measured parameters and there are only 5 unknowns for the track as it leave station one.