1 package org.lcsim.recon.tracking.trfxyp;
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 * Propagates tracks from one XYPlane to another in a constant field.
17 *<p>
18 * Propagation will fail if either the origin or destination is
19 * not a XYPlane.
20 * Propagator works incorrectly for tracks with very small curvatures
21 *
22 *
23 *@author Norman A. Graf
24 *@version 1.0
25 *
26 */
27 public class PropXYXY extends PropDirected
28 {
29
30 // attributes
31
32 private double _bfac;
33
34 private static final int IV = SurfXYPlane.IV;
35 private static final int IZ = SurfXYPlane.IZ;
36 private static final int IDVDU = SurfXYPlane.IDVDU;
37 private static final int IDZDU = SurfXYPlane.IDZDU;
38 private static final int IQP = SurfXYPlane.IQP;
39
40 // static methods
41
42 /**
43 *Return a String representation of the class' type name.
44 *Included for completeness with the C++ version.
45 *
46 * @return A String representation of the class' type name.
47 */
48 public static String typeName()
49 { return "PropXYXY"; }
50
51 /**
52 *Return a String representation of the class' type name.
53 *Included for completeness with the C++ version.
54 *
55 * @return A String representation of the class' type name.
56 */
57 public static String staticType()
58 { return typeName(); }
59
60
61 /**
62 *Construct an instance from a constant solenoidal magnetic field in Tesla.
63 *
64 * @param bfield The magnetic field strength in Tesla.
65 */
66 public PropXYXY(double bfield)
67 {
68 _bfac = TRFMath.BFAC*bfield;
69 }
70
71 /**
72 *Clone an instance.
73 *
74 * @return A Clone of this instance.
75 */
76 public Propagator newPropagator( )
77 {
78 return new PropXYXY( bField() );
79 }
80
81 /**
82 *Propagate a track without error in the specified direction
83 *and return the derivative matrix in deriv.
84 *
85 * @param trv The VTrack to propagate.
86 * @param srf The Surface to which to propagate.
87 * @param dir The direction in which to propagate.
88 * @param deriv The track derivatives to update at the surface srf.
89 * @return The propagation status.
90 **/
91 public PropStat vecDirProp( VTrack trv, Surface srf,
92 PropDir dir,TrackDerivative deriv )
93 {
94 PropStat pstat = vec_propagatexyxy_( _bfac, trv, srf, dir, deriv );
95 return pstat;
96 }
97
98 /**
99 *Propagate a track without error in the specified direction.
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 *Return the strength of the magnetic field in Tesla.
115 *
116 * @return The strength of the magnetic field in Tesla.
117 */
118 public double bField()
119 {
120 return _bfac/TRFMath.BFAC;
121 }
122
123
124 /**
125 *Return a String representation of the class' type name.
126 *Included for completeness with the C++ version.
127 *
128 * @return A String representation of the class' type name.
129 */
130 public String type()
131 { return staticType(); }
132
133
134 /**
135 *output stream
136 *
137 * @return A String representation of this instance.
138 */
139 public String toString()
140 {
141 return "XYPlane-XYPlane propagation with constant "
142 + bField() + " Tesla field";
143 }
144
145
146 // Private function to determine dphi of the propagation
147 double direction(int flag_forward, int sign_dphi,
148 double du,
149 double norm, double cat,
150 double sinphi, double cosphi)
151 {
152
153 int sign_cat = 0;
154 if(du*flag_forward>0.) sign_cat = 1;
155 if(du*flag_forward<0.) sign_cat = -1;
156 if(du == 0.)
157 {
158 if(sinphi >=0. ) sign_cat = 1;
159 if(sinphi < 0. ) sign_cat = -1;
160 }
161
162 double sin_dphi = norm*sinphi + cat*cosphi*sign_cat;
163 double cos_dphi = -norm*cosphi + cat*sinphi*sign_cat;
164
165 int sign_sindphi = 0;
166 if(sin_dphi> 0.) sign_sindphi = 1;
167 if(sin_dphi< 0.) sign_sindphi = -1;
168 if(sin_dphi == 0.) sign_sindphi = sign_dphi;
169
170 if( Math.abs(cos_dphi)>1. )
171 {
172 double sn= -1.;
173 if ( cos_dphi > 0 ) sn=1.;
174 cos_dphi= sn*Math.sqrt(1.-sin_dphi*sin_dphi);
175 }
176 if( Math.abs(sin_dphi)>1. )
177 {
178 double sn= -1.;
179 if ( sin_dphi > 0 ) sn=1.;
180 sin_dphi= sn*Math.sqrt(1.-cos_dphi*cos_dphi);
181 }
182
183 double dphi = Math.PI*(sign_dphi - sign_sindphi) + sign_sindphi*Math.acos(cos_dphi);
184 return dphi;
185
186 }
187 //**********************************************************************
188
189 // Private function to propagate a track without error
190 // The corresponding track parameters are:
191 // u (cm) is fixed
192 // 0 - v (cm)
193 // 1 - z (cm)
194 // 2 - dv/du
195 // 3 - dz/du
196 // 4 - q/p p is momentum of a track, q is its charge
197 // If deriv is nonzero, return the derivative matrix there.
198
199 PropStat
200 vec_propagatexyxy_( double B, VTrack trv, Surface srf,
201 PropDir dir1,
202 TrackDerivative deriv )
203 {
204
205 // construct return status
206 PropStat pstat = new PropStat();
207
208 PropDir dir = dir1; //need to check constness of this
209 boolean move = Propagator.reduceDirection(dir);
210 if(move) dir = reduce(dir);
211 // fetch the originating surface and vector
212 Surface srf1 = trv.surface();
213 // TrackVector vec1 = trv.get_vector();
214
215 // Check origin is a XYPlane.
216 Assert.assertTrue( srf1.pureType().equals(SurfXYPlane.staticType()) );
217 if ( !srf1.pureType( ).equals(SurfXYPlane.staticType()) )
218 return pstat;
219 SurfXYPlane sxyp1 = (SurfXYPlane) srf1;
220
221 // Check destination is an XYPlane.
222 Assert.assertTrue( srf.pureType().equals(SurfXYPlane.staticType()) );
223 if ( !srf.pureType( ).equals(SurfXYPlane.staticType()) )
224 return pstat;
225 SurfXYPlane sxyp2 = (SurfXYPlane) srf;
226
227 // If surfaces are the same and no XXX_MOVE is requested we can return now.
228 boolean same = srf.pureEqual(srf1);
229 if ( same && ! move )
230 {
231 if(deriv!=null)
232 {
233 deriv.setIdentity();
234 }
235 pstat.setSame();
236 return pstat;
237 }
238
239
240 if ( same && dir1.equals(PropDir.NEAREST_MOVE) ) dir = PropDir.FORWARD;
241
242 //v00.63.04 cng 01/12/01
243 if( Math.abs(B) < 1e-7 ) return zeroBField(trv,sxyp1,sxyp2,dir1,deriv);
244
245 // Fetch the u's and phi's of the planes and the starting track vector.
246 TrackVector vec = trv.vector();
247 int iphi = SurfXYPlane.NORMPHI;
248 int idist = SurfXYPlane.DISTNORM;
249 double phi_0 = sxyp1.parameter(iphi);
250 double u_0 = sxyp1.parameter(idist);
251 double phi_n = sxyp2.parameter(iphi);
252 double u_n = sxyp2.parameter(idist);
253 if ( same && move )
254 {
255 int forw= trv.isForward()? 1 : -1 ;
256 int forw_dir = (dir == PropDir.FORWARD) ? 1 : -1 ;
257 int sn = vec.get(IDVDU) > 0. ? 1 : -1 ;
258 u_0 += 1.e-7 * (double)(forw*forw_dir*sn);
259 }
260
261 double b1_0 = vec.get(IV); // v
262 double b2_0 = vec.get(IZ); // z
263 double b3_0 = vec.get(IDVDU); // dv/du
264 double b4_0 = vec.get(IDZDU); // dz/du
265 double b5_0 = vec.get(IQP); // q/p
266
267 double phi_u = phi_0 - phi_n;
268
269 double cosphi_u = Math.cos(phi_u);
270 double sinphi_u = Math.sin(phi_u);
271
272 // check if du == 0 ( that is track moves parallel to the destination plane )
273 double du_du_0 = cosphi_u - b3_0*sinphi_u;
274 if(du_du_0==0.) return pstat;
275
276 double a_hat_u = 1./du_du_0;
277 double a_hat_u2 = a_hat_u*a_hat_u;
278
279 double u = u_0*cosphi_u - b1_0*sinphi_u;
280 double b1 = b1_0*cosphi_u + u_0*sinphi_u;
281 double b2 = b2_0;
282 double b3 = (b3_0*cosphi_u + sinphi_u)*a_hat_u;
283 double b4 = b4_0*a_hat_u;
284 double b5 = b5_0;
285
286 int sign_du0 = 0;
287 if(trv.isForward()) sign_du0 = 1;
288 if(trv.isBackward()) sign_du0 = -1;
289 if(sign_du0 == 0)
290 {
291 System.out.println( "PropXYXY._vec_propagate: Unknown direction of a track ");
292 System.exit(1);
293 }
294 int sign_du = 0;
295 if(du_du_0*sign_du0 > 0) sign_du = 1;
296 if(du_du_0*sign_du0 < 0) sign_du = -1;
297
298 // check that q/p != 0
299 Assert.assertTrue(b5 !=0. );
300
301 // 1/curv of the track is r
302 double r = 1/(b5*B)*Math.sqrt(1 + b3_0*b3_0)/Math.sqrt(1 + b3_0*b3_0 + b4_0*b4_0);
303 double b3_hat = Math.sqrt(1 + b3*b3);
304 double b34_hat = Math.sqrt(1 + b3*b3 + b4*b4);
305 double b3_hat2 = b3_hat*b3_hat;
306 double b34_hat2 = b34_hat*b34_hat;
307 double cosphi = -b3*sign_du/b3_hat;
308 double sinphi = sign_du/b3_hat;
309 double rsinphi = 1./(b5*B)*sign_du/b34_hat;
310 double rcosphi = -b3/(b5*B)*sign_du/b34_hat;
311
312 double du = u_n - u;
313 double norm = du/r - cosphi;
314
315 // xy-xy propagation failed : noway to the new plane
316 if(Math.abs(norm)>1.) return pstat;
317
318 double cat = Math.sqrt(1.-norm*norm);
319 int flag_forward = 0;
320 int sign_dphi = 0;
321
322 if( dir.equals(PropDir.NEAREST))
323 {
324 // try forward propagation
325 flag_forward = 1;
326 if(b5*B>0.) sign_dphi = 1;
327 if(b5*B<0.) sign_dphi = -1;
328 double dphi1=
329 direction(flag_forward,sign_dphi,du,norm, cat, sinphi, cosphi);
330 // try backward propagation
331 flag_forward = -1;
332 if(b5*B>0.) sign_dphi = -1;
333 if(b5*B<0.) sign_dphi = 1;
334 double dphi2=
335 direction(flag_forward,sign_dphi,du,norm, cat, sinphi, cosphi);
336 if( Math.abs(dphi2)>Math.abs(dphi1) )
337 {
338 flag_forward = -flag_forward;
339 sign_dphi = - sign_dphi;
340 }
341 }
342
343 else if( dir.equals(PropDir.FORWARD))
344 {
345 flag_forward = 1;
346 if(b5*B>0.) sign_dphi = 1;
347 if(b5*B<0.) sign_dphi = -1;
348 }
349 else if( dir.equals(PropDir.BACKWARD))
350 {
351 flag_forward = -1;
352 if(b5*B>0.) sign_dphi = -1;
353 if(b5*B<0.) sign_dphi = 1;
354 }
355 else
356 {
357 System.out.println("PropXYXY._vec_propagate: Unknown direction.");
358 System.exit(1);
359 }
360
361 int sign_cat = 0;
362 if(du*sign_dphi*b5*B>0.) sign_cat = 1;
363 if(du*sign_dphi*b5*B<0.) sign_cat = -1;
364 if(du == 0.)
365 {
366 if(sinphi >=0. ) sign_cat = 1;
367 if(sinphi < 0. ) sign_cat = -1;
368 }
369
370 double sin_dphi = norm*sinphi + cat*cosphi*sign_cat;
371 double cos_dphi = -norm*cosphi + cat*sinphi*sign_cat;
372 if(cos_dphi>1.) cos_dphi=1.;
373 if(cos_dphi<-1.) cos_dphi= -1.;
374
375 int sign_sindphi = 0;
376 if(sin_dphi> 0.) sign_sindphi = 1;
377 if(sin_dphi< 0.) sign_sindphi = -1;
378 if(sin_dphi == 0.) sign_sindphi = sign_dphi;
379
380 double dphi = Math.PI*(sign_dphi - sign_sindphi) + sign_sindphi*Math.acos(cos_dphi);
381
382 // check that I didn't make any mistakes
383
384 Assert.assertTrue(Math.abs(Math.sin(dphi) - sin_dphi)<1.e-5);
385
386 // check if dun == 0 ( that is track moves parallel to the destination plane)
387 double du_n_du = cos_dphi - b3*sin_dphi;
388 if(du_n_du==0.) return pstat;
389
390 double a_hat_dphi = 1./du_n_du;
391 double a_hat_dphi2 = a_hat_dphi*a_hat_dphi;
392 double c_hat_dphi = sin_dphi + b3*cos_dphi ;
393
394 double b1_n = b1 + rsinphi*(1-cos_dphi) - rcosphi*sin_dphi;
395 double b2_n = b2 + b4/(b5*B)*sign_du/b34_hat*dphi;
396 double b3_n = c_hat_dphi*a_hat_dphi;
397 double b4_n = b4*a_hat_dphi;
398 double b5_n = b5;
399
400 // double u_n_0 = u_n*cosphi_u + b1_n*sinphi_u;
401
402 // check if track crossed original plane during the propagation
403 // switch (dir) {
404 // if( dir.equals(PropDir.FORWARD:
405 // if((u_n_0 - u_0)*sign_du0<0) return pstat;
406 // break;
407 // if( dir.equals(PropDir.BACKWARD:
408 // if((u_n_0 - u_0)*sign_du0>0) return pstat;
409 // break;
410 // }
411
412 int sign_dun = 0;
413 if(du_n_du*sign_du > 0) sign_dun = 1;
414 if(du_n_du*sign_du < 0) sign_dun = -1;
415
416 vec.set(IV ,b1_n);
417 vec.set(IZ ,b2_n);
418 vec.set(IDVDU ,b3_n);
419 vec.set(IDZDU ,b4_n);
420 vec.set(IQP ,b5_n);
421
422 // Update trv
423 trv.setSurface(srf.newPureSurface());
424 trv.setVector(vec);
425
426 // set new direction of the track
427 if(sign_dun == 1) trv.setForward();
428 if(sign_dun == -1) trv.setBackward();
429
430 // Calculate sT.
431 double crv = B*b5;
432 double dv = b1_n - b1;
433 double dxy = Math.sqrt( du*du + dv*dv );
434 double arg = 0.5*crv*dxy;
435 double dst = dxy*TRFMath.asinrat(arg);
436
437 // Calculate s.
438 double dz = b2_n - b2;
439 Assert.assertTrue( flag_forward==1 || flag_forward==-1 );
440 double ds = ((double)(flag_forward))*Math.sqrt(dst*dst+dz*dz);
441
442 // Set the return status.
443 pstat.setPathDistance(ds);
444 //(flag_forward==1)?pstat.set_forward():pstat.set_backward();
445
446 // exit now if user did not ask for error matrix.
447 if ( deriv == null ) return pstat;
448
449 // du_d0
450
451 double du_db1_0 = - sinphi_u;
452
453 // db1_d0
454
455 double db1_db1_0 = cosphi_u;
456
457 // db3_d0
458
459 double db3_db3_0 = a_hat_u2;
460
461 // db4_d0
462
463 double db4_db3_0 = b4_0*sinphi_u*a_hat_u2;
464 double db4_db4_0 = a_hat_u;
465
466 // dr_d
467
468 double dr_db3 = r*b3*b4*b4/(b3_hat2*b34_hat2);
469 double dr_db4 = -r*b4/b34_hat2;
470 double dr_db5 = -r/b5;
471
472 // dcosphi_d
473
474 double dcosphi_db3 = - sign_du/b3_hat - cosphi*b3/b3_hat2;
475
476 // dsinphi_d
477
478 double dsinphi_db3 = - sinphi*b3/b3_hat2;
479
480 // dcat_d
481
482 double dcat_db3 = norm/cat*(du/(r*r)*dr_db3 + dcosphi_db3 );
483 double dcat_db4 = norm/cat* du/(r*r)*dr_db4;
484 double dcat_db5 = norm/cat* du/(r*r)*dr_db5;
485 double dcat_du = norm/(cat*r);
486
487 // dnorm_d
488
489 double dnorm_db3 = - du/(r*r)*dr_db3 - dcosphi_db3;
490 double dnorm_db4 = - du/(r*r)*dr_db4;
491 double dnorm_db5 = - du/(r*r)*dr_db5;
492 double dnorm_du = - 1./r;
493
494 // dcos_dphi_d
495
496 double dcos_dphi_db3 = - cosphi*dnorm_db3 - norm*dcosphi_db3 +
497 sign_cat*(sinphi*dcat_db3 + cat*dsinphi_db3);
498 double dcos_dphi_db4 = - cosphi*dnorm_db4 + sign_cat*sinphi*dcat_db4;
499 double dcos_dphi_db5 = - cosphi*dnorm_db5 + sign_cat*sinphi*dcat_db5;
500 double dcos_dphi_du = - cosphi*dnorm_du + sign_cat*sinphi*dcat_du;
501
502 // dsin_dphi_d
503
504 double dsin_dphi_db3 = sinphi*dnorm_db3 + norm*dsinphi_db3 +
505 sign_cat*(cosphi*dcat_db3 + cat*dcosphi_db3);
506 double dsin_dphi_db4 = sinphi*dnorm_db4 + sign_cat*cosphi*dcat_db4;
507 double dsin_dphi_db5 = sinphi*dnorm_db5 + sign_cat*cosphi*dcat_db5;
508 double dsin_dphi_du = sinphi*dnorm_du + sign_cat*cosphi*dcat_du;
509
510 // ddphi_d
511
512 double ddphi_db3;
513 double ddphi_db4;
514 double ddphi_db5;
515 double ddphi_du;
516 if(Math.abs(sin_dphi)>0.5)
517 {
518 ddphi_db3 = - dcos_dphi_db3/sin_dphi;
519 ddphi_db4 = - dcos_dphi_db4/sin_dphi;
520 ddphi_db5 = - dcos_dphi_db5/sin_dphi;
521 ddphi_du = - dcos_dphi_du /sin_dphi;
522 }
523 else
524 {
525 ddphi_db3 = dsin_dphi_db3/cos_dphi;
526 ddphi_db4 = dsin_dphi_db4/cos_dphi;
527 ddphi_db5 = dsin_dphi_db5/cos_dphi;
528 ddphi_du = dsin_dphi_du /cos_dphi;
529 }
530
531 // da_hat_dphi_d
532
533 double da_hat_dphi_db3 = - a_hat_dphi2*
534 (dcos_dphi_db3 - sin_dphi - b3*dsin_dphi_db3);
535 double da_hat_dphi_db4 = - a_hat_dphi2*(dcos_dphi_db4 - b3*dsin_dphi_db4);
536 double da_hat_dphi_db5 = - a_hat_dphi2*(dcos_dphi_db5 - b3*dsin_dphi_db5);
537 double da_hat_dphi_du = - a_hat_dphi2*(dcos_dphi_du - b3*dsin_dphi_du );
538
539 // dc_hat_dphi_d
540
541 double dc_hat_dphi_db3 = b3*dcos_dphi_db3 + dsin_dphi_db3 + cos_dphi;
542 double dc_hat_dphi_db4 = b3*dcos_dphi_db4 + dsin_dphi_db4;
543 double dc_hat_dphi_db5 = b3*dcos_dphi_db5 + dsin_dphi_db5;
544 double dc_hat_dphi_du = b3*dcos_dphi_du + dsin_dphi_du ;
545
546 // db1_n_d
547
548 double db1_n_db1 = 1;
549 double db1_n_db3 = (dr_db3*sinphi+r*dsinphi_db3)*(1-cos_dphi)
550 - rsinphi*dcos_dphi_db3
551 - dr_db3*cosphi*sin_dphi - r*dcosphi_db3*sin_dphi
552 - rcosphi*dsin_dphi_db3;
553 double db1_n_db4 = dr_db4*sinphi*(1-cos_dphi) - rsinphi*dcos_dphi_db4
554 - dr_db4*cosphi*sin_dphi - rcosphi*dsin_dphi_db4;
555 double db1_n_db5 = dr_db5*sinphi*(1-cos_dphi) - rsinphi*dcos_dphi_db5
556 - dr_db5*cosphi*sin_dphi - rcosphi*dsin_dphi_db5;
557 double db1_n_du = - rsinphi*dcos_dphi_du - rcosphi*dsin_dphi_du;
558
559 // db2_n_d
560
561 double db2_n_db2 = 1.;
562 double db2_n_db3 = 1./(b5*B)*b4*sign_du/b34_hat*
563 ( - dphi*b3/b34_hat2 + ddphi_db3);
564 double db2_n_db4 = 1./(b5*B)*sign_du/b34_hat*
565 ( - dphi*b4*b4/b34_hat2 + b4*ddphi_db4 + dphi );
566 double db2_n_db5 = 1./(b5*B)*b4*sign_du/b34_hat*( ddphi_db5 - dphi/b5);
567 double db2_n_du = 1./(b5*B)*b4*sign_du/b34_hat*ddphi_du;
568
569 // db3_n_d
570
571 double db3_n_db3 = a_hat_dphi*dc_hat_dphi_db3 + da_hat_dphi_db3*c_hat_dphi;
572 double db3_n_db4 = a_hat_dphi*dc_hat_dphi_db4 + da_hat_dphi_db4*c_hat_dphi;
573 double db3_n_db5 = a_hat_dphi*dc_hat_dphi_db5 + da_hat_dphi_db5*c_hat_dphi;
574 double db3_n_du = a_hat_dphi*dc_hat_dphi_du + da_hat_dphi_du *c_hat_dphi;
575
576 // db4_n_d
577
578 double db4_n_db3 = b4*da_hat_dphi_db3;
579 double db4_n_db4 = b4*da_hat_dphi_db4 + a_hat_dphi;
580 double db4_n_db5 = b4*da_hat_dphi_db5;
581 double db4_n_du = b4*da_hat_dphi_du;
582
583 // db5_n_d
584
585 // db1_n_d0
586
587 double db1_n_db1_0 = db1_n_du * du_db1_0 + db1_n_db1*db1_db1_0;
588 double db1_n_db2_0 = 0.;
589 double db1_n_db3_0 = db1_n_db3*db3_db3_0 + db1_n_db4*db4_db3_0;
590 double db1_n_db4_0 = db1_n_db4*db4_db4_0;
591 double db1_n_db5_0 = db1_n_db5;
592
593 // db2_n_d0
594
595 double db2_n_db1_0 = db2_n_du *du_db1_0;
596 double db2_n_db2_0 = db2_n_db2;
597 double db2_n_db3_0 = db2_n_db3*db3_db3_0 + db2_n_db4*db4_db3_0;
598 double db2_n_db4_0 = db2_n_db4*db4_db4_0;
599 double db2_n_db5_0 = db2_n_db5;
600
601 // db3_n_d0
602
603 double db3_n_db1_0 = db3_n_du*du_db1_0;
604 double db3_n_db2_0 = 0.;
605 double db3_n_db3_0 = db3_n_db3*db3_db3_0 + db3_n_db4*db4_db3_0;
606 double db3_n_db4_0 = db3_n_db4*db4_db4_0;
607 double db3_n_db5_0 = db3_n_db5;
608
609 // db4_n_d0
610
611 double db4_n_db1_0 = db4_n_du*du_db1_0;
612 double db4_n_db2_0 = 0.;
613 double db4_n_db3_0 = db4_n_db3*db3_db3_0 + db4_n_db4*db4_db3_0;
614 double db4_n_db4_0 = db4_n_db4*db4_db4_0;
615 double db4_n_db5_0 = db4_n_db5;
616
617 // db5_n_d0
618
619 double db5_n_db1_0 = 0.;
620 double db5_n_db2_0 = 0.;
621 double db5_n_db3_0 = 0.;
622 double db5_n_db4_0 = 0.;
623 double db5_n_db5_0 = 1.;
624
625
626 deriv.set(IV,IV ,db1_n_db1_0);
627 deriv.set(IV,IZ ,db1_n_db2_0);
628 deriv.set(IV,IDVDU ,db1_n_db3_0);
629 deriv.set(IV,IDZDU ,db1_n_db4_0);
630 deriv.set(IV,IQP ,db1_n_db5_0);
631 deriv.set(IZ,IV ,db2_n_db1_0);
632 deriv.set(IZ,IZ ,db2_n_db2_0);
633 deriv.set(IZ,IDVDU ,db2_n_db3_0);
634 deriv.set(IZ,IDZDU ,db2_n_db4_0);
635 deriv.set(IZ,IQP ,db2_n_db5_0);
636 deriv.set(IDVDU,IV ,db3_n_db1_0);
637 deriv.set(IDVDU,IZ ,db3_n_db2_0);
638 deriv.set(IDVDU,IDVDU ,db3_n_db3_0);
639 deriv.set(IDVDU,IDZDU ,db3_n_db4_0);
640 deriv.set(IDVDU,IQP ,db3_n_db5_0);
641 deriv.set(IDZDU,IV ,db4_n_db1_0);
642 deriv.set(IDZDU,IZ ,db4_n_db2_0);
643 deriv.set(IDZDU,IDVDU ,db4_n_db3_0);
644 deriv.set(IDZDU,IDZDU ,db4_n_db4_0);
645 deriv.set(IDZDU,IQP ,db4_n_db5_0);
646 deriv.set(IQP,IV ,db5_n_db1_0);
647 deriv.set(IQP,IZ ,db5_n_db2_0);
648 deriv.set(IQP,IDVDU ,db5_n_db3_0);
649 deriv.set(IQP,IDZDU ,db5_n_db4_0);
650 deriv.set(IQP,IQP ,db5_n_db5_0);
651
652 return pstat;
653
654 }
655
656 //cng
657 PropStat zeroBField( VTrack trv, SurfXYPlane sxyp1,
658 SurfXYPlane sxyp2, PropDir dir1,
659 TrackDerivative deriv)
660 {
661 // construct return status
662 PropStat pstat = new PropStat();
663
664 PropDir dir = dir1; //need to check constness of this
665 boolean move = Propagator.reduceDirection(dir);
666 boolean same = sxyp2.pureEqual(sxyp1);
667
668 // There is only one solution. Can't XXX_MOVE
669 if ( same && move ) return pstat;
670
671 if ( same )
672 {
673 if(deriv!=null)
674 {
675 deriv.setIdentity();
676 }
677 pstat.setSame();
678 return pstat;
679 }
680
681 TrackVector vec = trv.vector();
682 double v0 = vec.get(IV);
683 double z0 = vec.get(IZ);
684 double dvdu0 = vec.get(IDVDU);
685 double dzdu0 = vec.get(IDZDU);
686
687 double du0 =1.;
688 if( trv.isBackward() ) du0 = -1.;
689
690 double phi0 = sxyp1.parameter(SurfXYPlane.NORMPHI);
691 double cphi0 = Math.cos(phi0);
692 double sphi0 = Math.sin(phi0);
693 double u0 = sxyp1.parameter(SurfXYPlane.DISTNORM);
694
695 double phi1 = sxyp2.parameter(SurfXYPlane.NORMPHI);
696 double cphi1 = Math.cos(phi1);
697 double sphi1 = Math.sin(phi1);
698 double u1 = sxyp2.parameter(SurfXYPlane.DISTNORM);
699
700 double a = (cphi0 - dvdu0*sphi0)*du0;
701 double b = (sphi0 + dvdu0*cphi0)*du0;
702 double c = dzdu0*du0;
703
704 double x0 = u0*cphi0 - v0*sphi0;
705 double y0 = u0*sphi0 + v0*cphi0;
706
707 double ap = u1*cphi1;
708 double bp = u1*sphi1;
709 double cp = 0;
710
711 double xp = ap;
712 double yp = bp;
713 double zp = 0.0;
714
715 double denom = a*ap + b*bp + c*cp;
716
717 if( denom == 0. ) return pstat;
718
719 double S = ( (xp-x0)*ap + (yp-y0)*bp + (zp-z0)*cp )/denom;
720
721 double x1 = x0 + S*a;
722 double y1 = y0 + S*b;
723 double z1 = z0 + S*c;
724
725 boolean forward = S > 0. ? true : false;
726
727 if( dir == PropDir.FORWARD && !forward ) return pstat;
728 if( dir == PropDir.BACKWARD && forward ) return pstat;
729
730 double v1 = y1*cphi1 - x1*sphi1;
731
732 double v01 = y0*cphi1 - x0*sphi1;
733 double u01 = y0*sphi1 + x0*cphi1;
734 double z01 = z0;
735
736 if( u01 == u1 ) return pstat;
737
738 double dvdu1 = (v1-v01)/(u1-u01);
739 double dzdu1 = (z1-z01)/(u1-u01);
740
741 vec.set(IV, v1);
742 vec.set(IZ, z1);
743 vec.set(IDVDU, dvdu1);
744 vec.set(IDZDU, dzdu1);
745
746 // Update trv
747 trv.setSurface(sxyp2.newPureSurface());
748 trv.setVector(vec);
749 // set new direction of the track
750 if( b*sphi1 + a*cphi1 >0 ) trv.setForward();
751 else trv.setBackward();
752
753 // Calculate s.
754 double ds = S*Math.sqrt(a*a+b*b+c*c);
755
756 // Set the return status.
757 pstat.setPathDistance(ds);
758
759 if( deriv == null ) return pstat;
760
761 double dx0dv0 = -sphi0;
762 double dy0dv0 = cphi0;
763
764 double dadv_du = -sphi0*du0;
765 double dbdv_du = cphi0*du0;
766 double dcdz_du = du0;
767
768 double ddenomdv_du = dadv_du*ap + dbdv_du*bp;
769 double ddenomdz_du = dcdz_du*cp;
770
771
772 double dSdv0 = -(dx0dv0*ap + dy0dv0*bp)/denom;
773 double dSdz0 = -cp/denom;
774 double dSdv_du = -S/denom*ddenomdv_du;
775 double dSdz_du = -S/denom*ddenomdz_du;
776
777 double dy1dv0 = dy0dv0 + dSdv0*b;
778 double dx1dv0 = dx0dv0 + dSdv0*a;
779 double dx1dv_du = dSdv_du*a + S*dadv_du;
780 double dy1dv_du = dSdv_du*b + S*dbdv_du;
781
782 double du01dv0 = dy0dv0*sphi1 + dx0dv0*cphi1;
783 double dv01dv0 = dy0dv0*cphi1 - dx0dv0*sphi1;
784
785 double dv1dv0 = dy1dv0*cphi1 - dx1dv0*sphi1;
786 double dv1dv_du = dy1dv_du*cphi1 - dx1dv_du*sphi1;
787
788 double dz1dz0 = 1 + dSdz0*c;
789 double dz1dv0 = dSdv0*c;
790 double dz1dv_du = dSdv_du*c;
791 double dz1dz_du = dSdz_du*c + S*dcdz_du;
792
793 double dv_du1dv0 = ((dv1dv0-dv01dv0)*(u1-u01) + du01dv0*(v1-v01))/(u1-u01)/(u1-u01);
794 double dv_du1dv_du = dv1dv_du/(u1-u01);
795
796 double dz_du1dv0 = (dz1dv0*(u1-u01) + du01dv0*(z1-z01))/(u1-u01)/(u1-u01);
797 double dz_du1dz0 = (dz1dz0 - 1.)/(u1-u01);
798 double dz_du1dv_du = dz1dv_du/(u1-u01);
799 double dz_du1dz_du = dz1dz_du/(u1-u01);
800
801 //Set the track derivatives...
802 deriv.setIdentity();
803
804 deriv.set(IV,IV, dv1dv0);
805 deriv.set(IV,IDVDU, dv1dv_du);
806
807 deriv.set(IZ,IV, dz1dv0);
808 deriv.set(IZ,IZ, dz1dz0);
809 deriv.set(IZ,IDVDU, dz1dv_du);
810 deriv.set(IZ,IDZDU, dz1dz_du);
811
812 deriv.set(IDVDU,IV, dv_du1dv0);
813 deriv.set(IDVDU,IDVDU, dv_du1dv_du);
814
815 deriv.set(IDZDU,IV, dz_du1dv0);
816 deriv.set(IDZDU,IZ, dz_du1dz0);
817 deriv.set(IDZDU,IDVDU, dz_du1dv_du);
818 deriv.set(IDZDU,IDZDU, dz_du1dz_du);
819
820 deriv.set(IQP,IQP, 1.0);
821
822 return pstat;
823
824 }
825 //cng
826
827
828
829
830 }