1 /*
2 * PropZCyl_Test.java
3 *
4 * Created on July 24, 2007, 10:42 PM
5 *
6 * $Id: PropZCyl_Test.java,v 1.1.1.1 2010/04/08 20:38:00 jeremy Exp $
7 */
8
9 package org.lcsim.recon.tracking.trfcylplane;
10
11 import junit.framework.TestCase;
12 import org.lcsim.recon.tracking.spacegeom.SpacePath;
13 import org.lcsim.recon.tracking.spacegeom.SpacePoint;
14 import org.lcsim.recon.tracking.trfbase.ETrack;
15 import org.lcsim.recon.tracking.trfbase.PropDir;
16 import org.lcsim.recon.tracking.trfbase.PropStat;
17 import org.lcsim.recon.tracking.trfbase.Propagator;
18 import org.lcsim.recon.tracking.trfbase.Surface;
19 import org.lcsim.recon.tracking.trfbase.TrackDerivative;
20 import org.lcsim.recon.tracking.trfbase.TrackError;
21 import org.lcsim.recon.tracking.trfbase.TrackVector;
22 import org.lcsim.recon.tracking.trfbase.VTrack;
23 import org.lcsim.recon.tracking.trfcyl.SurfCylinder;
24 import org.lcsim.recon.tracking.trfutil.Assert;
25 import org.lcsim.recon.tracking.trfutil.TRFMath;
26 import org.lcsim.recon.tracking.trfzp.SurfZPlane;
27
28 /**
29 *
30 * @author Norman Graf
31 */
32 public class PropZCyl_Test extends TestCase
33 {
34 private boolean debug;
35
36 // Assign track parameter indices.
37 private static final int IX = SurfZPlane.IX;
38 private static final int IY = SurfZPlane.IY;
39 private static final int IDXDZ = SurfZPlane.IDXDZ;
40 private static final int IDYDZ = SurfZPlane.IDYDZ;
41 private static final int IQP_Z = SurfZPlane.IQP;
42
43 private static final int IPHI = SurfCylinder.IPHI;
44 private static final int IZ = SurfCylinder.IZ;
45 private static final int IALF = SurfCylinder.IALF;
46 private static final int ITLM = SurfCylinder.ITLM;
47 private static final int IQPT = SurfCylinder.IQPT;
48
49
50 // compare two tracks without errors
51
52 private double compare(VTrack trv1,VTrack trv2)
53 {
54 Surface srf = trv2.surface();
55
56 Assert.assertTrue(trv1.surface().equals(srf));
57
58 double diff = srf.vecDiff(trv2.vector(),trv1.vector()).amax();
59
60 return diff;
61 }
62
63 //**********************************************************************
64 // compare two tracks with errors
65
66 private double[] compare(ETrack trv1,ETrack trv2 )
67 {
68 double[] tmp = new double[2];
69 Surface srf = trv2.surface();
70
71 Assert.assertTrue(trv1.surface().equals(srf));
72
73 tmp[0] = srf.vecDiff(trv2.vector(),trv1.vector()).amax();
74
75 TrackError dfc = trv2.error().minus(trv1.error());
76 tmp[1] = dfc.amax();
77
78
79 return tmp;
80 }
81 /** Creates a new instance of PropZCyl_Test */
82 public void testPropZCyl()
83 {
84 String ok_prefix = "PropZCyl (I): ";
85 String error_prefix = "PropZCyl test (E): ";
86
87 if(debug) System.out.println( ok_prefix
88 + "-------- Testing component PropZCyl. --------" );
89 {
90 if(debug) System.out.println( ok_prefix + "Test constructor." );
91 double BFIELD = 2.0;
92 PropZCyl prop = new PropZCyl(BFIELD);
93 if(debug) System.out.println( prop );
94 PropZCyl_Test tst = new PropZCyl_Test();
95 //********************************************************************
96
97 // Here we propagate some tracks both forward and backward and then
98 // the same track forward and backward but using method
99 // that we checked very thoroughly before.
100
101 if(debug) System.out.println( ok_prefix + "Check against correct propagation." );
102
103 PropZCyl propzcyl = new PropZCyl(BFIELD/TRFMath.BFAC);
104
105 double z1[] ={10., 10., 10., 10., 10., 10. };
106
107 int sign_dz[] ={ -1, 1, -1, 1, 1, -1 };
108 double x[] ={ -2., 2., 5., 5., -2., -2. };
109 double y[] ={ 3., 3., 3., 3., 3., 3. };
110 double dxdz[] ={ 1.5, -2.3, 0., 1.5, -1.5, 1.5 };
111 double dydz[] ={ 2.3, -1.5, -2.3, 0., -2.3, 0. };
112 double qp[] ={ 0.01, -0.01, 0.01, -0.01, -0.01, 0.01 };
113
114 double rcyl2[] ={ 20., 20., 20., 20., 20., 20. };
115 double rcyl2b[] ={ 20., 20., 20., 20., 20., 20. };
116
117 double maxdiff = 1.e-10;
118 double diff;
119 int ntrk = 6;
120 int i;
121
122 for ( i=0; i<ntrk; ++i )
123 {
124 if(debug) System.out.println( "********** Propagate track " + i + ". **********" );
125
126 PropStat pstat = new PropStat();
127
128 SurfZPlane szp1 = new SurfZPlane(z1[i]);
129 SurfCylinder scy2 = new SurfCylinder(rcyl2[i]);
130 SurfCylinder scy2b = new SurfCylinder(rcyl2b[i]);
131
132 TrackVector vec1 = new TrackVector();
133
134 vec1.set(IX , x[i]); // x
135 vec1.set(IY , y[i]); // y
136 vec1.set(IDXDZ , dxdz[i]); // dx/dz
137 vec1.set(IDYDZ , dydz[i]); // dy/dz
138 vec1.set(IQP_Z , qp[i]); // q/p
139
140 VTrack trv1 = new VTrack(szp1.newPureSurface(),vec1);
141 if (sign_dz[i]==1) trv1.setForward();
142 else trv1.setBackward();
143
144 if(debug) System.out.println( "\n starting: " + trv1 );
145
146 VTrack trv2f = new VTrack(trv1);
147 pstat = propzcyl.vecDirProp(trv2f,scy2,PropDir.FORWARD);
148 Assert.assertTrue( pstat.forward() );
149 if(debug) System.out.println( "\n forward: " + trv2f );
150
151 VTrack trv2f_my = new VTrack(trv1);
152 pstat = tst.vec_propzcyl(BFIELD,trv2f_my,scy2,PropDir.FORWARD);
153 Assert.assertTrue( pstat.forward() );
154 if(debug) System.out.println( "\n forward my: " + trv2f_my );
155 diff = tst.compare(trv2f_my,trv2f);
156
157 if(debug) System.out.println( "\n diff: " + diff );
158 Assert.assertTrue( diff < maxdiff );
159
160 VTrack trv2b = new VTrack(trv1);
161 pstat = propzcyl.vecDirProp(trv2b,scy2b,PropDir.BACKWARD);
162 Assert.assertTrue( pstat.backward() );
163 if(debug) System.out.println( "\n backward: " + trv2b );
164
165 VTrack trv2b_my = new VTrack(trv1);
166 pstat = tst.vec_propzcyl(BFIELD,trv2b_my,scy2b,PropDir.BACKWARD);
167 Assert.assertTrue( pstat.backward() );
168 if(debug) System.out.println( "\n backward my: " + trv2b_my );
169 diff = tst.compare(trv2b_my,trv2b);
170
171 if(debug) System.out.println( "\n diff: " + diff );
172 Assert.assertTrue( diff < maxdiff );
173
174 }
175 //********************************************************************
176
177 // Repeat the above with errors.
178 if(debug) System.out.println( ok_prefix + "Check against correct propagation with errors." );
179
180 double exx[] = { 0.01, 0.01, 0.01, 0.01, 0.01, 0.01 };
181 double exy[] = { 0.01, -0.01, 0.01, -0.01, 0.01, -0.01 };
182 double eyy[] = { 0.25, 0.25, 0.25, 0.25, 0.25, 0.25 };
183 double exdx[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004 };
184 double eydx[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004 };
185 double edxdx[] = { 0.01, 0.01, 0.01, 0.01, 0.01, 0.01 };
186 double exdy[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004 };
187 double edxdy[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004 };
188 double eydy[] = { 0.04, -0.04, 0.04, -0.04, 0.04, -0.04 };
189 double edydy[] = { 0.02, 0.02, 0.02, 0.02, 0.02, 0.02 };
190 double exqp[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004 };
191 double eyqp[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004 };
192 double edxqp[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004 };
193 double edyqp[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004 };
194 double eqpqp[] = { 0.01, 0.01, 0.01, 0.01, 0.01, 0.01 };
195
196 maxdiff = 1.e-8;
197 double ediff;
198 ntrk = 6;
199
200 for ( i=0; i<ntrk; ++i )
201 {
202 if(debug) System.out.println( "********** Propagate track " + i + ". **********" );
203
204 PropStat pstat = new PropStat();
205
206 SurfZPlane szp1 = new SurfZPlane(z1[i]);
207 SurfCylinder scy2 = new SurfCylinder(rcyl2[i]);
208 SurfCylinder scy2b = new SurfCylinder(rcyl2b[i]);
209
210 TrackVector vec1 = new TrackVector();
211
212 vec1.set(IX , x[i]); // x
213 vec1.set(IY , y[i]); // y
214 vec1.set(IDXDZ , dxdz[i]); // dx/dz
215 vec1.set(IDYDZ , dydz[i]); // dy/dz
216 vec1.set(IQP_Z , qp[i]); // q/p
217
218 TrackError err1 = new TrackError();
219
220 err1.set(IX,IX , exx[i]);
221 err1.set(IX,IY , exy[i]);
222 err1.set(IY,IY , eyy[i]);
223 err1.set(IX,IDXDZ , exdx[i]);
224 err1.set(IY,IDXDZ , eydx[i]);
225 err1.set(IDXDZ,IDXDZ , edxdx[i]);
226 err1.set(IX,IDYDZ , exdy[i]);
227 err1.set(IY,IDYDZ , eydy[i]);
228 err1.set(IDXDZ,IDYDZ , edxdy[i]);
229 err1.set(IDYDZ,IDYDZ , edydy[i]);
230 err1.set(IX,IQP_Z , exqp[i]);
231 err1.set(IY,IQP_Z , eyqp[i]);
232 err1.set(IDXDZ,IQP_Z , edxqp[i]);
233 err1.set(IDYDZ,IQP_Z , edyqp[i]);
234 err1.set(IQP_Z,IQP_Z , eqpqp[i]);
235
236 ETrack trv1 = new ETrack(szp1.newPureSurface(),vec1,err1);
237 if (sign_dz[i]==1)trv1.setForward();
238 else trv1.setBackward();
239
240 if(debug) System.out.println( "\n starting: " + trv1 );
241
242 ETrack trv2f = new ETrack(trv1);
243 pstat = propzcyl.errDirProp(trv2f,scy2,PropDir.FORWARD);
244 Assert.assertTrue( pstat.forward() );
245 if(debug) System.out.println( "\n forward: " + trv2f );
246
247 ETrack trv2f_my = new ETrack(trv1);
248 TrackDerivative deriv = new TrackDerivative();
249 pstat = tst.vec_propzcyl(BFIELD,trv2f_my,scy2,PropDir.FORWARD,deriv);
250 Assert.assertTrue( pstat.forward() );
251 TrackError err = trv2f_my.error();
252 trv2f_my.setError( err.Xform(deriv) );
253 if(debug) System.out.println( "\n forward my: " + trv2f_my );
254 double[] diffs = tst.compare(trv2f_my,trv2f);
255
256 if(debug) System.out.println( "\n diff: " + diffs[0] + ' ' + "ediff: "+ diffs[1] );
257 Assert.assertTrue( diffs[0] < maxdiff );
258 Assert.assertTrue( diffs[1] < maxdiff );
259
260
261 ETrack trv2b = new ETrack(trv1);
262 pstat = propzcyl.errDirProp(trv2b,scy2b,PropDir.BACKWARD);
263 Assert.assertTrue( pstat.backward() );
264 if(debug) System.out.println( "\n backward: " + trv2b );
265
266
267 ETrack trv2b_my = new ETrack(trv1);
268 pstat = tst.vec_propzcyl(BFIELD,trv2b_my,scy2b,PropDir.BACKWARD,deriv);
269 Assert.assertTrue( pstat.backward() );
270 err = trv2b_my.error();
271 trv2b_my.setError( err.Xform(deriv) );
272 if(debug) System.out.println( "\n backward my: " + trv2b_my );
273 diffs = tst.compare(trv2b_my,trv2b);
274
275 if(debug) System.out.println( "\n diff: " + diffs[0] + ' ' + "ediff: "+ diffs[1] );
276 Assert.assertTrue( diffs[0] < maxdiff );
277 Assert.assertTrue( diffs[1] < maxdiff );
278
279 }
280
281 //********************************************************************
282
283 if(debug) System.out.println( ok_prefix + "Test cloning." );
284 Assert.assertTrue( prop.newPropagator() != null );
285
286 //********************************************************************
287
288 if(debug) System.out.println( ok_prefix + "Test the field." );
289 Assert.assertTrue( prop.bField() == 2.0 );
290
291
292 }
293 //********************************************************************
294
295 if(debug) System.out.println( ok_prefix + "Test Zero Field Propagation." );
296 {
297 PropZCyl prop0 = new PropZCyl(0.0);
298 if(debug) System.out.println( prop0 );
299 Assert.assertTrue( prop0.bField() == 0. );
300
301 double z=10.;
302 Surface srf = new SurfZPlane(z);
303 VTrack trv0 = new VTrack(srf);
304 TrackVector vec = new TrackVector();
305 vec.set(SurfZPlane.IX, 2.);
306 vec.set(SurfZPlane.IY, 10.);
307 vec.set(SurfZPlane.IDXDZ, 4.);
308 vec.set(SurfZPlane.IDYDZ, 2.);
309 trv0.setVector(vec);
310 trv0.setForward();
311 Surface srf_to = new SurfCylinder(13.0);
312
313 VTrack trv = new VTrack(trv0);
314 VTrack trv_der = new VTrack(trv);
315 PropStat pstat = prop0.vecDirProp(trv,srf_to,PropDir.NEAREST);
316 Assert.assertTrue( pstat.success() );
317
318 Assert.assertTrue( pstat.forward() );
319 Assert.assertTrue(trv.surface().pureEqual(srf_to));
320
321 check_zero_propagation(trv0,trv,pstat);
322 check_derivatives(prop0,trv_der,srf_to);
323 if(debug) System.out.println("\n**trv= \n"+trv+"\n trv_der= \n"+trv_der);
324
325 srf_to = new SurfCylinder(trv0.spacePoint().rxy());
326 if(debug) System.out.println("\n**trv0.spacePoint()= "+trv0.spacePoint() );
327 if(debug) System.out.println("\n**trv0.spacePoint().rxy()= "+trv0.spacePoint().rxy());
328
329 trv = new VTrack(trv0);
330 trv_der = new VTrack(trv);
331 pstat = prop0.vecDirProp(trv,srf_to,PropDir.NEAREST);
332 Assert.assertTrue( pstat.success() );
333
334 Assert.assertTrue( pstat.same() );
335 Assert.assertTrue(trv.surface().pureEqual(srf_to));
336 if(debug) System.out.println("\n**trv= \n"+trv+"\n trv_der= \n"+trv_der);
337
338 check_zero_propagation(trv0,trv,pstat);
339 check_derivatives(prop0,trv_der,srf_to);
340
341 trv0.setBackward();
342 trv = new VTrack(trv0);
343 trv_der = new VTrack(trv);
344 pstat = prop0.vecDirProp(trv,srf_to,PropDir.NEAREST);
345 Assert.assertTrue( pstat.success() );
346
347 Assert.assertTrue( pstat.same() );
348 Assert.assertTrue(trv.surface().pureEqual(srf_to));
349
350 check_zero_propagation(trv0,trv,pstat);
351 check_derivatives(prop0,trv_der,srf_to);
352 }
353
354 if(debug) System.out.println( ok_prefix
355 + "------------- All tests passed. -------------" );
356 }
357
358 //********************************************************************
359 // Very well tested Z-Cyl propagator. Each new one can be tested against it
360
361 //**********************************************************************
362 // helpers
363 //**********************************************************************
364
365 //**********************************************************************
366 // The track parameters for a cylinder are:
367 // phi z alpha tan(lambda) curvature
368 // (NOT [. . . lambda .] as in TRF and earlier version of TRF++.)
369 //
370 // If pderiv is nonzero, return the derivative matrix there.
371 // On Cylinder:
372 // r (cm) is fixed
373 // 0 - phi
374 // 1 - z (cm)
375 // 2 - alp
376 // 3 - tlam
377 // 4 - q/pt pt is transverse momentum of a track, q is its charge
378 // On ZPlane:
379 // 0 - x (cm)
380 // 1 - y (cm)
381 // 2 - dx/dz
382 // 3 - dy/dz
383 // 4 - q/p p is momentum of a track, q is its charge
384 // If pderiv is nonzero, return the derivative matrix there.
385
386 PropStat
387 vec_propzcyl( double B, VTrack trv, Surface srf,
388 PropDir dir)
389 {
390 TrackDerivative deriv = null;
391 return vec_propzcyl(B, trv, srf, dir, deriv);
392 }
393
394
395
396 PropStat
397 vec_propzcyl( double B, VTrack trv, Surface srf,
398 PropDir dir,
399 TrackDerivative deriv )
400 {
401
402 // construct return status
403 PropStat pstat = new PropStat();
404
405 // fetch the originating surface and vector
406 Surface srf1 = trv.surface();
407 //TrackVector vec1 = trv.vector();
408
409 // Check origin is a Zplane.
410 Assert.assertTrue( srf1.pureType().equals(SurfZPlane.staticType()) );
411 if ( !srf1.pureType( ).equals(SurfZPlane.staticType()) )
412 return pstat;
413 SurfZPlane szp1 = ( SurfZPlane) srf1;
414
415 // Check destination is a cylinder.
416 Assert.assertTrue( srf.pureType().equals(SurfCylinder.staticType()) );
417 if ( ! srf.pureType( ).equals(SurfCylinder.staticType()) )
418 return pstat;
419 SurfCylinder scy2 = ( SurfCylinder) srf;
420
421
422 // Fetch the z of the plane and the starting track vector.
423 int iz = SurfZPlane.ZPOS;
424 double z = szp1.parameter(iz);
425
426 TrackVector vec = trv.vector();
427 double x = vec.get(IX); // v
428 double y = vec.get(IY); // z
429 double b = vec.get(IDXDZ); // dv/du
430 double a = vec.get(IDYDZ); // dz/du
431 double e = vec.get(IQP_Z); // q/p
432
433 // Fetch the radii and the starting track vector.
434 int irad = SurfCylinder.RADIUS;
435 double r2 = scy2.parameter(irad);
436
437 int sign_dz = 0;
438 if(trv.isForward()) sign_dz = 1;
439 if(trv.isBackward()) sign_dz = -1;
440 if(sign_dz == 0)
441 {
442 if(debug) System.out.println("PropZCyl._vec_propagate: Unknown direction of a track ");
443 System.exit(1);
444 }
445
446 // Calculate cylindrical coordinates
447
448 double cnst1=x>0?0.:Math.PI;
449 double cnst2=Math.PI;
450 if(a>0.&&b>0.&&sign_dz>0.) cnst2=0.;
451 if(a<0.&&b>0.&&sign_dz>0.) cnst2=0.;
452 if(a>0.&&b<0.&&sign_dz<0.) cnst2=0.;
453 if(a<0.&&b<0.&&sign_dz<0.) cnst2=0.;
454
455 double sign_y= y>0 ? 1.:-1;
456 double sign_a= a>0 ? 1.:-1;
457 if(y==0) sign_y=0.;
458 if(a==0) sign_a=0.;
459 double atnxy=(x!=0.?Math.atan(y/x):sign_y*Math.PI/2.);
460 double atnab=(b!=0.?Math.atan(a/b):sign_a*Math.PI/2.);
461
462 double phi1= TRFMath.fmod2(atnxy+cnst1,TRFMath.TWOPI);
463 double r1=Math.sqrt(x*x+y*y);
464 double z1=z;
465 double alp1= TRFMath.fmod2(atnab-phi1+cnst2,TRFMath.TWOPI);
466 double tlm1= sign_dz/Math.sqrt(a*a+b*b);
467 double qpt1=e*Math.sqrt((1+a*a+b*b)/(a*a+b*b));
468
469
470
471 // Check alpha range.
472 alp1 = TRFMath.fmod2( alp1, TRFMath.TWOPI );
473 Assert.assertTrue( Math.abs(alp1) <= Math.PI );
474
475 //if ( trv.is_forward() ) Assert.assertTrue( Math.abs(alp1) <= TRFMath.PI2 );
476 //else Assert.assertTrue( Math.abs(alp1) > TRFMath.PI2 );
477
478 // Calculate the cosine of lambda.
479 //double clam1 = 1.0/Math.sqrt(1+tlm1*tlm1);
480
481 // Calculate curvature: C = _bfac*(q/p)/cos(lambda)
482 // and its derivatives
483 // Assert.assertTrue( clam1 != 0.0 );
484 // double dcrv1_dqp1 = B/clam1;
485 // double crv1 = dcrv1_dqp1*qp1;
486 // double dcrv1_dtlm1 = crv1*clam1*clam1*tlm1;
487
488 // Calculate the curvature = _bfac*(q/pt)
489 double dcrv1_dqpt1 = B;
490 double crv1 = dcrv1_dqpt1*qpt1;
491 //double dcrv1_dtlm1 = 0.0;
492
493 // Evaluate the new track vector.
494 // See dla log I-044
495
496 // lambda and curvature do not change
497 double tlm2 = tlm1;
498 double crv2 = crv1;
499 double qpt2 = qpt1;
500
501 // We can evaluate sin(alp2), leaving two possibilities for alp2
502 // 1st solution: alp21, phi21, phid21, tht21
503 // 2nd solution: alp22, phi22, phid22, tht22
504 // evaluate phi2 to choose
505 double salp1 = Math.sin( alp1 );
506 double calp1 = Math.cos( alp1 );
507 double salp2 = r1/r2*salp1 + 0.5*crv1/r2*(r2*r2-r1*r1);
508 // if salp2 > 1, track does not cross cylinder
509 if ( Math.abs(salp2) > 1.0 ) return pstat;
510 double alp21 = Math.asin( salp2 );
511 double alp22 = alp21>0 ? Math.PI-alp21 : -Math.PI-alp21;
512 double calp21 = Math.cos( alp21 );
513 double calp22 = Math.cos( alp22 );
514 double phi20 = phi1 + Math.atan2( salp1-r1*crv1, calp1 );
515 double phi21 = phi20 - Math.atan2( salp2-r2*crv2, calp21 ); // phi position
516 double phi22 = phi20 - Math.atan2( salp2-r2*crv2, calp22 );
517 // Construct an sT object for each solution.
518 STCalcZ_CHECK sto1 = new STCalcZ_CHECK(r1,phi1,alp1,crv1,r2,phi21,alp21);
519 STCalcZ_CHECK sto2 = new STCalcZ_CHECK(r1,phi1,alp1,crv1,r2,phi22,alp22);
520 // Check the two solutions are nonzero and have opposite sign
521 // or at least one is nonzero.
522
523 // Choose the correct solution
524 boolean use_first_solution = false;
525
526 if (dir.equals(PropDir.NEAREST))
527 {
528 use_first_solution = Math.abs(sto2.st()) > Math.abs(sto1.st());
529 }
530 else if (dir.equals(PropDir.FORWARD))
531 {
532 use_first_solution = sto1.st() > 0.0;
533 }
534 else if (dir.equals(PropDir.BACKWARD))
535 {
536 use_first_solution = sto1.st() < 0.0;
537 }
538 else
539 {
540 if(debug) System.out.println("PropCyl._vec_propagate: Unknown direction." );
541 System.exit(1);
542 }
543
544 // Assign phi2, alp2 and sto2 for the chosen solution.
545 double phi2, alp2;
546 STCalcZ_CHECK sto;
547 double calp2;
548 if ( use_first_solution )
549 {
550 sto = sto1;
551 phi2 = phi21;
552 alp2 = alp21;
553 calp2 = calp21;
554 }
555 else
556 {
557 sto = sto2;
558 phi2 = phi22;
559 alp2 = alp22;
560 calp2 = calp22;
561 }
562
563 // fetch sT.
564 double st = sto.st();
565
566 // use sT to evaluate z2
567 double z2 = z1 + tlm1*st;
568
569 // Check alpha range.
570 Assert.assertTrue( Math.abs(alp2) <= Math.PI );
571
572 // put new values in vec
573 vec.set(IPHI , phi2);
574 vec.set(IZ , z2);
575 vec.set(IALF , alp2);
576 vec.set(ITLM , tlm2);
577 vec.set(IQPT , qpt2);
578
579 // Update trv
580 trv.setSurface(srf.newPureSurface());
581 trv.setVector(vec);
582 if ( Math.abs(alp2) <= TRFMath.PI2 ) trv.setForward();
583 else trv.setBackward();
584
585 // Set the return status.
586 if (st > 0) pstat.setForward() ;
587 else pstat.setBackward();
588
589 // exit now if user did not ask for error matrix.
590 if ( deriv == null) return pstat;
591
592 // Calculate derivatives.
593 // dphi1
594
595 double dphi1_dx= -y/(x*x+y*y);
596 double dphi1_dy= x/(x*x+y*y);
597
598 // dz1
599
600 //double dz1_dz=0.;
601
602 // dalf1
603
604 double dalp1_da= b/(a*a+b*b);
605 double dalp1_db= -a/(a*a+b*b);
606 double dalp1_dy= -dphi1_dy;
607 double dalp1_dx= -dphi1_dx;
608
609 // dr1
610
611 double dr1_dx= x/Math.sqrt(x*x+y*y);
612 double dr1_dy= y/Math.sqrt(x*x+y*y);
613
614 // dtlm1
615
616 double dtlm1_da= -sign_dz*a/Math.sqrt(a*a+b*b)/(a*a+b*b);
617 double dtlm1_db= -sign_dz*b/(a*a+b*b)/Math.sqrt(a*a+b*b);
618
619 // dcrv1
620
621 double dcrv1_de= Math.sqrt((1+a*a+b*b)/(a*a+b*b))*B;
622 double dcrv1_da= -a*e/Math.sqrt((a*a+b*b)*(1+a*a+b*b))/(a*a+b*b)*B;
623 double dcrv1_db= -e*b/Math.sqrt(1+a*a+b*b)/Math.sqrt(a*a+b*b)/(a*a+b*b)*B;
624
625 // alpha_2
626 double da2da1 = r1*calp1/r2/calp2;
627 double da2dc1 = (r2*r2-r1*r1)*0.5/r2/calp2;
628 double da2dr1 = (salp1-crv2*r1)/r2/calp2;
629
630 // phi2
631 double rcsal1 = r1*crv1*salp1;
632 double rcsal2 = r2*crv2*salp2;
633 double den1 = 1.0 + r1*r1*crv1*crv1 - 2.0*rcsal1;
634 double den2 = 1.0 + r2*r2*crv2*crv2 - 2.0*rcsal2;
635 double dp2dp1 = 1.0;
636 double dp2da1 = (1.0-rcsal1)/den1 - (1.0-rcsal2)/den2*da2da1;
637 double dp2dc1 = -r1*calp1/den1 + r2*calp2/den2
638 - (1.0-rcsal2)/den2*da2dc1;
639 double dp2dr1= -crv1*calp1/den1-(1.0-rcsal2)*da2dr1/den2;
640
641 // z2
642 //double dz2dz1 = 1.0;
643 double dz2dl1 = st;
644 double dz2da1 = tlm1*sto.d_st_dalp1(dp2da1, da2da1);
645 double dz2dc1 = tlm1*sto.d_st_dcrv1(dp2dc1, da2dc1);
646 double dz2dr1 = tlm1*sto.d_st_dr1( dp2dr1, da2dr1);
647
648
649 // final derivatives
650
651 // phi2
652 double dphi2_dx=dp2dp1*dphi1_dx+dp2da1*dalp1_dx+dp2dr1*dr1_dx;
653 double dphi2_dy=dp2dp1*dphi1_dy+dp2da1*dalp1_dy+dp2dr1*dr1_dy;
654 double dphi2_db=dp2da1*dalp1_db+dp2dc1*dcrv1_db;
655 double dphi2_da=dp2da1*dalp1_da+dp2dc1*dcrv1_da;
656 double dphi2_de=dp2dc1*dcrv1_de;
657
658 // alp2
659 double dalp2_dx= da2da1*dalp1_dx+da2dr1*dr1_dx;
660 double dalp2_dy= da2da1*dalp1_dy+da2dr1*dr1_dy;
661 double dalp2_db= da2da1*dalp1_db+da2dc1*dcrv1_db;
662 double dalp2_da= da2da1*dalp1_da+da2dc1*dcrv1_da;
663 double dalp2_de= da2dc1*dcrv1_de;
664
665 // crv2
666 double dcrv2_da=dcrv1_da;
667 double dcrv2_db=dcrv1_db;
668 double dcrv2_de=dcrv1_de;
669
670 // tlm2
671 double dtlm2_da= dtlm1_da;
672 double dtlm2_db= dtlm1_db;
673
674 // z2
675 double dz2_dx= dz2dr1*dr1_dx+dz2da1*dalp1_dx;
676 double dz2_dy= dz2dr1*dr1_dy+dz2da1*dalp1_dy;
677 double dz2_db= dz2da1*dalp1_db+dz2dl1*dtlm1_db+dz2dc1*dcrv1_db;
678 double dz2_da= dz2da1*dalp1_da+dz2dl1*dtlm1_da+dz2dc1*dcrv1_da;
679 double dz2_de= dz2dc1*dcrv1_de;
680
681
682 // Build derivative matrix.
683
684 deriv.set(IPHI,IX , dphi2_dx);
685 deriv.set(IPHI,IY , dphi2_dy);
686 deriv.set(IPHI,IDXDZ, dphi2_db);
687 deriv.set(IPHI,IDYDZ, dphi2_da);
688 deriv.set(IPHI,IQP_Z, dphi2_de);
689 deriv.set(IZ,IX , dz2_dx);
690 deriv.set(IZ,IY , dz2_dy);
691 deriv.set(IZ,IDXDZ , dz2_db);
692 deriv.set(IZ,IDYDZ , dz2_da);
693 deriv.set(IZ,IQP_Z, dz2_de);
694 deriv.set(IALF,IX , dalp2_dx);
695 deriv.set(IALF,IY , dalp2_dy);
696 deriv.set(IALF,IDXDZ , dalp2_db);
697 deriv.set(IALF,IDYDZ , dalp2_da);
698 deriv.set(IALF,IQP_Z , dalp2_de);
699 deriv.set(ITLM,IDXDZ , dtlm2_db);
700 deriv.set(ITLM,IDYDZ , dtlm2_da);
701 deriv.set(IQPT,IDXDZ , dcrv2_db/B);
702 deriv.set(IQPT,IDYDZ , dcrv2_da/B);
703 deriv.set(IQPT,IQP_Z , dcrv2_de/B);
704
705 return pstat;
706
707 }
708
709
710 // Private class STCalcZ_CHECK.
711 //
712 // An STCalcZ_CHECK_ object calculates sT (the signed transverse path length)
713 // and its derivatives w.r.t. alf1 and crv1. It is constructed from
714 // the starting (r1, phi1, alf1, crv1) and final track parameters
715 // (r2, phi2, alf2) assuming these are consistent. Methods are
716 // provided to retrieve sT and the two derivatives.
717
718 private class STCalcZ_CHECK
719 {
720
721 private boolean _big_crv;
722 private double _st;
723 private double _dst_dphi21;
724 private double _dst_dcrv1;
725 private double _dst_dr1;
726 private double _cnst1,_cnst2;
727 public double _crv1;
728
729 // constructor
730 public STCalcZ_CHECK()
731 {
732 }
733 public STCalcZ_CHECK(double r1, double phi1, double alf1, double crv1,
734 double r2, double phi2, double alf2)
735 {
736 _crv1 = crv1;
737 Assert.assertTrue( r1 > 0.0 );
738 Assert.assertTrue( r2 > 0.0 );
739 double rmax = r1+r2;
740
741 // Calculate the change in xy direction
742 double phi_dir_diff = TRFMath.fmod2(phi2+alf2-phi1-alf1,TRFMath.TWOPI);
743 Assert.assertTrue( Math.abs(phi_dir_diff) <= Math.PI );
744
745 // Evaluate whether |C| is" big"
746 _big_crv = rmax*Math.abs(crv1) > 0.001;
747
748 // If the curvature is big we can use
749 // sT = (phi_dir2 - phi_dir1)/crv1
750 if ( _big_crv )
751 {
752 Assert.assertTrue( crv1 != 0.0 );
753 _st = phi_dir_diff/crv1;
754 }
755
756 // Otherwise, we calculate the straight-line distance
757 // between the points and use an approximate correction
758 // for the (small) curvature.
759 else
760 {
761
762 // evaluate the distance
763 double d = Math.sqrt( r1*r1 + r2*r2 - 2.0*r1*r2*Math.cos(phi2-phi1) );
764 double arg = 0.5*d*crv1;
765 double arg2 = arg*arg;
766 double st_minus_d = d*arg2*( 1.0/6.0 + 3.0/40.0*arg2 );
767 _st = d + st_minus_d;
768
769 // evaluate the sign
770 // We define a metric xsign = abs( (dphid-d*C)/(d*C) ).
771 // Because sT*C = dphid and d = abs(sT):
772 // xsign = 0 for sT > 0
773 // xsign = 2 for sT < 0
774 // Numerical roundoff will smear these predictions.
775 double xsign = Math.abs( (phi_dir_diff - _st*crv1) / (_st*crv1) );
776 double sign = 0.0;
777 if ( crv1 != 0 )
778 {
779 if ( xsign < 0.5 ) sign = 1.0;
780 if ( xsign > 1.5 && xsign < 3.0 ) sign = -1.0;
781 }
782 // If the above is indeterminate, assume zero curvature.
783 // In this case abs(alpha) decreases monotonically
784 // with sT. Track passing through origin has alpha = 0 on one
785 // side and alpha = +/-pi on the other. If both points are on
786 // the same side, we use dr/ds > 0 for |alpha|<pi/2.
787 if ( sign == 0)
788 {
789 sign = 1.0;
790 if ( Math.abs(alf2) > Math.abs(alf1) ) sign = -1.0;
791 if ( Math.abs(alf2) == Math.abs(alf1) )
792 {
793 if ( Math.abs(alf2) < TRFMath.PI2 )
794 {
795 if ( r2 < r1 ) sign = -1.0;
796 }
797 else
798 {
799 if ( r2 > r1 ) sign = -1.0;
800 }
801 }
802 }
803
804 // Correct _st using the above sign.
805 Assert.assertTrue( Math.abs(sign) == 1.0 );
806 _st = sign*_st;
807
808 // save derivatives
809 _dst_dcrv1 = sign*d*d*arg*( 1.0/6.0 + 3.0/20.0*arg2);
810 double root = (1.0 + 0.5*arg*arg + 3.0/8.0*arg*arg*arg*arg );
811 _dst_dphi21 = sign*(r1*r2*Math.sin(phi2-phi1))*root/d;
812 _dst_dr1= (1.+arg2/2.*(1+3./4.*arg2))/d*sign;
813 _cnst1=r1-r2*Math.cos(phi2-phi1);
814 _cnst2=r1*r2*Math.sin(phi2-phi1);
815 }
816
817 }
818 public double st()
819 { return _st; };
820 public double d_st_dalp1(double d_phi2_dalf1, double d_alf2_dalf1 )
821 {
822 if ( _big_crv ) return ( d_phi2_dalf1 + d_alf2_dalf1 - 1.0 ) / _crv1;
823 else return _dst_dphi21 * d_phi2_dalf1;
824 }
825 public double d_st_dcrv1(double d_phi2_dcrv1, double d_alf2_dcrv1 )
826 {
827 if ( _big_crv ) return ( d_phi2_dcrv1 + d_alf2_dcrv1 - _st ) / _crv1;
828 else return _dst_dcrv1 + _dst_dphi21*d_phi2_dcrv1;
829
830 }
831 public double d_st_dr1( double d_phi2_dr1, double d_alf2_dr1 )
832 {
833 if ( _big_crv ) return ( d_phi2_dr1 + d_alf2_dr1 ) / _crv1;
834 else return _dst_dr1*(_cnst1+_cnst2*d_phi2_dr1);
835
836 }
837
838 }
839
840 private static void check_zero_propagation( VTrack trv0, VTrack trv, PropStat pstat)
841 {
842
843 SpacePoint sp = trv.spacePoint();
844 SpacePoint sp0 = trv0.spacePoint();
845
846 SpacePath sv = trv.spacePath();
847 SpacePath sv0 = trv0.spacePath();
848
849 Assert.assertTrue( Math.abs(sv0.dx() - sv.dx())<1e-7 );
850 Assert.assertTrue( Math.abs(sv0.dy() - sv.dy())<1e-7 );
851 Assert.assertTrue( Math.abs(sv0.dz() - sv.dz())<1e-7 );
852
853 double x0 = sp0.x();
854 double y0 = sp0.y();
855 double z0 = sp0.z();
856 double x1 = sp.x();
857 double y1 = sp.y();
858 double z1 = sp.z();
859
860 double dx = sv.dx();
861 double dy = sv.dy();
862 double dz = sv.dz();
863
864 double prod = dx*(x1-x0)+dy*(y1-y0)+dz*(z1-z0);
865 double moda = Math.sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0) + (z1-z0)*(z1-z0));
866 double modb = Math.sqrt(dx*dx+dy*dy+dz*dz);
867 double st = pstat.pathDistance();
868 Assert.assertTrue( Math.abs(prod-st) < 1.e-7 );
869 Assert.assertTrue( Math.abs(Math.abs(prod) - moda*modb) < 1.e-7 );
870 }
871
872 private static void check_derivatives( Propagator prop, VTrack trv0, Surface srf)
873 {
874 for(int i=0;i<4;++i)
875 for(int j=0;j<4;++j)
876 check_derivative(prop,trv0,srf,i,j);
877 }
878
879 private static void check_derivative( Propagator prop, VTrack trv0, Surface srf,int i,int j)
880 {
881
882 double dx = 1.e-3;
883 VTrack trv = new VTrack(trv0);
884 TrackVector vec = trv.vector();
885 boolean forward = trv0.isForward();
886
887 VTrack trv_0 = new VTrack(trv0);
888 TrackDerivative der = new TrackDerivative();
889 PropStat pstat = prop.vecProp(trv_0,srf,der);
890 Assert.assertTrue(pstat.success());
891
892 TrackVector tmp = new TrackVector(vec);
893 tmp.set(j, tmp.get(j)+dx);
894 trv.setVector(tmp);
895 if(forward) trv.setForward();
896 else trv.setBackward();
897
898 VTrack trv_pl = new VTrack(trv);
899 pstat = prop.vecProp(trv_pl,srf);
900 Assert.assertTrue(pstat.success());
901
902 TrackVector vecpl = trv_pl.vector();
903
904 tmp = new TrackVector(vec);
905 tmp.set(j, tmp.get(j)-dx);
906 trv.setVector(tmp);
907 if(forward) trv.setForward();
908 else trv.setBackward();
909
910 VTrack trv_mn = new VTrack(trv);
911 pstat = prop.vecProp(trv_mn,srf);
912 Assert.assertTrue(pstat.success());
913
914 TrackVector vecmn = trv_mn.vector();
915
916 double dy = (vecpl.get(i)-vecmn.get(i))/2.;
917
918 double dydx = dy/dx;
919
920 double didj = der.get(i,j);
921
922 if( Math.abs(didj) > 1e-10 )
923 {
924 if( Math.abs((dydx - didj)/didj) >= 1e-4 )
925 System.out.println("i="+ i +" j="+ j + " "+dydx + " "+ didj+'\n');
926 Assert.assertTrue( Math.abs((dydx - didj)/didj) < 1e-4 );
927 }
928 else
929 {
930 if( Math.abs((dydx - didj)/didj) >= 1e-4 )
931 System.out.println("i="+ i +" j="+ j + " "+dydx + " "+ didj+'\n');
932
933 Assert.assertTrue( Math.abs(dydx) < 1e-4 ) ;
934 }
935 }
936
937 }