1 package org.lcsim.recon.tracking.trfcylplane;
2
3 import org.lcsim.recon.tracking.trfbase.Propagator;
4 import org.lcsim.recon.tracking.trfbase.PropDirected;
5 import org.lcsim.recon.tracking.trfbase.PropDir;
6 import org.lcsim.recon.tracking.trfbase.TrackVector;
7 import org.lcsim.recon.tracking.trfbase.TrackDerivative;
8 import org.lcsim.recon.tracking.trfbase.PropStat;
9 import org.lcsim.recon.tracking.trfutil.TRFMath;
10 import org.lcsim.recon.tracking.trfutil.Assert;
11 import org.lcsim.recon.tracking.trfcyl.SurfCylinder;
12 import org.lcsim.recon.tracking.trfxyp.SurfXYPlane;
13 import org.lcsim.recon.tracking.trfbase.Surface;
14 import org.lcsim.recon.tracking.trfbase.VTrack;
15 import org.lcsim.recon.tracking.trfbase.ETrack;
16
17 /**
18 * Propagates tracks from an XYPlane to a Cylinder in a constant field.
19 *<p>
20 * Propagation will fail if either the origin is not an XYPlane
21 * or destination is not a Cylinder.
22 * Propagator works incorrectly for tracks with very small curvatures.
23 *<p>
24 *@author Norman A. Graf
25 *@version 1.0
26 *
27 */
28 public class PropXYCyl extends PropDirected
29 {
30
31 // Assign track parameter indices.
32
33 private static final int IV = SurfXYPlane.IV;
34 private static final int IZC = SurfXYPlane.IZ;
35 private static final int IDVDU = SurfXYPlane.IDVDU;
36 private static final int IDZDU = SurfXYPlane.IDZDU;
37 private static final int IQP_XY = SurfXYPlane.IQP;
38
39 private static final int IPHI = SurfCylinder.IPHI;
40 private static final int IZ = SurfCylinder.IZ;
41 private static final int IALF = SurfCylinder.IALF;
42 private static final int ITLM = SurfCylinder.ITLM;
43 private static final int IQPT = SurfCylinder.IQPT;
44
45
46
47 // attributes
48
49 // BFAC * bfield
50 private double _bfac;
51
52 // static methods
53
54 /**
55 *Return a String representation of the class' type name.
56 *Included for completeness with the C++ version.
57 *
58 * @return A String representation of the class' type name.
59 */
60 public static String typeName()
61 { return "PropXYCyl";
62 }
63
64 /**
65 *Return a String representation of the class' type name.
66 *Included for completeness with the C++ version.
67 *
68 * @return A String representation of the class' type name.
69 */
70 public static String staticType()
71 { return typeName();
72 }
73
74
75 /**
76 *Construct an instance from a constant solenoidal magnetic field in Tesla.
77 *
78 * @param bfield The magnetic field strength in Tesla.
79 */
80 public PropXYCyl(double bfield)
81 {
82 _bfac = TRFMath.BFAC*bfield;
83 }
84
85 /**
86 *Clone an instance.
87 *
88 * @return A Clone of this instance.
89 */
90 public Propagator newPropagator( )
91 {
92 return new PropXYCyl( bField() );
93 }
94
95 /**
96 *Propagate a track without error in the specified direction.
97 *
98 * The track parameters for a cylinder are:
99 * phi z alpha tan(lambda) curvature
100 *
101 * @param trv The VTrack to propagate.
102 * @param srf The Surface to which to propagate.
103 * @param dir The direction in which to propagate.
104 * @return The propagation status.
105 */
106 public PropStat vecDirProp(VTrack trv, Surface srf,
107 PropDir dir )
108 {
109 TrackDerivative deriv = null;
110 return vecDirProp(trv, srf, dir, deriv);
111 }
112
113 /**
114 *Propagate a track without error in the specified direction
115 *and return the derivative matrix in deriv.
116 *
117 * The track parameters for a cylinder are:
118 * phi z alpha tan(lambda) curvature
119 *
120 * @param trv The VTrack to propagate.
121 * @param srf The Surface to which to propagate.
122 * @param dir The direction in which to propagate.
123 * @param deriv The track derivatives to update at the surface srf.
124 * @return The propagation status.
125 */
126 public PropStat vecDirProp(VTrack trv, Surface srf,
127 PropDir dir, TrackDerivative deriv )
128 {
129 return vecPropagateXYCyl( _bfac, trv, srf, dir, deriv);
130
131 }
132
133 /**
134 *Return the strength of the magnetic field in Tesla.
135 *
136 * @return The strength of the magnetic field in Tesla.
137 */
138 public double bField()
139 {
140 return _bfac/TRFMath.BFAC;
141 }
142
143 /**
144 *Return a String representation of the class' type name.
145 *Included for completeness with the C++ version.
146 *
147 * @return A String representation of the class' type name.
148 */
149 public String type()
150 { return staticType();
151 }
152
153 /**
154 *output stream
155 * @return A String representation of this instance.
156 */
157 public String toString()
158 {
159 return "XYPlane-Cylinder propagation with constant "
160 + bField() +" Tesla field";
161
162 }
163
164
165 //**********************************************************************
166
167 // Private function to propagate a track without error
168 //
169 // The track parameters for a cylinder are:
170 // phi z alpha tan(lambda) curvature
171 // (NOT [. . . lambda .] as in TRF and earlier version of TRF++.)
172 //
173 // If pderiv is nonzero, return the derivative matrix there.
174 // On Cylinder:
175 // r (cm) is fixed
176 // 0 - phi
177 // 1 - z (cm)
178 // 2 - alp
179 // 3 - tlam
180 // 4 - q/pt pt is transverse momentum of a track, q is its charge
181 // On XYPlane:
182 // u (cm) is fixed
183 // 0 - v (cm)
184 // 1 - z (cm)
185 // 2 - dv/du
186 // 3 - dz/du
187 // 4 - q/p p is momentum of a track, q is its charge
188
189 PropStat
190 vecPropagateXYCyl( double bfac, VTrack trv, Surface srf,
191 PropDir dir,
192 TrackDerivative deriv )
193 {
194
195 // construct return status
196 PropStat pstat = new PropStat();
197
198 // fetch the originating surface and vector
199 Surface srf1 = trv.surface();
200 TrackVector vec1 = trv.vector();
201
202 // Check origin is a XYPlane.
203 Assert.assertTrue( srf1.pureType().equals(SurfXYPlane.staticType()) );
204 if ( !srf1.pureType( ).equals(SurfXYPlane.staticType()) )
205 return pstat;
206 SurfXYPlane sxyp1 = ( SurfXYPlane) srf1;
207
208 // Check destination is a cylinder.
209 Assert.assertTrue( srf.pureType().equals(SurfCylinder.staticType()) );
210 if ( !srf.pureType( ).equals(SurfCylinder.staticType()) )
211 return pstat;
212 SurfCylinder scy2 = ( SurfCylinder) srf;
213
214
215 // Fetch the dist and phi of the plane and the starting track vector.
216 int iphi = SurfXYPlane.NORMPHI;
217 int idist = SurfXYPlane.DISTNORM;
218 double phi = sxyp1.parameter(iphi);
219 double r = sxyp1.parameter(idist);
220
221 TrackVector vec = trv.vector();
222 double v = vec.get(IV); // v
223 double z = vec.get(IZC); // z
224 double b = vec.get(IDVDU); // dv/du
225 double a = vec.get(IDZDU); // dz/du
226 double e = vec.get(IQP_XY); // q/p
227
228 // Fetch the radii and the starting track vector.
229 int irad = SurfCylinder.RADIUS;
230 double r2 = scy2.parameter(irad);
231
232 int sign_du = 0;
233 if(trv.isForward()) sign_du = 1;
234 if(trv.isBackward()) sign_du = -1;
235 if(sign_du == 0)
236 {
237 System.out.println("PropXYCyl._vec_propagate: Unknown direction of a track ");
238 return pstat;
239 }
240
241 // Calculate cylindrical coordinates
242
243 double cnst=sign_du>0?0.: Math.PI;
244 double atn=Math.atan(v/r);
245 Assert.assertTrue(r!=0.);
246
247 double phi1= TRFMath.fmod2(phi+atn,TRFMath.TWOPI);
248 double r1=Math.sqrt(r*r+v*v);
249 double z1=z;
250 double alp1= TRFMath.fmod2(Math.atan(b)-atn+cnst,TRFMath.TWOPI);
251 double tlm1= a*sign_du/Math.sqrt(1+b*b);
252 double qpt1=e*Math.sqrt((1+a*a+b*b)/(1+b*b));
253
254
255 // Check alpha range.
256 alp1 = TRFMath.fmod2( alp1, TRFMath.TWOPI );
257 Assert.assertTrue( Math.abs(alp1) <= Math.PI );
258 //if ( trv.is_forward() ) Assert.assertTrue( Math.abs(alp1) <= TRFMath.PI2 );
259 //else Assert.assertTrue( Math.abs(alp1) > TRFMath.PI2 );
260
261 // Calculate the cosine of lambda.
262 //double clam1 = 1.0/Math.sqrt(1+tlm1*tlm1);
263
264 // Calculate curvature: C = _bfac*(q/p)/Math.cos(lambda)
265 // and its derivatives
266 // Assert.assertTrue( clam1 != 0.0 );
267 // double dcrv1_dqp1 = bfac/clam1;
268 // double crv1 = dcrv1_dqp1*qp1;
269 // double dcrv1_dtlm1 = crv1*clam1*clam1*tlm1;
270
271 // Calculate the curvature = _bfac*(q/pt)
272 double dcrv1_dqpt1 = bfac;
273 double crv1 = dcrv1_dqpt1*qpt1;
274 //double dcrv1_dtlm1 = 0.0;
275
276 // Evaluate the new track vector.
277 // See dla log I-044
278
279 // lambda and curvature do not change
280 double tlm2 = tlm1;
281 double crv2 = crv1;
282 double qpt2 = qpt1;
283
284 // We can evaluate Math.sin(alp2), leaving two possibilities for alp2
285 // 1st solution: alp21, phi21, phid21, tht21
286 // 2nd solution: alp22, phi22, phid22, tht22
287 // evaluate phi2 to choose
288 double salp1 = Math.sin( alp1 );
289 double calp1 = Math.cos( alp1 );
290 double salp2 = r1/r2*salp1 + 0.5*crv1/r2*(r2*r2-r1*r1);
291 // if salp2 > 1, track does not cross cylinder
292 if ( Math.abs(salp2) > 1.0 ) return pstat;
293 double alp21 = Math.asin( salp2 );
294 double alp22 = alp21>0 ? Math.PI-alp21 : -Math.PI-alp21;
295 double calp21 = Math.cos( alp21 );
296 double calp22 = Math.cos( alp22 );
297 double phi20 = phi1 + Math.atan2( salp1-r1*crv1, calp1 );
298 double phi21 = phi20 - Math.atan2( salp2-r2*crv2, calp21 ); // phi position
299 double phi22 = phi20 - Math.atan2( salp2-r2*crv2, calp22 );
300
301 // Construct an sT object for each solution.
302 STCalcXY sto1 = new STCalcXY(r1,phi1,alp1,crv1,r2,phi21,alp21);
303 STCalcXY sto2 = new STCalcXY(r1,phi1,alp1,crv1,r2,phi22,alp22);
304 // Check the two solutions are nonzero and have opposite sign
305 // or at least one is nonzero.
306
307 // Choose the correct solution
308 boolean use_first_solution = false;
309
310 if (dir.equals(PropDir.NEAREST))
311 {
312 use_first_solution = Math.abs(sto2.st()) > Math.abs(sto1.st());
313 }
314 else if (dir.equals(PropDir.FORWARD))
315 {
316 use_first_solution = sto1.st() > 0.0;
317 }
318 else if (dir.equals(PropDir.BACKWARD))
319 {
320 use_first_solution = sto1.st() < 0.0;
321 }
322 else
323 {
324 System.out.println("PropCyl._vec_propagate: Unknown direction.");
325 System.exit(1);
326 }
327
328 // Assign phi2, alp2 and sto2 for the chosen solution.
329 double phi2, alp2;
330 STCalcXY sto;
331 double calp2;
332 if ( use_first_solution )
333 {
334 sto = sto1;
335 phi2 = phi21;
336 alp2 = alp21;
337 calp2 = calp21;
338 }
339 else
340 {
341 sto = sto2;
342 phi2 = phi22;
343 alp2 = alp22;
344 calp2 = calp22;
345 }
346
347 // fetch sT.
348 double st = sto.st();
349
350 // use sT to evaluate z2
351 double z2 = z1 + tlm1*st;
352
353 // Check alpha range.
354 Assert.assertTrue( Math.abs(alp2) <= Math.PI );
355
356 // put new values in vec
357 vec.set(IPHI , phi2);
358 vec.set(IZ , z2);
359 vec.set(IALF , alp2);
360 vec.set(ITLM , tlm2);
361 vec.set(IQPT , qpt2);
362
363 // Update trv
364 trv.setSurface(srf.newPureSurface());
365 trv.setVector(vec);
366 if ( Math.abs(alp2) <= TRFMath.PI2 ) trv.setForward();
367 else trv.setBackward();
368
369 // Set the return status.
370 double s = st*Math.sqrt(1.0+tlm2*tlm2);
371 pstat.setPathDistance(s);
372 //st > 0 ? pstat.set_forward() : pstat.set_backward();
373
374 // exit now if user did not ask for error matrix.
375 if ( deriv == null ) return pstat;
376
377 // Calculate derivatives.
378 // dphi1
379
380 double dphi1_dv= r/(r*r+v*v);
381
382 // dz1
383
384 double dz1_dz=1.;
385
386 // dalf1
387
388 double dalp1_db= 1./(1.+b*b);
389 double dalp1_dv= -r/(r*r+v*v);
390
391 // dr1
392
393 double dr1_dv= v/Math.sqrt(v*v+r*r);
394
395 // dtlm1
396
397 double dtlm1_da= sign_du/Math.sqrt(1+b*b);
398 double dtlm1_db= -a*sign_du*b/(1+b*b)/Math.sqrt(1+b*b);
399
400 // dcrv1
401
402 double dqpt1_de= Math.sqrt((1+a*a+b*b)/(1+b*b));
403 double dqpt1_da= a*e/Math.sqrt((1+b*b)*(1+a*a+b*b));
404 double dqpt1_db= -e*b*a*a/Math.sqrt(1+a*a+b*b)/Math.sqrt(1+b*b)/(1+b*b);
405
406 double dcrv1_de= dqpt1_de*bfac;
407 double dcrv1_da= dqpt1_da*bfac;
408 double dcrv1_db= dqpt1_db*bfac;
409
410 // alpha_2
411 double da2da1 = r1*calp1/r2/calp2;
412 double da2dc1 = (r2*r2-r1*r1)*0.5/r2/calp2;
413 double da2dr1 = (salp1-crv2*r1)/r2/calp2;
414
415 // phi2
416 double rcsal1 = r1*crv1*salp1;
417 double rcsal2 = r2*crv2*salp2;
418 double den1 = 1.0 + r1*r1*crv1*crv1 - 2.0*rcsal1;
419 double den2 = 1.0 + r2*r2*crv2*crv2 - 2.0*rcsal2;
420 double dp2dp1 = 1.0;
421 double dp2da1 = (1.0-rcsal1)/den1 - (1.0-rcsal2)/den2*da2da1;
422 double dp2dc1 = -r1*calp1/den1 + r2*calp2/den2
423 - (1.0-rcsal2)/den2*da2dc1;
424 double dp2dr1= -crv1*calp1/den1-(1.0-rcsal2)*da2dr1/den2;
425
426 // z2
427 double dz2dz1 = 1.0;
428 double dz2dl1 = st;
429 double dz2da1 = tlm1*sto.d_st_dalp1(dp2da1, da2da1);
430 double dz2dc1 = tlm1*sto.d_st_dcrv1(dp2dc1, da2dc1);
431 double dz2dr1 = tlm1*sto.d_st_dr1( dp2dr1, da2dr1);
432
433 // final derivatives
434
435 // phi2
436 double dphi2_dv=dp2dp1*dphi1_dv+dp2da1*dalp1_dv+dp2dr1*dr1_dv;
437 double dphi2_db=dp2da1*dalp1_db+dp2dc1*dcrv1_db;
438 double dphi2_da=dp2dc1*dcrv1_da;
439 double dphi2_de=dp2dc1*dcrv1_de;
440
441 // alp2
442 double dalp2_dv= da2da1*dalp1_dv+da2dr1*dr1_dv;
443 double dalp2_db= da2da1*dalp1_db+da2dc1*dcrv1_db;
444 double dalp2_da= da2dc1*dcrv1_da;
445 double dalp2_de= da2dc1*dcrv1_de;
446
447 // crv2
448 double dqpt2_da=dqpt1_da;
449 double dqpt2_db=dqpt1_db;
450 double dqpt2_de=dqpt1_de;
451
452
453 // tlm2
454 double dtlm2_da= dtlm1_da;
455 double dtlm2_db= dtlm1_db;
456
457 // z2
458 double dz2_dz= dz2dz1*dz1_dz;
459 double dz2_dv= dz2dr1*dr1_dv+dz2da1*dalp1_dv;
460 double dz2_db= dz2da1*dalp1_db+dz2dl1*dtlm1_db+dz2dc1*dcrv1_db;
461 double dz2_da= dz2dl1*dtlm1_da+dz2dc1*dcrv1_da;
462 double dz2_de= dz2dc1*dcrv1_de;
463
464
465 // Build derivative matrix.
466
467 deriv.set(IPHI,IV , dphi2_dv);
468 deriv.set(IPHI,IDVDU, dphi2_db);
469 deriv.set(IPHI,IDZDU, dphi2_da);
470 deriv.set(IPHI,IQP_XY, dphi2_de);
471 deriv.set(IZ,IV , dz2_dv);
472 deriv.set(IZ,IZC , dz2_dz);
473 deriv.set(IZ,IDVDU , dz2_db);
474 deriv.set(IZ,IDZDU , dz2_da);
475 deriv.set(IZ,IQP_XY, dz2_de);
476 deriv.set(IALF,IV , dalp2_dv);
477 deriv.set(IALF,IDVDU , dalp2_db);
478 deriv.set(IALF,IDZDU , dalp2_da);
479 deriv.set(IALF,IQP_XY , dalp2_de);
480 deriv.set(ITLM,IDVDU , dtlm2_db);
481 deriv.set(ITLM,IDZDU , dtlm2_da);
482 deriv.set(IQPT,IDVDU , dqpt2_db);
483 deriv.set(IQPT,IDZDU , dqpt2_da);
484 deriv.set(IQPT,IQP_XY , dqpt2_de);
485
486 return pstat;
487
488 }
489
490 //**********************************************************************
491 // helpers
492 //**********************************************************************
493
494 // Private class STCalcXY.
495 //
496 // An STCalcXY_ object calculates sT (the signed transverse path length)
497 // and its derivatives w.r.t. alf1 and crv1. It is constructed from
498 // the starting (r1, phi1, alf1, crv1) and final track parameters
499 // (r2, phi2, alf2) assuming these are consistent. Methods are
500 // provided to retrieve sT and the two derivatives.
501
502 class STCalcXY
503 {
504
505 private boolean _big_crv;
506 private double _st;
507 private double _dst_dphi21;
508 private double _dst_dcrv1;
509 private double _dst_dr1;
510 private double _cnst1,_cnst2;
511 public double _crv1;
512
513 // constructor
514 public STCalcXY()
515 {
516 }
517 public STCalcXY(double r1, double phi1, double alf1, double crv1,
518 double r2, double phi2, double alf2)
519 {
520
521 _crv1 = crv1;
522 Assert.assertTrue( r1 > 0.0 );
523 Assert.assertTrue( r2 > 0.0 );
524 double rmax = r1+r2;
525
526 // Calculate the change in xy direction
527 double phi_dir_diff = TRFMath.fmod2(phi2+alf2-phi1-alf1,TRFMath.TWOPI);
528 Assert.assertTrue( Math.abs(phi_dir_diff) <= Math.PI );
529
530 // Evaluate whether |C| is" big"
531 _big_crv = rmax*Math.abs(crv1) > 0.001 || Math.abs(r1-r2)<1.e-9;
532 if( Math.abs(crv1) < 1.e-10 ) _big_crv=false;
533
534 // If the curvature is big we can use
535 // sT = (phi_dir2 - phi_dir1)/crv1
536 if ( _big_crv )
537 {
538 Assert.assertTrue( crv1 != 0.0 );
539 _st = phi_dir_diff/crv1;
540 }
541
542 // Otherwise, we calculate the straight-line distance
543 // between the points and use an approximate correction
544 // for the (small) curvature.
545 else
546 {
547
548 // evaluate the distance
549 double d = Math.sqrt( r1*r1 + r2*r2 - 2.0*r1*r2*Math.cos(phi2-phi1) );
550 double arg = 0.5*d*crv1;
551 double arg2 = arg*arg;
552 double st_minus_d = d*arg2*( 1.0/6.0 + 3.0/40.0*arg2 );
553 _st = d + st_minus_d;
554
555 // evaluate the sign
556 // We define a metric xsign = abs( (dphid-d*C)/(d*C) ).
557 // Because sT*C = dphid and d = abs(sT):
558 // xsign = 0 for sT > 0
559 // xsign = 2 for sT < 0
560 // Numerical roundoff will smear these predictions.
561 double xsign = (crv1 == 0. ? 0.: Math.abs( (phi_dir_diff - _st*crv1) / (_st*crv1)) );
562 double sign = 0.0;
563 if ( crv1 != 0. )
564 {
565 if ( xsign < 0.5 ) sign = 1.0;
566 if ( xsign > 1.5 && xsign < 3.0 ) sign = -1.0;
567 }
568 // If the above is indeterminate, assume zero curvature.
569 // In this case abs(alpha) decreases monotonically
570 // with sT. Track passing through origin has alpha = 0 on one
571 // side and alpha = +/-pi on the other. If both points are on
572 // the same side, we use dr/ds > 0 for |alpha|<pi/2.
573 if ( sign == 0. )
574 {
575 sign = 1.0;
576 if ( Math.abs(alf2) > Math.abs(alf1) ) sign = -1.0;
577 if ( Math.abs(alf2) == Math.abs(alf1) )
578 {
579 if ( Math.abs(alf2) < TRFMath.PI2 )
580 {
581 if ( r2 < r1 ) sign = -1.0;
582 }
583 else
584 {
585 if ( r2 > r1 ) sign = -1.0;
586 }
587 }
588 }
589
590 // Correct _st using the above sign.
591 Assert.assertTrue( Math.abs(sign) == 1.0 );
592 _st = sign*_st;
593
594 // save derivatives
595 _dst_dcrv1 = sign*d*d*arg*( 1.0/6.0 + 3.0/20.0*arg2);
596 double root = (1.0 + 0.5*arg*arg + 3.0/8.0*arg*arg*arg*arg );
597 _dst_dphi21 = sign*(r1*r2*Math.sin(phi2-phi1))*root/d;
598 _dst_dr1= (1.+arg2/2.*(1+3./4.*arg2))/d*sign;
599 _cnst1=r1-r2*Math.cos(phi2-phi1);
600 _cnst2=r1*r2*Math.sin(phi2-phi1);
601 }
602 }
603
604 public double st()
605 { return _st;
606 }
607
608 public double d_st_dalp1(double d_phi2_dalf1, double d_alf2_dalf1 )
609 {
610 if ( _big_crv ) return ( d_phi2_dalf1 + d_alf2_dalf1 - 1.0 ) / _crv1;
611 else return _dst_dphi21 * d_phi2_dalf1;
612
613 }
614 public double d_st_dcrv1(double d_phi2_dcrv1, double d_alf2_dcrv1 )
615 {
616 if ( _big_crv ) return ( d_phi2_dcrv1 + d_alf2_dcrv1 - _st ) / _crv1;
617 else return _dst_dcrv1 + _dst_dphi21*d_phi2_dcrv1;
618
619 }
620 public double d_st_dr1( double d_phi2_dr1, double d_alf2_dr1 )
621 {
622 if ( _big_crv ) return ( d_phi2_dr1 + d_alf2_dr1 ) / _crv1;
623 else return _dst_dr1*(_cnst1+_cnst2*d_phi2_dr1);
624 }
625 }
626 }