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 }