This program finds the shower-detector plane (SDP). It works by identifying the great circle on the eye's celestial sphere which best contains the center directions of the pixels which received the light flux from the shower axis. The great circle is designated by the "plane normal vector" pln[3], which is a unit vector. A great circle on a sphere is the intersection of the sphere with a plane that contains the sphere's center. A great circle on the unit sphere of directions is conveniently represented by a unit vector perpendicular to its defining plane. The output of PlaneFit is therefore a unit "plane normal vector" denoted as pln[] in this module. The unit vector in the opposite direction, -pln[], defines the same great circle. We choose the sense of pln[] in accordance with (first)x(last), i.e. the cross product of the first pixel direction with the last pixel direction (ordered by time) should be (approximately) parallel (not anti-parallel) to pln[].
In the present version, the function we minimize is the "moment of inertia" about the shower track. If you think of each pixel's charge integral as a point mass m located at its angular position, then this moment of inertia about a trial track is the sum of terms m*r^2, where r is the shortest angular distance from the pixel center to the that track. The best track is the one which minimizes this moment-of-inertia sum. In practice, we define r = dot(pix[i].dir[],pln[]), where pix[i].dir[] is the direction unit vector of pixel number i, pln[] is the trial plane normal vector, and r is their scalar dot product. The charge integral pix[i].pint is used for each "mass" m in this moment-of-inertia function" sum{m*r^2}.
The subroutine PlaneGuess is used to start amoeba with one of its simplex vertices located close to the correct answer. The guess is obtained by examining n/2 cross products of pixel directions, where the first pixel direction in each product is taken from the pixels in the first half of the track (ordered by time) and the second direction in the cross product is that of a pixel in the second half. These candidate unit normal vectors are then tested against each other. The one that is "most like all the others" is taken as the initial plane normal guess. Explicitly, it is the candidate unit vector which gives the greatest sum of dot products with all of the candidates.
The set of possible plane normal vectors is a two-dimensional space. It could be parametrized by spherical coordinates (theta,phi), for example. We have adopted a different parametrization, however, in which the uncertainties in the two parameters have clear geometric meaning. One parameter (denoted by a) relates to shifts of the track parallel to itself near the centroid of the pixels. The other parameter (denoted by b) relates to variations of the track which keep the track centroid fixed. These tilts are rotations of the track about that centroid point. In general, the pixel directions fix the a-parameter more precisely than the b-parameter, so uncertainties in b are much greater than uncertainties in a. These two parameters (a and b) are defined with respect to the guessed plane normal vector. Their best-fit values are not meaningful without reference to that particular vector. Their uncertainties, however, are not sensitive to small changes in the choice of that starting vector. To be specific, let "cen[]" be the centroid vector obtained by averaging the pixel directions and normalizing the resultant to a unit vector. Let "pln[]" be the unit normal guess. Define "top[]" by the cross product (cen)x(pln). Any trial normal vector is then given by a pair of numbers (a,b) according to pln[]+a*cen[]+b*top[] followed by renormalization to a unit vector.
The chisquare-type function that is presently minimized to
determine the optimal parameters (a,b) is the moment of inertia
described above. A contour plot of this function typically
reveals a long valley extending in the b-direction, which
indicates that the pixel directions constrain sideways shifts of
the track better than tilts that leave the centroid fixed. The
long valley can have multiple minima, and the absolute minimum
need not occur at the (a,b) values of the true plane normal
vector. In the example plot below, the true (Monte Carlo) normal
vector is used for the guess, so the correct normal vector is at
(0,0). The fitted value is off in the a-parameter (shift) by
less than 0.01 radian, but the chisquare has minima that are
approximately 0.08 radian away from the correct
value in the b-parameter (tilt). The chisquare function can be
seen to have multiple local minima, and the absolute minimum
yields a best-fit plane normal vector that differs from the true
Monte-Carlo normal vector by 4.5 degrees. The track in this
example is only 10 degrees long and is not typical. The plane
normal can generally be determined more accurately than this.
The PlaneFit C code has the
following properties:
Input: