8 added files
lcsim/src/org/lcsim/contrib/RobKutschke/TRFSelfTest/util
diff -N CovSmear.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ CovSmear.java 8 Jun 2009 06:04:55 -0000 1.1
@@ -0,0 +1,74 @@
+package org.lcsim.contrib.RobKutschke.TRFSelfTest.util;
+
+import java.util.Random;
+
+import org.lcsim.recon.tracking.trfbase.TrackVector;
+import org.lcsim.recon.tracking.trfbase.TrackError;
+
+import Jama.Matrix;
+import Jama.EigenvalueDecomposition;
+
+/**
+ *
+ * Given a set of track parameters and a covariance matrix
+ * smear the track parameters by a multi-dimensional gaussian
+ * drawn from the covariance matrix.
+ *
+ * The method is to find the transformation that diagonalizes the
+ * covariance matrix, smear the diagonalized matrix and then transform
+ * the smearing back to the orignal basis.
+ *
+ *@author $Author: kutschke $
+ *@version $Id: CovSmear.java,v 1.1 2009/06/08 06:04:55 kutschke Exp $
+ *
+ * Date $Date: 2009/06/08 06:04:55 $
+ *
+ */
+
+public class CovSmear{
+
+ private Random ran = null;
+
+ public CovSmear( long seed ){
+ ran = new Random(seed);
+ }
+
+ public TrackVector Smear ( TrackError e, TrackVector v ){
+ Matrix m = new Matrix(e.matrix());
+
+ // Extract eigenvalues and the matrix of eigenvectors.
+ EigenvalueDecomposition eigen = new EigenvalueDecomposition(m);
+ double[] ev = eigen.getRealEigenvalues();
+ Matrix a = eigen.getV();
+
+ // Smear the diagonalized eigenvalue matrix.
+ int n = ev.length;
+ Matrix mout = new Matrix(n,1);
+ for ( int i=0; i<n; ++i ){
+ double x = ev[i];
+ if ( x < 0. ){
+ System.out.println("CovSmear Warning negative element on diag: "
+ + i + " "
+ + x
+ );
+ x = Math.abs(x);
+ }
+ double g = ran.nextGaussian();
+ x = Math.sqrt(x)*g;
+ mout.set(i,0,x);
+ }
+
+ // Rotate smeared vector to original basis.
+ Matrix m2 = a.times(mout);
+
+ // Apply smearing to the input vector.
+ TrackVector vout = new TrackVector();
+ for ( int i=0; i<n; ++i ){
+ vout.set(i,v.get(i)+m2.get(i,0));
+ }
+ return vout;
+
+ }
+
+
+}
lcsim/src/org/lcsim/contrib/RobKutschke/TRFSelfTest/util
diff -N PrintSymMatrix.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ PrintSymMatrix.java 8 Jun 2009 06:04:55 -0000 1.1
@@ -0,0 +1,79 @@
+package org.lcsim.contrib.RobKutschke.TRFSelfTest.util;
+
+import Jama.Matrix;
+import org.lcsim.recon.tracking.trfbase.TrackError;
+import org.lcsim.recon.tracking.trfbase.ETrack;
+
+/**
+ *
+ * Print a square matrix, assumed to be symmetric, as sqrt(diagonal elements)
+ * and correlation coefficients.
+ *
+ *@author $Author: kutschke $
+ *@version $Id: PrintSymMatrix.java,v 1.1 2009/06/08 06:04:55 kutschke Exp $
+ *
+ * Date $Date: 2009/06/08 06:04:55 $
+ *
+ */
+
+
+public class PrintSymMatrix{
+
+ public PrintSymMatrix(){
+ }
+
+ public void Print( TrackError e){
+ Print (e.getMatrix());
+ }
+
+ public void Print( ETrack tre){
+ System.out.printf ( "Params: %14.8e %14.8e %14.8e %14.8e %14.8e\n",
+ tre.vector(0),
+ tre.vector(1),
+ tre.vector(2),
+ tre.vector(3),
+ tre.vector(4)
+ );
+ Print (tre.error().getMatrix());
+ }
+
+ public void Print( Matrix eta, Matrix m){
+ System.out.printf ( "Params: %14.8e %14.8e %14.8e %14.8e %14.8e\n",
+ eta.get(0,0),
+ eta.get(1,0),
+ eta.get(2,0),
+ eta.get(3,0),
+ eta.get(4,0)
+ );
+ Print (m);
+ }
+
+
+ public void Print( Matrix m){
+ int dim = m.getColumnDimension();
+ if ( m.getRowDimension() != dim ) return;
+
+ double[] diag = new double[dim];
+ System.out.printf ( "Diag: " );
+ for ( int i=0; i<dim; ++i ){
+ if ( m.get(i,i) >= 0 ){
+ diag[i] = Math.sqrt(m.get(i,i));
+ }else{
+ diag[i] = -1.;
+ }
+ System.out.printf ( " %14.8e", diag[i] );
+ }
+ System.out.printf("\n");
+
+ for ( int i=0; i<dim; ++i ){
+ for ( int j=i+1; j<dim; ++j ){
+ double rho = -2.;
+ if ( diag[i] >0. && diag[j] >0. ){
+ rho = m.get(i,j)/diag[i]/diag[j];
+ }
+ System.out.printf ("%22.16f", rho );
+ }
+ System.out.printf("\n");
+ }
+ }
+}
lcsim/src/org/lcsim/contrib/RobKutschke/TRFSelfTest/util
diff -N RKDebug.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ RKDebug.java 8 Jun 2009 06:04:55 -0000 1.1
@@ -0,0 +1,108 @@
+package org.lcsim.contrib.RobKutschke.TRFSelfTest.util;
+
+import org.lcsim.recon.tracking.trfbase.PropDir;
+import org.lcsim.contrib.RobKutschke.TRFSelfTest.generator.RKTrack;
+
+/**
+ *
+ * A way to propagate debug info to low level routines.
+ *
+ *@author $Author: kutschke $
+ *@version $Id: RKDebug.java,v 1.1 2009/06/08 06:04:55 kutschke Exp $
+ *
+ * Date $Date: 2009/06/08 06:04:55 $
+ *
+ */
+
+public class RKDebug {
+
+ static private int track = -1;
+ static private PropDir pdir;
+ static private RKTrack rkt;
+ static private int start_type;
+
+ static private double msfac = 1.;
+ static private boolean doms = true;
+ static private boolean doeloss = true;
+
+ static private boolean printoutliers = false;
+
+ static private RKDebug instance = null;
+
+ // A "common block" to pass info around.
+ static private double degen = 0.;
+
+ private RKDebug(){
+ }
+
+ static public RKDebug Instance(){
+ if ( instance == null ){
+ instance = new RKDebug();
+ }
+ return instance;
+ }
+
+ static public int getTrack(){
+ return track;
+ }
+ static public void setTrack( int trk){
+ track = trk;
+ }
+
+ static public PropDir getPropDir(){
+ return pdir;
+ }
+ static public void setPropDir( PropDir dir){
+ pdir = dir;
+ }
+
+ static public RKTrack getRKTrack(){
+ return rkt;
+ }
+ static public void setRKTrack( RKTrack t){
+ rkt = t;
+ }
+
+ static public int getStartType(){
+ return start_type;
+ }
+ static public void setStartType( int type){
+ start_type = type;
+ }
+
+ static public double getMsFac(){
+ return msfac;
+ }
+ static public void setMsFac( double f){
+ msfac = f;
+ }
+
+ static public boolean getDoMs(){
+ return doms;
+ }
+ static public void setDoMs( boolean b){
+ doms = b;
+ }
+
+ static public boolean getDoEloss(){
+ return doeloss;
+ }
+ static public void setDoEloss( boolean b){
+ doeloss = b;
+ }
+
+ static public boolean getPrintOutliers(){
+ return printoutliers;
+ }
+ static public void setPrintOutliers( boolean b){
+ printoutliers = b;
+ }
+
+ static public double getDeGen(){
+ return degen;
+ }
+ static public void setDeGen( double b){
+ degen = b;
+ }
+
+}
lcsim/src/org/lcsim/contrib/RobKutschke/TRFSelfTest/util
diff -N RKZot.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ RKZot.java 8 Jun 2009 06:04:55 -0000 1.1
@@ -0,0 +1,96 @@
+package org.lcsim.contrib.RobKutschke.TRFSelfTest.util;
+
+import org.lcsim.recon.tracking.trfbase.ETrack;
+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.trffit.HTrack;
+
+import org.lcsim.recon.tracking.trfcyl.SurfCylinder;
+
+
+/**
+ *
+ * A utility used to debug precision problems with
+ * inwards going fits. It resets the z dependent parts
+ * of the covariance matrix of a track to the values
+ * that it had at the start of the fit.
+ *
+ * This is used to effectively do a 2-D circle fit instead
+ * of a 3D helix fit but throwing out the z information
+ * after each hit has been added.
+ *
+ * For this to work properly the track must start
+ * at the SurfCyl surface but I never put a check in the
+ * code for this.
+ *
+ * Usage:
+ * - at the start of a fit, instantiate a new RKZot object.
+ * - whenever you feel like throwing out the z information
+ * call the Zot method.
+ *
+ *@author $Author: kutschke $
+ *@version $Id: RKZot.java,v 1.1 2009/06/08 06:04:55 kutschke Exp $
+ *
+ * Date $Date: 2009/06/08 06:04:55 $
+ *
+ *
+ */
+
+
+public class RKZot{
+
+ static final int IPHI = SurfCylinder.IPHI;
+ static final int IZ = SurfCylinder.IZ;
+ static final int IALF = SurfCylinder.IALF;
+ static final int ITLM = SurfCylinder.ITLM;
+ static final int IQPT = SurfCylinder.IQPT;
+
+ TrackError e0;
+
+ public RKZot( HTrack h0 ){
+ ETrack et =h0.newTrack();
+
+ // Do I want to check that this is a SurfCyl?
+
+ e0 = et.error();
+
+ }
+
+ public void Zot( HTrack h ){
+ ETrack et = h.newTrack();
+
+ if ( !(et.surface() instanceof SurfCylinder) ){
+ System.out.println ("Error RKZot2: Not a cylinder ... " );
+ return;
+ }
+
+ TrackError err = et.error();
+
+ err.set( IZ, IZ, e0.get( IZ, IZ) );
+ err.set( ITLM, ITLM, e0.get(ITLM,ITLM) );
+
+ err.set( IZ, IPHI, 0. );
+ err.set( IZ, IALF, 0. );
+ err.set( IZ, ITLM, 0. );
+ err.set( IZ, IQPT, 0. );
+
+ err.set( IPHI, IZ, 0. );
+ err.set( IALF, IZ, 0. );
+ err.set( ITLM, IZ, 0. );
+ err.set( IQPT, IZ, 0. );
+
+ err.set( ITLM, IPHI, 0. );
+ err.set( ITLM, IALF, 0. );
+ err.set( ITLM, IQPT, 0. );
+
+ err.set( IPHI, ITLM, 0. );
+ err.set( IALF, ITLM, 0. );
+ err.set( IQPT, ITLM, 0. );
+
+ et.setError( err );
+ h.setFit( et, h.chisquared() );
+ }
+
+}
+
lcsim/src/org/lcsim/contrib/RobKutschke/TRFSelfTest/util
diff -N SurfaceCode.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ SurfaceCode.java 8 Jun 2009 06:04:55 -0000 1.1
@@ -0,0 +1,45 @@
+package org.lcsim.contrib.RobKutschke.TRFSelfTest.util;
+
+import org.lcsim.recon.tracking.trfbase.Surface;
+
+import org.lcsim.recon.tracking.trfcyl.SurfCylinder;
+import org.lcsim.recon.tracking.trfzp.SurfZPlane;
+import org.lcsim.recon.tracking.trfdca.SurfDCA;
+
+/**
+ *
+ * Define integer codes, and corresponding strings, to identify a type of surface.
+ *
+ *@author $Author: kutschke $
+ *@version $Id: SurfaceCode.java,v 1.1 2009/06/08 06:04:55 kutschke Exp $
+ *
+ * Date $Date: 2009/06/08 06:04:55 $
+ *
+ */
+
+
+public class SurfaceCode {
+
+ static final String[] namelist = { "Unknown", "Cylinder", "ZPlane", "DCA" };
+
+ private int code=0;
+
+ public SurfaceCode ( Surface s ){
+ if ( s.pureType().equals( SurfCylinder.staticType() ) ){
+ code = 1;
+ } else if ( s.pureType().equals( SurfZPlane.staticType() ) ){
+ code = 2;
+ } else if ( s.pureType().equals( SurfDCA.staticType() ) ){
+ code = 3;
+ }
+
+ }
+
+ public int getCode(){
+ return code;
+ }
+
+ public String getName(){
+ return namelist[code];
+ }
+}
lcsim/src/org/lcsim/contrib/RobKutschke/TRFSelfTest/util
diff -N VTUtil.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ VTUtil.java 8 Jun 2009 06:04:55 -0000 1.1
@@ -0,0 +1,166 @@
+package org.lcsim.contrib.RobKutschke.TRFSelfTest.util;
+
+import org.lcsim.recon.tracking.trfbase.ETrack;
+import org.lcsim.recon.tracking.trfbase.VTrack;
+import org.lcsim.recon.tracking.trfbase.Surface;
+
+import org.lcsim.recon.tracking.trfcyl.SurfCylinder;
+import org.lcsim.recon.tracking.trfzp.SurfZPlane;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.BasicHep3Vector;
+
+/**
+ *
+ * Compute various derived quantities from a VTrack or an ETrack.
+ *
+ *@author $Author: kutschke $
+ *@version $Id: VTUtil.java,v 1.1 2009/06/08 06:04:55 kutschke Exp $
+ *
+ * Date $Date: 2009/06/08 06:04:55 $
+ *
+ */
+
+
+public class VTUtil {
+
+ // Matches the precision in the TRF eloss code.
+ private double m_pi = 0.13957;
+
+ public VTUtil( VTrack v){
+ this.v = v;
+ Surface s = v.surface();
+
+ if ( s.pureType().equals( SurfCylinder.staticType()) ){
+ DoCylinder();
+ } else if ( s.pureType().equals( SurfZPlane.staticType()) ){
+ DoZPlane();
+ }
+ else{
+ DoNull();
+ }
+ }
+
+ public VTUtil( ETrack e){
+ v = new VTrack( e.surface(), e.vector() );
+ Surface s = v.surface();
+
+ if ( s.pureType().equals( SurfCylinder.staticType()) ){
+ DoCylinder();
+ } else if ( s.pureType().equals( SurfZPlane.staticType()) ){
+ DoZPlane();
+ }
+ else{
+ DoNull();
+ }
+ }
+
+
+ public double momentum(){
+ return p;
+ }
+
+ public double p(){
+ return p;
+ }
+
+ public double e(){
+ return e(m_pi);
+ }
+
+ public double e(double m){
+ return Math.sqrt( p()*p() +m*m );
+ }
+
+ public Hep3Vector asHep3Vector(){
+ double px = pt*Math.cos(phi);
+ double py = pt*Math.sin(phi);
+ double pz = p *costh;
+ return new BasicHep3Vector(px,py,pz);
+ }
+
+ public double costh(){
+ return costh;
+ }
+
+ public double sinth(){
+ return sinth;
+ }
+ public double q(){
+ return q;
+ }
+
+ public double z(){
+ return z;
+ }
+
+ public double r(){
+ return r;
+ }
+
+ private VTrack v = null;
+ private double costh;
+ private double sinth;
+ private double cotth;
+ private double pt;
+ private double p;
+ private double phi;
+ private double phi0;
+ private double q;
+
+ private double r;
+ private double z;
+
+ private void DoCylinder(){
+ cotth = v.vector(SurfCylinder.ITLM);
+ sinth = 1./Math.sqrt( 1. + cotth*cotth );
+ costh = cotth*sinth;
+ pt = 1./Math.abs(v.vector(SurfCylinder.IQPT));
+ p = pt/sinth;
+ phi = v.vector(SurfCylinder.IPHI)+v.vector(SurfCylinder.IALF);
+ q = Math.signum(v.vector(SurfCylinder.IQPT));
+
+ SurfCylinder cyl = (SurfCylinder) v.surface();
+ r = cyl.radius();
+ z = v.vector(SurfCylinder.IZ);
+ }
+
+ private void DoNull(){
+ cotth = 0.;
+ costh = 0.;
+ sinth = 0.;
+ pt = 0.;
+ p = 0.;
+ phi = 0.;
+ q = 1.;
+ r = 0.;
+ z = 0.;
+
+ }
+
+ private void DoZPlane(){
+ double xpr = v.vector(SurfZPlane.IDXDZ);
+ double ypr = v.vector(SurfZPlane.IDYDZ);
+
+ costh = 1./Math.sqrt( 1. + xpr*xpr + ypr*ypr);
+ SurfZPlane sz = (SurfZPlane) v.surface();
+
+ // This is a hack and needs to be done more generally.
+ if ( sz.z() < 0) costh = -costh;
+
+ sinth = Math.sqrt( 1. - costh*costh);
+ cotth = costh/sinth;
+ p = Math.abs(1./v.vector(SurfZPlane.IQP));
+ pt = p*sinth;
+ phi = Math.atan2( ypr, xpr );
+
+ q = Math.signum(v.vector(SurfZPlane.IQP));
+
+ double x = v.vector(SurfZPlane.IX);
+ double y = v.vector(SurfZPlane.IY);
+ r = Math.sqrt( x*x + y*y );
+ z = ((SurfZPlane)v.surface()).z();
+
+ }
+
+
+}
lcsim/src/org/lcsim/contrib/RobKutschke/TRFSelfTest/util
diff -N cyl_int.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ cyl_int.java 8 Jun 2009 06:04:55 -0000 1.1
@@ -0,0 +1,126 @@
+package org.lcsim.contrib.RobKutschke.TRFSelfTest.util;
+
+import java.lang.Math;
+
+/**
+ *
+ * Find the interesection of a track, specified by the CLEO paramters
+ * with the a cylinder of radius r, centered on the origin.
+ *
+ *@author $Author: kutschke $
+ *@version $Id: cyl_int.java,v 1.1 2009/06/08 06:04:55 kutschke Exp $
+ *
+ * Date $Date: 2009/06/08 06:04:55 $
+ *
+ */
+
+public class cyl_int{
+
+ public double x;
+ public double y;
+ public double psi;
+
+ public double k;
+ public double d0;
+ public double phi;
+ public double z;
+ public double t;
+
+ public double q;
+ public double rho;
+
+ public cyl_int( double k,
+ double d0,
+ double phi,
+ double z,
+ double t,
+ double r
+ ){
+
+ this.k = k;
+ this.d0 = d0;
+ this.phi = phi;
+ this.z = z;
+ this.t = t;
+
+ // Some other quantities derived from the track parameters.
+
+ q = (k > 0 ) ? 1.0 : -1.0;
+ rho = Math.abs( 0.5/k );
+
+
+ double u0 = Math.cos(phi);
+ double v0 = Math.sin(phi);
+ double x0 = -d0 * v0;
+ double y0 = d0 * u0;
+
+ double x0t = x0-q*rho*v0;
+ double y0t = y0+q*rho*u0;
+
+ double dpsi = Math.atan2(-y0t,x0t);
+
+ double sang = (r*r-x0t*x0t-y0t*y0t-rho*rho)/
+ (2.*q*rho*Math.sqrt(x0t*x0t+y0t*y0t));
+ if ( sang > 1.0 ){
+ sang = 1.0;
+ }else if ( sang < -1.0 ){
+ sang = -1.0;
+ }
+
+ double psi1=(Math.asin(sang)-phi-dpsi);
+ double psi2=(Math.PI-Math.asin(sang)-phi-dpsi);
+
+ double twopi = 2.*Math.PI;
+
+ // Put angles in the correct quadrant.
+ if (psi1 > twopi){
+ psi1 -= twopi;
+ }
+ if (psi1 < 0.0) {
+ psi1 += twopi;
+ }
+ if (psi1 < 0.0) {
+ psi1 += twopi;
+ }
+
+ if (psi2 > twopi) {
+ psi2 -= twopi;
+ }
+ if (psi2 < 0.0) {
+ psi2 += twopi;
+ }
+ if (psi2 < 0.0) {
+ psi2 += twopi;
+ }
+
+ psi = ( psi2 < psi1) ? psi2 : psi1;
+
+ x = x0t + q*rho*Math.sin(phi+q*psi);
+ y = y0t - q*rho*Math.cos(phi+q*psi);
+ }
+
+ public double getX(){
+ return x;
+ }
+
+ public double getY(){
+ return y;
+ }
+
+ public double getZ(){
+ double cz = t / Math.sqrt( 1. + t*t );
+ return z + rho*psi*t;
+ }
+
+
+ public double getPsi(){
+ return psi;
+ }
+
+
+}
+
+
+
+
+
lcsim/src/org/lcsim/contrib/RobKutschke/TRFSelfTest/util
diff -N zplane_int.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ zplane_int.java 8 Jun 2009 06:04:55 -0000 1.1
@@ -0,0 +1,88 @@
+package org.lcsim.contrib.RobKutschke.TRFSelfTest.util;
+
+import java.lang.Math;
+
+/**
+ *
+ * Find the interesection of a track, specified by the cleo paramters
+ * with a plane normal to the z axis and at the specfied z.
+ *
+ *@author $Author: kutschke $
+ *@version $Id: zplane_int.java,v 1.1 2009/06/08 06:04:55 kutschke Exp $
+ *
+ * Date $Date: 2009/06/08 06:04:55 $
+ *
+ */
+
+
+public class zplane_int{
+
+ public double x;
+ public double y;
+ public double z;
+ public double psi;
+
+ // Track parameters referenced to PCA (0,0,0);
+ public double k;
+ public double d0;
+ public double phi;
+ public double z0;
+ public double t;
+
+ public double q;
+ public double rho;
+
+ public zplane_int( double k,
+ double d0,
+ double phi,
+ double z0,
+ double t,
+ double zp
+ ){
+
+ this.k = k;
+ this.d0 = d0;
+ this.phi = phi;
+ this.z0 = z0;
+ this.t = t;
+
+ // Some other quantities derived from the track parameters.
+
+ q = (k > 0 ) ? 1.0 : -1.0;
+ rho = Math.abs( 0.5/k );
+
+ double rhopsi = (zp-z0)/t;
+ psi = rhopsi/rho;
+
+ double u0 = Math.cos(phi);
+ double v0 = Math.sin(phi);
+ double x0 = -d0 * v0;
+ double y0 = d0 * u0;
+
+ double xc = x0-q*rho*v0;
+ double yc = y0+q*rho*u0;
+
+ x = xc + q*rho*Math.sin(phi+q*psi);
+ y = yc - q*rho*Math.cos(phi+q*psi);
+ z = zp;
+ }
+
+ public double getX(){
+ return x;
+ }
+
+ public double getY(){
+ return y;
+ }
+
+ public double getZ(){
+ return z;
+ }
+
+
+ public double getPsi(){
+ return psi;
+ }
+
+
+}
CVSspam 0.2.8
