1 /*
2 * PropCylXY_Test.java
3 *
4 * Created on July 24, 2007, 10:53 PM
5 *
6 * $Id: PropCylXY_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 PropCylXY_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 IZ = 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_CYL = 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 /** Creates a new instance of PropCylXY_Test */
85 public void testPropCylXY()
86 {
87 String ok_prefix = "PropCylXY (I): ";
88 String error_prefix = "PropCylXY test (E): ";
89
90 if(debug) System.out.println( ok_prefix
91 + "-------- Testing component PropCylXY. --------" );
92
93
94 if(debug) System.out.println( ok_prefix + "Test constructor." );
95 double BFIELD = 2.0;
96 PropCylXY prop = new PropCylXY(BFIELD);
97 if(debug) System.out.println( prop );
98
99 PropCylXY_Test tst = new PropCylXY_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 PropCylXY propcylxy = new PropCylXY(BFIELD/TRFMath.BFAC);
110
111 double phi[] ={ Math.PI/5, Math.PI/6, Math.PI/6, 4/3*Math.PI, 5/3*Math.PI, 5/3*Math.PI };
112 double z[] ={ 1.5, -2.3, 0., 1.5, -1.5, 1.5 };
113 double alpha[] ={ 0.1, -0.1, 0., 0.2, -0.2, 0. };
114 double tlm[] ={ 2.3, -1.5, -2.3, 2.3, -2.3, 0. };
115 double qpt[] ={ 0.01, -0.01, 0.01, -0.01, -0.01, 0.01 };
116
117 double u2b[] ={ 6., 6., 6., 6., 6., 6. };
118 double phic2b[] ={ Math.PI/6., Math.PI/5., Math.PI/7., 5/3*Math.PI, 4/3*Math.PI, 7/4*Math.PI };
119 double u2[] ={ 15., 15., 15., 15., 15., 15. };
120 double phic2[] ={ Math.PI/6., Math.PI/5., Math.PI/7., 5/3*Math.PI, 4/3*Math.PI, 7/4*Math.PI };
121
122 double maxdiff = 1.e-10;
123 double diff;
124 int ntrk = 6;
125 int i;
126
127 for ( i=0; i<ntrk; ++i )
128 {
129 if(debug) System.out.println( "********** Propagate track " + i + ". **********" );
130 PropStat pstat = new PropStat();
131 SurfCylinder scy1 = new SurfCylinder(10.);
132 SurfXYPlane sxyp2 = new SurfXYPlane(u2[i],phic2[i]);
133 SurfXYPlane sxyp2b= new SurfXYPlane(u2b[i],phic2b[i]);
134
135 TrackVector vec1 = new TrackVector();
136
137 vec1.set(IPHI , phi[i]); // phi
138 vec1.set(IZ_CYL , z[i]); // z
139 vec1.set(IALF , alpha[i]); // alpha
140 vec1.set(ITLM , tlm[i]); // tan(lambda)
141 vec1.set(IQPT , qpt[i]); // q/pt
142
143 VTrack trv1 = new VTrack(scy1.newPureSurface(),vec1);
144
145 if(debug) System.out.println( "\n starting: " + trv1 );
146
147
148 VTrack trv2f = new VTrack(trv1);
149 pstat = propcylxy.vecDirProp(trv2f,sxyp2,PropDir.FORWARD);
150 Assert.assertTrue( pstat.forward() );
151 if(debug) System.out.println( "\n forward: " + trv2f );
152
153 VTrack trv2f_my = new VTrack(trv1);
154 pstat = tst.vec_propcylxy(BFIELD,trv2f_my,sxyp2,PropDir.FORWARD);
155 Assert.assertTrue( pstat.forward() );
156 if(debug) System.out.println( "\n forward my: " + trv2f_my );
157 diff = tst.compare(trv2f_my,trv2f);
158
159 if(debug) System.out.println( "\n diff: " + diff );
160 Assert.assertTrue( diff < maxdiff );
161
162
163 VTrack trv2b = new VTrack(trv1);
164 pstat = propcylxy.vecDirProp(trv2b,sxyp2b,PropDir.BACKWARD_MOVE);
165 Assert.assertTrue( pstat.backward() );
166 if(debug) System.out.println( "\n backward: " + trv2b );
167
168 VTrack trv2b_my = new VTrack(trv1);
169 pstat = tst.vec_propcylxy(BFIELD,trv2b_my,sxyp2b,PropDir.BACKWARD);
170 Assert.assertTrue( pstat.backward() );
171 if(debug) System.out.println( "\n backward my: " + trv2b_my );
172 diff = tst.compare(trv2b_my,trv2b);
173
174 if(debug) System.out.println( "\n diff: " + diff );
175 Assert.assertTrue( diff < maxdiff );
176
177
178 }
179 //********************************************************************
180
181 // Repeat the above with errors.
182 if(debug) System.out.println( ok_prefix + "Check against correct propagation with errors." );
183
184 double epp[] = { 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01 };
185 double epz[] = { 0.01, -0.01, 0.01, -0.01, 0.01, -0.01, 0.01 };
186 double ezz[] = { 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25 };
187 double epa[] = { 0.003, -0.003, 0.003, -0.003, 0.003, -0.003, 0.003 };
188 double eza[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004 };
189 double eaa[] = { 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01 };
190 double epl[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004 };
191 double eal[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004 };
192 double ezl[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004 };
193 double ell[] = { 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02 };
194 double epc[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004 };
195 double ezc[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004 };
196 double eac[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004 };
197 double elc[] = { 0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004 };
198 double ecc[] = { 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01 };
199
200 maxdiff = 1.e-8;
201 double ediff;
202 ntrk = 6;
203
204 for ( i=0; i<ntrk; ++i )
205 {
206 if(debug) System.out.println( "********** Propagate track " + i + ". **********" );
207 PropStat pstat = new PropStat();
208 SurfCylinder scy1 = new SurfCylinder(10.);
209 SurfXYPlane sxyp2 = new SurfXYPlane(u2[i],phic2[i]);
210 SurfXYPlane sxyp2b = new SurfXYPlane(u2b[i],phic2b[i]);
211
212 TrackVector vec1 = new TrackVector();
213
214 vec1.set(IPHI , phi[i]); // phi
215 vec1.set(IZ_CYL , z[i]); // z
216 vec1.set(IALF , alpha[i]); // alpha
217 vec1.set(ITLM , tlm[i]); // tan(lambda)
218 vec1.set(IQPT , qpt[i]); // q/pt
219
220 TrackError err1 = new TrackError();
221
222 err1.set(IPHI,IPHI , epp[i]);
223 err1.set(IPHI,IZ_CYL , epz[i]);
224 err1.set(IZ_CYL,IZ_CYL , ezz[i]);
225 err1.set(IPHI,IALF , epa[i]);
226 err1.set(IZ_CYL,IALF , eza[i]);
227 err1.set(IALF,IALF , eaa[i]);
228 err1.set(IPHI,ITLM , epl[i]);
229 err1.set(IZ_CYL,ITLM , ezl[i]);
230 err1.set(IALF,ITLM , eal[i]);
231 err1.set(ITLM,ITLM , ell[i]);
232 err1.set(IPHI,IQPT , epc[i]);
233 err1.set(IZ_CYL,IQPT , ezc[i]);
234 err1.set(IALF,IQPT , eac[i]);
235 err1.set(ITLM,IQPT , elc[i]);
236 err1.set(IQPT,IQPT , ecc[i]);
237
238 ETrack trv1 = new ETrack(scy1.newPureSurface(),vec1,err1);
239
240 if(debug) System.out.println( "\n starting: " + trv1 );
241
242 ETrack trv2f = new ETrack(trv1);
243 pstat = propcylxy.errDirProp(trv2f,sxyp2,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_propcylxy(BFIELD,trv2f_my,sxyp2,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 = propcylxy.errDirProp(trv2b,sxyp2b,PropDir.BACKWARD);
263 Assert.assertTrue( pstat.backward() );
264 if(debug) System.out.println( "\n backward: " + trv2b );
265
266 ETrack trv2b_my = new ETrack(trv1);
267 pstat = tst.vec_propcylxy(BFIELD,trv2b_my,sxyp2b,PropDir.BACKWARD,deriv);
268 Assert.assertTrue( pstat.backward() );
269 err = trv2b_my.error();
270 trv2b_my.setError( err.Xform(deriv) );
271 if(debug) System.out.println( "\n backward my: " + trv2b_my );
272 diffs = tst.compare(trv2b_my,trv2b);
273
274 if(debug) System.out.println( "\n diff: " + diffs[0] + ' ' + "ediff: "+ diffs[1] );
275 Assert.assertTrue( diffs[0] < maxdiff );
276 Assert.assertTrue( diffs[1] < maxdiff );
277
278
279 }
280 //********************************************************************
281
282 if(debug) System.out.println( ok_prefix + "Test cloning." );
283 Assert.assertTrue( prop.newPropagator() != null );
284
285 //********************************************************************
286
287 if(debug) System.out.println( ok_prefix + "Test the field." );
288 Assert.assertTrue( prop.bField() == 2.0 );
289
290 //********************************************************************
291
292 if(debug) System.out.println( ok_prefix + "Test Zero Field Propagation." );
293 {
294 PropCylXY prop0 = new PropCylXY(0.0);
295 if(debug) System.out.println( prop0 );
296 Assert.assertTrue( prop0.bField() == 0. );
297
298 double u = 10.;
299 double phi2 = 0.1;
300 Surface srf = new SurfCylinder(13.0);
301 VTrack trv0 = new VTrack(srf);
302 TrackVector vec = new TrackVector();
303 vec.set(SurfCylinder.IPHI, phi2+0.1);
304 vec.set(SurfCylinder.IZ, 10.);
305 vec.set(SurfCylinder.IALF, 0.1);
306 vec.set(SurfCylinder.ITLM, 2.);
307 vec.set(SurfCylinder.IQPT, 0.);
308 trv0.setVector(vec);
309 trv0.setForward();
310 Surface srf_to = new SurfXYPlane(u,phi2);
311
312 VTrack trv = new VTrack(trv0);
313 VTrack trv_der = new VTrack(trv);
314 PropStat pstat = prop0.vecDirProp(trv,srf_to,PropDir.NEAREST);
315 if(debug) System.out.println("trv= \n"+trv+'\n');
316 Assert.assertTrue( pstat.success() );
317
318 Assert.assertTrue( pstat.backward() );
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 }
324
325 //********************************************************************
326
327 if(debug) System.out.println( ok_prefix
328 + "------------- All tests passed. -------------" );
329
330 //********************************************************************
331 }
332
333 // Very well tested Cyl-XY propagator. Each new one can be tested against it
334 PropStat
335 vec_propcylxy( double B, VTrack trv, Surface srf,
336 PropDir dir
337 )
338 {
339 TrackDerivative deriv = null;
340 return vec_propcylxy(B, trv, srf, dir, deriv);
341 }
342
343
344 PropStat
345 vec_propcylxy( double B, VTrack trv, Surface srf,
346 PropDir dir,
347 TrackDerivative deriv )
348 {
349
350
351 // construct return status
352 PropStat pstat = new PropStat();
353
354 // fetch the originating surface and vector
355 Surface srf1 = trv.surface();
356 // TrackVector vec1 = trv.vector();
357
358 // Check origin is a Cylinder.
359 Assert.assertTrue( srf1.pureType().equals(SurfCylinder.staticType()) );
360 if (! srf1.pureType( ).equals(SurfCylinder.staticType()) )
361 return pstat;
362 SurfCylinder scy1 = ( SurfCylinder ) srf1;
363
364 // Fetch the R of the cylinder and the starting track vector.
365 int ir = SurfCylinder.RADIUS;
366 double Rcyl = scy1.parameter(ir);
367
368 TrackVector vec = trv.vector();
369 double c1 = vec.get(IPHI); // phi
370 double c2 = vec.get(IZ); // z
371 double c3 = vec.get(IALF); // alpha
372 double c4 = vec.get(ITLM); // tan(lambda)
373 double c5 = vec.get(IQPT); // q/pt
374
375
376 // Check destination is a XYPlane.
377 Assert.assertTrue( srf.pureType().equals(SurfXYPlane.staticType()) );
378 if ( !srf.pureType( ).equals(SurfXYPlane.staticType()) )
379 return pstat;
380 SurfXYPlane sxyp2 = ( SurfXYPlane ) srf;
381
382 // Fetch the u and phi of the plane.
383 int iphi = SurfXYPlane.NORMPHI;
384 int idist = SurfXYPlane.DISTNORM;
385
386 double phi_n = sxyp2.parameter(iphi);
387 double u_n = sxyp2.parameter(idist);
388
389 // rotate coordinate system on phi_n
390
391 c1 -= phi_n;
392 if(c1 < 0.) c1 += TRFMath.TWOPI;
393
394 double cos_c1 = Math.cos(c1);
395
396 double u = Rcyl*cos_c1;
397
398 double sin_c1 = Math.sin(c1);
399 double cos_dir = Math.cos(c1+c3);
400 double sin_dir = Math.sin(c1+c3);
401 double c4_hat2 = 1+c4*c4;
402 double c4_hat = Math.sqrt(c4_hat2);
403
404 // check if du == 0 ( that is track moves parallel to the destination plane )
405 // du = pt*cos_dir
406 if(cos_dir/c5 == 0.) return pstat;
407
408 double tan_dir = sin_dir/cos_dir;
409
410 double b1 = Rcyl*sin_c1;
411 double b2 = c2;
412 double b3 = tan_dir;
413 double b4 = c4/cos_dir;
414 double b5 = c5/c4_hat;
415
416 int sign_du = 0;
417 if(cos_dir > 0) sign_du = 1;
418 if(cos_dir < 0) sign_du = -1;
419
420 double b3_hat = Math.sqrt(1 + b3*b3);
421 double b34_hat = Math.sqrt(1 + b3*b3 + b4*b4);
422 double b3_hat2 = b3_hat*b3_hat;
423 double b34_hat2 = b34_hat*b34_hat;
424 double r = 1/(b5*B)*b3_hat/b34_hat;
425 double cosphi = -b3*sign_du/b3_hat;
426 double sinphi = sign_du/b3_hat;
427 double rsinphi = 1./(b5*B)*sign_du/b34_hat;
428 double rcosphi = -b3/(b5*B)*sign_du/b34_hat;
429
430 double du = u_n - u;
431 double norm = du/r - cosphi;
432
433 // cyl-xy propagation failed : noway to the new plane
434 if(Math.abs(norm)>1.) return pstat;
435
436 double cat = Math.sqrt(1.-norm*norm);
437 int sign_dphi = 0;
438
439 if (dir.equals(PropDir.NEAREST))
440 {
441 if(debug) System.out.println("PropCylXY._vec_prop: Propagation in NEAREST direction is undefined");
442 System.exit(1);
443 }
444 else if (dir.equals(PropDir.FORWARD))
445 {
446 if(b5>0.) sign_dphi = 1;
447 if(b5<0.) sign_dphi = -1;
448 }
449 else if (dir.equals(PropDir.BACKWARD))
450 {
451 if(b5>0.) sign_dphi = -1;
452 if(b5<0.) sign_dphi = 1;
453 }
454 else
455 {
456 if(debug) System.out.println("PropCylXY._vec_prop: Unknown direction." );
457 System.exit(1);
458 }
459
460 int sign_cat = 0;
461 if(du*sign_dphi*b5>0.) sign_cat = 1;
462 if(du*sign_dphi*b5<0.) sign_cat = -1;
463 if(du == 0.)
464 {
465 if(sinphi >=0. ) sign_cat = 1;
466 if(sinphi < 0. ) sign_cat = -1;
467 }
468
469 double sin_dphi = norm*sinphi + cat*cosphi*sign_cat;
470 double cos_dphi = -norm*cosphi + cat*sinphi*sign_cat;
471
472 int sign_sindphi = 0;
473 if(sin_dphi> 0.) sign_sindphi = 1;
474 if(sin_dphi< 0.) sign_sindphi = -1;
475 if(sin_dphi == 0.) sign_sindphi = sign_dphi;
476
477 double dphi = Math.PI*(sign_dphi - sign_sindphi) + sign_sindphi*Math.acos(cos_dphi);
478
479 // check if dun == 0 ( that is track moves parallel to the destination plane)
480 double du_n_du = cos_dphi - b3*sin_dphi;
481 if(du_n_du==0.) return pstat;
482
483 double a_hat_dphi = 1./du_n_du;
484 double a_hat_dphi2 = a_hat_dphi*a_hat_dphi;
485 double c_hat_dphi = sin_dphi + b3*cos_dphi ;
486
487 double b1_n = b1 + rsinphi*(1-cos_dphi) - rcosphi*sin_dphi;
488 double b2_n = b2 + b4/(b5*B)*sign_du/b34_hat*dphi;
489 double b3_n = c_hat_dphi*a_hat_dphi;
490 double b4_n = b4*a_hat_dphi;
491 double b5_n = b5;
492
493 int sign_dun = 0;
494 if(du_n_du*sign_du > 0) sign_dun = 1;
495 if(du_n_du*sign_du < 0) sign_dun = -1;
496
497 vec.set(IV , b1_n);
498 vec.set(IZ , b2_n);
499 vec.set(IDVDU , b3_n);
500 vec.set(IDZDU , b4_n);
501 vec.set(IQP_XY , b5_n);
502
503 // Set the return status.
504
505 if (dir.equals(PropDir.FORWARD))
506 {
507 pstat.setForward();
508 }
509 else if (dir.equals(PropDir.BACKWARD))
510 {
511 pstat.setBackward();
512 }
513
514 // Update trv
515 trv.setSurface(srf.newPureSurface());
516 trv.setVector(vec);
517 // set new direction of the track
518 if(sign_dun == 1) trv.setForward();
519 if(sign_dun == -1) trv.setBackward();
520
521 // exit now if user did not ask for error matrix.
522 if (deriv == null ) return pstat;
523
524 // du_dc
525
526 double du_dc1 = -Rcyl*sin_c1;
527
528 // db1_dc
529
530 double db1_dc1 = Rcyl*cos_c1;
531
532 // db2_dc
533
534 double db2_dc2 = 1.;
535
536 // db3_dc
537
538 double db3_dc1 = 1./(cos_dir*cos_dir);
539 double db3_dc3 = 1./(cos_dir*cos_dir);
540
541 // db4_dc
542
543 double db4_dc1 = b4*tan_dir;
544 double db4_dc3 = b4*tan_dir;
545 double db4_dc4 = 1/cos_dir;
546
547 // db5_dc
548
549 double db5_dc4 = -c4*c5/(c4_hat*c4_hat2);
550 double db5_dc5 = 1./c4_hat;
551
552
553 // dr_db
554
555 double dr_db3 = r*b3*b4*b4/(b3_hat2*b34_hat2);
556 double dr_db4 = -r*b4/b34_hat2;
557 double dr_db5 = -r/b5;
558
559 // dcosphi_db
560
561 double dcosphi_db3 = - sign_du/b3_hat - cosphi*b3/b3_hat2;
562
563 // dsinphi_db
564
565 double dsinphi_db3 = - sinphi*b3/b3_hat2;
566
567 // dcat_db
568
569 double dcat_db3 = norm/cat*(du/(r*r)*dr_db3 + dcosphi_db3 );
570 double dcat_db4 = norm/cat* du/(r*r)*dr_db4;
571 double dcat_db5 = norm/cat* du/(r*r)*dr_db5;
572 double dcat_du = norm/(cat*r);
573
574 // dnorm_db
575
576 double dnorm_db3 = - du/(r*r)*dr_db3 - dcosphi_db3;
577 double dnorm_db4 = - du/(r*r)*dr_db4;
578 double dnorm_db5 = - du/(r*r)*dr_db5;
579 double dnorm_du = - 1./r;
580
581 // dcos_dphi_db
582
583 double dcos_dphi_db3 = - cosphi*dnorm_db3 - norm*dcosphi_db3 +
584 sign_cat*(sinphi*dcat_db3 + cat*dsinphi_db3);
585 double dcos_dphi_db4 = - cosphi*dnorm_db4 + sign_cat*sinphi*dcat_db4;
586 double dcos_dphi_db5 = - cosphi*dnorm_db5 + sign_cat*sinphi*dcat_db5;
587 double dcos_dphi_du = - cosphi*dnorm_du + sign_cat*sinphi*dcat_du;
588
589 // dsin_dphi_db
590
591 double dsin_dphi_db3 = sinphi*dnorm_db3 + norm*dsinphi_db3 +
592 sign_cat*(cosphi*dcat_db3 + cat*dcosphi_db3);
593 double dsin_dphi_db4 = sinphi*dnorm_db4 + sign_cat*cosphi*dcat_db4;
594 double dsin_dphi_db5 = sinphi*dnorm_db5 + sign_cat*cosphi*dcat_db5;
595 double dsin_dphi_du = sinphi*dnorm_du + sign_cat*cosphi*dcat_du;
596
597 // ddphi_db
598
599 double ddphi_db3;
600 double ddphi_db4;
601 double ddphi_db5;
602 double ddphi_du;
603 if(Math.abs(sin_dphi)>0.5)
604 {
605 ddphi_db3 = - dcos_dphi_db3/sin_dphi;
606 ddphi_db4 = - dcos_dphi_db4/sin_dphi;
607 ddphi_db5 = - dcos_dphi_db5/sin_dphi;
608 ddphi_du = - dcos_dphi_du /sin_dphi;
609 }
610 else
611 {
612 ddphi_db3 = dsin_dphi_db3/cos_dphi;
613 ddphi_db4 = dsin_dphi_db4/cos_dphi;
614 ddphi_db5 = dsin_dphi_db5/cos_dphi;
615 ddphi_du = dsin_dphi_du /cos_dphi;
616 }
617
618 // da_hat_dphi_db
619
620 double da_hat_dphi_db3 = - a_hat_dphi2*
621 (dcos_dphi_db3 - sin_dphi - b3*dsin_dphi_db3);
622 double da_hat_dphi_db4 = - a_hat_dphi2*(dcos_dphi_db4 - b3*dsin_dphi_db4);
623 double da_hat_dphi_db5 = - a_hat_dphi2*(dcos_dphi_db5 - b3*dsin_dphi_db5);
624 double da_hat_dphi_du = - a_hat_dphi2*(dcos_dphi_du - b3*dsin_dphi_du );
625
626 // dc_hat_dphi_db
627
628 double dc_hat_dphi_db3 = b3*dcos_dphi_db3 + dsin_dphi_db3 + cos_dphi;
629 double dc_hat_dphi_db4 = b3*dcos_dphi_db4 + dsin_dphi_db4;
630 double dc_hat_dphi_db5 = b3*dcos_dphi_db5 + dsin_dphi_db5;
631 double dc_hat_dphi_du = b3*dcos_dphi_du + dsin_dphi_du ;
632
633 // db1_n_db
634
635 double db1_n_db1 = 1;
636 double db1_n_db3 = (dr_db3*sinphi+r*dsinphi_db3)*(1-cos_dphi)
637 - rsinphi*dcos_dphi_db3
638 - dr_db3*cosphi*sin_dphi - r*dcosphi_db3*sin_dphi
639 - rcosphi*dsin_dphi_db3;
640 double db1_n_db4 = dr_db4*sinphi*(1-cos_dphi) - rsinphi*dcos_dphi_db4
641 - dr_db4*cosphi*sin_dphi - rcosphi*dsin_dphi_db4;
642 double db1_n_db5 = dr_db5*sinphi*(1-cos_dphi) - rsinphi*dcos_dphi_db5
643 - dr_db5*cosphi*sin_dphi - rcosphi*dsin_dphi_db5;
644 double db1_n_du = - rsinphi*dcos_dphi_du - rcosphi*dsin_dphi_du;
645
646 // db2_n_db
647
648 double db2_n_db2 = 1.;
649 double db2_n_db3 = 1./(b5*B)*b4*sign_du/b34_hat*
650 ( - dphi*b3/b34_hat2 + ddphi_db3);
651 double db2_n_db4 = 1./(b5*B)*sign_du/b34_hat*
652 ( - dphi*b4*b4/b34_hat2 + b4*ddphi_db4 + dphi );
653 double db2_n_db5 = 1./(b5*B)*b4*sign_du/b34_hat*( ddphi_db5 - dphi/b5);
654 double db2_n_du = 1./(b5*B)*b4*sign_du/b34_hat*ddphi_du;
655
656 // db3_n_db
657
658 double db3_n_db3 = a_hat_dphi*dc_hat_dphi_db3 + da_hat_dphi_db3*c_hat_dphi;
659 double db3_n_db4 = a_hat_dphi*dc_hat_dphi_db4 + da_hat_dphi_db4*c_hat_dphi;
660 double db3_n_db5 = a_hat_dphi*dc_hat_dphi_db5 + da_hat_dphi_db5*c_hat_dphi;
661 double db3_n_du = a_hat_dphi*dc_hat_dphi_du + da_hat_dphi_du *c_hat_dphi;
662
663 // db4_n_db
664
665 double db4_n_db3 = b4*da_hat_dphi_db3;
666 double db4_n_db4 = b4*da_hat_dphi_db4 + a_hat_dphi;
667 double db4_n_db5 = b4*da_hat_dphi_db5;
668 double db4_n_du = b4*da_hat_dphi_du;
669
670 // db5_n_db
671
672 // double db5_n_db5 = 1.;
673
674 // db1_n_dc
675
676 double db1_n_dc1 = db1_n_du * du_dc1 + db1_n_db1*db1_dc1 +
677 db1_n_db3*db3_dc1 + db1_n_db4*db4_dc1;
678 double db1_n_dc2 = 0.;
679 double db1_n_dc3 = db1_n_db3*db3_dc3 + db1_n_db4*db4_dc3;
680 double db1_n_dc4 = db1_n_db4*db4_dc4 + db1_n_db5*db5_dc4;
681 double db1_n_dc5 = db1_n_db5*db5_dc5;
682
683
684 // db2_n_dc
685
686 double db2_n_dc1 = db2_n_du * du_dc1 + db2_n_db3*db3_dc1 + db2_n_db4*db4_dc1;
687 double db2_n_dc2 = db2_n_db2 * db2_dc2;
688 double db2_n_dc3 = db2_n_db3*db3_dc3 + db2_n_db4*db4_dc3;
689 double db2_n_dc4 = db2_n_db4*db4_dc4 + db2_n_db5*db5_dc4;
690 double db2_n_dc5 = db2_n_db5*db5_dc5;
691
692 // db3_n_dc
693
694 double db3_n_dc1 = db3_n_du * du_dc1 + db3_n_db3*db3_dc1 + db3_n_db4*db4_dc1;
695 double db3_n_dc2 = 0.;
696 double db3_n_dc3 = db3_n_db3*db3_dc3 + db3_n_db4*db4_dc3;
697 double db3_n_dc4 = db3_n_db4*db4_dc4 + db3_n_db5*db5_dc4;
698 double db3_n_dc5 = db3_n_db5*db5_dc5;
699
700 // db4_n_dc
701
702 double db4_n_dc1 = db4_n_du * du_dc1 + db4_n_db3*db3_dc1 + db4_n_db4*db4_dc1;
703 double db4_n_dc2 = 0.;
704 double db4_n_dc3 = db4_n_db3*db3_dc3 + db4_n_db4*db4_dc3;
705 double db4_n_dc4 = db4_n_db4*db4_dc4 + db4_n_db5*db5_dc4;
706 double db4_n_dc5 = db4_n_db5*db5_dc5;
707
708 // db5_n_dc
709
710 double db5_n_dc1 = 0.;
711 double db5_n_dc2 = 0.;
712 double db5_n_dc3 = 0.;
713 double db5_n_dc4 = db5_dc4;
714 double db5_n_dc5 = db5_dc5;
715
716 deriv.set(IV,IPHI , db1_n_dc1);
717 deriv.set(IV,IZ_CYL , db1_n_dc2);
718 deriv.set(IV,IALF , db1_n_dc3);
719 deriv.set(IV,ITLM , db1_n_dc4);
720 deriv.set(IV,IQPT , db1_n_dc5);
721 deriv.set(IZ,IPHI , db2_n_dc1);
722 deriv.set(IZ,IZ_CYL , db2_n_dc2);
723 deriv.set(IZ,IALF , db2_n_dc3);
724 deriv.set(IZ,ITLM , db2_n_dc4);
725 deriv.set(IZ,IQPT , db2_n_dc5);
726 deriv.set(IDVDU,IPHI , db3_n_dc1);
727 deriv.set(IDVDU,IZ_CYL , db3_n_dc2);
728 deriv.set(IDVDU,IALF , db3_n_dc3);
729 deriv.set(IDVDU,ITLM , db3_n_dc4);
730 deriv.set(IDVDU,IQPT , db3_n_dc5);
731 deriv.set(IDZDU,IPHI , db4_n_dc1);
732 deriv.set(IDZDU,IZ_CYL , db4_n_dc2);
733 deriv.set(IDZDU,IALF , db4_n_dc3);
734 deriv.set(IDZDU,ITLM , db4_n_dc4);
735 deriv.set(IDZDU,IQPT , db4_n_dc5);
736 deriv.set(IQP_XY,IPHI , db5_n_dc1);
737 deriv.set(IQP_XY,IZ_CYL , db5_n_dc2);
738 deriv.set(IQP_XY,IALF , db5_n_dc3);
739 deriv.set(IQP_XY,ITLM , db5_n_dc4);
740 deriv.set(IQP_XY,IQPT , db5_n_dc5);
741
742 return pstat;
743 }
744
745 static void check_zero_propagation( VTrack trv0, VTrack trv, PropStat pstat)
746 {
747
748 SpacePoint sp = trv.spacePoint();
749 SpacePoint sp0 = trv0.spacePoint();
750
751 SpacePath sv = trv.spacePath();
752 SpacePath sv0 = trv0.spacePath();
753
754 Assert.assertTrue( Math.abs(sv0.dx() - sv.dx())<1e-7 );
755 Assert.assertTrue( Math.abs(sv0.dy() - sv.dy())<1e-7 );
756 Assert.assertTrue( Math.abs(sv0.dz() - sv.dz())<1e-7 );
757
758 double x0 = sp0.x();
759 double y0 = sp0.y();
760 double z0 = sp0.z();
761 double x1 = sp.x();
762 double y1 = sp.y();
763 double z1 = sp.z();
764
765 double dx = sv.dx();
766 double dy = sv.dy();
767 double dz = sv.dz();
768
769 double prod = dx*(x1-x0)+dy*(y1-y0)+dz*(z1-z0);
770 double moda = Math.sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0) + (z1-z0)*(z1-z0));
771 double modb = Math.sqrt(dx*dx+dy*dy+dz*dz);
772 double st = pstat.pathDistance();
773 // Assert.assertTrue( Math.abs(prod-st) < 1.e-7 );
774 // Assert.assertTrue( Math.abs(Math.abs(prod) - moda*modb) < 1.e-7 );
775 }
776
777 static void check_derivatives( Propagator prop, VTrack trv0, Surface srf)
778 {
779 for(int i=0;i<4;++i)
780 for(int j=0;j<4;++j)
781 check_derivative(prop,trv0,srf,i,j);
782 }
783
784 static void check_derivative( Propagator prop, VTrack trv0, Surface srf,int i,int j)
785 {
786
787 double dx = 1.e-3;
788 VTrack trv = new VTrack(trv0);
789 TrackVector vec = trv.vector();
790 boolean forward = trv0.isForward();
791
792 VTrack trv_0 = new VTrack(trv0);
793 TrackDerivative der = new TrackDerivative();
794 PropStat pstat = prop.vecProp(trv_0,srf,der);
795 Assert.assertTrue(pstat.success());
796
797 TrackVector tmp = new TrackVector(vec);
798 tmp.set(j, tmp.get(j)+dx);
799 trv.setVector(tmp);
800 if(forward) trv.setForward();
801 else trv.setBackward();
802
803 VTrack trv_pl = new VTrack(trv);
804 pstat = prop.vecProp(trv_pl,srf);
805 Assert.assertTrue(pstat.success());
806
807 TrackVector vecpl = trv_pl.vector();
808
809 tmp = new TrackVector(vec);
810 tmp.set(j,tmp.get(j)-dx);
811 trv.setVector(tmp);
812 if(forward) trv.setForward();
813 else trv.setBackward();
814
815 VTrack trv_mn = new VTrack(trv);
816 pstat = prop.vecProp(trv_mn,srf);
817 Assert.assertTrue(pstat.success());
818
819 TrackVector vecmn = trv_mn.vector();
820
821 double dy = (vecpl.get(i)-vecmn.get(i))/2.;
822
823 double dydx = dy/dx;
824
825 double didj = der.get(i,j);
826
827 if( Math.abs(didj) > 1e-10 )
828 Assert.assertTrue( Math.abs((dydx - didj)/didj) < 1e-4 );
829 else
830 Assert.assertTrue( Math.abs(dydx) < 1e-4 );
831 }
832 }