View Javadoc

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 }