1 /*
2 * PropCylZ_Test.java
3 *
4 * Created on July 24, 2007, 10:51 PM
5 *
6 * $Id: PropCylZ_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 PropCylZ_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_CYL = 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
51 // compare two tracks without errors
52
53 private double compare(VTrack trv1,VTrack trv2)
54 {
55 Surface srf = trv2.surface();
56
57 Assert.assertTrue(trv1.surface().equals(srf));
58
59 double diff = srf.vecDiff(trv2.vector(),trv1.vector()).amax();
60
61 return diff;
62 }
63
64 //**********************************************************************
65 // compare two tracks with errors
66
67 private double[] compare(ETrack trv1,ETrack trv2 )
68 {
69 double[] tmp = new double[2];
70 Surface srf = trv2.surface();
71
72 Assert.assertTrue(trv1.surface().equals(srf));
73
74 tmp[0] = srf.vecDiff(trv2.vector(),trv1.vector()).amax();
75
76 TrackError dfc = trv2.error().minus(trv1.error());
77 tmp[1] = dfc.amax();
78
79
80 return tmp;
81 }
82 /** Creates a new instance of PropCylZ_Test */
83 public void testPropCylZ()
84 {
85 String ok_prefix = "PropCylZ (I): ";
86 String error_prefix = "PropCylZ test (E): ";
87
88 if(debug) System.out.println( ok_prefix
89 + "-------- Testing component PropCylZ. --------" );
90 {
91 if(debug) System.out.println( ok_prefix + "Test constructor." );
92 double BFIELD = 2.0;
93 PropCylZ prop = new PropCylZ(BFIELD);
94 if(debug) System.out.println( prop );
95
96 PropCylZ_Test tst = new PropCylZ_Test();
97
98 //********************************************************************
99
100 // Here we propagate some tracks both forward and backward and then
101 // the same track forward and backward but using method
102 // that we checked very thoroughly before.
103
104 if(debug) System.out.println( ok_prefix + "Check against correct propagation." );
105
106 PropCylZ propcylz = new PropCylZ(BFIELD/TRFMath.BFAC);
107
108 double phi[] ={ Math.PI/5, Math.PI/6, Math.PI/6, 4/3*Math.PI, 5/3*Math.PI, 5/3*Math.PI };
109 double z[] ={ 1.5, -2.3, 0., 1.5, -1.5, -1.5 };
110 double alpha[] ={ 0.1, -0.1, 0., 0.2, -0.2, 0. };
111 double tlm[] ={ 2.3, -1.5, -2.3, 2.3, -2.3, 2.3 };
112 double qpt[] ={ 0.01, -0.01, 0.01, -0.01, -0.01, 0.01 };
113
114 double z2[] ={ 5.5, -6.0, -5.0, 6.0, -5.7, 4.0 };
115 double z2b[] ={ -4.0, 4.0, 5.0, -3.0, 3.0, -6.0 };
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 PropStat pstat = new PropStat();
126 SurfCylinder scy1 = new SurfCylinder(10.);
127 SurfZPlane szp2 = new SurfZPlane(z2[i]);
128 SurfZPlane szp2b = new SurfZPlane(z2b[i]);
129
130 TrackVector vec1 = new TrackVector();
131
132 vec1.set(IPHI , phi[i]); // phi
133 vec1.set(IZ_CYL , z[i]); // z
134 vec1.set(IALF , alpha[i]); // alpha
135 vec1.set(ITLM , tlm[i]); // tan(lambda)
136 vec1.set(IQPT , qpt[i]); // q/pt
137
138 VTrack trv1 = new VTrack(scy1.newPureSurface(),vec1);
139
140 if(debug) System.out.println( "\n starting: " + trv1 );
141
142 VTrack trv2f = new VTrack(trv1);
143 pstat = propcylz.vecDirProp(trv2f,szp2,PropDir.FORWARD);
144 Assert.assertTrue( pstat.forward() );
145 if(debug) System.out.println( "\n forward: " + trv2f );
146
147 VTrack trv2f_my = new VTrack(trv1);
148 pstat = tst.vec_propcylz(BFIELD,trv2f_my,szp2,PropDir.FORWARD);
149 Assert.assertTrue( pstat.forward() );
150 if(debug) System.out.println( "\n forward my: " + trv2f_my );
151 diff = tst.compare(trv2f_my,trv2f);
152
153 if(debug) System.out.println( "\n diff: " + diff );
154 Assert.assertTrue( diff < maxdiff );
155
156 VTrack trv2b = new VTrack(trv1);
157 pstat = propcylz.vecDirProp(trv2b,szp2b,PropDir.BACKWARD);
158 Assert.assertTrue( pstat.backward() );
159 if(debug) System.out.println( "\n backward: " + trv2b );
160
161 VTrack trv2b_my = new VTrack(trv1);
162 pstat = tst.vec_propcylz(BFIELD,trv2b_my,szp2b,PropDir.BACKWARD);
163 Assert.assertTrue( pstat.backward() );
164 if(debug) System.out.println( "\n backward my: " + trv2b_my );
165 diff = tst.compare(trv2b_my,trv2b);
166
167 if(debug) System.out.println( "\n diff: " + diff );
168 Assert.assertTrue( diff < maxdiff );
169
170 }
171 //********************************************************************
172
173 // Repeat the above with errors.
174 if(debug) System.out.println( ok_prefix + "Check against correct propagation with errors." );
175
176 double epp[] = { 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01 };
177 double epz[] = { 0.01, -0.01, 0.01, -0.01, 0.01, -0.01, 0.01 };
178 double ezz[] = { 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25 };
179 double epa[] = { 0.003, -0.003, 0.003, -0.003, 0.003, -0.003, 0.003 };
180 double eza[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004 };
181 double eaa[] = { 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01 };
182 double epl[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004 };
183 double eal[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004 };
184 double ezl[] = { 0.04, -0.04, 0.04, -0.04, 0.04, -0.04, 0.04 };
185 double ell[] = { 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02 };
186 double epc[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004 };
187 double ezc[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004 };
188 double eac[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004 };
189 double elc[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004 };
190 double ecc[] = { 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01 };
191
192 maxdiff = 1.e-8;
193 double ediff;
194 ntrk = 6;
195
196 for ( i=0; i<ntrk; ++i )
197 {
198 if(debug) System.out.println( "********** Propagate track " + i + ". **********" );
199 PropStat pstat = new PropStat();
200 SurfCylinder scy1 = new SurfCylinder(10.);
201 SurfZPlane szp2 = new SurfZPlane(z2[i]);
202 SurfZPlane szp2b = new SurfZPlane(z2b[i]);
203
204 TrackVector vec1 = new TrackVector();
205
206 vec1.set(IPHI , phi[i]); // phi
207 vec1.set(IZ_CYL, z[i]); // z
208 vec1.set(IALF , alpha[i]); // alpha
209 vec1.set(ITLM , tlm[i]); // tan(lambda)
210 vec1.set(IQPT , qpt[i]); // q/pt
211
212 TrackError err1 = new TrackError();
213
214 err1.set(IPHI,IPHI , epp[i]);
215 err1.set(IPHI,IZ_CYL , epz[i]);
216 err1.set(IZ_CYL,IZ_CYL, ezz[i]);
217 err1.set(IPHI,IALF , epa[i]);
218 err1.set(IZ_CYL,IALF , eza[i]);
219 err1.set(IALF,IALF , eaa[i]);
220 err1.set(IPHI,ITLM , epl[i]);
221 err1.set(IZ_CYL,ITLM , ezl[i]);
222 err1.set(IALF,ITLM , eal[i]);
223 err1.set(ITLM,ITLM , ell[i]);
224 err1.set(IPHI,IQPT , epc[i]);
225 err1.set(IZ_CYL,IQPT , ezc[i]);
226 err1.set(IALF,IQPT , eac[i]);
227 err1.set(ITLM,IQPT , elc[i]);
228 err1.set(IQPT,IQPT , ecc[i]);
229
230 ETrack trv1 = new ETrack(scy1.newPureSurface(),vec1,err1);
231
232 if(debug) System.out.println( "\n starting: " + trv1 );
233
234 ETrack trv2f = new ETrack(trv1);
235 pstat = propcylz.errDirProp(trv2f,szp2,PropDir.FORWARD);
236 Assert.assertTrue( pstat.forward() );
237 if(debug) System.out.println( "\n forward: " + trv2f );
238
239 ETrack trv2f_my = new ETrack(trv1);
240 TrackDerivative deriv = new TrackDerivative();
241 pstat = tst.vec_propcylz(BFIELD,trv2f_my,szp2,PropDir.FORWARD,deriv);
242 Assert.assertTrue( pstat.forward() );
243 TrackError err = trv2f_my.error();
244 trv2f_my.setError( err.Xform(deriv) );
245 if(debug) System.out.println( "\n forward my: " + trv2f_my );
246 double[] diffs = tst.compare(trv2f_my,trv2f);
247
248 if(debug) System.out.println( "\n diff: " + diffs[0] + ' ' + "ediff: "+ diffs[1] );
249 Assert.assertTrue( diffs[0] < maxdiff );
250 Assert.assertTrue( diffs[1] < maxdiff );
251
252
253 ETrack trv2b = new ETrack(trv1);
254 pstat = propcylz.errDirProp(trv2b,szp2b,PropDir.BACKWARD);
255 Assert.assertTrue( pstat.backward() );
256 if(debug) System.out.println( "\n backward: " + trv2b );
257
258 ETrack trv2b_my = new ETrack(trv1);
259 pstat = tst.vec_propcylz(BFIELD,trv2b_my,szp2b,PropDir.BACKWARD,deriv);
260 Assert.assertTrue( pstat.backward() );
261 err = trv2b_my.error();
262 trv2b_my.setError( err.Xform(deriv) );
263 if(debug) System.out.println( "\n backward my: " + trv2b_my );
264 diffs = tst.compare(trv2b_my,trv2b);
265
266 if(debug) System.out.println( "\n diff: " + diffs[0] + ' ' + "ediff: "+ diffs[1] );
267 Assert.assertTrue( diffs[0] < maxdiff );
268 Assert.assertTrue( diffs[1] < maxdiff );
269
270 }
271
272 //********************************************************************
273
274 if(debug) System.out.println( ok_prefix + "Test cloning." );
275 Assert.assertTrue( prop.newPropagator() != null );
276
277 //********************************************************************
278
279 if(debug) System.out.println( ok_prefix + "Test the field." );
280 Assert.assertTrue( prop.bField() == 2.0 );
281
282 }
283
284
285 //********************************************************************
286
287 if(debug) System.out.println( ok_prefix + "Test Zero Field Propagation." );
288 {
289 PropCylZ prop0 = new PropCylZ(0.0);
290 if(debug) System.out.println( prop0 );
291 Assert.assertTrue( prop0.bField() == 0. );
292
293 double z=6.;
294 Surface srf = new SurfCylinder(13.0);
295 VTrack trv0 = new VTrack(srf);
296 TrackVector vec = new TrackVector();
297 vec.set(SurfCylinder.IPHI, 0.1);
298 vec.set(SurfCylinder.IZ, 10.);
299 vec.set(SurfCylinder.IALF, 0.1);
300 vec.set(SurfCylinder.ITLM, 2.);
301 vec.set(SurfCylinder.IQPT, 0.);
302 trv0.setVector(vec);
303 trv0.setForward();
304 Surface srf_to = new SurfZPlane(z);
305
306 VTrack trv = new VTrack(trv0);
307 VTrack trv_der = new VTrack(trv);
308 PropStat pstat = prop0.vecDirProp(trv,srf_to,PropDir.NEAREST);
309 if(debug) System.out.println(trv);
310 Assert.assertTrue( pstat.success() );
311
312 Assert.assertTrue( pstat.backward() );
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 srf_to = new SurfZPlane(10.);
319 trv = new VTrack(trv0);
320 trv_der = new VTrack(trv);
321 pstat = prop0.vecDirProp(trv,srf_to,PropDir.NEAREST);
322 if(debug) System.out.println(trv);
323 Assert.assertTrue( pstat.success() );
324
325 Assert.assertTrue( pstat.same() );
326 Assert.assertTrue(trv.surface().pureEqual(srf_to));
327
328 check_zero_propagation(trv0,trv,pstat);
329 check_derivatives(prop0,trv_der,srf_to);
330
331 }
332
333 //********************************************************************
334
335 if(debug) System.out.println( ok_prefix
336 + "------------- All tests passed. -------------" );
337 //********************************************************************
338 }
339
340
341 // Very well tested Cyl-Z propagator. Each new one can be tested against it
342 private PropStat
343 vec_propcylz( double B, VTrack trv, Surface srf,
344 PropDir dir )
345 {
346 TrackDerivative deriv = null;
347 return vec_propcylz(B, trv, srf, dir, deriv);
348 }
349 private PropStat
350 vec_propcylz( double B, VTrack trv, Surface srf,
351 PropDir dir,
352 TrackDerivative deriv )
353 {
354
355 // phi, z, alpha=phi_dir-phi, tan(lambda)=pz/pt, q/pt
356 // c1 c2 c3 c4 c5
357 //
358 // x y dx/dz dy/dz q/p
359 // a1 a2 a3 a4 a5
360 //
361 // phi is polar angle of track momentum, R radius of the track in field B
362 // dphi is angle of rotation of momentum
363 // Rcyl is cylinder radius
364 //
365 // a1n = a1 + Rsin_phi*(cos_dphi - 1) + Rcos_phi*sin_dphi;
366 // a2n = a2 - Rcos_phi*(cos_dphi - 1) + Rsin_phi*sin_dphi;
367 // a3n = a3*cos_dphi - a4*sin_dphi;
368 // a4n = a3*sin_dphi + a4*cos_dphi;
369 // a5n = a5;
370
371 // a1 = Rcyl*cos_c1;
372 // a2 = Rcyl*sin_c1;
373 // a3 = cos(c1+c3)/c4;
374 // a4 = sin(c1+c3)/c4;
375 // a5 = c5/sqrt(1+c4*c4);
376 //
377 // dz = zpos - c2
378
379
380 // construct return status
381 PropStat pstat = new PropStat();;
382 // fetch the originating surface and vector
383 Surface srf1 = trv.surface();
384 // TrackVector vec1 = trv.vector();
385
386 // Check origin is a Cylinder.
387 Assert.assertTrue( srf1.pureType().equals(SurfCylinder.staticType()) );
388 if ( ! srf1.pureType( ).equals(SurfCylinder.staticType()) )
389 return pstat;
390 SurfCylinder scy1 = (SurfCylinder) srf1;
391
392
393 // Fetch the R of the cylinder and the starting track vector.
394 int ir = SurfCylinder.RADIUS;
395 double Rcyl = scy1.parameter(ir);
396
397 TrackVector vec = trv.vector();
398 double c1 = vec.get(IPHI); // phi
399 double c2 = vec.get(IZ_CYL); // z
400 double c3 = vec.get(IALF); // alpha
401 double c4 = vec.get(ITLM); // tan(lambda)
402 double c5 = vec.get(IQPT); // q/pt
403
404 // check that dz != 0
405 if(c4 == 0.) return pstat;
406
407 double cos_c1 = Math.cos(c1);
408 double sin_c1 = Math.sin(c1);
409 double cos_dir = Math.cos(c1+c3);
410 double sin_dir = Math.sin(c1+c3);
411 double c4_hat2 = 1+c4*c4;
412 double c4_hat = Math.sqrt(c4_hat2);
413
414 double a1 = Rcyl*cos_c1;
415 double a2 = Rcyl*sin_c1;
416 double a3 = cos_dir/c4;
417 double a4 = sin_dir/c4;
418 double a5 = c5/c4_hat;
419
420
421 // if tan(lambda)=dz/dst > 0 : dz>0; dzdst < 0 dz<0
422
423 int sign_dz = 0;
424 if(c4 > 0) sign_dz = 1;
425 if(c4 < 0) sign_dz = -1;
426
427
428 // Check destination is a ZPlane.
429 Assert.assertTrue( srf.pureType().equals(SurfZPlane.staticType()) );
430 if ( ! srf.pureType( ).equals(SurfZPlane.staticType()) )
431 return pstat;
432 SurfZPlane szp2 = (SurfZPlane) srf;
433
434 // Fetch the zpos's of the planes.
435 int izpos = SurfZPlane.ZPOS;
436
437 double zpos = szp2.parameter(izpos);
438
439 double dz = zpos - c2;
440
441
442 if (dir.equals(PropDir.NEAREST))
443 {
444
445 if(debug) System.out.println("vec_propcylz: Propagation in NEAREST direction is undefined");
446 System.exit(1);
447 }
448 else if (dir.equals(PropDir.FORWARD))
449 {
450 // cyl-z propagation failed
451 if(sign_dz*dz<0.) return pstat;
452 }
453 else if (dir.equals(PropDir.BACKWARD))
454 {
455 // cyl-z propagation failed
456 if(sign_dz*dz>0.) return pstat;
457 }
458 else
459 {
460 if(debug) System.out.println("vec_propcylz: Unknown direction." );
461 System.exit(1);
462 }
463
464 double a34_hat = sign_dz*Math.sqrt(a3*a3 + a4*a4);
465 double a34_hat2 = a34_hat*a34_hat;
466 double dphi = B*dz*a5*Math.sqrt(1.+a34_hat2)*sign_dz;
467 double cos_dphi = Math.cos(dphi);
468 double sin_dphi = Math.sin(dphi);
469
470 double Rcos_phi = 1./(B*a5)*a3/Math.sqrt(1.+a34_hat2)*sign_dz;
471 double Rsin_phi = 1./(B*a5)*a4/Math.sqrt(1.+a34_hat2)*sign_dz;
472
473 double a1n = a1 + Rsin_phi*(cos_dphi - 1) + Rcos_phi*sin_dphi;
474 double a2n = a2 - Rcos_phi*(cos_dphi - 1) + Rsin_phi*sin_dphi;
475 double a3n = a3*cos_dphi - a4*sin_dphi;
476 double a4n = a3*sin_dphi + a4*cos_dphi;
477 double a5n = a5;
478
479 vec.set(IX , a1n);
480 vec.set(IY , a2n);
481 vec.set(IDXDZ , a3n);
482 vec.set(IDYDZ , a4n);
483 vec.set(IQP_Z , a5n);
484
485 // Update trv
486 trv.setSurface(srf.newPureSurface());
487 trv.setVector(vec);
488
489 // set new direction of the track
490
491 if(sign_dz == 1) trv.setForward();
492 if(sign_dz == -1) trv.setBackward();
493
494 // Set the return status.
495
496 if (dir.equals(PropDir.FORWARD))
497 {
498 pstat.setForward();
499 }
500 else if (dir.equals(PropDir.BACKWARD))
501 {
502 pstat.setBackward();
503 }
504
505 // exit now if user did not ask for error matrix.
506 if ( deriv == null ) return pstat;
507
508 //da34_hat
509
510 double da34_hat_da3 = 0.;
511 double da34_hat_da4 = 0.;
512 if(Math.abs(a3)>=Math.abs(a4))
513 {
514 int sign3=1;
515 if(a3<0) sign3 = -1;
516 if(a3 !=0.)
517 {
518 da34_hat_da3 = sign_dz*sign3/Math.sqrt(1.+(a4/a3)*(a4/a3) );
519 da34_hat_da4 = sign_dz*(a4/Math.abs(a3))/Math.sqrt(1.+(a4/a3)*(a4/a3) );
520 }
521 else
522 {
523 da34_hat_da3 = 0.;
524 da34_hat_da4 = 0.;
525 }
526 }
527 if(Math.abs(a4)>Math.abs(a3))
528 {
529 int sign4=1;
530 if(a4<0) sign4 = -1;
531 if(a4 !=0.)
532 {
533 da34_hat_da3 = sign_dz*(a3/Math.abs(a4))/Math.sqrt(1.+(a3/a4)*(a3/a4) );
534 da34_hat_da4 = sign_dz*sign4/Math.sqrt(1.+(a3/a4)*(a3/a4) );
535 }
536 else
537 {
538 da34_hat_da3 = 0.;
539 da34_hat_da4 = 0.;
540 }
541 }
542
543 // ddphi / da
544
545 double ddphi_da3 = B*dz*a5*a34_hat*sign_dz/Math.sqrt(1.+a34_hat2)*da34_hat_da3;
546 double ddphi_da4 = B*dz*a5*a34_hat*sign_dz/Math.sqrt(1.+a34_hat2)*da34_hat_da4;
547 double ddphi_da5 = B*dz*Math.sqrt(a34_hat2+1.)*sign_dz;
548
549 // dRsin_phi
550
551 double dRsin_phi_da3 = -Rsin_phi*a34_hat/(1.+a34_hat2)*da34_hat_da3;
552 double dRsin_phi_da4 = -Rsin_phi*a34_hat/(1.+a34_hat2)*da34_hat_da4 +
553 sign_dz/(B*a5)/Math.sqrt(1.+a34_hat2);
554 double dRsin_phi_da5 = -Rsin_phi/a5;
555
556 // dRcos_phi
557
558 double dRcos_phi_da3 = -Rcos_phi*a34_hat/(1.+a34_hat2)*da34_hat_da3 +
559 sign_dz/(B*a5)/Math.sqrt(1.+a34_hat2);
560 double dRcos_phi_da4 = -Rcos_phi*a34_hat/(1.+a34_hat2)*da34_hat_da4;
561 double dRcos_phi_da5 = -Rcos_phi/a5;
562
563 // da1n first two are simple because dR,dphi _da1,_da2 = 0.
564
565 // double da1n_da1 = 1.;
566 double da1n_da3 = dRsin_phi_da3*(cos_dphi-1.) + dRcos_phi_da3*sin_dphi
567 - Rsin_phi*sin_dphi*ddphi_da3 + Rcos_phi*cos_dphi*ddphi_da3;
568 double da1n_da4 = dRsin_phi_da4*(cos_dphi-1.) + dRcos_phi_da4*sin_dphi
569 - Rsin_phi*sin_dphi*ddphi_da4 + Rcos_phi*cos_dphi*ddphi_da4;
570 double da1n_da5 = dRsin_phi_da5*(cos_dphi-1.) + dRcos_phi_da5*sin_dphi
571 - Rsin_phi*sin_dphi*ddphi_da5 + Rcos_phi*cos_dphi*ddphi_da5;
572
573 // da2n first two are simple because dR,dphi _da1,_da2 = 0.
574
575 // double da2n_da2 = 1.;
576 double da2n_da3 = -dRcos_phi_da3*(cos_dphi-1.) + dRsin_phi_da3*sin_dphi
577 + Rcos_phi*sin_dphi*ddphi_da3 + Rsin_phi*cos_dphi*ddphi_da3;
578 double da2n_da4 = -dRcos_phi_da4*(cos_dphi-1.) + dRsin_phi_da4*sin_dphi
579 + Rcos_phi*sin_dphi*ddphi_da4 + Rsin_phi*cos_dphi*ddphi_da4;
580 double da2n_da5 = -dRcos_phi_da5*(cos_dphi-1.) + dRsin_phi_da5*sin_dphi
581 + Rcos_phi*sin_dphi*ddphi_da5 + Rsin_phi*cos_dphi*ddphi_da5;
582
583 // da3n first two are simple because dphi _da1,_da2 = 0.
584
585 double da3n_da3 = cos_dphi - a3*sin_dphi*ddphi_da3 - a4*cos_dphi*ddphi_da3;
586 double da3n_da4 = - sin_dphi - a3*sin_dphi*ddphi_da4 - a4*cos_dphi*ddphi_da4;
587 double da3n_da5 = - a3*sin_dphi*ddphi_da5 - a4*cos_dphi*ddphi_da5;
588
589 // da4n first two are simple because dphi _da1,_da2 = 0.
590
591 double da4n_da3 = sin_dphi + a3*cos_dphi*ddphi_da3 - a4*sin_dphi*ddphi_da3;
592 double da4n_da4 = cos_dphi + a3*cos_dphi*ddphi_da4 - a4*sin_dphi*ddphi_da4;
593 double da4n_da5 = a3*cos_dphi*ddphi_da5 - a4*sin_dphi*ddphi_da5;
594
595 // da5n
596
597 // double da5n_da5 = 1.;
598
599 // ddphi / dc
600
601 double ddphi_dc2 = -B*a5*Math.sqrt(a34_hat2+1.)*sign_dz;
602
603 // da1 / dc
604
605 double da1_dc1 = -Rcyl*sin_c1;
606
607 // da2 / dc
608
609 double da2_dc1 = Rcyl*cos_c1;
610
611 // da3 / dc
612
613 double da3_dc1 = -sin_dir/c4;
614 double da3_dc3 = -sin_dir/c4;
615 double da3_dc4 = -cos_dir/(c4*c4);
616
617 // da4 / dc
618
619 double da4_dc1 = cos_dir/c4;
620 double da4_dc3 = cos_dir/c4;
621 double da4_dc4 = -sin_dir/(c4*c4);
622
623 // da5 / dc
624
625 double da5_dc4 = -c5*c4/(c4_hat*c4_hat2);
626 double da5_dc5 = 1./c4_hat;
627
628 // da1n/ dc
629
630 double da1n_dc1 = da1_dc1 + da1n_da3*da3_dc1 + da1n_da4*da4_dc1;
631 double da1n_dc2 = -Rsin_phi*sin_dphi*ddphi_dc2 + Rcos_phi*cos_dphi*ddphi_dc2;
632 double da1n_dc3 = da1n_da3*da3_dc3 + da1n_da4*da4_dc3;
633 double da1n_dc4 = da1n_da3*da3_dc4 + da1n_da4*da4_dc4 + da1n_da5*da5_dc4;
634 double da1n_dc5 = da1n_da5*da5_dc5;
635
636 // da2n/ dc
637
638 double da2n_dc1 = da2_dc1 + da2n_da3*da3_dc1 + da2n_da4*da4_dc1;
639 double da2n_dc2 = Rcos_phi*sin_dphi*ddphi_dc2 + Rsin_phi*cos_dphi*ddphi_dc2;
640 double da2n_dc3 = da2n_da3*da3_dc3 + da2n_da4*da4_dc3;
641 double da2n_dc4 = da2n_da3*da3_dc4 + da2n_da4*da4_dc4 + da2n_da5*da5_dc4;
642 double da2n_dc5 = da2n_da5*da5_dc5;
643
644 // da3n/ dc
645
646 double da3n_dc1 = da3n_da3*da3_dc1 + da3n_da4*da4_dc1;
647 double da3n_dc2 = -a3*sin_dphi*ddphi_dc2 - a4*cos_dphi*ddphi_dc2;
648 double da3n_dc3 = da3n_da3*da3_dc3 + da3n_da4*da4_dc3;
649 double da3n_dc4 = da3n_da3*da3_dc4 + da3n_da4*da4_dc4 + da3n_da5*da5_dc4;
650 double da3n_dc5 = da3n_da5*da5_dc5;
651
652 // da4n/ dc
653
654 double da4n_dc1 = da4n_da3*da3_dc1 + da4n_da4*da4_dc1;
655 double da4n_dc2 = a3*cos_dphi*ddphi_dc2 - a4*sin_dphi*ddphi_dc2;
656 double da4n_dc3 = da4n_da3*da3_dc3 + da4n_da4*da4_dc3;
657 double da4n_dc4 = da4n_da3*da3_dc4 + da4n_da4*da4_dc4 + da4n_da5*da5_dc4;
658 double da4n_dc5 = da4n_da5*da5_dc5;
659
660 // da5n/ dc
661
662 double da5n_dc1 = 0.;
663 double da5n_dc2 = 0.;
664 double da5n_dc3 = 0.;
665 double da5n_dc4 = da5_dc4;
666 double da5n_dc5 = da5_dc5;
667
668 deriv.set(IX,IPHI , da1n_dc1);
669 deriv.set(IX,IZ_CYL , da1n_dc2);
670 deriv.set(IX,IALF , da1n_dc3);
671 deriv.set(IX,ITLM , da1n_dc4);
672 deriv.set(IX,IQPT , da1n_dc5);
673 deriv.set(IY,IPHI , da2n_dc1);
674 deriv.set(IY,IZ_CYL , da2n_dc2);
675 deriv.set(IY,IALF , da2n_dc3);
676 deriv.set(IY,ITLM , da2n_dc4);
677 deriv.set(IY,IQPT , da2n_dc5);
678 deriv.set(IDXDZ,IPHI , da3n_dc1);
679 deriv.set(IDXDZ,IZ_CYL , da3n_dc2);
680 deriv.set(IDXDZ,IALF , da3n_dc3);
681 deriv.set(IDXDZ,ITLM , da3n_dc4);
682 deriv.set(IDXDZ,IQPT , da3n_dc5);
683 deriv.set(IDYDZ,IPHI , da4n_dc1);
684 deriv.set(IDYDZ,IZ_CYL , da4n_dc2);
685 deriv.set(IDYDZ,IALF , da4n_dc3);
686 deriv.set(IDYDZ,ITLM , da4n_dc4);
687 deriv.set(IDYDZ,IQPT , da4n_dc5);
688 deriv.set(IQP_Z,IPHI , da5n_dc1);
689 deriv.set(IQP_Z,IZ_CYL , da5n_dc2);
690 deriv.set(IQP_Z,IALF , da5n_dc3);
691 deriv.set(IQP_Z,ITLM , da5n_dc4);
692 deriv.set(IQP_Z,IQPT , da5n_dc5);
693
694 return pstat;
695 }
696
697 private static void check_zero_propagation( VTrack trv0, VTrack trv, PropStat pstat)
698 {
699
700 SpacePoint sp = trv.spacePoint();
701 SpacePoint sp0 = trv0.spacePoint();
702
703 SpacePath sv = trv.spacePath();
704 SpacePath sv0 = trv0.spacePath();
705
706 Assert.assertTrue( Math.abs(sv0.dx() - sv.dx())<1e-7 );
707 Assert.assertTrue( Math.abs(sv0.dy() - sv.dy())<1e-7 );
708 Assert.assertTrue( Math.abs(sv0.dz() - sv.dz())<1e-7 );
709
710 double x0 = sp0.x();
711 double y0 = sp0.y();
712 double z0 = sp0.z();
713 double x1 = sp.x();
714 double y1 = sp.y();
715 double z1 = sp.z();
716
717 double dx = sv.dx();
718 double dy = sv.dy();
719 double dz = sv.dz();
720
721 double prod = dx*(x1-x0)+dy*(y1-y0)+dz*(z1-z0);
722 double moda = Math.sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0) + (z1-z0)*(z1-z0));
723 double modb = Math.sqrt(dx*dx+dy*dy+dz*dz);
724 double st = pstat.pathDistance();
725 Assert.assertTrue( Math.abs(prod-st) < 1.e-7 );
726 Assert.assertTrue( Math.abs(Math.abs(prod) - moda*modb) < 1.e-7 );
727 }
728
729 private static void check_derivatives( Propagator prop, VTrack trv0, Surface srf)
730 {
731 for(int i=0;i<4;++i)
732 for(int j=0;j<4;++j)
733 check_derivative(prop,trv0,srf,i,j);
734 }
735
736 private static void check_derivative( Propagator prop, VTrack trv0, Surface srf,int i,int j)
737 {
738
739 double dx = 1.e-3;
740 VTrack trv = new VTrack(trv0);
741 TrackVector vec = trv.vector();
742 boolean forward = trv0.isForward();
743
744 VTrack trv_0 = new VTrack(trv0);
745 TrackDerivative der = new TrackDerivative();
746 PropStat pstat = prop.vecProp(trv_0,srf,der);
747 Assert.assertTrue(pstat.success());
748
749 TrackVector tmp = new TrackVector(vec);
750 tmp.set(j, tmp.get(j)+dx);
751 trv.setVector(tmp);
752 if(forward) trv.setForward();
753 else trv.setBackward();
754
755 VTrack trv_pl = new VTrack(trv);
756 pstat = prop.vecProp(trv_pl,srf);
757 Assert.assertTrue(pstat.success());
758
759 TrackVector vecpl = trv_pl.vector();
760
761 tmp = new TrackVector(vec);
762 tmp.set(j, tmp.get(j)-dx);
763 trv.setVector(tmp);
764 if(forward) trv.setForward();
765 else trv.setBackward();
766
767 VTrack trv_mn = new VTrack(trv);
768 pstat = prop.vecProp(trv_mn,srf);
769 Assert.assertTrue(pstat.success());
770
771 TrackVector vecmn = trv_mn.vector();
772
773 double dy = (vecpl.get(i)-vecmn.get(i))/2.;
774
775 double dydx = dy/dx;
776
777 double didj = der.get(i,j);
778
779 if( Math.abs(didj) > 1e-10 )
780 Assert.assertTrue( Math.abs((dydx - didj)/didj) < 1e-4 );
781 else
782 Assert.assertTrue( Math.abs(dydx) < 1e-4 );
783 }
784
785 }