lcsim/src/org/lcsim/contrib/RobKutschke/TRF/trffit
diff -N AddFitKalman.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ AddFitKalman.java 8 Jun 2009 06:02:25 -0000 1.1
@@ -0,0 +1,706 @@
+package org.lcsim.contrib.RobKutschke.TRF.trffit;
+
+import org.lcsim.util.aida.AIDA;
+
+import java.util.List;
+import java.util.ArrayList;
+import org.lcsim.recon.tracking.trfutil.Assert;
+import org.lcsim.recon.tracking.trfbase.ETrack;
+import org.lcsim.recon.tracking.trfbase.Hit;
+import org.lcsim.recon.tracking.trfbase.TrackVector;
+import org.lcsim.recon.tracking.trfbase.TrackError;
+import org.lcsim.recon.tracking.trfbase.Surface;
+import org.lcsim.recon.tracking.trfbase.PropDir;
+import Jama.Matrix;
+
+import java.math.BigDecimal;
+import java.math.MathContext;
+
+import org.lcsim.recon.tracking.trffit.AddFitter;
+
+import org.lcsim.contrib.RobKutschke.TRFSelfTest.util.PrintSymMatrix;
+import org.lcsim.contrib.RobKutschke.TRFSelfTest.util.VTUtil;
+import org.lcsim.contrib.RobKutschke.TRFSelfTest.util.RKDebug;
+
+import org.lcsim.contrib.RobKutschke.ToyConfig.*;
+
+
+// Fit tracks using Kalman filter.
+
+/**
+ * Copied from org.lcsim.recon.tracking.trffit.AddFitKalman
+ * and added debug printout.
+ *
+ * Original author was Norman Graf.
+ *
+ * AddFitKalman uses a Kalman filter to update the track fit.
+ *
+ *
+ *@author $Author: kutschke $
+ *@version $Id: AddFitKalman.java,v 1.1 2009/06/08 06:02:25 kutschke Exp $
+ *
+ * Date $Date: 2009/06/08 06:02:25 $
+ */
+
+public class AddFitKalman extends AddFitter
+{
+
+ private AIDA aida = AIDA.defaultInstance();
+
+
+ // Maximum allowed hit dimension.
+ private static final int MAXDIM = 3;
+
+ // Maximum number of track vectors.
+ private static final int MAXVECTOR = 1;
+
+ // Maximum number of track errors.
+ private static final int MAXERROR = 1;
+
+ //private: // nested classes
+
+ // The nested class Box holds the vectors, matrices and symmetric
+ // matrices needed for adding a hit in the main class.
+ class Box
+ {
+ /*
+ private: // typedefs
+ typedef Ptr<TrfVector,AutoPolicy> VectorPtr;
+ typedef Ptr<TrfSMatrix,AutoPolicy> SMatrixPtr;
+ typedef Ptr<TrfMatrix,AutoPolicy> MatrixPtr;
+ typedef vector<VectorPtr> VectorList;
+ typedef vector<SMatrixPtr> SMatrixList;
+ typedef vector<MatrixPtr> MatrixList;
+ */
+
+ // enums
+ // number of vectors
+ private static final int NVECTOR = 2 ;
+ // number of errors
+ private static final int NERROR = 3 ;
+ // number of vectors
+ private static final int NDERIV = 2 ;
+ // number of gains
+ private static final int NGAIN = 2 ;
+
+ // attributes
+ // dimension of the vector, matrix, etc
+ private int _size;
+ // array of vectors
+ List _vectors;
+ // array of error matrices
+ List _errors;
+ // array of derivatives (Nx5 matrices)
+ List _derivs;
+ // array of gains (5xN matrices)
+ List _gains;
+
+ // methods
+ // constructor
+ public Box(int size)
+ {
+ _size= size;
+ _vectors = new ArrayList();
+ _errors = new ArrayList();
+ _derivs = new ArrayList();
+ _gains = new ArrayList();
+
+ int icnt;
+ // Track vectors
+ for ( icnt=0; icnt<NVECTOR; ++icnt )
+ _vectors.add( new Matrix(size,1) );
+ // Track errors
+ for ( icnt=0; icnt<NERROR; ++icnt )
+ _errors.add( new Matrix(size,size) );
+ // Track derivatives
+ for ( icnt=0; icnt<NDERIV; ++icnt )
+ _derivs.add( new Matrix(size,5) );
+ // Gains
+ for ( icnt=0; icnt<NDERIV; ++icnt )
+ _gains.add( new Matrix(5,size) );
+
+ }
+
+ // return the dimension
+ public int get_size()
+ { return _size;
+ }
+ // fetch a vector
+ public Matrix get_vector(int ivec)
+ {
+ return (Matrix) _vectors.get(ivec);
+ }
+ // fetch an error
+ public Matrix get_error(int ierr)
+ {
+ return (Matrix) _errors.get(ierr);
+ }
+ // fetch a derivative
+ public Matrix get_deriv(int ider)
+ {
+ return (Matrix) _derivs.get(ider);
+ }
+ // fetch a gain
+ public Matrix get_gain(int igai)
+ {
+ return (Matrix) _gains.get(igai);
+ }
+ } // end of Box inner class
+
+ // static methods
+
+ //
+
+ /**
+ *Return a String representation of the class' the type name.
+ *Included for completeness with the C++ version.
+ *
+ * @return A String representation of the class' the type name.
+ */
+ public static String typeName()
+ { return "AddFitKalman";
+ }
+
+ //
+
+ /**
+ *Return a String representation of the class' the type name.
+ *Included for completeness with the C++ version.
+ *
+ * @return A String representation of the class' the type name.
+ */
+ public static String staticType()
+ { return typeName();
+ }
+
+ // attributes
+
+ // Array of boxes for each supported size.
+ private List _boxes;
+
+ // Track vectors.
+ private List _tvectors;
+
+ // Track errors.
+ private List _terrors;
+
+ //methods
+
+ //
+
+ int AddFitKalmanDebugLevel = 0;
+
+
+ /**
+ *Construct a default instance.
+ * Allocate space needed for hits of dimension from 1 to MAXDIM.
+ *
+ */
+ public AddFitKalman()
+ {
+ _boxes = new ArrayList();
+ _tvectors = new ArrayList();
+ _terrors = new ArrayList();
+
+ // Create boxes for hit containers.
+ for ( int dim=1; dim<MAXDIM; ++dim )
+ _boxes.add( new Box(dim) );
+
+ int icnt;
+
+ for ( icnt=0; icnt<MAXVECTOR; ++icnt )
+ _tvectors.add( new Matrix(5, 1) );
+
+ for ( icnt=0; icnt<MAXERROR; ++icnt )
+ _terrors.add( new Matrix(5,5) );
+
+ try{
+ ToyConfig config = ToyConfig.getInstance();
+ AddFitKalmanDebugLevel = config.getInt( "AddFitKalmanDebugLevel",
+ AddFitKalmanDebugLevel );
+
+
+ } catch (ToyConfigException e){
+ System.out.println (e.getMessage() );
+ System.out.println ("Stopping now." );
+ System.exit(-1);
+ }
+
+
+ }
+
+ //
+
+ /**
+ *Return a String representation of the class' the type name.
+ *Included for completeness with the C++ version.
+ *
+ * @return A String representation of the class' the type name.
+ */
+ public String type()
+ { return staticType();
+ }
+
+ //
+
+
+ /**
+ *Add a hit and fit with the new hit.
+ * Use a Kalman filter to add a hit to a track.
+ * The hit is updated with the input track.
+ * Note: We make direct use of the underlying vector and matrix
+ * classes here. It will probably be neccessary to modify
+ * this routine if these are changed.
+ *
+ * @param tre The ETrack to update.
+ * @param chsq The chi-square for the fit.
+ * @param hit The Hit to add to the track.
+ * @return 0 if successful.
+ */
+ public int addHitFit(ETrack tre, double chsq, Hit hit)
+ {
+
+ // Update the hit with the input track.
+ hit.update(tre);
+
+ // Fetch hit size.
+ int dim = hit.size();
+ Assert.assertTrue( dim <= MAXDIM );
+
+ // Fetch the box holding the needed hit containers.
+ // The chice of boxes depends on the size of the hit.
+ Box box = (Box)_boxes.get(dim-1);
+ Assert.assertTrue( box.get_size() == dim );
+
+ // Fetch the hit containers.
+ Matrix diff = box.get_vector(0);
+ Matrix hit_res = box.get_vector(1);
+ Matrix hit_err = box.get_error(0);
+ Matrix hit_err_tot = box.get_error(1);
+ Matrix hit_res_err = box.get_error(2);
+ Matrix dhit_dtrk = box.get_deriv(0);
+ Matrix new_dhit_dtrk = box.get_deriv(1);
+ Matrix trk_err_dhit_dtrk = box.get_gain(0);
+ Matrix gain = box.get_gain(1);
+ Matrix gain2 = gain.copy();
+
+ // Fetch the track containers.
+ Matrix new_vec = (Matrix) _tvectors.get(0);
+ Matrix new_err = (Matrix) _terrors.get(0);
+
+ hit_err = hit.measuredError().matrix();
+ Matrix Vm = hit_err.copy();
+ Matrix VmInv = Vm.inverse();
+ //System.out.println("hit_err= \n"+hit_err);
+
+ // Fetch track prediction of hit.
+ diff = hit.differenceVector().matrix();
+ //System.out.println("diff= \n"+diff);
+ dhit_dtrk = hit.dHitdTrack().matrix();
+ Matrix D = dhit_dtrk.copy();
+ //System.out.println("dhit_dtrk= \n"+dhit_dtrk);
+ // Fetch track info.
+ Matrix trk_vec = tre.vector().matrix();
+ Matrix trk_err = tre.error().getMatrix(); //need to fix this!
+ Matrix V = tre.error().getMatrix();
+ VTUtil vtu = new VTUtil( tre );
+
+ //System.out.println("trk_vec= \n"+trk_vec);
+ //System.out.println("trk_err= \n"+trk_err);
+
+ // Build gain matrix.
+ hit_err_tot = hit.predictedError().matrix().plus(
+ hit_err);
+ //System.out.println("hit_err_tot= \n"+hit_err_tot);
+ //if ( invert(hit_err_tot)!=0 ) return 3;
+ Matrix hetot = new Matrix(2,2);
+ Matrix hetoti = new Matrix(2,2);
+ hetot = hit.predictedError().matrix().plus(hit_err);
+ hetoti = hetot.inverse();
+ hit_err_tot= hit_err_tot.inverse();
+
+ //System.out.println("hit_err_tot inverse= \n"+hit_err_tot);
+ trk_err_dhit_dtrk = trk_err.times(dhit_dtrk.transpose());
+ gain = trk_err_dhit_dtrk.times(hit_err_tot);
+ //System.out.println("trk_err_dhit_dtrk= \n"+trk_err_dhit_dtrk);
+ //System.out.println("gain= \n"+gain);
+
+ // if ( get_debug() ) {
+ //System.out.println("\n");
+ //System.out.println(" trk_vec: " + "\n"+ trk_vec + "\n");
+ //System.out.println(" trk_err: " + "\n"+ trk_err + "\n");
+ //System.out.println(" dhit_dtrk: " + "\n"+ dhit_dtrk + "\n");
+ // }
+
+ // We need to return dhit_dtrk to its original state for the
+ // next call.
+ // dhit_dtrk = dhit_dtrk.transpose(); //need to check this!
+ dhit_dtrk.transpose();
+ // Build new track vector.
+ new_vec = trk_vec.minus(gain.times(diff));
+ //System.out.println("new_vec= \n"+new_vec);
+
+ // Build new error;
+ new_err = trk_err.minus(trk_err_dhit_dtrk.times(hit_err_tot.times(trk_err_dhit_dtrk.transpose())));
+
+ Matrix Vnew = new_err.copy();
+ D = D.transpose();
+ gain2 = Vnew.times(D.times(VmInv));
+ Matrix etap = trk_vec.minus(gain2.times(diff));
+
+ //System.out.println("new_err= \n"+new_err);
+ // Check the error.
+ if ( AddFitKalmanDebugLevel > 0 ) {
+ int nbad = 0;
+ for ( int i=0; i<5; ++i )
+ {
+ if ( new_err.get(i,i) < 0.0 ) {
+ if ( nbad == 0 ) System.out.println("");
+ System.out.println ( "Bad on: " + i + " " + new_err.get(i,i) );
+ ++nbad;
+ }
+ double eii = new_err.get(i,i);
+ for ( int j=0; j<i; ++j )
+ {
+ double ejj = new_err.get(j,j);
+ double eij = new_err.get(j,i);
+ if ( Math.abs(eij*eij) >= eii*ejj ) {
+ double delta = -Math.abs(eij*eij) + eii*ejj;
+ if ( nbad == 0 ) System.out.println("");
+ System.out.println( "Bad off: " + i + " "
+ + j + " "
+ + eii + " "
+ + ejj + " "
+ + eij + " "
+ + delta
+ );
+ ++nbad;
+ }
+ }
+ }
+ if ( nbad > 0 ){
+ PrintSymMatrix psm = new PrintSymMatrix();
+ System.out.println ("Illegal cov in addfitkalman: "
+ + RKDebug.Instance().getTrack() + " "
+ + RKDebug.Instance().getPropDir() + " "
+ );
+ System.out.println ("Surface: " + tre.surface() );
+
+ System.out.println ("old_err \n" );
+ psm.Print(tre.error());
+
+ System.out.println ("Predicted error: " );
+ psm.Print(hit.predictedError().matrix() );
+
+ System.out.println ("Measured error: " );
+ psm.Print(hit.measuredError().matrix() );
+
+ System.out.println ("Total error: " );
+ psm.Print( hetot );
+
+ System.out.println ("Inverse: " );
+ psm.Print( hetoti );
+
+ psm.Print( hetoti.times(hetot) );
+
+
+ System.out.println ("Derivs: \n" + hit.dHitdTrack() );
+
+ System.out.println ("new_err \n" );
+ psm.Print(new_err);
+ }
+ if ( nbad > 0 ) return 5;
+ }
+
+ // Create track vector with new values.
+ tre.setVectorAndKeepDirection(new TrackVector(new_vec));
+ tre.setError(new TrackError(new_err));
+
+ // Calculate residual vector.
+
+ // Update the hit with the new track.
+ //System.out.println("update the hit");
+ hit.update(tre);
+
+ hit_res = hit.differenceVector().matrix();
+ new_dhit_dtrk = hit.dHitdTrack().matrix();
+ //System.out.println("new_dhit_dtrk= \n"+new_dhit_dtrk);
+
+ // Calculate residual covariance and invert.
+ hit_res_err = hit_err.minus(dhit_dtrk.times(new_err.times(dhit_dtrk.transpose())));
+ // System.out.println("hit_res_err= \n"+hit_res_err);
+ hit_res_err = hit_res_err.inverse();
+ //System.out.println("hit_res_err inverse= \n"+hit_res_err);
+ // Update chi-square.
+ // result should be 1x1 matrix, so should be able to do the following
+ //System.out.println( " hr*hre*hrT=\n"+(hit_res.transpose()).times(hit_res_err.times(hit_res)) );
+ double dchsq = (hit_res.transpose()).times(hit_res_err.times(hit_res)).get(0,0);
+ //System.out.println("chsq= "+chsq+", dchsq= "+dchsq);
+
+ Matrix dEta = etap.minus(trk_vec);
+ Matrix dM = hit.differenceVector().matrix();
+ Matrix VInv = V.inverse();
+
+ /*
+ System.out.printf ( "Dim1: (%d,%d) (%d,%d) (%d,%d)\n",
+ dEta.getRowDimension(), dEta.getColumnDimension(),
+ dM.getRowDimension(), dM.getColumnDimension(),
+ VmInv.getRowDimension(), VmInv.getColumnDimension()
+ );
+ */
+ double dch1 = ((dM.transpose()).times(VmInv.times(dM))).get(0,0);
+ double dch2 = ((dEta.transpose()).times(VInv.times(dEta))).get(0,0);
+ double dch = dch1 + dch2;
+
+ double ddd = dch-dchsq;
+
+ if ( Math.abs(ddd) > 0.01 ){
+ if ( RKDebug.Instance().getPropDir() == PropDir.BACKWARD ){
+ aida.cloud1D("/Bugs/Chisq/Back: Delta Chisq: ").fill(dch-dchsq);
+ aida.cloud2D("/Bugs/Chisq/Back: Chisq TRF vs CLEO").fill(dch,dchsq);
+ aida.cloud1D("/Bugs/Chisq/radius of bad chisq").fill( vtu.r() );
+ aida.cloud1D("/Bugs/Chisq/z of bad chisq").fill( vtu.z() );
+ aida.cloud2D("/Bugs/Chisq/Back: Chisq TRF vs r").fill( vtu.r(), dchsq );
+ aida.cloud2D("/Bugs/Chisq/Back: Chisq TRF vs z").fill( vtu.z(), dchsq );
+ } else {
+ aida.cloud1D("/Bugs/Chisq/Forw: Delta Chisq: ").fill(dch-dchsq);
+ aida.cloud2D("/Bugs/Chisq/Forw: Chisq TRF vs CLEO").fill(dch,dchsq);
+ }
+ } else {
+ if ( RKDebug.Instance().getPropDir() == PropDir.BACKWARD ){
+ aida.cloud1D("/Bugs/Chisq/Back: Small Delta Chisq: ").fill(dch-dchsq);
+ } else {
+ aida.cloud1D("/Bugs/Chisq/Forw: Small Delta Chisq: ").fill(dch-dchsq);
+ }
+ }
+
+ dchsq = dch;
+
+
+ /*
+ System.out.println ("Chitest: " + dchsq + " " + dch + " | " + dch1 + " " + dch2 );
+ for ( int i=0; i<5; ++i ){
+ System.out.printf ("Eta %3d %12.7f %12.7f %12.7f %12.7f %12.7f | %12.7f \n", i,
+ trk_vec.get(i,0),
+ new_vec.get(i,0),
+ etap.get(i,0),
+ (new_vec.minus(etap)).get(i,0),
+ Math.sqrt(new_err.get(i,i)),
+ dEta.get(i,0)
+ );
+
+ }
+
+ if ( Vm.getRowDimension() == 1){
+ Hack1D ( trk_vec, V, Vm, new_vec, new_err );
+ } else if ( Vm.getRowDimension() == 2){
+ Hack2D ( trk_vec, V, Vm, diff, new_vec, new_err );
+ }
+ */
+
+ // Warn about bad chisquared contribution, ignoring small problems that might
+ // well be just round off error.
+
+ if ( dchsq < -0.001 ){
+ /*
+ System.out.println ("\ndchisq: "
+ + RKDebug.Instance().getTrack() + " "
+ + RKDebug.Instance().getPropDir() + " "
+ + dchsq + " "
+ + tre.surface() + " "
+ + RKDebug.Instance().getRKTrack().cz() + " "
+ );
+ */
+ if ( RKDebug.Instance().getPropDir() == PropDir.BACKWARD ){
+ aida.histogram1D( "/Bugs/Backward: cz for tracks with bad dchisquared",100,-1.,1.)
+ .fill( RKDebug.Instance().getRKTrack().cz() );
+ aida.histogram1D( "/Bugs/Backward: bad dchisquared",100,-20.,0.).fill( dchsq );
+ aida.cloud2D( "/Bugs/Backward: dch vs cz for tracks with bad dchisquared",-1)
+ .fill( RKDebug.Instance().getRKTrack().cz(), dchsq );
+ }else{
+ aida.histogram1D( "/Bugs/Forward: cz for tracks with bad dchisquared",100,-1.,1.)
+ .fill( RKDebug.Instance().getRKTrack().cz() );
+ aida.histogram1D( "/BugsForward: bad dchisquared",100,-20.,0.).fill( dchsq );
+ aida.cloud2D( "/Bugs/Forward: dch vs cz for tracks with bad dchisquared",-1)
+ .fill( RKDebug.Instance().getRKTrack().cz(), dchsq );
+ }
+ }
+
+ chsq = chsq + dchsq;
+ setChisquared(chsq);
+
+ // if ( get_debug() ) {
+ //System.out.println(" gain: " + "\n"+ gain + "\n");
+ //System.out.println(" new_vec: " + "\n"+ new_vec + "\n");
+ //System.out.println(" new_err: " + "\n"+ new_err + "\n");
+ //System.out.println("new_dhit_dtrk: " + "\n"+ new_dhit_dtrk + "\n");
+ //System.out.println(" hit_res: " + "\n"+ hit_res + "\n");
+ //System.out.println(" hit_res_err: " + "\n"+ hit_res_err + "\n");
+ //System.out.println(" dchsq & chsq: " + "\n"+ dchsq + " " + chsq + "\n");
+ // }
+
+ return 0;
+
+ }
+
+
+ /**
+ *output stream
+ *
+ * @return The String representation of this instance.
+ */
+ public String toString()
+ {
+ return getClass().getName();
+ }
+
+ /**
+ * Perform matrix operations as BigDecimal, specialized to a 1D measurement.
+ *
+ *
+ */
+ private void Hack1D ( Matrix etain,
+ Matrix Vin,
+ Matrix Vm,
+ Matrix etat,
+ Matrix Vt ){
+
+ BigDecimal sigsq = new BigDecimal( Vm.get(0,0), MathContext.DECIMAL128 );
+ BigDecimal V00 = new BigDecimal( Vin.get(0,0), MathContext.DECIMAL128 );
+ BigDecimal denom = sigsq.add(V00);
+
+ BigDecimal[] Vi0 = { new BigDecimal( Vin.get(0,0), MathContext.DECIMAL128 ),
+ new BigDecimal( Vin.get(1,0), MathContext.DECIMAL128 ),
+ new BigDecimal( Vin.get(2,0), MathContext.DECIMAL128 ),
+ new BigDecimal( Vin.get(3,0), MathContext.DECIMAL128 ),
+ new BigDecimal( Vin.get(4,0), MathContext.DECIMAL128 ) };
+
+ Matrix Vp = new Matrix(5,5);
+ for ( int i=0; i<5; ++i ){
+ for ( int j=i; j<5; ++j ){
+ BigDecimal Vij = new BigDecimal( Vin.get(i,j), MathContext.DECIMAL128 );
+ BigDecimal numer = Vi0[i].multiply(Vi0[j]);
+ BigDecimal ratio = numer.divide(denom,MathContext.DECIMAL128);
+ BigDecimal Vijnew = Vij.subtract(ratio);
+ Vp.set(i,j,Vijnew.doubleValue() );
+ }
+ }
+ System.out.println ("Hack 1D");
+ PrintSymMatrix psm = new PrintSymMatrix();
+ psm.Print(Vp);
+ }
+
+ /**
+ * Perform matrix operations as BigDecimal, specialized to a 2D measurement.
+ *
+ *
+ */
+ private void Hack2D ( Matrix etain,
+ Matrix Vin,
+ Matrix Vm,
+ Matrix dM,
+ Matrix etat,
+ Matrix Vt ){
+
+ BigDecimal DTVD00 = new BigDecimal( Vin.get(0,0), MathContext.DECIMAL128 );
+ BigDecimal DTVD11 = new BigDecimal( Vin.get(1,1), MathContext.DECIMAL128 );
+ BigDecimal DTVD01 = new BigDecimal( Vin.get(0,1), MathContext.DECIMAL128 );
+
+ BigDecimal sigsqx = new BigDecimal( Vm.get(0,0), MathContext.DECIMAL128 );
+ BigDecimal sigsqy = new BigDecimal( Vm.get(1,1), MathContext.DECIMAL128 );
+
+ BigDecimal Sum00 = DTVD00.add(sigsqx);
+ BigDecimal Sum11 = DTVD11.add(sigsqy);
+ BigDecimal Sum01 = DTVD01;
+
+ BigDecimal Disc = (Sum00.multiply(Sum11)).subtract(Sum01.multiply(Sum01));
+
+ BigDecimal U00 = Sum11.divide(Disc,MathContext.DECIMAL128);
+ BigDecimal U11 = Sum00.divide(Disc,MathContext.DECIMAL128);
+ BigDecimal U01 = (Sum01.negate()).divide(Disc,MathContext.DECIMAL128);
+ BigDecimal U10 = U01;
+
+
+
+ BigDecimal[] Vi0 = { new BigDecimal( Vin.get(0,0), MathContext.DECIMAL128 ),
+ new BigDecimal( Vin.get(1,0), MathContext.DECIMAL128 ),
+ new BigDecimal( Vin.get(2,0), MathContext.DECIMAL128 ),
+ new BigDecimal( Vin.get(3,0), MathContext.DECIMAL128 ),
+ new BigDecimal( Vin.get(4,0), MathContext.DECIMAL128 ) };
+
+ BigDecimal[] Vi1 = { new BigDecimal( Vin.get(0,1), MathContext.DECIMAL128 ),
+ new BigDecimal( Vin.get(1,1), MathContext.DECIMAL128 ),
+ new BigDecimal( Vin.get(2,1), MathContext.DECIMAL128 ),
+ new BigDecimal( Vin.get(3,1), MathContext.DECIMAL128 ),
+ new BigDecimal( Vin.get(4,1), MathContext.DECIMAL128 ) };
+
+ // New covariance matrix, both as Matrix and as 2D array of BigDecimal.
+ Matrix Vp = new Matrix(5,5);
+ BigDecimal BigZero = new BigDecimal ( 0., MathContext.DECIMAL128 );
+ BigDecimal[][] VpBig = { { BigZero, BigZero, BigZero, BigZero, BigZero},
+ { BigZero, BigZero, BigZero, BigZero, BigZero},
+ { BigZero, BigZero, BigZero, BigZero, BigZero},
+ { BigZero, BigZero, BigZero, BigZero, BigZero},
+ { BigZero, BigZero, BigZero, BigZero, BigZero} };
+
+
+ for ( int i=0; i<5; ++i ){
+ for ( int j=i; j<5; ++j ){
+
+ BigDecimal a = U00.multiply(Vi0[j]);
+ BigDecimal b = U01.multiply(Vi1[j]);
+ BigDecimal c = U01.multiply(Vi0[j]);
+ BigDecimal d = U11.multiply(Vi1[j]);
+
+ BigDecimal s1 = Vi0[i].multiply( a.add(b));
+ BigDecimal s2 = Vi1[i].multiply( c.add(d));
+
+ BigDecimal Vij = new BigDecimal( Vin.get(i,j), MathContext.DECIMAL128 );
+ BigDecimal Vijnew = Vij.subtract( s1.add(s2) );
+
+ Vp.set(i,j,Vijnew.doubleValue() );
+ VpBig[i][j] = Vijnew;
+ double denom = (Vp.get(i,j)+Vt.get(i,j))/2.;
+ if ( denom != 0. ){
+ double frac = (Vp.get(i,j)-Vt.get(i,j))/denom;
+ aida.cloud1D( "/Bugs/AddFit/Frac Vdiff " + i + " " + j ).fill(frac);
+ }
+ if ( i != j ){
+ Vp.set(j,i,Vp.get(i,j));
+ VpBig[j][i] = Vijnew;
+ }
+ }
+ }
+
+ // Predicted-measured vector as BigDecimal.
+ BigDecimal dM0 = new BigDecimal( dM.get(0,0), MathContext.DECIMAL128);
+ BigDecimal dM1 = new BigDecimal( dM.get(1,0), MathContext.DECIMAL128);
+
+ // New track parameters, both as Matrix and as 2D array of BigDecimal.
+ BigDecimal[] etanew = { BigZero, BigZero, BigZero, BigZero, BigZero };
+ Matrix etap = new Matrix(5,1);
+ for ( int i=0; i<5; ++i ){
+ BigDecimal rx = VpBig[i][0].divide(sigsqx,MathContext.DECIMAL128);
+ BigDecimal ry = VpBig[i][1].divide(sigsqy,MathContext.DECIMAL128);
+ BigDecimal deta = (rx.multiply(dM0)).add(ry.multiply(dM1));
+ BigDecimal etai = new BigDecimal( etain.get(i,0),MathContext.DECIMAL128);
+ etanew[i] = etai.subtract(deta);
+ etap.set(i,0,etanew[i].doubleValue());
+ aida.cloud1D("/Bugs/AddFit/Diff eta " + i).fill(etap.get(i,0)-etat.get(i,0));
+ double denom = (etap.get(i,0)+etat.get(i,0))/2.;
+ if ( denom != 0. ){
+ double frac = (etap.get(i,0)-etat.get(i,0))/denom;
+ aida.cloud1D("/Bugs/AddFit/Frac Diff eta " + i).fill(frac);
+ }
+
+ }
+
+ System.out.println ("Hack 2D " + dM0 + " " + dM1 );
+ PrintSymMatrix psm = new PrintSymMatrix();
+ psm.Print(etap, Vp);
+
+
+ }
+
+}
+
lcsim/src/org/lcsim/contrib/RobKutschke/TRF/trffit
diff -N FullFitKalman.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ FullFitKalman.java 8 Jun 2009 06:02:25 -0000 1.1
@@ -0,0 +1,610 @@
+package org.lcsim.contrib.RobKutschke.TRF.trffit;
+
+import org.lcsim.recon.tracking.trfbase.Propagator;
+import org.lcsim.recon.tracking.trfbase.PropDir;
+import org.lcsim.recon.tracking.trfbase.PropStat;
+import org.lcsim.recon.tracking.trfbase.Hit;
+import org.lcsim.recon.tracking.trfbase.ETrack;
+import org.lcsim.recon.tracking.trfbase.TrackError;
+import org.lcsim.recon.tracking.trfbase.Surface;
+import java.util.*;
+
+import org.lcsim.recon.tracking.trffit.HTrack;
+import org.lcsim.recon.tracking.trffit.FullFitter;
+
+import org.lcsim.recon.tracking.trfdca.SurfDCA;
+
+import org.lcsim.recon.tracking.trfcyl.SurfCylinder;
+import org.lcsim.recon.tracking.trfcyl.ThinCylMs;
+
+import org.lcsim.recon.tracking.trfzp.SurfZPlane;
+import org.lcsim.recon.tracking.trfzp.ThinZPlaneMs;
+
+import org.lcsim.recon.tracking.trfcyl.HitCylPhi;
+import org.lcsim.recon.tracking.trfcyl.HitCylPhiZ2D;
+import org.lcsim.recon.tracking.trfzp.HitZPlane1;
+import org.lcsim.recon.tracking.trfzp.HitZPlane2;
+
+import org.lcsim.contrib.RobKutschke.TRFSelfTest.util.SurfaceCode;
+import org.lcsim.contrib.RobKutschke.TRFSelfTest.util.VTUtil;
+import org.lcsim.contrib.RobKutschke.TRFSelfTest.util.RKDebug;
+import org.lcsim.contrib.RobKutschke.TRFSelfTest.util.RKZot;
+
+
+import org.lcsim.contrib.RobKutschke.ToyConfig.*;
+
+import org.lcsim.util.aida.AIDA;
+
+/**
+ *
+ * Copied from org.lcsim.recon.tracking.trffit.FullFitKalman
+ * - added fitForward() and fitBackward() methods.
+ * - added debug printout
+ *
+ * Full track fit using Kalman filter. The propagator is specified
+ * when the fitter is constructed. The starting surface, vector and
+ * error matrix are taken from the input track. Errors should be
+ * increased appropriately if the fitter is applied repeatedly to
+ * a single track.
+ *
+ *@author $Author: kutschke $
+ *@version $Id: FullFitKalman.java,v 1.1 2009/06/08 06:02:25 kutschke Exp $
+ *
+ * Date $Date: 2009/06/08 06:02:25 $
+ *
+ */
+
+public class FullFitKalman extends FullFitter
+{
+
+ private AIDA aida = AIDA.defaultInstance();
+
+ // Flags to control: multiple scattering, energy loss and adding the hit.
+ private boolean doMs = true;
+ private boolean doEloss = true;
+ private boolean doMeas = true;
+
+ private double dedxscale = 1.;
+ private double dedxsigma = 0.0;
+
+ // static methods
+
+ //
+
+ /**
+ *Return a String representation of the class' the type name.
+ *Included for completeness with the C++ version.
+ *
+ * @return A String representation of the class' the type name.
+ */
+ public static String typeName()
+ { return "FullFitKalman"; }
+
+ //
+
+ /**
+ *Return a String representation of the class' the type name.
+ *Included for completeness with the C++ version.
+ *
+ * @return A String representation of the class' the type name.
+ */
+ public static String staticType()
+ { return typeName(); }
+
+ // The propagator.
+ private Propagator _pprop;
+
+ // The add fitter.
+ private AddFitKalman _addfit;
+
+ int AddFitKalmanDebugLevel = 0;
+ //
+
+ /**
+ *Construct an instance specifying a propagator.
+ *
+ * @param prop The Propagator to be used during the fit.
+ */
+ public FullFitKalman(Propagator prop)
+ {
+ _pprop = prop;
+ _addfit = new AddFitKalman();
+
+ try{
+ ToyConfig config = ToyConfig.getInstance();
+ AddFitKalmanDebugLevel = config.getInt( "AddFitKalmanDebugLevel",
+ AddFitKalmanDebugLevel );
+ dedxscale = config.getDouble("dEdXScale");
+ dedxsigma = config.getDouble("dEdXSigma");
+
+
+ } catch (ToyConfigException e){
+ System.out.println (e.getMessage() );
+ System.out.println ("Stopping now." );
+ System.exit(-1);
+ }
+ System.out.println ("FullfitKalman dedxscale: " + dedxscale );
+
+ }
+
+ //
+
+ /**
+ *Return a String representation of the class' type name.
+ *Included for completeness with the C++ version.
+ *
+ * @return A String representation of the class' the type name.
+ */
+ public String type()
+ { return staticType(); }
+
+ //
+
+ /**
+ *Return the propagator.
+ *
+ * @return The Propagator used in the fit.
+ */
+ public Propagator propagator()
+ { return _pprop; }
+
+ //
+
+ public void setDoMs( boolean b){
+ doMs = b;
+ }
+
+ public void setDoEloss( boolean b){
+ doEloss = b;
+ }
+
+ /**
+ *Fit the specified track.
+ *
+ * @param trh The HTrack to fit.
+ * @return 0 if successful.
+ */
+ public int fit(HTrack trh)
+ {
+ // Copy the hits from the track.
+ List hits = trh.hits();
+ //System.out.println("Hits has "+hits.size()+" elements");
+ // Delete the list of hits from the track.
+ while ( trh.hits().size()>0 ) trh.dropHit();
+ //System.out.println("Hits has "+hits.size()+" elements");
+
+ // Set direction to be nearest.
+ //PropDir dir = PropDir.NEAREST;
+ PropDir dir = PropDir.FORWARD;
+ RKDebug.Instance().setPropDir(dir);
+
+ // Loop over hits and fit.
+ int icount = 0;
+ for ( Iterator ihit=hits.iterator(); ihit.hasNext(); )
+ {
+
+ // Extract the next hit pointer.
+ Hit hit = (Hit)ihit.next();
+ //System.out.println("Hit "+icount+" is: \n"+hit);
+ // propagate to the surface
+ //System.out.println ("Before prop: " + trh.newTrack() );
+ PropStat pstat = trh.propagate(_pprop,hit.surface(),dir);
+
+ //System.out.println ("After prop: " + trh.newTrack() );
+ // fit track
+ //System.out.println("trh= \n"+trh+", hit= \n"+hit);
+ //System.out.println("_addfit= "+_addfit);
+ int fstat = _addfit.addHit(trh,hit);
+ //System.out.println ("After addhit: " + fstat + " "
+ // + trh.chisquared() + "\n" + trh.newTrack() );
+ if ( fstat>0 ) return 10000 + 1000*fstat + icount;
+
+ }
+ return 0;
+
+ }
+
+ public int fitForward(HTrack trh)
+ {
+
+ PropDir dir = PropDir.FORWARD;
+ RKDebug.Instance().setPropDir(dir);
+
+ // Copy the hits from the track.
+ List hits = trh.hits();
+
+ // Delete the list of hits from the track.
+ while ( trh.hits().size()>0 ) trh.dropHit();
+
+ double sumde = 0.;
+
+ // Loop over hits and fit.
+ int icount = 0;
+ for ( Iterator ihit=hits.iterator(); ihit.hasNext(); )
+ {
+ Surface s_save = trh.newTrack().surface().newPureSurface();
+ ETrack e_save = trh.newTrack();
+
+ // Extract the next hit pointer.
+ Hit hit = (Hit)ihit.next();
+
+ int from = (new SurfaceCode(s_save)).getCode();
+ int to = (new SurfaceCode(hit.surface())).getCode();
+
+ // Propagate to the next surface.
+ PropStat pstat = trh.propagate(_pprop,hit.surface(),dir);
+ if ( ! pstat.success() ) {
+ if ( AddFitKalmanDebugLevel > 0 ) {
+ System.out.println ("Error: "
+ + RKDebug.Instance().getTrack() + " "
+ + RKDebug.Instance().getPropDir() + " "
+ + icount
+ );
+ System.out.println ("From surface 5: " + s_save );
+ System.out.println ("To surface 5: " + hit.surface());
+ System.out.println ("Params: " + e_save.vector() );
+ }
+ aida.histogram1D("/Bugs/Fit/Failed Fwd prop from Surface",5,0,5).fill( from );
+ aida.histogram1D("/Bugs/Fit/Failed Fwd prop to Surface",5,0,5).fill( to );
+ aida.cloud2D("/Bugs/Fit/Failed Fwd prop to vs from Surface").fill( from, to );
+
+ return icount+1;
+ }
+
+ //if ( icount != 0 ) {
+ // int istat = interact( trh, hit, dir );
+ //}
+
+
+ // Add the hit.
+ int fstat = _addfit.addHit(trh,hit);
+ if ( fstat>0 ){
+ if ( AddFitKalmanDebugLevel > 0){
+
+ System.out.println ("Error: "
+ + RKDebug.Instance().getTrack() + " "
+ + RKDebug.Instance().getPropDir() + " "
+ );
+ System.out.println ("From surface 4: " + s_save );
+ System.out.println ("To surface 4: " + hit.surface());
+ }
+ aida.histogram1D("/Bugs/Fit/Failed Fwd addhit from Surface",5,5,5).fill( from );
+ aida.histogram1D("/Bugs/Fit/Failed Fwd addhit to Surface",5,0,5).fill( to );
+ aida.cloud2D("/Bugs/Fit/Failed Fwd addhit to vs from Surface").fill( from, to );
+
+ }
+ if ( fstat>0 ) return 10000 + 1000*fstat + icount;
+
+ VTUtil before = new VTUtil( trh.newTrack() );
+ int istat = interact( trh, hit, dir );
+ VTUtil after = new VTUtil( trh.newTrack() );
+
+ double de = before.e() - after.e();
+ sumde += de;
+
+ //SurfCylinder ss = (SurfCylinder)trh.newTrack().surface();
+ //System.out.printf ("Forw dedx: %10.4f %12.8f %12.8f\n", ss.radius(), de, sumde );
+
+
+ ++icount;
+ }
+
+
+ //System.out.println ("Forward fit sumde: " + sumde );
+ aida.cloud1D("Forward dedx check:").fill(sumde-RKDebug.Instance().getDeGen());
+
+ return 0;
+
+ }
+
+ public int fitBackward(HTrack trh)
+ {
+ PropDir dir = PropDir.BACKWARD;
+ RKDebug.Instance().setPropDir(dir);
+
+ //RKPrintSymMatrix psm = new RKPrintSymMatrix();
+
+ // Copy the hits from the track.
+ List hits = trh.hits();
+
+ // Delete the list of hits from the track.
+ while ( trh.hits().size()>0 ) trh.dropHit();
+
+ double chold = 0.;
+
+ double sumde = 0.;
+
+ RKZot zot = new RKZot(trh);
+ boolean zottable = true;
+
+ int nc = 0;
+ int nz = 0;
+ int nu = 0;
+ String thishit;
+
+
+ // Loop over hits and fit.
+ int icount = 0;
+ for ( Iterator ihit=hits.iterator(); ihit.hasNext(); )
+ {
+ // Extract the next hit pointer.
+ Hit hit = (Hit)ihit.next();
+
+ Surface s_save = trh.newTrack().surface().newPureSurface();
+
+ int from = (new SurfaceCode(s_save)).getCode();
+ int to = (new SurfaceCode(hit.surface())).getCode();
+
+ /*
+ Surface s_new = hit.surface();
+ if ( s_new instanceof SurfCylinder ){
+ System.out.printf ("Next: %3d Cylinder\n",icount);
+ ++nc;
+ thishit = "Cylinder";
+ } else if ( s_new instanceof SurfZPlane ){
+ System.out.printf ("Next: %3d ZPlane\n",icount);
+ ++nz;
+ thishit = "ZPlane";
+ zottable = false;
+ } else {
+ System.out.printf ("Next: %3d Unknown\n",icount);
+ ++nu;
+ thishit = "Unknown";
+ }
+ */
+ if ( hit instanceof HitCylPhi ){
+ //System.out.printf ("Next: %3d Cylinder 1D\n",icount);
+ ++nc;
+ thishit = "Cylinder1D";
+ } else if ( hit instanceof HitCylPhiZ2D ){
+ //System.out.printf ("Next: %3d Cylinder 2D\n",icount);
+ ++nc;
+ ++nz;
+ zottable = false;
+ thishit = "Cylinder2D";
+ } else if ( hit instanceof HitZPlane1 ){
+ //System.out.printf ("Next: %3d ZPlane 1D\n",icount);
+ ++nc;
+ ++nz;
+ zottable = false;
+ thishit = "ZPlane1D";
+ } else if ( hit instanceof HitZPlane2 ){
+ //System.out.printf ("Next: %3d ZPlane 2D\n",icount);
+ ++nc;
+ ++nz;
+ zottable = false;
+ thishit = "ZPlane2D";
+ } else{
+ System.out.printf ("Next: %3d Unknown\n",icount);
+ thishit = "Unknown";
+ }
+ //System.out.println("Error before prop: ");
+ //psm.Print(trh.newTrack());
+
+ PropStat pstat = trh.propagate(_pprop,hit.surface(),dir);
+
+ if ( ! pstat.success() ) {
+ if ( AddFitKalmanDebugLevel > 0 ){
+ System.out.println ("Error: "
+ + RKDebug.Instance().getTrack() + " "
+ + RKDebug.Instance().getPropDir() + " "
+ + icount
+ );
+ System.out.println ("From surface 1: " + s_save );
+ System.out.println ("To surface 1: " + hit.surface());
+ System.out.println ("Failed prop: " + nc + " " + nz + " " + nu + " " + thishit );
+ }
+
+ aida.histogram1D("/Bugs/Fit/Failed prop from Surface",5,0,5).fill( from );
+ aida.histogram1D("/Bugs/Fit/Failed prop to Surface",5,0,5).fill( to );
+ aida.cloud2D("/Bugs/Fit/Failed prop to vs from Surface").fill( from, to );
+
+ return icount+1;
+ }
+ //System.out.println("Error after prop: ");
+ //psm.Print(trh.newTrack());
+
+ if ( zottable ){
+ zot.Zot(trh);
+ //System.out.println("Error after zot: ");
+ //psm.Print(trh.newTrack());
+ }
+
+
+ //if ( icount != 0 ) {
+ //int istat = interact( trh, hit, dir );
+ //}
+ //VTUtil before = new VTUtil( trh.newTrack() );
+ //int istat = interact( trh, hit, dir );
+ //VTUtil after = new VTUtil( trh.newTrack() );
+
+
+ // Add the hit.
+ int fstat = _addfit.addHit(trh,hit);
+
+ //System.out.println ("Hit info: " + hit );
+
+
+ if ( fstat>0 ){
+ if ( AddFitKalmanDebugLevel > 0){
+ System.out.println ("Error: "
+ + RKDebug.Instance().getTrack() + " "
+ + RKDebug.Instance().getPropDir() + " "
+ );
+ System.out.println ("From surface 2: " + s_save );
+ System.out.println ("To surface 2: " + hit.surface());
+ System.out.println ("Failed addhit: " + nc + " " + nz + " " + nu + " " + thishit );
+ }
+ aida.histogram1D("/Bugs/Fit/Failed addhit from Surface",5,0,5).fill( from );
+ aida.histogram1D("/Bugs/Fit/Failed addhit to Surface",5,0,5).fill( to );
+ aida.cloud2D("/Bugs/Fit/Failed addhit to vs from Surface").fill( from, to );
+
+ }
+ if ( fstat>0 ) return 10000 + 1000*fstat + icount;
+
+
+ double chnew = trh.chisquared();
+ double dch = chnew - chold;
+
+ //System.out.printf("Error after addhit: %15.5f\n", dch);
+ //psm.Print(trh.newTrack());
+
+ /*
+ if( dch < -0.001 ){
+ System.out.println ("From surface 3: " + s_save );
+ System.out.println ("To surface 3: " + hit.surface());
+ }
+ */
+ chold = chnew;
+
+ // Save track for diagnostics.
+ VTUtil before = new VTUtil( trh.newTrack() );
+
+ // Apply multiple scattering and energy loss
+ int istat = interact( trh, hit, dir );
+
+ // Compute change in energy.
+ VTUtil after = new VTUtil( trh.newTrack() );
+ double de = after.e() - before.e();
+
+ sumde += de;
+
+ //SurfCylinder ss = (SurfCylinder)trh.newTrack().surface();
+ //System.out.printf ("Back dedx: %10.4f %12.8f %12.8f\n", ss.radius(), de, sumde );
+ ++icount;
+ }
+
+ // Propagate to the beampipe and interact.
+ if ( doMs ){
+
+ // Beampipe parameters for sid02
+ double radius = 1.22;
+ double l_over_radl = 0.006136;
+
+ SurfCylinder sbp = new SurfCylinder(radius);
+ PropStat pstat = trh.propagate( _pprop, sbp, dir );
+ if ( ! pstat.success() ) return -2;
+
+ int istat = interactonly( trh, radius, l_over_radl );
+
+ }
+
+ // Propagate to the PCAO.
+ SurfDCA sdca = new SurfDCA( 0., 0. );
+ PropStat pstat = trh.propagate( _pprop, sdca, dir );
+ if ( ! pstat.success() ) return -1;
+
+ //System.out.println ("Backwards fit sumde: " + sumde );
+ aida.cloud1D("Backward dedx check:").fill(sumde-RKDebug.Instance().getDeGen());
+
+ //System.out.println("Error at PCAO: ");
+ //psm.Print(trh.newTrack());
+
+ return 0;
+
+ }
+
+ private int interact ( HTrack trh, Hit hit, PropDir dir ){
+
+ if ( hit.surface().pureType().equals(SurfCylinder.staticType()) ){
+
+ SurfCylinder s = (SurfCylinder) hit.surface();
+ double r = s.radius();
+ if ( doMs ){
+ TrackError eold = trh.newTrack().error();
+
+ aida.histogram1D("/Bugs/Fit/Fit scat radius:",300,0.,150.).fill(r);
+
+ double l_over_radl = 0.;
+ if ( r < 1.3 ) {
+ l_over_radl = 0.006136;
+ }else if ( r < 10. ){
+ l_over_radl = 0.000916;
+ }else {
+ l_over_radl = 0.013747;
+ }
+
+ ThinCylMs scat = new ThinCylMs(l_over_radl*RKDebug.Instance().getMsFac());
+ ETrack et = trh.newTrack();
+ ETrack ets = new ETrack(et);
+ double chnew = trh.chisquared();
+
+ scat.interact(et);
+ hit.update(et);
+ trh.setFit(et,chnew);
+
+ /*
+ for ( int i=0; i<5; ++i ){
+ double ex1 = et.error().get(i,i);
+ double ex2 = ets.error().get(i,i);
+ double sigsq = ex1-ex2;
+ //double pull = -9.;
+ double sig=-1.;
+ if ( sigsq > 0. ){
+ sig = Math.sqrt(sigsq);
+ //pull = (et.vector(i)-ets.vector(i))/sig;
+ }
+ aida.cloud1D("/Bugs/Fit/Forward Fit Delta param:"+i).fill(et.vector(i)-ets.vector(i));
+ if ( sig > 0. ) aida.cloud1D("/Bugs/Fit/Forward Fit Delta error:"+i).fill(sig);
+ }
+ //System.out.println( "Error after: " + trh.newTrack().error().minus(eold) );
+ */
+ } // end if MS enabled
+
+
+ } // end CYlinder MS
+
+ if ( hit.surface().pureType().equals(SurfZPlane.staticType()) ){
+
+ SurfZPlane s = (SurfZPlane) hit.surface();
+ double z = s.z();
+ if ( doMs ){
+ TrackError eold = trh.newTrack().error();
+
+ aida.histogram1D("/Bugs/Fit/Fit scat z forward:",300,-150,150.).fill(z);
+
+ double l_over_radl = ( Math.abs(z)< 25) ? 0.000916 : 0.013747;
+ ThinZPlaneMs scat = new ThinZPlaneMs(l_over_radl*RKDebug.Instance().getMsFac());
+ ETrack et = trh.newTrack();
+ // ETrack ets = new ETrack(et);
+ double chnew = trh.chisquared();
+
+ scat.interact(et);
+ hit.update(et);
+ trh.setFit(et,chnew);
+
+ } // end if MS enabled.
+
+ } // end ZPlane MS
+
+ // Successful return;
+ return 0;
+ }
+
+
+ // There is no hit at this site so we only need to do the interaction, not update hits.
+ private int interactonly ( HTrack trh, double r, double l_over_radl ){
+
+ ThinCylMs scat = new ThinCylMs(l_over_radl*RKDebug.Instance().getMsFac());
+ ETrack et = trh.newTrack();
+ double chnew = trh.chisquared();
+
+ scat.interact(et);
+ trh.setFit(et,chnew);
+
+ // Successful return;
+ return 0;
+ }
+
+ /**
+ *output stream
+ *
+ * @return The String representation of this instance.
+ */
+ public String toString()
+ {
+ return getClass().getName();
+ }
+
+}
+