1 /* 2 * PropDCACyl_Test.java 3 * 4 * Created on July 24, 2007, 9:57 PM 5 * 6 * $Id: PropDCACyl_Test.java,v 1.1.1.1 2010/04/08 20:38:00 jeremy Exp $ 7 */ 8 9 package org.lcsim.recon.tracking.trfdca; 10 11 import junit.framework.TestCase; 12 import org.lcsim.recon.tracking.spacegeom.SpacePoint; 13 import org.lcsim.recon.tracking.trfbase.ETrack; 14 import org.lcsim.recon.tracking.trfbase.PropDir; 15 import org.lcsim.recon.tracking.trfbase.PropStat; 16 import org.lcsim.recon.tracking.trfbase.Propagator; 17 import org.lcsim.recon.tracking.trfbase.Surface; 18 import org.lcsim.recon.tracking.trfbase.TrackDerivative; 19 import org.lcsim.recon.tracking.trfbase.TrackError; 20 import org.lcsim.recon.tracking.trfbase.TrackVector; 21 import org.lcsim.recon.tracking.trfbase.VTrack; 22 import org.lcsim.recon.tracking.trfcyl.SurfCylinder; 23 import org.lcsim.recon.tracking.trfutil.Assert; 24 import org.lcsim.recon.tracking.trfutil.TRFMath; 25 26 /** 27 * 28 * @author Norman Graf 29 */ 30 public class PropDCACyl_Test extends TestCase 31 { 32 private boolean debug; 33 /** Creates a new instance of PropDCACyl_Test */ 34 public void testPropDCACyl() 35 { 36 String ok_prefix = "PropDCACyl (I): "; 37 String error_prefix = "PropDCACyl test (E): "; 38 39 if(debug) System.out.println(ok_prefix 40 + "-------- Testing component PropDCACyl. --------" ); 41 42 //******************************************************************** 43 44 if(debug) System.out.println(ok_prefix + "Test constructor." ); 45 double bfield = 2.0; 46 PropDCACyl prop = new PropDCACyl(bfield); 47 if(debug) System.out.println(prop ); 48 49 //******************************************************************** 50 51 if(debug) System.out.println(ok_prefix + "Test cloning." ); 52 Assert.assertTrue( prop.newPropagator() != null); 53 54 //******************************************************************** 55 56 if(debug) System.out.println(ok_prefix 57 + "Test the magnetic field." ); 58 Assert.assertTrue( prop.bField() == bfield ); 59 60 //******************************************************************** 61 62 if(debug) System.out.println(ok_prefix 63 + "Test propagation of a track - without error." ); 64 65 PropStat pstat = new PropStat(); 66 67 TrackVector tvec = new TrackVector(); 68 tvec.set(0, 10.0); // r_signed 69 tvec.set(1, 2.0); // z 70 tvec.set(2, 3.0); // phi_direction 71 tvec.set(3, 4.0); // tlm 72 tvec.set(4, 5.0); // qpt 73 SurfDCA sdca = new SurfDCA(); 74 VTrack trv = new VTrack( sdca.newPureSurface(),tvec ); 75 SpacePoint spt = trv.spacePoint(); 76 77 if(debug) System.out.println(" *** Propagation to Cylinder of same r *** " ); 78 SurfCylinder scys = new SurfCylinder(10.0); 79 VTrack trvsf = new VTrack(trv); 80 pstat = prop.vecDirProp(trvsf,scys,PropDir.FORWARD); 81 Assert.assertTrue( pstat.success() ); 82 SpacePoint sptsf = trvsf.spacePoint(); 83 TrackVector tvec_scys = trvsf.vector(); 84 if(debug) System.out.println(" sptsf = " ); 85 if(debug) System.out.println( sptsf ); 86 if(debug) System.out.println(" tvec_scys = " + tvec_scys ); 87 Assert.assertTrue( sptsf.rxy() == tvec.get(SurfDCA.IRSIGNED) ); 88 Assert.assertTrue( tvec_scys.get(SurfCylinder.IZ) == tvec.get(SurfDCA.IZ) ); 89 Assert.assertTrue( tvec_scys.get(SurfCylinder.IPHI) == spt.phi() ); 90 Assert.assertTrue( tvec_scys.get(SurfCylinder.ITLM) == tvec.get(SurfDCA.ITLM) ); 91 Assert.assertTrue( tvec_scys.get(SurfCylinder.IQPT) == tvec.get(SurfDCA.IQPT) ); 92 93 SurfCylinder scyl = new SurfCylinder(15.); 94 95 if(debug) System.out.println(" *** Forward propagation *** " ); 96 VTrack trv2f = new VTrack(trv); 97 if(debug) System.out.println(" before propagation: trv2f = " + trv2f ); 98 SpacePoint spt2f = trv2f.spacePoint(); 99 if(debug) System.out.println(" before propagation: spt2f = " ); 100 if(debug) System.out.println(spt2f ); 101 Assert.assertTrue( spt2f.rxy() == 10.0 ); 102 Assert.assertTrue( Math.abs(spt2f.phi()-spt.phi()) < 1.e-4 ); 103 Assert.assertTrue( spt2f.z() == 2.0 ); 104 pstat = prop.vecDirProp(trv2f,scyl,PropDir.FORWARD); 105 Assert.assertTrue( pstat.forward() ); 106 if(debug) System.out.println(" after propagation: trv2f = " + trv2f ); 107 spt2f = trv2f.spacePoint(); 108 if(debug) System.out.println(" after propagation: spt2f = " ); 109 if(debug) System.out.println(spt2f ); 110 Assert.assertTrue( spt2f.rxy() == 15.0 ); 111 112 if(debug) System.out.println(" *** Backward propagation *** " ); 113 VTrack trv2b = new VTrack(trv); 114 if(debug) System.out.println(" before propagation: trv2b = " + trv2b ); 115 SpacePoint spt2b = trv2b.spacePoint(); 116 if(debug) System.out.println(" before propagation: spt2b = " ); 117 if(debug) System.out.println(spt2b ); 118 Assert.assertTrue( spt2b.rxy() == 10.0 ); 119 Assert.assertTrue( Math.abs(spt2b.phi()-spt.phi()) < 1.e-4 ); 120 Assert.assertTrue( spt2b.z() == 2.0 ); 121 pstat = prop.vecDirProp(trv2b,scyl,PropDir.BACKWARD); 122 Assert.assertTrue( pstat.backward() ); 123 if(debug) System.out.println(" after propagation: trv2b = " + trv2b ); 124 spt2b = trv2b.spacePoint(); 125 if(debug) System.out.println(" after propagation: spt2b = " ); 126 if(debug) System.out.println(spt2b ); 127 Assert.assertTrue( spt2b.rxy() == 15.0 ); 128 129 130 //******************************************************************** 131 132 if(debug) System.out.println(ok_prefix 133 + "Test propagation of a track - with error." ); 134 135 // PropStat pstat = new PropStat(); 136 137 // TrackVector tvec = new TrackVector(); 138 // tvec.set(0, 10.0); // r_signed 139 // tvec.set(1, 2.0); // z 140 // tvec.set(2, 3.0); // phi_direction 141 // tvec.set(3, 4.0); // tlm 142 // tvec.set(4, 5.0); // qpt 143 TrackError err = new TrackError(); 144 err.set(0,0, 2.0); 145 err.set(0,1, 1.0); 146 err.set(0,2, 1.0); 147 err.set(0,3, 1.0); 148 err.set(0,4, 1.0); 149 err.set(1,1, 3.0); 150 err.set(1,2, 1.0); 151 err.set(1,3, 1.0); 152 err.set(1,4, 1.0); 153 err.set(2,2, 4.0); 154 err.set(2,3, 1.0); 155 err.set(2,4, 1.0); 156 err.set(3,3, 5.0); 157 err.set(3,4, 1.0); 158 err.set(4,4, 6.0); 159 // SurfDCA sdca; 160 ETrack etrv = new ETrack( sdca.newPureSurface(),tvec,err ); 161 SpacePoint espt = etrv.spacePoint(); 162 163 // SurfCylinder scyl = new SurfCylinder(15.); 164 165 if(debug) System.out.println(" *** Forward propagation *** " ); 166 ETrack etrv2f = new ETrack(etrv); 167 if(debug) System.out.println(" before propagation: etrv2f = " + etrv2f ); 168 SpacePoint espt2f = etrv2f.spacePoint(); 169 if(debug) System.out.println(" before propagation: espt2f = " ); 170 if(debug) System.out.println(espt2f ); 171 Assert.assertTrue( espt2f.rxy() == 10.0 ); 172 Assert.assertTrue( Math.abs(espt2f.phi()-espt.phi()) < 1.e-4 ); 173 Assert.assertTrue( espt2f.z() == 2.0 ); 174 if(debug) System.out.println("scyl= "+scyl); 175 if(debug) System.out.println("\n* etrv2f before: "+etrv2f); 176 pstat = prop.errDirProp(etrv2f,scyl,PropDir.FORWARD); 177 if(debug) System.out.println("\n* etrv2f after: "+etrv2f); 178 if(debug) System.out.println("pstat after propagation= "+pstat); 179 Assert.assertTrue( pstat.forward() ); 180 if(debug) System.out.println(" after propagation: etrv2f = " + etrv2f ); 181 espt2f = etrv2f.spacePoint(); 182 if(debug) System.out.println(" after propagation: espt2f = " ); 183 if(debug) System.out.println(espt2f ); 184 Assert.assertTrue( espt2f.rxy() == 15.0 ); 185 //******************************************************************** 186 187 // Test non 0.0 DCA 188 if(debug) System.out.println("\n\n ***Testing non zero DCA \n\n"); 189 int NUM = 10; 190 double[] rad = 191 {20.,2.0,Math.sqrt(5.)}; 192 // double[] prec = { 1e-5,1e-5, 1e-5}; 193 double[] xv = 194 { 0.5 , 0., 2. , 2., 2., 2., 2. ,2., 2. ,2.}; 195 double[] yv = 196 { 0.4 , 0., -1., -1., -1., -1., -1. ,-1., -1. ,-1.}; 197 double[] r = 198 { 1. , 0., 0.1, -0.1, -0.1, 0.1 , 0. , 0., -0.1 , 0.1}; 199 double[] phid = 200 { 0.1, 0., 0.2, 0.2, 0.2 , 0.2 , 0. , 0., 0. , 0.}; 201 double[] z = 202 { 2., 2., 2. , 2., 2., 2., 2. , 3., 2. , 3.}; 203 double[] tlm = 204 { -4., 4., 4. , 4., 4., -4., 4. , 2., 4. , 2.}; 205 double[] qpt = 206 { 0.1, 0., 0.01, 0.01,-0.01,-0.01,-0.01 , 0.01,0.0 ,0.}; 207 for( int j=0; j <3 ;++j) 208 { 209 for(int i=0 ; i< NUM; ++i ) 210 { 211 if(debug) System.out.println("\n\n\n\n "+j+" "+i); 212 SurfDCA srf = new SurfDCA(xv[i],yv[i]); 213 TrackVector vec = new TrackVector(); 214 TrackError errnzdca = new TrackError(); 215 for(int k=1;k<=5;++k) errnzdca.set(k-1,k-1,k); 216 vec.set(SurfDCA.IRSIGNED, r[i]); 217 vec.set(SurfDCA.IZ, z[i]); 218 vec.set(SurfDCA.IPHID, phid[i]); 219 vec.set(SurfDCA.ITLM, tlm[i]); 220 vec.set(SurfDCA.IQPT, qpt[i]); 221 ETrack tre = new ETrack( srf.newPureSurface()); 222 tre.setVector(vec); 223 tre.setError(errnzdca); 224 tre.setForward(); 225 SurfCylinder cyl = new SurfCylinder(rad[j]); 226 PropDCACyl propdca_cyl = new PropDCACyl( -2.); 227 PropCylDCA propcyl_dca = new PropCylDCA(-2.); 228 ETrack tre0 = new ETrack(tre); 229 ETrack treb = new ETrack(tre); 230 PropStat pstatnzdca = propdca_cyl.errDirProp(tre,cyl,PropDir.FORWARD); 231 Assert.assertTrue( pstatnzdca.success()); 232 pstatnzdca = propdca_cyl.errDirProp(treb,cyl,PropDir.BACKWARD); 233 Assert.assertTrue( pstatnzdca.success()); 234 check_derivatives(propdca_cyl,tre0,cyl); 235 ETrack tre1 = new ETrack(tre); 236 ETrack tre1b = new ETrack(treb); 237 pstatnzdca = propcyl_dca.errProp(tre,srf); 238 Assert.assertTrue( pstatnzdca.success()); 239 pstatnzdca = propcyl_dca.errProp(treb,srf); 240 Assert.assertTrue( pstatnzdca.success()); 241 check_derivatives(propcyl_dca,tre1,srf); 242 if(debug) System.out.println("checking equality \n"+tre0+" \n"+tre); 243 if(debug) System.out.println("about_equal: "+about_equal(tre0,tre)); 244 //Assert.assertTrue(about_equal(tre0,tre)); 245 check_derivatives(propcyl_dca,tre1b,srf); 246 // Assert.assertTrue(about_equal(tre0,treb)); 247 248 } 249 } 250 251 //******************************************************************** 252 if(debug) System.out.println(ok_prefix 253 + "------------- All tests passed. -------------" ); 254 255 //******************************************************************** 256 257 } 258 public static boolean about_equal(ETrack tre1, ETrack tre2) 259 { 260 boolean equal = true; 261 double maxdif = 1.e-6; 262 // Check vector. 263 TrackVector vdiff = tre1.surface().vecDiff( tre1.vector(), tre2.vector() ); 264 for (int i=0; i<5; ++i) 265 { 266 double adif =Math.abs( vdiff.get(i) ); 267 if ( ! (adif < maxdif) ) 268 { 269 System.out.println("vecdif(" + i + ") = " + adif ); 270 equal = false; 271 } 272 } 273 // Check errors. 274 TrackError ediff = tre1.error().minus(tre2.error()); 275 TrackError eavg = tre1.error().plus(tre2.error()); 276 for (int i=0; i<5; ++i) 277 { 278 for (int j=0; j<=i; ++j) 279 { 280 double adif = Math.abs( ediff.get(i,j) ); 281 double afrac = ediff.get(i,j)/Math.sqrt(eavg.get(i,i)*eavg.get(j,j)); 282 if ( !(adif < maxdif) ) 283 { 284 System.out.println("errdif(" + i + "," + j + ") = " + adif + " frac = " + afrac ); 285 equal = false; 286 } 287 } 288 } 289 return equal; 290 } 291 292 public static void check_derivatives( Propagator prop, VTrack trv0, Surface srf) 293 { 294 check_derivatives( prop, trv0, srf, 1.e-5); 295 } 296 297 public static void check_derivatives( Propagator prop, VTrack trv0, Surface srf,double prec) 298 { 299 for(int i=0;i<4;++i) 300 { 301 for(int j=0;j<4;++j) 302 { 303 check_derivative(prop,trv0,srf,i,j,prec); 304 } 305 } 306 } 307 308 309 public static void check_derivative( Propagator prop, VTrack trv0, Surface srf,int i,int j,double prec) 310 { 311 312 double dx = 1.e-3; 313 VTrack trv = new VTrack(trv0); 314 TrackVector vec = trv.vector(); 315 boolean forward = trv0.isForward(); 316 317 VTrack trv_0 = new VTrack(trv0); 318 TrackDerivative der = new TrackDerivative(); 319 PropStat pstat = prop.vecProp(trv_0,srf,der); 320 Assert.assertTrue(pstat.success()); 321 322 TrackVector tmp=vec; 323 double tmpValue = tmp.get(j); 324 tmp.set(j, tmpValue+dx); 325 trv.setVector(tmp); 326 if(forward) trv.setForward(); 327 else trv.setBackward(); 328 329 VTrack trv_pl = new VTrack(trv); 330 pstat = prop.vecProp(trv_pl,srf); 331 Assert.assertTrue(pstat.success()); 332 333 TrackVector vecpl = trv_pl.vector(); 334 335 tmp=vec; 336 tmpValue = tmp.get(j); 337 tmp.set(j, tmpValue-dx); 338 trv.setVector(tmp); 339 if(forward) trv.setForward(); 340 else trv.setBackward(); 341 342 VTrack trv_mn = trv; 343 pstat = prop.vecProp(trv_mn,srf); 344 Assert.assertTrue(pstat.success()); 345 346 TrackVector vecmn = trv_mn.vector(); 347 348 double dy = (vecpl.get(i)-vecmn.get(i))/2.; 349 350 double dy1 = vecpl.get(i)-trv_0.vector(i); 351 352 if( !TRFMath.isZero(dy) ) 353 { 354 if( Math.abs((dy1-dy)/dy) > 0.1 ) 355 return; 356 } 357 else 358 { 359 if( Math.abs(dy1) > 1e-5 ) 360 return; 361 } 362 if( Math.abs(Math.abs(dy1)-Math.abs(dy)) > 0.01 ) 363 dy=dy1; 364 365 double dydx = dy/dx; 366 367 double didj = der.get(i,j); 368 369 if( Math.abs(didj) > 1e-10 ) 370 Assert.assertTrue( Math.abs((dydx - didj)/didj) < prec ); 371 else 372 Assert.assertTrue( Math.abs(dydx) < prec ); 373 } 374 }