View Javadoc

1   /*
2    * To change this template, choose Tools | Templates
3    * and open the template in the editor.
4    */
5   package org.lcsim.recon.tracking.trfzp;
6   
7   import junit.framework.TestCase;
8   import org.lcsim.recon.tracking.magfield.ConstantMagneticField;
9   import org.lcsim.recon.tracking.trfbase.*;
10  
11  import static java.lang.Math.abs;
12  import static java.lang.Math.max;
13  
14  /**
15   *
16   * @author ngraf
17   */
18  public class PropZZRKTest extends TestCase
19  {
20  
21      private boolean debug = false;
22  
23      public void testPropZZRK()
24      {
25          ConstantMagneticField field = new ConstantMagneticField(0., 0., 2.);
26          PropZZRK prop = new PropZZRK(field);
27          if(debug) System.out.println(prop);
28  
29          // Construct equivalent PropZZ propagator.
30          PropZZ prop0 = new PropZZ(2.0);
31          if (debug) {
32              System.out.println(prop0);
33          }
34  
35          // Here we propagate some tracks both forward and backward and then
36          // each back to the original track.  We check that the returned
37          // track parameters match those of the original track.
38          if(debug) System.out.println("Check reversibility.");
39          double z1[] = {100.0, -100.0, 100.0, -100.0, 100.0, -100.0, 100.0, 100.0};
40          double z2[] = {150.0, -50.0, 50.0, -150.0, 50.0, -50.0, 150.0, 50.0};
41          double x[] = {10.0, 1.0, -1.0, 2.0, 2.0, -2.0, 3.0, 0.0};
42          double xv[] = {0.5, 0.03, -0.03, 0.5, -0.5, 1.0, 1.0, 2.0};
43          double crv[] = {0.1, -0.1, 0.1, 0.01, 0.01, -1.0, 1.0, 0.02};
44          double y[] = {20.0, 0.0, 0.0, 0.0, 0.0, 0.0, 15.0, 0.0};
45          double yv[] = {-0.5, 0.01, -0.02, 0.0, 0.0, 0.5, -0.5, 0.0};
46          double fbdf[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
47          double bfdf[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
48          boolean forw[] = {true, false, false, true, false, true, true, true};
49          double maxdiff = 1.e-7;
50          int ntrk = 8;
51          int i;
52  
53  
54          for (i = 0; i < ntrk; ++i) {
55              if(debug) System.out.println("********** Propagate track " + i + ". **********");
56              PropStat pstat = new PropStat();
57              SurfZPlane sz1 = new SurfZPlane(z1[i]);
58              SurfZPlane sz2 = new SurfZPlane(z2[i]);
59              TrackVector vec1 = new TrackVector();
60              vec1.set(0, x[i]);    // x
61              vec1.set(1, y[i]);    // y
62              vec1.set(2, xv[i]);   // dx/dz
63              vec1.set(3, yv[i]);   // dy/dz
64              vec1.set(4, crv[i]);  // q/p
65              TrackSurfaceDirection tdir;
66              if (forw[i]) {
67                  tdir = TrackSurfaceDirection.TSD_FORWARD;
68              } else {
69                  tdir = TrackSurfaceDirection.TSD_BACKWARD;
70              }
71              VTrack trv1 = new VTrack(sz1.newPureSurface(), vec1, tdir);
72              if(debug) System.out.println(" starting: " + trv1);
73              //
74              // Find the direction that will propagate this track from z1 to z2.
75              //
76              PropDir dir = PropDir.FORWARD;
77              PropDir rdir = PropDir.BACKWARD;
78              if (z2[i] > z1[i] && tdir == TrackSurfaceDirection.TSD_BACKWARD
79                      || z2[i] < z1[i] && tdir == TrackSurfaceDirection.TSD_FORWARD) {
80                  dir = PropDir.BACKWARD;
81                  rdir = PropDir.FORWARD;
82              }
83  
84              //
85              // Propagate.
86              VTrack trv2f = trv1;
87              pstat = prop.vecDirProp(trv2f, sz2, dir);
88              assert (pstat.success());
89              if(debug) System.out.println("  forward: " + trv2f);
90              if(debug) System.out.println(pstat);
91              if (dir == PropDir.FORWARD) {
92                  assert (pstat.forward());
93              }
94              if (dir == PropDir.BACKWARD) {
95                  assert (pstat.backward());
96              }
97  
98              //
99              // Propagate using PropZZ and check difference.
100             VTrack trv2f0 = trv1;
101             pstat = prop0.vecDirProp(trv2f0, sz2, dir);
102             assert (pstat.success());
103             if(debug) System.out.println("  forward: " + trv2f0);
104             if(debug) System.out.println(pstat);
105             double diff0 =
106                     sz2.vecDiff(trv2f.vector(), trv2f0.vector()).amax();
107             if(debug) System.out.println("diff: " + diff0);
108             assert (diff0 < maxdiff);
109 
110             //
111             // Propagate in reverse direction.
112             VTrack trv2fb = trv2f;
113             pstat = prop.vecDirProp(trv2fb, sz1, rdir);
114             assert (pstat.success());
115             if(debug) System.out.println(" f return: " + trv2fb);
116             if(debug) System.out.println(pstat);
117             if (rdir == PropDir.FORWARD) {
118                 assert (pstat.forward());
119             }
120             if (rdir == PropDir.BACKWARD) {
121                 assert (pstat.backward());
122             }
123             // Check the return differences.
124             double difff =
125                     sz1.vecDiff(trv2fb.vector(), trv1.vector()).amax();
126             if(debug) System.out.println("diff: " + difff);
127             assert (difff < maxdiff);
128             //
129             // Check no-move forward propagation to the same surface.
130             VTrack trv1s = trv1;
131             pstat = prop.vecDirProp(trv1s, sz1, PropDir.FORWARD);
132             assert (pstat.success());
133             if(debug) System.out.println(" same f: " + trv1s);
134             if(debug) System.out.println(pstat);
135             assert (pstat.same());
136             assert (pstat.pathDistance() == 0);
137             //
138             // Check no-move backward propagation to the same surface.
139             trv1s = trv1;
140             pstat = prop.vecDirProp(trv1s, sz1, PropDir.BACKWARD);
141             assert (pstat.success());
142             if(debug) System.out.println(" same b: " + trv1s);
143             if(debug) System.out.println(pstat);
144             assert (pstat.same());
145             assert (pstat.pathDistance() == 0);
146             //
147             // Check move propagation to the same surface.
148             //
149             // forward
150             int successes = 0;
151             trv1s = trv1;
152             pstat = prop.vecDirProp(trv1s, sz1, PropDir.FORWARD_MOVE);
153             if(debug) System.out.println(" forward move: " + trv1s);
154             if(debug) System.out.println(pstat);
155             if (pstat.forward()) {
156                 ++successes;
157             }
158             // backward
159             trv1s = trv1;
160             pstat = prop.vecDirProp(trv1s, sz1, PropDir.BACKWARD_MOVE);
161             if(debug) System.out.println(" backward move: " + trv1s);
162             if(debug) System.out.println(pstat);
163             if (pstat.backward()) {
164                 ++successes;
165             }
166             // Neither of these should have succeeded.
167             assert (successes == 0);
168             //
169             // nearest
170             trv1s = trv1;
171             pstat = prop.vecDirProp(trv1s, sz1, PropDir.NEAREST_MOVE);
172             if(debug) System.out.println(" nearest move: " + trv1s);
173             if(debug) System.out.println(pstat);
174             assert (!pstat.success());
175 
176 
177 //               cng 120905 problems here... 
178 //            // Test derivatives numerically using uniform field. 
179 //            System.out.println("Testing derivatives with uniform field");
180 //            VTrack trv1a = trv1;
181 //            TrackDerivative deriv = new TrackDerivative();
182 //            pstat = prop.vecDirProp(trv1a, sz2, PropDir.NEAREST, deriv);
183 //            assert (pstat.success());
184 //            VTrack trv1a0 = trv1;
185 //            TrackDerivative deriv0 = new TrackDerivative();
186 //            pstat = prop0.vecDirProp(trv1a0, sz2, PropDir.NEAREST, deriv0);
187 //            assert (pstat.success());
188 //            for (int j = 0; j < 5; ++j) {
189 //                double d;
190 //                if (j < 2) {
191 //                    d = 0.25;
192 //                } else if (j < 4) {
193 //                    d = 0.01;
194 //                } else {
195 //                    d = 0.001;
196 //                }
197 //                TrackVector vec1b = vec1;
198 //                TrackVector vec1c = vec1;
199 //                vec1b.set(j, vec1.get(j) + d);
200 //                vec1c.set(j, vec1.get(j) - d);
201 //                VTrack trv1b = new VTrack(sz1.newPureSurface(), vec1b, tdir);
202 //                VTrack trv1c = new VTrack(sz1.newPureSurface(), vec1c, tdir);
203 //                pstat = prop.vecDirProp(trv1b, sz2, PropDir.NEAREST);
204 //                assert (pstat.success());
205 //                pstat = prop.vecDirProp(trv1c, sz2, PropDir.NEAREST);
206 //                assert (pstat.success());
207 //                System.out.println("Testing diffs");
208 //                System.out.println("ii, j \t dij deriv(ii,j) deriv0(ii,j)");
209 //                for (int ii = 0; ii < 5; ++ii) {
210 //                    double dij = (trv1b.vector(ii) - trv1c.vector(ii)) / (2. * d);
211 //                    System.out.println(ii + ", " + j + '\t' + dij + '\t' + deriv.get(ii, j) + '\t' + deriv0.get(ii, j));
212 //                    double scale = max(abs(deriv.get(ii, j)), 1.);
213 //                    assert (abs(deriv.get(ii, j) - deriv0.get(ii, j)) / scale < maxdiff);
214 //                    double tmp = abs(dij - deriv.get(ii, j)) / scale;
215 //                    System.out.println(tmp + " " + scale);
216 //                    assert (abs(dij - deriv.get(ii, j)) / scale < 1.e-3);
217 //                }
218 //            }
219 
220 // following was commented out originally
221 //    /*
222 //    // Test derivatives numerically using non-uniform field.
223 //    System.out.println("Testing derivatives with non-uniform field");
224 //    VTrack trv1an = trv1;
225 //    TrackDerivative derivn;
226 //    pstat = propmc.vecDirProp(trv1an,sz2,PropDir.NEAREST, &derivn);
227 //    assert(pstat.success());
228 //    for(int j=0; j<5; ++j) {
229 //      double d;
230 //      if(j < 2)
231 //	d = 0.1;
232 //      else if(j < 4)
233 //	d = 0.01;
234 //      else
235 //	d = 0.001;
236 //      TrackVector vec1bn = vec1;
237 //      TrackVector vec1cn = vec1;
238 //      vec1bn(j, vec1(j) + d;
239 //      vec1cn(j, vec1(j) - d;
240 //      VTrack trv1bn(SurfacePtr(sz1.newPureSurface()), vec1bn, tdir);
241 //      VTrack trv1cn(SurfacePtr(sz1.newPureSurface()), vec1cn, tdir);
242 //      pstat = propmc.vecDirProp(trv1bn,sz2,PropDir.NEAREST);
243 //      assert(pstat.success());
244 //      pstat = propmc.vecDirProp(trv1cn,sz2,PropDir.NEAREST);
245 //      assert(pstat.success());
246 //      for(int i=0; i<5; ++i) {
247 //	double dijn = (trv1bn.vector(i) - trv1cn.vector(i))/(2.*d);
248 //	cout + i + ", " + j + '\t' +  dijn + '\t' + derivn(i,j) );
249 //	double scale = max(abs(derivn(i,j)), 1.);
250 //	assert(abs(dijn-derivn(i,j))/scale < 0.01);
251 //      }
252 //    }
253 //    */
254         }
255 
256         //********************************************************************
257 
258         // Repeat the above with errors.
259         if(debug) System.out.println("Check reversibility with errors.");
260         double exx[] = {0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1, 0.1};
261         double exy[] = {0.05, -0.05, 0.05, -0.05, 0.05, -0.05, 0.05, 0.05};
262         double eyy[] = {0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25};
263         double exxv[] = {0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004, 0.004};
264         double eyxv[] = {0.02, -0.02, 0.02, -0.02, 0.02, -0.02, 0.02, 0.02};
265         double exvxv[] = {0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01};
266         double exyv[] = {0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004, 0.004};
267         double eyyv[] = {0.04, -0.04, 0.04, -0.04, 0.04, -0.04, 0.04, 0.04};
268         double exvyv[] = {0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004, 0.004};
269         double eyvyv[] = {0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02};
270         double exc[] = {0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004, 0.004};
271         double eyc[] = {0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004, 0.004};
272         double exvc[] = {0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004, 0.004};
273         double eyvc[] = {0.004, -0.004, 0.004, -0.004, 0.004, -0.004, 0.004, 0.004};
274         double ecc[] = {0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01};
275         for (i = 0; i < ntrk; ++i) {
276             if(debug) System.out.println("********** Propagate track " + i + ". **********");
277             PropStat pstat = new PropStat();
278             SurfZPlane dz1 = new SurfZPlane(z1[i]);
279             SurfZPlane dz2 = new SurfZPlane(z2[i]);
280             TrackVector vec1 = new TrackVector();
281             vec1.set(0, x[i]);    // x
282             vec1.set(1, y[i]);    // y
283             vec1.set(2, xv[i]);   // dx/dz
284             vec1.set(3, yv[i]);   // dy/dz
285             vec1.set(4, crv[i]);  // curvature
286             TrackError err1 = new TrackError();
287             err1.set(0, 0, exx[i]);
288             err1.set(0, 1, exy[i]);
289             err1.set(1, 1, eyy[i]);
290             err1.set(0, 2, exxv[i]);
291             err1.set(1, 2, eyxv[i]);
292             err1.set(2, 2, exvxv[i]);
293             err1.set(0, 3, exyv[i]);
294             err1.set(1, 3, eyyv[i]);
295             err1.set(2, 3, exvyv[i]);
296             err1.set(3, 3, eyvyv[i]);
297             err1.set(0, 4, exc[i]);
298             err1.set(1, 4, eyc[i]);
299             err1.set(2, 4, exvc[i]);
300             err1.set(3, 4, eyvc[i]);
301             err1.set(4, 4, ecc[i]);
302             TrackSurfaceDirection tdir;
303             if (forw[i]) {
304                 tdir = TrackSurfaceDirection.TSD_FORWARD;
305             } else {
306                 tdir = TrackSurfaceDirection.TSD_BACKWARD;
307             }
308             ETrack tre1 = new ETrack((dz1.newPureSurface()), vec1, err1, tdir);
309             if(debug) System.out.println(" starting: " + tre1);
310             //
311             // Find the direction that will propagate this track from r1 to r2.
312             //
313             PropDir dir = PropDir.FORWARD;
314             PropDir rdir = PropDir.BACKWARD;
315             if (z2[i] > z1[i] && tdir == TrackSurfaceDirection.TSD_BACKWARD
316                     || z2[i] < z1[i] && tdir == TrackSurfaceDirection.TSD_FORWARD) {
317                 dir = PropDir.BACKWARD;
318                 rdir = PropDir.FORWARD;
319             }
320 
321             //
322             // Propagate.
323             ETrack tre2f = tre1;
324             pstat = prop.errDirProp(tre2f, dz2, dir);
325             assert (pstat.success());
326             if(debug) System.out.println("  forward: " + tre2f);
327             if(debug) System.out.println(pstat);
328             if (dir == PropDir.FORWARD) {
329                 assert (pstat.forward());
330             }
331             if (dir == PropDir.BACKWARD) {
332                 assert (pstat.backward());
333             }
334 
335             //
336             // Propagate using PropZZ and check difference.
337             ETrack tre2f0 = tre1;
338             pstat = prop0.errDirProp(tre2f0, dz2, dir);
339             assert (pstat.success());
340             if(debug) System.out.println("  forward: " + tre2f0);
341             if(debug) System.out.println(pstat);
342             double vdiff0 =
343                     dz2.vecDiff(tre2f.vector(), tre2f0.vector()).amax();
344             if(debug) System.out.println("vec diff: " + vdiff0);
345             assert (vdiff0 < maxdiff);
346             TrackError df0 = tre2f.error().minus(tre2f0.error());
347             double ediff0 = df0.amax();
348             if(debug) System.out.println("err diff: " + ediff0);
349             assert (ediff0 < maxdiff);
350 
351             //
352             // Propagate backward
353             ETrack tre2fb = tre2f;
354             pstat = prop.errDirProp(tre2fb, dz1, rdir);
355             assert (pstat.success());
356             if(debug) System.out.println(" f return: " + tre2fb);
357             if(debug) System.out.println(pstat);
358             if (rdir == PropDir.FORWARD) {
359                 assert (pstat.forward());
360             }
361             if (rdir == PropDir.BACKWARD) {
362                 assert (pstat.backward());
363             }
364             double difff =
365                     dz1.vecDiff(tre2fb.vector(), tre1.vector()).amax();
366             if(debug) System.out.println("vec diff: " + difff);
367             assert (difff < maxdiff);
368             TrackError dfb = tre2fb.error().minus(tre1.error());
369             double edifff = dfb.amax();
370             if(debug) System.out.println("err diff: " + edifff);
371             assert (edifff < maxdiff);
372         }
373 
374     }
375 }