1 package org.lcsim.recon.tracking.trfzp;
2
3 import org.lcsim.recon.tracking.trfutil.TRFMath;
4 import org.lcsim.recon.tracking.trfutil.Assert;
5
6 import org.lcsim.recon.tracking.trfbase.Propagator;
7 import org.lcsim.recon.tracking.trfbase.PropDirected;
8 import org.lcsim.recon.tracking.trfbase.PropDir;
9 import org.lcsim.recon.tracking.trfbase.PropStat;
10 import org.lcsim.recon.tracking.trfbase.VTrack;
11 import org.lcsim.recon.tracking.trfbase.Surface;
12 import org.lcsim.recon.tracking.trfbase.TrackVector;
13 import org.lcsim.recon.tracking.trfbase.TrackDerivative;
14
15
16
17
18
19
20
21
22
23
24
25
26 public class PropZZ extends PropDirected
27 {
28
29
30
31 private double _bfac;
32
33
34
35 private static final int IX = SurfZPlane.IX;
36 private static final int IY = SurfZPlane.IY;
37 private static final int IDXDZ = SurfZPlane.IDXDZ;
38 private static final int IDYDZ = SurfZPlane.IDYDZ;
39 private static final int IQP = SurfZPlane.IQP;
40
41
42
43
44
45
46
47
48
49
50
51
52 public static String typeName()
53 { return "PropZZ"; }
54
55
56
57
58
59
60
61
62
63 public static String staticType()
64 { return typeName(); }
65
66
67
68
69
70
71
72
73
74
75 public PropZZ(double bfield)
76 {
77 _bfac = TRFMath.BFAC*bfield;
78 }
79
80
81
82
83
84
85
86
87 public Propagator newPropagator( )
88 {
89 return new PropZZ( bField() );
90 }
91
92
93
94
95
96
97
98
99
100
101
102 public PropStat vecDirProp( VTrack trv, Surface srf,
103 PropDir dir)
104 {
105 TrackDerivative pder = null;
106 return vecDirProp(trv, srf, dir, pder);
107
108 }
109
110
111
112
113
114
115
116
117
118
119
120
121 public PropStat vecDirProp( VTrack trv, Surface srf,
122 PropDir dir, TrackDerivative pder )
123 {
124 return vec_propagatezz_( _bfac, trv, srf, dir , pder );
125 }
126
127
128
129
130
131
132
133
134
135 public double bField()
136 {
137 return _bfac/TRFMath.BFAC;
138 }
139
140
141
142
143
144
145
146
147
148 public String type()
149 { return staticType(); }
150
151
152
153
154
155
156
157 public String toString()
158 {
159 return "ZPlane-ZPlane propagation with constant "
160 +bField() +" Tesla field";
161
162 }
163
164
165
166
167
168
169
170
171
172
173
174
175 PropStat
176 vec_propagatezz_( double B, VTrack trv, Surface srf,
177 PropDir dir1,
178 TrackDerivative deriv )
179 {
180
181
182 PropStat pstat = new PropStat();
183 PropDir dir = dir1;
184 boolean move = Propagator.reduceDirection(dir);
185 if(move) dir = reduce(dir);
186
187 Surface srf1 = trv.surface();
188
189
190
191 Assert.assertTrue( srf1.pureType().equals(SurfZPlane.staticType()) );
192 if ( !srf1.pureType( ).equals(SurfZPlane.staticType()) )
193 return pstat;
194
195 SurfZPlane szp1 = (SurfZPlane) srf1;
196
197
198 Assert.assertTrue( srf.pureType().equals(SurfZPlane.staticType()) );
199 if ( !srf.pureType( ).equals(SurfZPlane.staticType()) )
200 return pstat;
201 SurfZPlane szp2 = (SurfZPlane) srf;
202
203
204 if ( srf.pureEqual(srf1) )
205 {
206
207
208 if ( move ) return new PropStat();
209 if(deriv!=null)
210 {
211 deriv.setIdentity();
212 }
213 pstat.setSame();
214 return pstat;
215 }
216
217 if( Math.abs(B) < 1e-7 ) return _zeroBField(trv,szp1,szp2,dir1,deriv);
218
219
220 int izpos = SurfZPlane.ZPOS;
221
222 double zpos_0 = szp1.parameter(izpos);
223 double zpos_n = szp2.parameter(izpos);
224
225
226
227 double dz = zpos_n - zpos_0;
228
229 TrackVector vec = trv.vector();
230 double a1 = vec.get(IX);
231 double a2 = vec.get(IY);
232 double a3 = vec.get(IDXDZ);
233 double a4 = vec.get(IDYDZ);
234 double a5 = vec.get(IQP);
235
236 int sign_dz = 0;
237 if(trv.isForward()) sign_dz = 1;
238 if(trv.isBackward()) sign_dz = -1;
239 if(sign_dz == 0)
240 {
241 System.out.println("PropZZ._vec_propagate: Unknown direction of a track ");
242 System.exit(1);
243 }
244
245
246
247
248 int flag_forward = 0;
249
250 if(dir.equals(PropDir.NEAREST))
251 {
252 if(sign_dz*dz>0.) flag_forward = 1;
253 else flag_forward = -1;
254 }
255
256 else if (dir.equals(PropDir.FORWARD) )
257 {
258
259
260 if(sign_dz*dz<0.) return pstat;
261 flag_forward = 1;
262 }
263 else if (dir.equals(PropDir.BACKWARD) )
264 {
265
266
267 if(sign_dz*dz>0.) return pstat;
268 flag_forward = -1;
269 }
270 else
271 {
272 System.out.println( "PropZZ._vec_propagate: Unknown direction.");
273 System.exit(1);
274 }
275
276 double a34_hat = sign_dz*Math.sqrt(a3*a3 + a4*a4);
277 double a34_hat2 = a34_hat*a34_hat;
278 double dphi = B*dz*a5*Math.sqrt(a34_hat2+1.)*sign_dz;
279 double cos_dphi = Math.cos(dphi);
280 double sin_dphi = Math.sin(dphi);
281
282
283 double Rcos_phi = 1./(B*a5)*a3/Math.sqrt(1.+a34_hat2)*sign_dz;
284 double Rsin_phi = 1./(B*a5)*a4/Math.sqrt(1.+a34_hat2)*sign_dz;
285
286 double a1n = a1 + Rsin_phi*(cos_dphi - 1) + Rcos_phi*sin_dphi;
287 double a2n = a2 - Rcos_phi*(cos_dphi - 1) + Rsin_phi*sin_dphi;
288 double a3n = a3*cos_dphi - a4*sin_dphi;
289 double a4n = a3*sin_dphi + a4*cos_dphi;
290 double a5n = a5;
291
292 vec.set(IX, a1n);
293 vec.set(IY, a2n);
294 vec.set(IDXDZ, a3n);
295 vec.set(IDYDZ, a4n);
296 vec.set(IQP, a5n);
297
298
299 trv.setSurface(srf.newPureSurface());
300 trv.setVector(vec);
301
302
303
304 if(sign_dz == 1) trv.setForward();
305 if(sign_dz == -1) trv.setBackward();
306
307
308 double ctlamsq = a3*a3 + a4*a4;
309 double slam_inv = sign_dz*Math.sqrt(1.0+ctlamsq);
310 double ds = dz*slam_inv;
311
312
313
314 pstat.setPathDistance(ds);
315
316
317
318 if ( flag_forward == 1 ) Assert.assertTrue( pstat.forward() );
319 else Assert.assertTrue( pstat.backward() );
320
321
322 if ( deriv == null ) return pstat;
323
324
325
326 double da34_hat_da3 = 0.;
327 double da34_hat_da4 = 0.;
328 if(Math.abs(a3)>=Math.abs(a4))
329 {
330 int sign3=1;
331 if(a3<0) sign3 = -1;
332 if(a3 !=0.)
333 {
334 da34_hat_da3 = sign_dz*sign3/Math.sqrt(1.+(a4/a3)*(a4/a3) );
335 da34_hat_da4 = sign_dz*(a4/Math.abs(a3))/Math.sqrt(1.+(a4/a3)*(a4/a3) );
336 }
337 else
338 {
339 da34_hat_da3 = 0.;
340 da34_hat_da4 = 0.;
341 }
342 }
343 if(Math.abs(a4)>Math.abs(a3))
344 {
345 int sign4=1;
346 if(a4<0) sign4 = -1;
347 if(a4 !=0.)
348 {
349 da34_hat_da3 = sign_dz*(a3/Math.abs(a4))/Math.sqrt(1.+(a3/a4)*(a3/a4) );
350 da34_hat_da4 = sign_dz*sign4/Math.sqrt(1.+(a3/a4)*(a3/a4) );
351 }
352 else
353 {
354 da34_hat_da3 = 0.;
355 da34_hat_da4 = 0.;
356 }
357 }
358
359
360
361 double ddphi_da3 = B*dz*a5*a34_hat*sign_dz/Math.sqrt(1.+a34_hat2)*da34_hat_da3;
362 double ddphi_da4 = B*dz*a5*a34_hat*sign_dz/Math.sqrt(1.+a34_hat2)*da34_hat_da4;
363 double ddphi_da5 = B*dz*Math.sqrt(a34_hat2+1.)*sign_dz;
364
365
366
367
368 double dRsin_phi_da3 = -Rsin_phi*a34_hat/(1.+a34_hat2)*da34_hat_da3;
369 double dRsin_phi_da4 = -Rsin_phi*a34_hat/(1.+a34_hat2)*da34_hat_da4 +
370 sign_dz/(B*a5)/Math.sqrt(1.+a34_hat2);
371 double dRsin_phi_da5 = -Rsin_phi/a5;
372
373
374
375
376 double dRcos_phi_da3 = -Rcos_phi*a34_hat/(1.+a34_hat2)*da34_hat_da3 +
377 sign_dz/(B*a5)/Math.sqrt(1.+a34_hat2);
378 double dRcos_phi_da4 = -Rcos_phi*a34_hat/(1.+a34_hat2)*da34_hat_da4;
379 double dRcos_phi_da5 = -Rcos_phi/a5;
380
381
382
383 double da1n_da1 = 1.;
384 double da1n_da2 = 0.;
385 double da1n_da3 = dRsin_phi_da3*(cos_dphi-1.) + dRcos_phi_da3*sin_dphi
386 - Rsin_phi*sin_dphi*ddphi_da3 + Rcos_phi*cos_dphi*ddphi_da3;
387 double da1n_da4 = dRsin_phi_da4*(cos_dphi-1.) + dRcos_phi_da4*sin_dphi
388 - Rsin_phi*sin_dphi*ddphi_da4 + Rcos_phi*cos_dphi*ddphi_da4;
389 double da1n_da5 = dRsin_phi_da5*(cos_dphi-1.) + dRcos_phi_da5*sin_dphi
390 - Rsin_phi*sin_dphi*ddphi_da5 + Rcos_phi*cos_dphi*ddphi_da5;
391
392
393
394 double da2n_da1 = 0.;
395 double da2n_da2 = 1.;
396 double da2n_da3 = -dRcos_phi_da3*(cos_dphi-1.) + dRsin_phi_da3*sin_dphi
397 + Rcos_phi*sin_dphi*ddphi_da3 + Rsin_phi*cos_dphi*ddphi_da3;
398 double da2n_da4 = -dRcos_phi_da4*(cos_dphi-1.) + dRsin_phi_da4*sin_dphi
399 + Rcos_phi*sin_dphi*ddphi_da4 + Rsin_phi*cos_dphi*ddphi_da4;
400 double da2n_da5 = -dRcos_phi_da5*(cos_dphi-1.) + dRsin_phi_da5*sin_dphi
401 + Rcos_phi*sin_dphi*ddphi_da5 + Rsin_phi*cos_dphi*ddphi_da5;
402
403
404
405 double da3n_da1 = 0.;
406 double da3n_da2 = 0.;
407 double da3n_da3 = cos_dphi - a3*sin_dphi*ddphi_da3 - a4*cos_dphi*ddphi_da3;
408 double da3n_da4 = - sin_dphi - a3*sin_dphi*ddphi_da4 - a4*cos_dphi*ddphi_da4;
409 double da3n_da5 = - a3*sin_dphi*ddphi_da5 - a4*cos_dphi*ddphi_da5;
410
411
412
413 double da4n_da1 = 0.;
414 double da4n_da2 = 0.;
415 double da4n_da3 = sin_dphi + a3*cos_dphi*ddphi_da3 - a4*sin_dphi*ddphi_da3;
416 double da4n_da4 = cos_dphi + a3*cos_dphi*ddphi_da4 - a4*sin_dphi*ddphi_da4;
417 double da4n_da5 = a3*cos_dphi*ddphi_da5 - a4*sin_dphi*ddphi_da5;
418
419
420
421 double da5n_da1 = 0.;
422 double da5n_da2 = 0.;
423 double da5n_da3 = 0.;
424 double da5n_da4 = 0.;
425 double da5n_da5 = 1.;
426
427 deriv.set(IX,IX ,da1n_da1);
428 deriv.set(IX,IY ,da1n_da2);
429 deriv.set(IX,IDXDZ ,da1n_da3);
430 deriv.set(IX,IDYDZ ,da1n_da4);
431 deriv.set(IX,IQP ,da1n_da5);
432 deriv.set(IY,IX ,da2n_da1);
433 deriv.set(IY,IY ,da2n_da2);
434 deriv.set(IY,IDXDZ ,da2n_da3);
435 deriv.set(IY,IDYDZ ,da2n_da4);
436 deriv.set(IY,IQP ,da2n_da5);
437 deriv.set(IDXDZ,IX ,da3n_da1);
438 deriv.set(IDXDZ,IY ,da3n_da2);
439 deriv.set(IDXDZ,IDXDZ ,da3n_da3);
440 deriv.set(IDXDZ,IDYDZ ,da3n_da4);
441 deriv.set(IDXDZ,IQP ,da3n_da5);
442 deriv.set(IDYDZ,IX ,da4n_da1);
443 deriv.set(IDYDZ,IY ,da4n_da2);
444 deriv.set(IDYDZ,IDXDZ ,da4n_da3);
445 deriv.set(IDYDZ,IDYDZ ,da4n_da4);
446 deriv.set(IDYDZ,IQP ,da4n_da5);
447 deriv.set(IQP,IX ,da5n_da1);
448 deriv.set(IQP,IY ,da5n_da2);
449 deriv.set(IQP,IDXDZ ,da5n_da3);
450 deriv.set(IQP,IDYDZ ,da5n_da4);
451 deriv.set(IQP,IQP ,da5n_da5);
452
453 return pstat;
454
455 }
456
457 private PropStat _zeroBField( VTrack trv, SurfZPlane szp1,
458 SurfZPlane szp2,PropDir dir1,
459 TrackDerivative deriv)
460 {
461 PropStat pstat = new PropStat();
462 PropDir dir = dir1;
463 boolean move = Propagator.reduceDirection(dir);
464 boolean same = szp2.pureEqual(szp1);
465
466
467 if ( same && move ) return pstat;
468
469 if ( same )
470 {
471 if(deriv != null)
472 {
473
474 deriv.setIdentity();
475 }
476 pstat.setSame();
477 return pstat;
478 }
479
480 TrackVector vec = trv.vector();
481 double x0 = vec.get(IX);
482 double y0 = vec.get(IY);
483 double dxdz0 = vec.get(IDXDZ);
484 double dydz0 = vec.get(IDYDZ);
485
486 double dz0 =1.;
487 if( trv.isBackward() ) dz0 = -1.;
488
489 double z0 = szp1.parameter(SurfZPlane.ZPOS);
490
491 double z1 = szp2.parameter(SurfZPlane.ZPOS);
492
493 double a = dxdz0*dz0;
494 double b = dydz0*dz0;
495 double c = dz0;
496
497 double ap = 0.;
498 double bp = 0.;
499 double cp = 1.0;
500
501 double xp = 0.;
502 double yp = 0.;
503 double zp = z1;
504
505 double denom = a*ap + b*bp + c*cp;
506
507 if( denom == 0. ) return pstat;
508
509 double S = ( (xp-x0)*ap + (yp-y0)*bp + (zp-z0)*cp )/denom;
510
511 double x1 = x0 + S*a;
512 double y1 = y0 + S*b;
513
514 boolean forward = S > 0. ? true : false;
515
516 if( dir == PropDir.FORWARD && !forward ) return pstat;
517 if( dir == PropDir.BACKWARD && forward ) return pstat;
518
519
520 double dxdz1 = (x1-x0)/(z1-z0);
521 double dydz1 = (y1-y0)/(z1-z0);
522
523 vec.set(IX, x1);
524 vec.set(IY, y1);
525 vec.set(IDXDZ, dxdz1);
526 vec.set(IDYDZ, dydz1);
527
528 trv.setSurface( szp2.newPureSurface());
529 trv.setVector(vec);
530 if( dz0 >0 ) trv.setForward();
531 else trv.setBackward();
532
533 double ds = S*Math.sqrt(a*a+b*b+c*c);
534
535 pstat.setPathDistance(ds);
536
537 if( deriv == null ) return pstat;
538
539
540
541
542
543
544 double da_dxdz0 = dz0;
545 double db_dydz0 = dz0;
546
547 double dx1dx0 =1.;
548 double dy1dy0 =1.;
549
550 double dx1_dxdz0 = S*da_dxdz0;
551 double dy1_dydz0 = S*db_dydz0;
552
553 double dxdz1_dx0 = (dx1dx0 - 1.)/(z1-z0);
554 double dydz1_dy0 = (dy1dy0 - 1.)/(z1-z0);
555
556 double dxdz1_dxdz0 = dx1_dxdz0/(z1-z0);
557 double dydz1_dydz0 = dy1_dydz0/(z1-z0);
558
559
560 deriv.setIdentity();
561
562 deriv.set(IX,IX, dx1dx0);
563 deriv.set(IX,IDXDZ, dx1_dxdz0);
564
565 deriv.set(IY,IY, dy1dy0);
566 deriv.set(IY,IDYDZ, dy1_dydz0);
567
568 deriv.set(IDXDZ,IX, dxdz1_dx0);
569 deriv.set(IDXDZ,IDXDZ, dxdz1_dxdz0);
570
571 deriv.set(IDYDZ,IY, dydz1_dy0);
572 deriv.set(IDYDZ,IDYDZ, dydz1_dydz0);
573
574 deriv.set(IQP,IQP, 1.0);
575
576 return pstat;
577
578 }
579
580
581
582 }