View Javadoc

1   
2   /*
3    * HitUtils.java
4    *
5    * Created on July 2, 2008, 9:31 AM
6    *
7    */
8   
9   package org.lcsim.fit.helicaltrack;
10  
11  import hep.physics.matrix.BasicMatrix;
12  import hep.physics.matrix.Matrix;
13  import hep.physics.matrix.MatrixOp;
14  import hep.physics.matrix.MutableMatrix;
15  import hep.physics.matrix.SymmetricMatrix;
16  import hep.physics.vec.Hep3Vector;
17  import hep.physics.vec.VecOp;
18  
19  import java.util.Map;
20  
21  import org.lcsim.geometry.subdetector.BarrelEndcapFlag;
22  
23  /**
24   *
25   * @author Richard Partridge
26   * @version 1.0
27   */
28  public class HitUtils {
29      private static double _eps = 1.0e-6;
30      private  boolean _debug = false;
31  
32      /** Creates a new instance of HitUtils */
33      public HitUtils() {
34      }
35  
36      /**
37       * Turn a pixel hit into a pseudo-strip hit.  Used for ZSegment hits that
38       * include a pixel hit.  If the pixel is in an endcap disk, we need to
39       * calculate an effective z uncertainty by calling the zres method.
40       * @param hit pixel hit
41       * @param smap map of x-y path lengths
42       * @param msmap map of multiple scatterings
43       * @param helix approximate helix
44       * @return strip hit
45       */
46      public static HelicalTrack2DHit PixelToStrip(HelicalTrackHit hit, Map<HelicalTrackHit, Double> smap,
47              Map<HelicalTrackHit, MultipleScatter> msmap, HelicalTrackFit helix, double tolerance) {
48  
49          //  Take the strip to be span the interval +- tolerance * dz about the pixel hit position
50          double dz = zres(hit, msmap, helix);
51          double zmin = hit.z() - tolerance * dz;
52          double zmax = hit.z() + tolerance * dz;
53  
54          //  Make a new strip hit using the pixel's corrected position and covariance matrix
55          HelicalTrack2DHit striphit = new HelicalTrack2DHit(hit.getCorrectedPosition(),
56                  hit.getCorrectedCovMatrix(), hit.getdEdx(), hit.getTime(), hit.getRawHits(),
57                  hit.Detector(), hit.Layer(), hit.BarrelEndcapFlag(), zmin, zmax);
58          //  Save the path length for the strip hit
59          smap.put(striphit, smap.get(hit));
60  
61          return striphit;
62      }
63  
64      public static Hep3Vector StripCenter(HelicalTrackStrip strip) {
65          return  VecOp.add(strip.origin(), VecOp.mult(strip.umeas(), strip.u()));
66      }
67  
68      public static SymmetricMatrix StripCov(HelicalTrackStrip strip) {
69          SymmetricMatrix cov = new SymmetricMatrix(3);
70          Hep3Vector pos = StripCenter(strip);
71          double x = pos.x();
72          double y = pos.y();
73          double r2 = x*x + y*y;
74          double du = strip.du();
75          cov.setElement(0, 0, y*y * du*du / r2);
76          cov.setElement(0, 1, -x * y * du*du / r2);
77          cov.setElement(1, 1, x*x * du*du / r2);
78          return cov;
79      }
80  
81      public static SymmetricMatrix PixelCov(double x, double y, double drphi, double dz) {
82          SymmetricMatrix cov = new SymmetricMatrix(3);
83          double r2 = x*x + y*y;
84          cov.setElement(0, 0, y*y * drphi*drphi / r2);
85          cov.setElement(0, 1, -x * y * drphi*drphi / r2);
86          cov.setElement(1, 1, x*x * drphi*drphi / r2);
87          cov.setElement(2, 2, dz);
88          return cov;
89      }
90  
91      /**
92       * Find the effective z uncertainty to use in the s-z line fit.  Include
93       * both resolution and multiple scattering contributions.  Endcap disk
94       * hits require converting an uncertainty in r to an effective uncertainty
95       * in z using dz = dr * slope.  If a helix is supplied, it is used to
96       * calculate the slope.  Otherwise, the track is assumed to travel along
97       * a straight line from the origin.
98       * @param hit the hit we want dz for
99       * @param msmap map of the multiple scattering uncertainties
100      * @param helix approximate helix for the track (or null)
101      * @return effective z uncertainty
102      */
103     public static double zres(HelicalTrackHit hit, Map<HelicalTrackHit, MultipleScatter> msmap,
104             HelicalTrackFit helix) {
105 
106         //  Get the multiple scattering uncertainty (if any)
107         double dz_ms = 0.;
108         if (msmap.containsKey(hit)) dz_ms = msmap.get(hit).dz();
109 
110         //  Barrels and disks are treated differently
111         if (hit.BarrelEndcapFlag() == BarrelEndcapFlag.BARREL) {
112             //  For barrel hits, take the resolution uncertainty from the hit covariance matrix
113             double dz_res2 = hit.getCorrectedCovMatrix().diagonal(2);
114             //  Combine resolution and multiple scattering uncertainties in quadrature
115             return Math.sqrt(dz_res2 + dz_ms*dz_ms);
116         } else {
117             //  For endcap disks, take dz = dr * |slope|
118             //  First find the slope - default to the slope for a straight-line track from the origin
119             double slope = hit.z() / hit.r();
120             //  If we have a helix, see if we can use the slope from the helix
121             if (helix != null) {
122                 //  Don't use the helix slope if the magnitude of the slope is smaller than its uncertainty
123                 if (Math.abs(helix.slope()) > helix.getSlopeError()) slope = helix.slope();
124 
125             }
126             //  Take the resolution uncertainty to be dr * |slope|
127             double dzres = hit.dr() * Math.abs(slope);
128             //  Combine resolution and multiple scattering uncertainties in quadrature
129             return Math.sqrt(dzres*dzres + dz_ms*dz_ms);
130         }
131     }
132 
133     /**
134      * Return the hit position assuming the track originated at the origin.
135      * @return hit position
136      */
137     public static Hep3Vector PositionFromOrigin(HelicalTrackStrip strip1, HelicalTrackStrip strip2) {
138 
139         //  Assume the particle is coming from the origin, so x2 = gamma * x1
140         //  gamma = Origin2 . w_hat / Origin1 . w_hat
141         double gamma = VecOp.dot(strip2.origin(), strip2.w()) / NonZeroDotProduct(strip1.origin(), strip1.w());
142 
143         //  Calculate v1hat . u2hat, which is equivalent to sin(alpha) where alpha is the stereo angle
144         double salpha = getSinAlpha(strip1, strip2);
145         //  Calculate the hit locations for v = 0:  p = Origin + u * u_hat
146         Hep3Vector p1 = StripCenter(strip1);
147         Hep3Vector p2 = StripCenter(strip2);
148 
149         //  dp = (p2 - gamma * p1)
150         Hep3Vector dp = VecOp.sub(p2, VecOp.mult(gamma, p1));
151 
152         //  Unmeasured coordinate v1:  v1 = (dp . u2_hat) / (gamma * sin(alpha))
153         double v1 = VecOp.dot(dp, strip2.u()) / (gamma * salpha);
154         if (v1 < strip1.vmin()) v1 = strip1.vmin();
155         if (v1 > strip1.vmax()) v1 = strip1.vmax();
156 
157         //  Position of hit on strip 1:  r1 = p1 + v1 * v1_hat
158         Hep3Vector r1 = VecOp.add(p1, VecOp.mult(v1, strip1.v()));
159         //  Take position to be the midpoint of r1 and r2: r = 0.5 * (1 + gamma) * r1
160      return VecOp.mult(0.5 * (1 + gamma), r1);
161     }
162 
163     /**
164      * Return the covariance matrix assuming the track originated from the
165      * origin.
166      * @return covariance matrix
167      */
168     public static SymmetricMatrix CovarianceFromOrigin(HelicalTrackStrip strip1, HelicalTrackStrip strip2) {
169         //  Assume the particle is coming from the origin, so x2 = gamma * x1
170         //  gamma = Origin2 . w_hat / Origin1 . w_hat
171         double gamma = VecOp.dot(strip2.origin(), strip2.w()) / NonZeroDotProduct(strip1.origin(), strip1.w());
172         //  Calculate the seperation between the sensor planes
173         double separation = SensorSeperation(strip1, strip2);
174         //  Calculate v1hat . u2hat, which is equivalent to sin(alpha) where alpha is the stereo angle
175         double salpha = getSinAlpha(strip1, strip2);
176         //  Calculate the scale factor:  factor = (1 + gamma)^2 / 4 * sin^2(alpha)
177         double factor = (1 + gamma)*(1 + gamma) / (4. * salpha*salpha);
178         //  Calculate v * v^T for strips 1 and 2
179         Matrix v1 = v2m(strip1.v());
180         Matrix v2 = v2m(strip2.v());
181         Matrix v1v1t = MatrixOp.mult(v1, MatrixOp.transposed(v1));
182         Matrix v2v2t = MatrixOp.mult(v2, MatrixOp.transposed(v2));
183         //  Find measurement uncertainties for strips 1 and 2
184         double du1 = strip1.du();
185         double du2 = strip2.du();
186         //  Calculate the uncertainty in the unmeasured coordinate due to not knowing the track direction
187         //  by assuming phat . u has an uncertainty of 2/sqrt(12) so dv = 2 / sqrt(12) * seperation / sin(alpha)
188         double dv = Math.abs(2. * separation / (Math.sqrt(12) * salpha));
189         //  Don't let dv by greater than the strip length / sqrt(12)
190         double dv1 = Math.min(dv, (strip1.vmax()-strip1.vmin()) / Math.sqrt(12.));
191         double dv2 = Math.min(dv, (strip2.vmax()-strip2.vmin()) / Math.sqrt(12.));
192         //  Calculate the covariance matrix.       
193         //    From resolution:  cov = factor * (v2 * v2^T * du1^2 + v1 * v1^T * du2^2)
194         //    From direction:                + (v1 * v1^T * (dv1/2)^2 + v2 * v2^T * (dv2/2)^2
195         Matrix cov1 = MatrixOp.mult(factor * du2*du2 + 0.25 * dv1*dv1, v1v1t);
196         Matrix cov2 = MatrixOp.mult(factor * du1*du1 + 0.25 * dv2*dv2, v2v2t);
197         Matrix cov = MatrixOp.add(cov1, cov2);
198         return new SymmetricMatrix(cov);
199     }
200 
201     /**
202      * Return the hit position given the track direction.
203      * @param trkdir TrackDirection object containing direction and derivatives
204      * @return Corrected hit position
205      */
206     public static Hep3Vector PositionOnHelix(TrackDirection trkdir, HelicalTrackStrip strip1,
207             HelicalTrackStrip strip2) {
208         //  Get the track direction unit vector
209         Hep3Vector dir = trkdir.Direction();
210         double salpha = getSinAlpha(strip1, strip2);
211         //  Gamma is the distance the particle travels between the two sensors:  gamma = separation / (what . dir)
212         double gamma = SensorSeperation(strip1, strip2) / NonZeroDotProduct(strip1.w(), dir);
213         Hep3Vector p1 = StripCenter(strip1);
214         Hep3Vector p2 = StripCenter(strip2);
215         //  dp = p2 - (p1 + gamma * dir)
216         Hep3Vector dp = VecOp.sub(p2, VecOp.add(p1, VecOp.mult(gamma, dir)));
217         //  Unmeasured coordinate v1: v1 = (dp . u2hat) / sin(alpha)
218         double v1 = VecOp.dot(dp, strip2.u()) / salpha;
219         //  Position of hit on strip 1:  r1 = p1 + v1 * v1_hat
220         Hep3Vector r1 = VecOp.add(p1, VecOp.mult(v1, strip1.v()));
221         //  Take position to be the midpoint of x1 and x2: r = r1 + 0.5 * gamma * dir
222         return VecOp.add(r1, VecOp.mult(0.5 * gamma, dir));
223     }
224 
225     
226     /**
227      * Return the covariance matrix given a track direction and helix
228      * covariance matrix.
229      * @param trkdir TrackDirection object containing direction and derivatives
230      * @param hcov covariance matrix for helix parameters
231      * @return corrected covariance matrix
232      */
233     public static SymmetricMatrix CovarianceOnHelix(TrackDirection trkdir, SymmetricMatrix hcov,
234             HelicalTrackStrip strip1, HelicalTrackStrip strip2) {
235         //  Get the track direction and direction derivatives with respect to the helix parameters
236         Hep3Vector dir = trkdir.Direction();
237         Matrix dirderiv = trkdir.Derivatives();
238         //  Calculate phat x v1
239         Matrix pcrossv1 = v2m(VecOp.cross(dir, strip1.v()));
240         //  Calculate phat x v2
241         Matrix pcrossv2 = v2m(VecOp.cross(dir, strip2.v()));
242         //  Calculate phat . w
243         double pdotw = NonZeroDotProduct(dir, strip1.w());
244         //  salpha is the sin of the stereo angle
245         double salpha = getSinAlpha(strip1, strip2);
246         //  Calculate the scale factor:  _separation / 2 * sin(alpha) * (phat . w)^2
247         double factor = SensorSeperation(strip1, strip2) / (2 * salpha * pdotw*pdotw);
248         //  Create the matrix of position derivatives:  d_i,j = dr_i / dphat_j
249         //  The matrix d is given by factor * (v1 * (phat x v2)^T + v2 * (phat x v1)^T  where ^T means transpose
250         Matrix v1 = v2m(strip1.v());
251         Matrix v2 = v2m(strip2.v());
252         Matrix d = MatrixOp.mult(factor, MatrixOp.add(MatrixOp.mult(v1, MatrixOp.transposed(pcrossv2)),
253                 MatrixOp.mult(v2, MatrixOp.transposed(pcrossv1))));
254         Matrix dh = MatrixOp.mult(d, dirderiv);
255         //  Construct the transpose of dh
256         Matrix dht = MatrixOp.transposed(dh);
257         //  Calculate the covariance contributions from the direction uncertainty:  cov = dh * hcov * dh^T
258         Matrix cov_dir = MatrixOp.mult(dh, MatrixOp.mult(hcov, dht));
259         //  Calculate the contributions from measurement errors:  cov += (v2 * v2^T * du1^2 + v1 * v1^T * du2^2) / sin(alpha)^2
260         double du1 = strip1.du();
261         double du2 = strip2.du();
262         Matrix cov1 = MatrixOp.mult(du1*du1 / (salpha*salpha), MatrixOp.mult(v2, MatrixOp.transposed(v2)));
263         Matrix cov2 = MatrixOp.mult(du2*du2 / (salpha*salpha), MatrixOp.mult(v1, MatrixOp.transposed(v1)));
264         //  Sum all the contributions
265         Matrix cov = MatrixOp.add(cov_dir, MatrixOp.add(cov1, cov2));
266 
267         //  Convert to a symmetric matrix
268         return new SymmetricMatrix(cov);
269     }
270 
271     public static double UnmeasuredCoordinate(TrackDirection trkdir, HelicalTrackStrip strip1,
272             HelicalTrackStrip strip2) {
273         //  Get the track direction unit vector
274         Hep3Vector dir = trkdir.Direction();
275         //  Gamma is the distance the particle travels between the two sensors:  gamma = separation / (what . dir)
276         double gamma = SensorSeperation(strip1, strip2) / NonZeroDotProduct(strip1.w(), dir);
277         double salpha = getSinAlpha(strip1, strip2);
278         Hep3Vector p1 = StripCenter(strip1);
279         Hep3Vector p2 = StripCenter(strip2);
280         //  dp = p2 - (p1 + gamma * dir)
281         Hep3Vector dp = VecOp.sub(p2, VecOp.add(p1, VecOp.mult(gamma, dir)));
282         //  Unmeasured coordinate v1: v1 = (dp . u2hat) / sin(alpha)
283         return VecOp.dot(dp, strip2.u()) / salpha;
284     }
285 
286     /**
287      * Calculate the uncertainty in the unmeasured coordinate v1.
288      * @param trkdir track direction
289      * @param hcov helix covariance matrix
290      * @return uncertainty in v1
291      */
292     public static double dv(TrackDirection trkdir, SymmetricMatrix hcov, HelicalTrackStrip strip1,
293             HelicalTrackStrip strip2) {
294         //  Get the track direction and the direction derivatives with respect to the helix parameters
295         Hep3Vector dir = trkdir.Direction();
296         Matrix dirderiv = trkdir.Derivatives();
297         //  Calculate u1 . u2
298         double u1dotu2 = VecOp.dot(strip1.u(), strip2.u());
299         //  Calculate phat . w
300         double pdotw = NonZeroDotProduct(dir, strip1.w());
301         //  Calculate phat x v2
302         Hep3Vector pcrossv2 = VecOp.cross(dir, strip2.v());
303         double salpha = getSinAlpha(strip1, strip2);
304         //  Calculate the scale factor:  separation / sin(alpha) * (phat . w)^2
305         double factor = SensorSeperation(strip1, strip2) / (salpha * pdotw*pdotw);
306         //  The matrix d^T is a row vector given by factor * (phat x v2)
307         MutableMatrix dT = new BasicMatrix(1, 3);
308         for (int i = 0; i < 3; i++) {
309             dT.setElement(0, i, factor * pcrossv2.v()[i]);
310         }
311         //  Construct the matrix dh = d^T * dirderiv
312         Matrix dh = MatrixOp.mult(dT, dirderiv);
313         //  Construct the transpose of dh
314         Matrix dhT = MatrixOp.transposed(dh);
315         //  Calculate the uncertainty in v1 from the direction uncertainty:  cov = dh * hcov * dh^T
316         MutableMatrix cov = (MutableMatrix) MatrixOp.mult(dh, MatrixOp.mult(hcov, dhT));
317         //  Return the uncertainty in v1: dv1^2 = ((u1 . u2)^2 * du1^2 + du2^2) / sin^2(alpha) + cov
318         double dvsq = (Math.pow(u1dotu2 * strip1.du(), 2) + Math.pow(strip2.du(), 2))/ (salpha*salpha);
319         return Math.sqrt(dvsq + cov.e(0, 0));
320     }
321 
322     public static double SensorSeperation(HelicalTrackStrip strip1, HelicalTrackStrip strip2) {
323         //  Calculate the seperation between the sensor planes
324         return VecOp.dot(strip1.w(), VecOp.sub(strip2.origin(), strip1.origin()));
325     }
326 
327     public static double v1Dotu2(HelicalTrackStrip strip1, HelicalTrackStrip strip2) {
328         //  Calculate v1hat . u2hat, which is equivalent to sin(alpha) where alpha is the stereo angle
329         return VecOp.dot(strip1.v(), strip2.u());
330     }
331 
332     public static double getSinAlpha(HelicalTrackStrip strip1, HelicalTrackStrip strip2) {
333         //  Calculate v1hat . u2hat, which is equivalent to sin(alpha) where alpha is the stereo angle
334         double salpha = v1Dotu2(strip1, strip2);
335 
336 //get cos(alpha) and check if the meaurement directions are ~parallel or flipped
337 //mg..5/23/2012...this is wrong for some reason...go back to original
338         //        double calpha = VecOp.dot(strip1.u(), strip2.u());
339 //        if (calpha < 0)
340 //            salpha = -salpha;
341         return salpha;
342     }
343 
344     private static double NonZeroDotProduct(Hep3Vector v1, Hep3Vector v2) {
345         double cth = VecOp.dot(v1, v2);
346         if (Math.abs(cth) < _eps) {
347           if (cth < 0.) cth = -_eps;
348           else cth = _eps;
349         }
350         return cth;
351     }
352 
353     private static Matrix v2m(Hep3Vector v) {
354         BasicMatrix mat = new BasicMatrix(3, 1);
355         mat.setElement(0, 0, v.x());
356         mat.setElement(1, 0, v.y());
357         mat.setElement(2, 0, v.z());
358         return mat;
359     }
360 }