lcsim/test/org/lcsim/fit/circle
diff -N CircleFitterTest.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ CircleFitterTest.java 21 Mar 2006 19:52:25 -0000 1.1
@@ -0,0 +1,313 @@
+package org.lcsim.fit.circle;
+// Circle Fitting Classes based on Karimaki's fast circle fit.
+/*
+$! **************************************************
+$! * CIRCLE FIT SAMPLE PROGRAM *
+$! **************************************************
+$ CREATE CIRCLEFIT.FOR
+ ******************************************************************
+ PROGRAM CIRCLEFIT
+ * ******* ********* 4.6.1991 V. Karimaki *
+ * *
+ * Main program to test the circle fitting routine CIRCLF *
+ * *
+ * Circle data points are Monte Carlo generated and then *
+ * used to test CIRCLF NTRGEN times. The pull values are *
+ * calculated for the three fitted parameters. *
+ * *
+ * The test is passed, if the statistical distributions of *
+ * the pull values follow the normal distribution. The means *
+ * and standard deviations of the three pull values are *
+ * printed in the test run output. *
+ * *
+ * Constants (units are supposed to be in metres): *
+ * AVEDEV = average residual *
+ * STEP = spacing between points along circle *
+ * RHOMAX = maximal (absolute) value of curvature *
+ * DISRMX = maximal ratio |DCA|/track length *
+ *----------------------------------------------------------------*
+ */
+import java.util.Random;
+import junit.framework.TestCase;
+
+public class CircleFitterTest extends TestCase
+{
+
+ public void testCircleFitter()
+ {
+ Random _r = new Random();
+ CircleFitter _fitter = new CircleFitter();
+ int _ncmax = 500;
+ double[] _xx = new double[_ncmax];
+ double[] _yy = new double[_ncmax];
+ double[] _wxy = new double[_ncmax];
+ double[] _xyr = new double[_ncmax];
+ double[] _pulls = new double[3];
+ double[][] _pulsum = new double[2][3];
+ double[] _pullerr = new double[3];
+
+ //
+ //--- Number of generations
+ //
+ int NTRGEN=2000;
+ int NPRINT=10;
+ int NMODUL=NTRGEN/NPRINT;
+
+ double step = 0.01;
+ double AVEDEV = 0.00030;
+ double RHOMAX = 1.0;
+ double DISRMX=0.40;
+
+// System.out.println("TEST RUN OUTPUT");
+// System.out.println("Circle fit routine is tested with "+NTRGEN+
+// " randomly generated sets of circle data.");
+ //
+ //--- Loop to generate circle points and test the routine CIRCLF
+ //
+ //DO 100 ITR=1,NTRGEN
+ for(int ITR=0; ITR<NTRGEN; ++ITR)
+ {
+ //--- Generate random curvature, direction and dca
+
+ double RND1=_r.nextDouble();
+ double RND2=_r.nextDouble();
+ double RND3=_r.nextDouble();
+ double RND4=_r.nextDouble();
+
+ double RHOGEN=2.*RHOMAX*(RND1-0.5);
+ double PHIGEN=2.*Math.PI*RND2;
+ double AROGEN=Math.abs(RHOGEN);
+ //
+ //--- Number of points on circle 4<=NPO<=300
+ //
+ int NPO=4+(int)(296*RND4);
+ double XLENG=(NPO-1)*step;
+ double DISGEN=DISRMX*2.*XLENG*(RND3-0.5);
+ if (AROGEN>0.)
+ {
+ double RADIUS=1./AROGEN;
+ double DISLIM=0.75*RADIUS;
+ if (Math.abs(DISGEN)>DISLIM)
+ {
+ DISGEN = DISLIM;
+ if(DISGEN<0) DISGEN = -DISLIM;
+ }
+ int NPFULL=(int)(6.2832*RADIUS/step-2);
+ NPO=Math.min(NPO,NPFULL);
+ }
+ double XREFE=RND1-0.5;
+ double YREFE=RND2-0.5;
+ //CNG
+ // if I fix these parameters, and calculate pulls
+ // with respect to these fixed parameters, everything
+ // is fine. When I don't everything goes to hell.
+ // Error must be in coding of modifications, or in propagate
+ // method.
+ // May, 18, 2000
+ // I have now fixed the propagate method.
+ // Everything works well if I fix these starting parameters
+ // and propagate to a random XTEST, YTEST
+ /*
+ NPO = 40;
+ RHOGEN = 0.01;
+ PHIGEN = 1.23;
+ DISGEN = 0.1;
+ XREFE = 0.;
+ YREFE = 0.;
+ */
+ //CNG
+ double SINGEN=Math.sin(PHIGEN);
+ double COSGEN=Math.cos(PHIGEN);
+ double XSTART= XREFE+DISGEN*SINGEN;
+ double YSTART= YREFE-DISGEN*COSGEN;
+ //
+ //--- Generate a circle with random parameters
+ //
+ simCirclePoints(RHOGEN,PHIGEN,XSTART,YSTART,NPO,AVEDEV,step, _xx, _yy, _wxy, _xyr);
+ //Check the generation here...
+ //System.out.println("phigen= "+PHIGEN+", rhogen= "+RHOGEN+", dca= "+DISGEN);
+ /*for(int i = 0; i<NPO; ++i)
+ {
+ System.out.println("_xx["+i+"]= "+_xx[i]+", _yy["+i+"]= "+_yy[i]+", _wxy["+i+"]= "+_wxy[i]);
+ }
+ */
+
+ //
+ //--- Perform circle fitting with error estimation
+ //
+ // CALL CIRCLF(XREFE,YREFE,XX,YY,WXY,NPO,1,IER)
+ _fitter.setreferenceposition( XREFE, YREFE);
+ //System.out.println("xref= "+XREFE+", yref= "+YREFE);
+ boolean OK = _fitter.fit( _xx, _yy, _wxy,NPO);
+ assertTrue(OK);
+ //if(OK) System.out.println("Fit OK");
+ CircleFit cf = _fitter.getfit();
+ //System.out.println(cf);
+
+ //
+ //--- Define also another reference point to test propagation
+ double XTEST=XSTART+DISGEN*(RND3-0.5);
+ double YTEST=YSTART+DISGEN*(RND4-0.5);
+ //
+ //--- Calculate the true parameters at the propagation test point
+ //
+ double XMOVE=XREFE-XTEST;
+ double YMOVE=YREFE-YTEST;
+ double ROD1=1.+RHOGEN*DISGEN;
+ double DPERP=XMOVE*SINGEN-YMOVE*COSGEN + DISGEN;
+ double DPARA=XMOVE*COSGEN+YMOVE*SINGEN;
+ double AA=2.*DPERP+RHOGEN*(DPERP*DPERP+DPARA*DPARA);
+ double UU=Math.sqrt(1.+RHOGEN*AA);
+ double SQ1AI=1./(1.+UU);
+ double BB= RHOGEN*XMOVE+ROD1*SINGEN;
+ double CC=-RHOGEN*YMOVE+ROD1*COSGEN;
+ double RHOTES=RHOGEN;
+ double PHITES=Math.atan2(BB,CC);
+ double DISTES=AA*SQ1AI;
+ //
+ //--- Propagate to a test point
+ //
+ _fitter.propagatefit(XTEST, YTEST);
+ cf = _fitter.getfit();
+
+ double FIC=cf.phi();
+ if(PHITES-cf.phi()>Math.PI) FIC+=2.*Math.PI;
+ if(cf.phi()-PHITES>Math.PI) FIC-=2.*Math.PI;
+ double[] covmat = cf.cov();
+ double ROERR=Math.sqrt(covmat[0]);
+ double FIERR=Math.sqrt(covmat[2]);
+ double DCERR=Math.sqrt(covmat[5]);
+ //
+ //--- Print a few fit results
+ //
+ /*
+ IF (MOD(ITR,NMODUL).EQ.1) THEN
+ WRITE(IUNIT,1005) RHO,ROERR,FIC,FIERR,DCA,DCERR,CHICIR,NPO
+ WRITE(IUNIT,1006) RHOTES,PHITES,DISTES
+ 1005 FORMAT(' Fitted',3(2X,F7.4,'+-',F6.4),4X,F5.1,I6)
+ 1006 FORMAT(' True ',3(2X,F7.4,8X))
+ ENDIF
+ */
+ //
+ //--- Calculate the pull values
+ //
+ // use when I have implemented propagate...
+ _pulls[0]=(RHOTES-cf.curvature())/ROERR;
+ _pulls[1]=(PHITES-FIC)/FIERR;
+ _pulls[2]=(DISTES-cf.dca())/DCERR;
+
+ for(int i =0; i<3;++i)
+ {
+ // System.out.println("pulls["+i+"]= "+_pulls[i]);
+ _pulsum[0][i]+=_pulls[i];
+ _pulsum[1][i]+=_pulls[i]*_pulls[i];
+ }
+
+ }
+
+ //
+ //--- Test run output
+ //
+ double FACTO=1./Math.sqrt((double)(NTRGEN));
+ //
+ //--- Calculate means and std's of the pull values
+ //
+// System.out.println("Mean +/- error and standard deviation of pulls");
+ for (int i=0; i<3; ++i)
+ {
+ _pulsum[0][i]=_pulsum[0][i]/(double)NTRGEN;
+ _pulsum[1][i]=Math.sqrt(_pulsum[1][i]/((double)NTRGEN)-_pulsum[0][i]*_pulsum[0][i]);
+ _pullerr[i]=FACTO*_pulsum[1][i];
+// System.out.println("Mean: "+_pulsum[0][i]+" +/- "+_pullerr[i]+" sigma= "+_pulsum[1][i]);
+ assertEquals(_pulsum[0][i],0., .05); // mean should be 0.
+ assertEquals(_pulsum[1][i], 1., 0.5); // width of pull should be 1.
+ }
+
+ //
+ //--- Timing tests
+ //
+ //
+ //--- Generate a circle with random parameters with 500 points
+ //
+ simCirclePoints(.01,1.23,0.009,0.0033,500,0.0003,0.01,_xx, _yy, _wxy, _xyr);
+ int[] nptest = {4, 20, 100, 500};
+ long[] times = new long[nptest.length];
+ _fitter.setreferenceposition( 0., 0.);
+ //Fit this circle a number of times and save times as
+ //a function of the number of points fit.
+ for(int i=0; i<nptest.length; ++i)
+ {
+ long t1 = System.currentTimeMillis();
+ int ndo = 5000;
+ for(int j=0; j<ndo; ++j)
+ {
+ boolean OK = _fitter.fit( _xx, _yy, _wxy,nptest[i]);
+ }
+ times[i] = System.currentTimeMillis()-t1;
+// System.out.println("It took "+(double)times[i]/(double)(ndo*nptest[i])+" milliseconds for circles with "+nptest[i]+" hits");
+ }
+
+ }
+
+
+ private void simCirclePoints(double rho, double phi, double x, double y,
+ int npoints, double avedev, double step, double[] xx, double[] yy, double[] wxy, double[] xyr)
+ {
+ /*
+ *****************************************************************
+ SUBROUTINE CIRGEN(RHO,PHI,X,Y,NPO,AVEDEV,STEP)
+ * *
+ * SUBROUTINE TO GENERATE MEASURED POINTS ALONG CIRCLE *
+ * POINTS ARE GAUSSIAN FLUCTUATED ABOUT THE TRUE CIRCLE. *
+ * *
+ * THE ROUTINE IS PART OF THE CIRCLE FIT TEST PROGRAM *
+ * *
+ * INPUT PARAMETERS: *
+ * RHO = SIGNED CURVATURE (+VE BENDING CLOCKWISE) *
+ * PHI = DIRECTION ANGLE IN XY-PROJECTION *
+ * X,Y = STARTING POINT COORDINATES *
+ * NPO = NUMBER OF POINTS TO GENERATE *
+ * AVEDEV = MEAN DEVIATION NORMAL TO TRACK *
+ * STEP = STEP LENGTH ALONG TRACK *
+ * *
+ * OUTPUT PARAMETERS: *
+ * OUTPUT COORDINATES ARE IN ARRAYS XX,YY *
+ * --------------------------------------------------------------*
+ */
+ /*
+ PARAMETER (NRAND=1000)
+ COMMON/RANTAB/ LUNI,LGAU,UNITAB(NRAND),GAUTAB(NRAND)
+ PARAMETER (NCMAX=500)
+ COMMON /CIRPOI/ XX(NCMAX),YY(NCMAX),WXY(NCMAX),XYR(NCMAX)
+ REAL*8 PHIS,HDPHI,CORD,XS,YS
+ C-
+ */
+ double HDPHI = -rho*step/2.0;
+ double CORD = step;
+ // if(Math.abs(HDPHI)<1.0e-5f) double CORD = step;
+ if(Math.abs(HDPHI)>=1.0e-5f) CORD = 2.*Math.sin(-HDPHI)/rho;
+ double XS = x; // x start
+ double YS = y; // y start
+ double PHIS = phi; // phi start
+ Random r = new Random();
+ // GENERATE POINTS ALONG CIRCLE
+ for(int NP=0; NP<npoints; ++NP)
+ {
+ double AVERES = (0.5+r.nextDouble())*avedev;
+ double RESIXY = r.nextGaussian()*AVERES;
+
+ // FLUCTUATE NORMAL TO TRAJECTORY AND STORE POINT
+ xx[NP] = XS + Math.sin(PHIS)*RESIXY;
+ yy[NP] = YS - Math.cos(PHIS)*RESIXY;
+ wxy[NP]= 1./(AVERES*AVERES);
+ xyr[NP]= RESIXY;
+
+ // STEP TO THE NEXT POINT ON CIRCLE
+ XS += CORD*Math.cos(PHIS+HDPHI);
+ YS += CORD*Math.sin(PHIS+HDPHI);
+ PHIS += 2.*HDPHI;
+ }
+
+ }
+
+}
