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