View Javadoc

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   * Propagates tracks from one ZPlane to another in a constant field.
16   *<p>
17   * Propagation will fail if either the origin or destination is
18   * not a ZPlane.
19   * Propagator works incorrectly for tracks with very small curvatures
20   *
21   *
22   *@author Norman A. Graf
23   *@version 1.0
24   *
25   */
26  public class PropZZ extends PropDirected
27  {
28      
29      // attributes
30      
31      private double _bfac;
32      
33      // Assign track parameter indices.
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      // static methods
43      
44      //
45      
46      /**
47       *Return a String representation of the class' type name.
48       *Included for completeness with the C++ version.
49       *
50       * @return   A String representation of the class' type name.
51       */
52      public static String typeName()
53      { return "PropZZ"; }
54      
55      //
56      
57      /**
58       *Return a String representation of the class' type name.
59       *Included for completeness with the C++ version.
60       *
61       * @return   A String representation of the class' type name.
62       */
63      public static String staticType()
64      { return typeName(); }
65      
66      
67      
68      //
69      
70      /**
71       *Construct an instance from a constant solenoidal magnetic field in Tesla.
72       *
73       * @param   bfield The magnetic field strength in Tesla.
74       */
75      public PropZZ(double bfield)
76      {
77          _bfac = TRFMath.BFAC*bfield;
78      }
79      
80      //
81      
82       /**
83       *Clone an instance.
84       *
85       * @return A Clone of this instance.
86       */
87      public Propagator newPropagator( )
88      {
89          return new PropZZ( bField() );
90      }
91      
92      //
93      
94      /**
95       *Propagate a track without error in the specified direction.
96       *
97       * @param   trv The VTrack to propagate.
98       * @param   srf The Surface to which to propagate.
99       * @param   dir The direction in which to propagate.
100      * @return The propagation status.
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      *Propagate a track without error in the specified direction
113      *and return the derivative matrix in deriv.
114      *
115      * @param   trv The VTrack to propagate.
116      * @param   srf The Surface to which to propagate.
117      * @param   dir The direction in which to propagate.
118      * @param   pder The track derivatives to update at the surface srf.
119      * @return The propagation status.
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      *Return the strength of the magnetic field in Tesla.
132      *
133      * @return The strength of the magnetic field in Tesla.
134      */
135     public double bField()
136     {
137         return _bfac/TRFMath.BFAC;
138     }
139     
140     //
141     
142     /**
143      *Return a String representation of the class' type name.
144      *Included for completeness with the C++ version.
145      *
146      * @return   A String representation of the class' type name.
147      */
148     public String type()
149     { return staticType(); }
150     
151     
152     /**
153      *output stream
154      *
155      * @return  A String representation of this instance.
156      */
157     public String toString()
158     {
159         return "ZPlane-ZPlane propagation with constant "
160         +bField() +" Tesla field";
161         
162     }
163     
164     
165     // Private function to propagate a track without error
166     // The corresponding track parameters are:
167     // z (cm) is fixed
168     // 0 - x (cm)
169     // 1 - y (cm)
170     // 2 - dx/dz
171     // 3 - dy/dz
172     // 4 - q/p   p is momentum of a track, q is its charge
173     // If pderiv is nonzero, return the derivative matrix there.
174     
175     PropStat
176     vec_propagatezz_( double B, VTrack trv, Surface srf,
177     PropDir dir1,
178     TrackDerivative deriv )
179     {
180         
181         // construct return status
182         PropStat pstat = new PropStat();
183         PropDir dir = dir1; //need to check this
184         boolean move = Propagator.reduceDirection(dir);
185         if(move) dir = reduce(dir);
186         // fetch the originating surface and vector
187         Surface srf1 = trv.surface();
188         // TrackVector vec1 = trv.get_vector();
189         
190         // Check origin is a ZPlane.
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         // Check destination is a ZPlane.
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         // If surfaces are the same, we can return now.
204         if ( srf.pureEqual(srf1) )
205         {
206             // if XXX_MOVE was requested , return fail because track doesn't cross
207             // Z place in more than one place
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         // Fetch the zpos's of the planes and the starting track vector.
220         int izpos  = SurfZPlane.ZPOS;
221         
222         double zpos_0 = szp1.parameter(izpos);
223         double zpos_n = szp2.parameter(izpos);
224         
225         // Calculate difference in zpos
226         
227         double dz = zpos_n - zpos_0;
228         
229         TrackVector vec = trv.vector();
230         double a1 = vec.get(IX);                  // x
231         double a2 = vec.get(IY);                  // y
232         double a3 = vec.get(IDXDZ);               // dx/dz
233         double a4 = vec.get(IDYDZ);               // dy/dz
234         double a5 = vec.get(IQP);                 // q/p
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         // if delta_z*dz/dt >0 and backward - fail
246         // if delta_z*dz/dt <0 and forward  - fail
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             //   z-z propagation failed
260             if(sign_dz*dz<0.) return pstat;
261             flag_forward =  1;
262         }
263         else if (dir.equals(PropDir.BACKWARD) )
264         {
265             
266             //   z-z propagation failed
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         // Update trv
299         trv.setSurface(srf.newPureSurface());
300         trv.setVector(vec);
301         
302         // set new direction of the track
303         
304         if(sign_dz ==  1) trv.setForward();
305         if(sign_dz == -1) trv.setBackward();
306         
307         // Calculate the path distance s.
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         // Set the return status.
313         //(flag_forward==1)?pstat.set_forward():pstat.set_backward();
314         pstat.setPathDistance(ds);
315         
316         // Check the direction of propgation is consistent with the old
317         // calculation.
318         if ( flag_forward == 1 ) Assert.assertTrue( pstat.forward() );
319         else Assert.assertTrue( pstat.backward() );
320         
321         // exit now if user did not ask for error matrix.
322         if ( deriv == null ) return pstat;
323         
324         //da34_hat
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         // ddphi
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         // dRsin_phi
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         // dRcos_phi
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         // da1n first two are simple because dR,dphi _da1,_da2 = 0.
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         // da2n first two are simple because dR,dphi _da1,_da2 = 0.
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         // da3n first two are simple because dphi _da1,_da2 = 0.
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         // da4n first two are simple because dphi _da1,_da2 = 0.
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         // da5n
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; //need to check constness
463         boolean move = Propagator.reduceDirection(dir);
464         boolean same = szp2.pureEqual(szp1);
465         
466         // There is only one solution. Can't XXX_MOVE
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         // S= (zp-z0)/dz0
541         // a(dxdz dz0) b(dydz dz0) c(dz0)
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 }