1 package org.lcsim.recon.tracking.trfdca;
2
3
4 import org.lcsim.recon.tracking.trfutil.Assert;
5 import org.lcsim.recon.tracking.trfutil.TRFMath;
6 import org.lcsim.recon.tracking.trfbase.PropDirected;
7 import org.lcsim.recon.tracking.trfbase.PropStat;
8 import org.lcsim.recon.tracking.trfbase.Surface;
9 import org.lcsim.recon.tracking.trfbase.VTrack;
10 import org.lcsim.recon.tracking.trfbase.TrackDerivative;
11 import org.lcsim.recon.tracking.trfbase.TrackVector;
12
13 import org.lcsim.recon.tracking.trfcyl.SurfCylinder;
14
15 import org.lcsim.recon.tracking.trfbase.Propagator;
16 import org.lcsim.recon.tracking.trfbase.PropDirected;
17 import org.lcsim.recon.tracking.trfbase.PropDir;
18
19 /**
20 * Propagates tracks from a Cylinder to a DCA surface.
21 *<p>
22 * Propagation will fail if either the origin is not a cylinder,
23 * or the destination is not a DCA surface.
24 *
25 *
26 *@author Norman A. Graf
27 *@version 1.0
28 *
29 */
30 public class PropCylDCA extends PropDirected
31 {
32
33 // static Attributes
34
35 // Assign track parameter indices
36
37
38 private static final int IPHI = SurfCylinder.IPHI;
39 private static final int IZ_CYL = SurfCylinder.IZ;
40 private static final int IALF = SurfCylinder.IALF;
41 private static final int ITLM_CYL = SurfCylinder.ITLM;
42 private static final int IQPT_CYL = SurfCylinder.IQPT;
43
44 private static final int IRSIGNED = SurfDCA.IRSIGNED;
45 private static final int IZ_DCA = SurfDCA.IZ;
46 private static final int IPHID = SurfDCA.IPHID;
47 private static final int ITLM_DCA = SurfDCA.ITLM;
48 private static final int IQPT_DCA = SurfDCA.IQPT;
49
50
51
52
53
54 // Attributes ***************************************************
55
56 // bfield * BFAC
57 private double _bfac;
58
59
60 // Methods ******************************************************
61
62 // static methods
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 typeName()
71 { return "PropCylDCA";
72 }
73
74 /**
75 *Return a String representation of the class' type name.
76 *Included for completeness with the C++ version.
77 *
78 * @return A String representation of the class' type name.
79 */
80 public String staticType()
81 { return typeName();
82 }
83
84 /**
85 *Construct an instance from a constant solenoidal magnetic field in Tesla.
86 *
87 * @param bfield The magnetic field strength in Tesla.
88 */
89 public PropCylDCA(double bfield)
90 {
91 _bfac = bfield * TRFMath.BFAC;
92 }
93
94 /**
95 *Clone an instance.
96 *
97 * @return A Clone of this instance.
98 */
99 public Propagator newPropagator()
100 {
101 return new PropCylDCA( bField() );
102 }
103
104 /**
105 *Return a String representation of the class' type name.
106 *Included for completeness with the C++ version.
107 *
108 * @return A String representation of the class' type name.
109 */
110 public String type()
111 { return staticType();
112 }
113
114
115 /**
116 *Propagate a track without error in the specified direction.
117 *
118 * The track parameters for a cylinder are:
119 * phi z alpha tan(lambda) curvature
120 *
121 * @param trv The VTrack to propagate.
122 * @param srf The Surface to which to propagate.
123 * @param dir The direction in which to propagate.
124 * @return The propagation status.
125 */
126 public PropStat vecDirProp( VTrack trv, Surface srf,
127 PropDir dir)
128 {
129 TrackDerivative deriv = null;
130 return cylDcaPropagate( _bfac, trv, srf, dir, deriv );
131 }
132
133
134 /**
135 *Propagate a track without error in the specified direction.
136 *and return the derivative matrix in deriv.
137 *
138 * The track parameters for a cylinder are:
139 * phi z alpha tan(lambda) curvature
140 *
141 * @param trv The VTrack to propagate.
142 * @param srf The Surface to which to propagate.
143 * @param dir The direction in which to propagate.
144 * @param deriv The track derivatives to update at the surface srf.
145 * @return The propagation status.
146 */
147 public PropStat vecDirProp( VTrack trv, Surface srf,
148 PropDir dir, TrackDerivative deriv)
149 {
150 return cylDcaPropagate( _bfac, trv, srf, dir, deriv );
151 }
152 //
153 // PropStat err_dir_prop( ETrack& tre, const Surface& srf,
154 // PropDir dir ) const;
155
156
157
158 /**
159 * Propagate a track with error in the specified direction.
160 *
161 * @param _bfac Numerical factor times the magnetic field strength.
162 * @param trv The VTrack to propagate.
163 * @param srf The Surface to which to propagate.
164 * @param dir The direction in which to propagate.
165 * @param deriv The track derivatives to update at the surface srf.
166 * @return The propagation status.
167 */
168 public PropStat cylDcaPropagate( double _bfac,
169 VTrack trv,
170 Surface srf,
171 PropDir dir,
172 TrackDerivative deriv )
173 {
174
175 // construct return status
176 PropStat pstat = new PropStat();
177
178 // fetch the originating Surface and check it is a Cylinder
179 Surface srf_scy = trv.surface();
180 Assert.assertTrue( srf_scy.pureType().equals(SurfCylinder.staticType()) );
181 if ( !srf_scy.pureType( ).equals(SurfCylinder.staticType()) )
182 return pstat;
183
184 // check that the destination surface is a DCA surface
185 Assert.assertTrue( srf.pureType().equals(SurfDCA.staticType()) );
186 if ( !srf.pureType( ).equals(SurfDCA.staticType()) )
187 return pstat;
188 SurfDCA srf_dca = ( SurfDCA ) srf;
189
190 // Check that dca surface has zero tilt.
191
192 boolean tilted = srf_dca.dXdZ() != 0 || srf_dca.dYdZ() != 0;
193 Assert.assertTrue(!tilted);
194 if(tilted) return pstat;
195
196 // original coordinates are marked with prime ( phip,alfp,r1p)
197 // coordinates centered at (xv,yv) are phi,alf,r1
198
199
200 // fetch the originating radius
201 double r1p = srf_scy.parameter(SurfCylinder.RADIUS);
202 double xv = srf_dca.x();
203 double yv = srf_dca.y();
204
205 // fetch the originating TrackVector
206 TrackVector vec_scy = trv.vector();
207 double phi1p = vec_scy.get(SurfCylinder.IPHI); // phi_position
208 double z1 = vec_scy.get(SurfCylinder.IZ); // z
209 double alf1p = vec_scy.get(SurfCylinder.IALF); // phi_dir - phi_pos
210 double tlam1 = vec_scy.get(SurfCylinder.ITLM); // tan(lambda)
211 double qpt1 = vec_scy.get(SurfCylinder.IQPT); // q/pT
212
213 double cosphi1p=Math.cos(phi1p);
214 double sinphi1p=Math.sin(phi1p);
215
216 double x1=r1p*cosphi1p-xv;
217 double y1=r1p*sinphi1p-yv;
218
219 double phi1 = Math.atan2(y1,x1); // phi position in (xv,yv) coord system
220 double r1 = Math.sqrt(x1*x1+y1*y1); // r1 in (xv,yv) coord system
221 double alf1 = alf1p + phi1p - phi1; // alf1 in (xv,yv) coord system
222
223
224 // calculate salf1, calf1 and crv1
225 double salf1 = Math.sin( alf1 );
226 double calf1 = Math.cos( alf1 );
227 double crv1 = _bfac * qpt1;
228
229 // calculate r2 and alf2 of the destination DCA Surface
230
231 double r2;
232 final double sign_alf1 = alf1 > 0.0 ? 1.0 : -1.0;
233 double sign_alf2 = 1.0;
234
235 double del = (r1*crv1 - 2.0*salf1);
236 double rcdel = r1*crv1*del;
237 double sroot = Math.sqrt(1.0 + rcdel);
238 double sroot_minus_one = sroot - 1.0;
239 if ( rcdel < 1.e-6 )
240 sroot_minus_one = 0.5*rcdel - 0.125*rcdel*rcdel + 0.0625*rcdel*rcdel*rcdel;
241
242 if ( Math.abs(del)<1e-14 )
243 {
244 sign_alf2 = sign_alf1;
245 r2 = 0.0;
246 }
247 else if ( TRFMath.isZero(crv1) )
248 {
249 sign_alf2 = sign_alf1;
250 r2 = r1 * Math.abs(salf1);
251 }
252 else
253 {
254 sign_alf2 = crv1*r1 > 2.0*salf1 ? -1.0 : 1.0;
255 r2 = -sign_alf2*sroot_minus_one/crv1;
256 Assert.assertTrue( r2 > 0.0 );
257 }
258 /*
259 if ( del == 0. )
260 {
261 sign_alf2 = sign_alf1;
262 r2 = 0.0;
263 }
264 else if ( crv1 == 0. )
265 {
266 sign_alf2 = sign_alf1;
267 r2 = r1p * Math.abs(salf1);
268 }
269 else
270 {
271 sign_alf2 = crv1*r1p > 2.0*salf1 ? -1.0 : 1.0;
272 r2 = -sign_alf2*sroot_minus_one/crv1;
273 Assert.assertTrue( r2 > 0.0 );
274 }
275 */
276
277 double alf2 = sign_alf2 * TRFMath.PI2;
278
279 // calculate phi2 of the destination DCA Surface
280 double sign_crv = 1.;
281
282 // calculate tlam2, qpt2 and crv2 of the destination DCA Surface
283 double tlam2 = tlam1;
284 double qpt2 = qpt1;
285 double crv2 = crv1;
286
287
288 double salf2 = sign_alf2 > 0 ? 1.: (sign_alf2 < 0 ? -1.: 0. );
289 //double calf2 = 0.;
290 double cnst = salf2-r2*crv2 > 0. ? TRFMath.PI2: (sign_alf2 == 0. ? 0. : -TRFMath.PI2);
291 double phi2 = phi1p + Math.atan2( salf1-r1p*crv1, calf1 )
292 - cnst;
293 if(crv1 == 0. ) phi2= phi1p + alf1 - cnst;
294
295 double phid2 = phi2 + sign_alf2 * TRFMath.PI2;
296 phid2 = TRFMath.fmod2( phid2, TRFMath.TWOPI );
297
298 // construct an object of ST_CylDCA
299 ST_CylDCA sto = new ST_CylDCA(r1,phi1,alf1,crv1,r2,phi2,alf2);
300 // fetch sT
301 double st = sto.st();
302 double s = st*Math.sqrt(1+tlam1*tlam1);
303
304 // calculate z2 of the destination DCA Surface
305 double z2 = z1 + st * tlam1;
306
307 // construct the destination TrackVector
308 TrackVector vec_dca = new TrackVector();
309 vec_dca.set(IRSIGNED , sign_alf2 * r2);
310 vec_dca.set(IZ_DCA , z2);
311 /*
312 #ifdef TRF_PHIRANGE
313 vec_dca(IPHID) = fmod1(phid2, TWOPI);
314 #else
315 vec_dca(IPHID) = phid2;
316 #endif
317 */
318 vec_dca.set(IPHID , phid2);
319 vec_dca.set(ITLM_DCA , tlam2);
320 vec_dca.set(IQPT_DCA , qpt2);
321
322 /*
323 // For axial tracks, zero z and tan(lambda).
324
325 if(trv.is_axial()) {
326 vec_dca(SurfDCA::IZ) = 0.;
327 vec_dca(SurfDCA::ITLM) = 0.;
328 }
329 */
330
331 // set the surface of trv to the destination DCA surface
332 trv.setSurface( srf_dca.newPureSurface() );
333
334 // set the vector of trv to the destination TrackVector (DCA coord.)
335 trv.setVector(vec_dca);
336
337 // set the direction of trv to the default value for VTrack on DCA
338 trv.setForward();
339
340 // set the return status
341 pstat.setPathDistance(s);
342
343 // exit now if user did not ask for error matrix
344 if ( deriv == null ) return pstat;
345
346
347 // calculate derivatives
348
349 double dr2_dalf1_or = calf1/(sign_alf2-crv1*r2);
350 double dr2_dalf1 = r1*dr2_dalf1_or;
351 double dr2_dcrv1 = (r2*r2-r1*r1)/(2.0*(sign_alf2-crv1*r2));
352 double dr2_dr1 = -sign_alf2/sroot*(r1*crv1 - salf1);
353
354 double dphi2_dalf1 = sign_crv*(1.0-r1*crv1*salf1)/(sroot*sroot);
355 double dphi2_dalf1_m1_or = sign_crv*(-crv1*salf1-crv1*del)/(sroot*sroot);
356 double dphi2_dcrv1 = -sign_crv*(r1*calf1)/(sroot*sroot);
357 double dphi2_dr1 = -crv1*calf1/(1.+r1*crv1*r1*crv1-2.*salf1*r1*crv1);
358
359
360 double dst_dalf1 = sto.d_st_dalf1(dr2_dalf1, dphi2_dalf1);
361 double dst_dalf1_or = sto.d_st_dalf1_or(r1,dr2_dalf1_or, dphi2_dalf1_m1_or);
362 double dst_dcrv1 = sto.d_st_dcrv1(dr2_dcrv1, dphi2_dcrv1);
363 double dst_dr1 = sto.d_st_dr1();
364
365
366 // derivatives between phi,alf and phip (prime) alfp
367 double dphi1_dphi1p_tr = r1p*Math.cos(phi1p-phi1);
368 double dalf1_dphi1p_tr = r1 - dphi1_dphi1p_tr;
369 double dalf1_dalf1p = 1.;
370 double dr1_dphi1p = r1p*Math.sin(phi1-phi1p);
371
372
373 // derivatives to original coordinates
374
375 double dr2_dphi1p = dr2_dalf1_or*dalf1_dphi1p_tr + dr2_dr1*dr1_dphi1p;
376 double dr2_dalf1p = dr2_dalf1 * dalf1_dalf1p;
377
378 double dst_dphi1p = dst_dalf1_or*dalf1_dphi1p_tr + dst_dr1*dr1_dphi1p;
379 double dst_dalf1p = dst_dalf1*dalf1_dalf1p;
380
381 double dphi2_dphi1p = -dphi1_dphi1p_tr*dphi2_dalf1_m1_or + dphi2_dr1*dr1_dphi1p + dphi2_dalf1;
382 double dphi2_dalf1p = dphi2_dalf1*dalf1_dalf1p;
383
384
385 // build the derivative matrix
386
387 // matrix<double>& deriv = *deriv;
388 // deriv.fill( 0.0 );
389
390 deriv.set(IRSIGNED, IPHI , sign_alf2 * dr2_dphi1p);
391 deriv.set(IRSIGNED, IZ_CYL , 0.0);
392 deriv.set(IRSIGNED, IALF , sign_alf2 * dr2_dalf1p);
393 deriv.set(IRSIGNED, ITLM_CYL , 0.0);
394 deriv.set(IRSIGNED, IQPT_CYL , sign_alf2 * _bfac * dr2_dcrv1);
395
396 deriv.set(IZ_DCA, IPHI , tlam1 * dst_dphi1p);
397 deriv.set(IZ_DCA, IZ_CYL , 1.0);
398 deriv.set(IZ_DCA, IALF , tlam1 * dst_dalf1p);
399 deriv.set(IZ_DCA, ITLM_CYL , st);
400 deriv.set(IZ_DCA, IQPT_CYL , tlam1 * _bfac * dst_dcrv1);
401
402 deriv.set(IPHID, IPHI , dphi2_dphi1p);
403 deriv.set(IPHID, IZ_CYL , 0.0);
404 deriv.set(IPHID, IALF , dphi2_dalf1p);
405 deriv.set(IPHID, ITLM_CYL , 0.0);
406 deriv.set(IPHID, IQPT_CYL , _bfac * dphi2_dcrv1);
407
408 deriv.set(ITLM_DCA, IPHI , 0.0);
409 deriv.set(ITLM_DCA, IZ_CYL , 0.0);
410 deriv.set(ITLM_DCA, IALF , 0.0);
411 deriv.set(ITLM_DCA, ITLM_CYL , 1.0);
412 deriv.set(ITLM_DCA, IQPT_CYL , 0.0);
413
414 deriv.set(IQPT_DCA, IPHI , 0.0);
415 deriv.set(IQPT_DCA, IZ_CYL , 0.0);
416 deriv.set(IQPT_DCA, IALF , 0.0);
417 deriv.set(IQPT_DCA, ITLM_CYL , 0.0);
418 deriv.set(IQPT_DCA, IQPT_CYL , 1.0);
419
420 // For axial tracks, zero all derivatives of or with respect to z or
421 // tan(lambda), that are not already zero. This will force the error
422 // matrix to have zero errors for z and tan(lambda).
423
424 /*
425 if(trv.is_axial()) {
426 deriv.set(IZ_DCA, IPHI, 0.);
427 deriv.set(IZ_DCA, IZ_CYL, 0.);
428 deriv.set(IZ_DCA, IALF, 0.);
429 deriv.set(IZ_DCA, ITLM_CYL, 0.);
430 deriv.set(IZ_DCA, IQPT_CYL, 0.);
431
432 deriv.set(ITLM_DCA, ITLM_CYL, 0.);
433 }
434 */
435
436 return pstat;
437
438 }
439
440
441 /**
442 *Return the strength of the magnetic field in Tesla.
443 *
444 * @return The strength of the magnetic field in Tesla.
445 */
446 public double bField()
447 {
448 return _bfac/TRFMath.BFAC;
449 }
450
451
452 /**
453 *output stream
454 * @return A String representation of this instance.
455 */
456 public String toString()
457 {
458 return "Propagation from Cylinder to DCA with constant "
459 + bField() + " Tesla field";
460
461 }
462
463
464 //**********************************************************************
465
466 // Private class ST_CylDCA.
467 //
468 // An ST_CylDCA_ object calculates sT (the signed transverse path length)
469 // and its derivatives w.r.t. alf1 and crv1. It is constructed from
470 // the starting (r1, phi1, alf1, crv1) and final track parameters
471 // (r2, phi2, alf2) assuming these are consistent. Methods are
472 // provided to retrieve sT and the two derivatives.
473
474 class ST_CylDCA
475 {
476
477 private boolean _big_crv;
478 private double _st;
479 private double _dst_dr2;
480 private double _dst_dcrv1;
481 private double _dst_dphi2;
482 private double _dst_dr1;
483 private double _dst_dphi2_or;
484 private double _crv1; // was public? need to check
485
486 // constructor
487 public ST_CylDCA()
488 {
489 }
490 public ST_CylDCA(double r1, double phi1, double alf1, double crv1,
491 double r2, double phi2, double alf2)
492 {
493
494 _crv1 = crv1;
495 Assert.assertTrue( r1 >= 0.0 );
496 Assert.assertTrue( r2 >= 0.0 );
497 double rmax = r1+r2;
498
499 // Calculate the change in xy direction
500 double phi_dir_diff = TRFMath.fmod2(phi2+alf2-phi1-alf1,TRFMath.TWOPI);
501 Assert.assertTrue( Math.abs(phi_dir_diff) <= Math.PI );
502
503 // Evaluate whether |C| is" big"
504 _big_crv = rmax*Math.abs(crv1) > 0.001;
505
506 // If the curvature is big we can use
507 // sT = (phi_dir2 - phi_dir1)/crv1
508 if ( _big_crv )
509 {
510 Assert.assertTrue( crv1 != 0.0 );
511 _st = phi_dir_diff/crv1;
512 }
513
514 // Otherwise, we calculate the straight-line distance
515 // between the points and use an approximate correction
516 // for the (small) curvature.
517 else
518 {
519
520 // evaluate the distance
521 double d = Math.sqrt( r1*r1 + r2*r2 - 2.0*r1*r2*Math.cos(phi2-phi1) );
522 double arg = 0.5*d*crv1;
523 double arg2 = arg*arg;
524 double st_minus_d = d*arg2*( 1.0/6.0 + 3.0/40.0*arg2 );
525 _st = d + st_minus_d;
526
527 // evaluate the sign
528 // We define a metric xsign = abs( (dphid-d*C)/(d*C) ).
529 // Because sT*C = dphid and d = abs(sT):
530 // xsign = 0 for sT > 0
531 // xsign = 2 for sT < 0
532 // Numerical roundoff will smear these predictions.
533 double sign = 0.0;
534 if ( crv1 != 0. )
535 {
536 double xsign = Math.abs( (phi_dir_diff - _st*crv1) / (_st*crv1) );
537 if ( xsign < 0.5 ) sign = 1.0;
538 if ( xsign > 1.5 && xsign < 3.0 ) sign = -1.0;
539 }
540 // If the above is indeterminate, assume zero curvature.
541 // In this case abs(alpha) decreases monotonically
542 // with sT. Track passing through origin has alpha = 0 on one
543 // side and alpha = +/-pi on the other. If both points are on
544 // the same side, we use dr/ds > 0 for |alpha|<pi/2.
545 if ( sign == 0. )
546 {
547 sign = 1.0;
548 if ( Math.abs(alf2) > Math.abs(alf1) ) sign = -1.0;
549 if ( Math.abs(alf2) == Math.abs(alf1) )
550 {
551 if ( Math.abs(alf2) < TRFMath.PI2 )
552 {
553 if ( r2 < r1 ) sign = -1.0;
554 }
555 else
556 {
557 if ( r2 > r1 ) sign = -1.0;
558 }
559 }
560 }
561
562 // Correct _st using the above sign.
563 Assert.assertTrue( Math.abs(sign) == 1.0 );
564 _st = sign*_st;
565
566 // save derivatives
567 if(TRFMath.isZero(d))
568 {
569 _dst_dcrv1 = 0.0;
570 _dst_dphi2 = sign*r1;
571 _dst_dphi2_or = sign;
572 _dst_dr2 = 0.0;
573 }
574 else
575 {
576 _dst_dcrv1 = sign*d*d*arg*( 1.0/6.0 + 3.0/20.0*arg2);
577 double root = (1.0 + 0.5*arg*arg + 3.0/8.0*arg*arg*arg*arg );
578 _dst_dphi2_or = sign*(r2*Math.sin(phi2-phi1))*root/d;
579 _dst_dphi2 = r1*_dst_dphi2_or;
580 _dst_dr2 = sign*(r2-r1*Math.cos(phi2-phi1))*root/d;
581 }
582
583 }
584 _dst_dr1 = -Math.cos(alf1)/(1.+r1*crv1*r1*crv1-2.*Math.sin(alf1)*r1*crv1);
585 }
586
587 public double st()
588 {
589 return _st;
590 }
591
592 public double d_st_dr1()
593 {
594 return _dst_dr1;
595 }
596
597 public double d_st_dalf1_or(double r1,double dr2_dalf1_or, double dphi2_dalf1_m1_or)
598 {
599 if ( _big_crv ) return ( dphi2_dalf1_m1_or ) / _crv1;
600 else return _dst_dr2 * dr2_dalf1_or + _dst_dphi2_or * (dphi2_dalf1_m1_or*r1+1.);
601 }
602
603 public double d_st_dalf1( double d_r2_dalf1, double d_phi2_dalf1 )
604 {
605 if ( _big_crv ) return ( d_phi2_dalf1 - 1.0 ) / _crv1;
606 else return _dst_dr2 * d_r2_dalf1 + _dst_dphi2 * d_phi2_dalf1;
607 }
608
609 public double d_st_dcrv1( double d_r2_dcrv1, double d_phi2_dcrv1 )
610 {
611 if ( _big_crv ) return ( d_phi2_dcrv1 - _st ) / _crv1;
612 else return _dst_dr2 * d_r2_dcrv1 + _dst_dphi2 * d_phi2_dcrv1+ _dst_dcrv1;
613
614 }
615
616
617 }
618
619 }