1
2
3
4
5
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
26
27
28 public class HitUtils {
29 private static double _eps = 1.0e-6;
30 private boolean _debug = false;
31
32
33 public HitUtils() {
34 }
35
36
37
38
39
40
41
42
43
44
45
46 public static HelicalTrack2DHit PixelToStrip(HelicalTrackHit hit, Map<HelicalTrackHit, Double> smap,
47 Map<HelicalTrackHit, MultipleScatter> msmap, HelicalTrackFit helix, double tolerance) {
48
49
50 double dz = zres(hit, msmap, helix);
51 double zmin = hit.z() - tolerance * dz;
52 double zmax = hit.z() + tolerance * dz;
53
54
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
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
93
94
95
96
97
98
99
100
101
102
103 public static double zres(HelicalTrackHit hit, Map<HelicalTrackHit, MultipleScatter> msmap,
104 HelicalTrackFit helix) {
105
106
107 double dz_ms = 0.;
108 if (msmap.containsKey(hit)) dz_ms = msmap.get(hit).dz();
109
110
111 if (hit.BarrelEndcapFlag() == BarrelEndcapFlag.BARREL) {
112
113 double dz_res2 = hit.getCorrectedCovMatrix().diagonal(2);
114
115 return Math.sqrt(dz_res2 + dz_ms*dz_ms);
116 } else {
117
118
119 double slope = hit.z() / hit.r();
120
121 if (helix != null) {
122
123 if (Math.abs(helix.slope()) > helix.getSlopeError()) slope = helix.slope();
124
125 }
126
127 double dzres = hit.dr() * Math.abs(slope);
128
129 return Math.sqrt(dzres*dzres + dz_ms*dz_ms);
130 }
131 }
132
133
134
135
136
137 public static Hep3Vector PositionFromOrigin(HelicalTrackStrip strip1, HelicalTrackStrip strip2) {
138
139
140
141 double gamma = VecOp.dot(strip2.origin(), strip2.w()) / NonZeroDotProduct(strip1.origin(), strip1.w());
142
143
144 double salpha = getSinAlpha(strip1, strip2);
145
146 Hep3Vector p1 = StripCenter(strip1);
147 Hep3Vector p2 = StripCenter(strip2);
148
149
150 Hep3Vector dp = VecOp.sub(p2, VecOp.mult(gamma, p1));
151
152
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
158 Hep3Vector r1 = VecOp.add(p1, VecOp.mult(v1, strip1.v()));
159
160 return VecOp.mult(0.5 * (1 + gamma), r1);
161 }
162
163
164
165
166
167
168 public static SymmetricMatrix CovarianceFromOrigin(HelicalTrackStrip strip1, HelicalTrackStrip strip2) {
169
170
171 double gamma = VecOp.dot(strip2.origin(), strip2.w()) / NonZeroDotProduct(strip1.origin(), strip1.w());
172
173 double separation = SensorSeperation(strip1, strip2);
174
175 double salpha = getSinAlpha(strip1, strip2);
176
177 double factor = (1 + gamma)*(1 + gamma) / (4. * salpha*salpha);
178
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
184 double du1 = strip1.du();
185 double du2 = strip2.du();
186
187
188 double dv = Math.abs(2. * separation / (Math.sqrt(12) * salpha));
189
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
193
194
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
203
204
205
206 public static Hep3Vector PositionOnHelix(TrackDirection trkdir, HelicalTrackStrip strip1,
207 HelicalTrackStrip strip2) {
208
209 Hep3Vector dir = trkdir.Direction();
210 double salpha = getSinAlpha(strip1, strip2);
211
212 double gamma = SensorSeperation(strip1, strip2) / NonZeroDotProduct(strip1.w(), dir);
213 Hep3Vector p1 = StripCenter(strip1);
214 Hep3Vector p2 = StripCenter(strip2);
215
216 Hep3Vector dp = VecOp.sub(p2, VecOp.add(p1, VecOp.mult(gamma, dir)));
217
218 double v1 = VecOp.dot(dp, strip2.u()) / salpha;
219
220 Hep3Vector r1 = VecOp.add(p1, VecOp.mult(v1, strip1.v()));
221
222 return VecOp.add(r1, VecOp.mult(0.5 * gamma, dir));
223 }
224
225
226
227
228
229
230
231
232
233 public static SymmetricMatrix CovarianceOnHelix(TrackDirection trkdir, SymmetricMatrix hcov,
234 HelicalTrackStrip strip1, HelicalTrackStrip strip2) {
235
236 Hep3Vector dir = trkdir.Direction();
237 Matrix dirderiv = trkdir.Derivatives();
238
239 Matrix pcrossv1 = v2m(VecOp.cross(dir, strip1.v()));
240
241 Matrix pcrossv2 = v2m(VecOp.cross(dir, strip2.v()));
242
243 double pdotw = NonZeroDotProduct(dir, strip1.w());
244
245 double salpha = getSinAlpha(strip1, strip2);
246
247 double factor = SensorSeperation(strip1, strip2) / (2 * salpha * pdotw*pdotw);
248
249
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
256 Matrix dht = MatrixOp.transposed(dh);
257
258 Matrix cov_dir = MatrixOp.mult(dh, MatrixOp.mult(hcov, dht));
259
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
265 Matrix cov = MatrixOp.add(cov_dir, MatrixOp.add(cov1, cov2));
266
267
268 return new SymmetricMatrix(cov);
269 }
270
271 public static double UnmeasuredCoordinate(TrackDirection trkdir, HelicalTrackStrip strip1,
272 HelicalTrackStrip strip2) {
273
274 Hep3Vector dir = trkdir.Direction();
275
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
281 Hep3Vector dp = VecOp.sub(p2, VecOp.add(p1, VecOp.mult(gamma, dir)));
282
283 return VecOp.dot(dp, strip2.u()) / salpha;
284 }
285
286
287
288
289
290
291
292 public static double dv(TrackDirection trkdir, SymmetricMatrix hcov, HelicalTrackStrip strip1,
293 HelicalTrackStrip strip2) {
294
295 Hep3Vector dir = trkdir.Direction();
296 Matrix dirderiv = trkdir.Derivatives();
297
298 double u1dotu2 = VecOp.dot(strip1.u(), strip2.u());
299
300 double pdotw = NonZeroDotProduct(dir, strip1.w());
301
302 Hep3Vector pcrossv2 = VecOp.cross(dir, strip2.v());
303 double salpha = getSinAlpha(strip1, strip2);
304
305 double factor = SensorSeperation(strip1, strip2) / (salpha * pdotw*pdotw);
306
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
312 Matrix dh = MatrixOp.mult(dT, dirderiv);
313
314 Matrix dhT = MatrixOp.transposed(dh);
315
316 MutableMatrix cov = (MutableMatrix) MatrixOp.mult(dh, MatrixOp.mult(hcov, dhT));
317
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
324 return VecOp.dot(strip1.w(), VecOp.sub(strip2.origin(), strip1.origin()));
325 }
326
327 public static double v1Dotu2(HelicalTrackStrip strip1, HelicalTrackStrip strip2) {
328
329 return VecOp.dot(strip1.v(), strip2.u());
330 }
331
332 public static double getSinAlpha(HelicalTrackStrip strip1, HelicalTrackStrip strip2) {
333
334 double salpha = v1Dotu2(strip1, strip2);
335
336
337
338
339
340
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 }