8 added files
lcsim/src/org/lcsim/contrib/garfield/tester
diff -N DocaTrackParameters.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ DocaTrackParameters.java 23 Feb 2006 19:10:09 -0000 1.1
@@ -0,0 +1,765 @@
+package org.lcsim.contrib.garfield.tester;
+
+import Jama.Matrix;
+
+import hep.physics.particle.Particle;
+
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+import org.lcsim.constants.Constants;
+
+
+/**
+ * MODIFIED COPY OF org.lcsim.mc.fast.tracking.DocaTrackParameters.
+ * (wish it was public)
+ *
+ * Holds DOCA parameters and error matrix of track. Can be initialized
+ * with a MC truth particle. <br>
+ *
+ * OriginalAuthor Tony Johnson, Wolfgang Walkowiak
+ * OriginalVersion DocaTrackParameters.java,v 1.5 2005/08/20 23:24:14 tonyj
+ *
+ * @version $Id: DocaTrackParameters.java,v 1.1 2006/02/23 19:10:09 onoprien Exp $
+ * InternalRelease 2.02.01 Feb 23, 2006
+ */
+class DocaTrackParameters
+{
+ private Hep3Vector m_pdoca_ref = null;
+ private Hep3Vector m_xdoca_ref = null;
+ private Hep3Vector m_xref = null;
+ private double[][] m_err = new double[5][5];
+ private double[] m_parm = new double[5];
+ private double m_Bz = 0.;
+ private double m_chi2 = -1.;
+ private double m_l0 = 0.;
+ private double m_l_ref = 0.;
+ private int m_ndf = 5;
+
+ //====================================================
+ //
+ // Constructors
+ //
+ //====================================================
+ DocaTrackParameters(double bField)
+ {
+ reset();
+ m_Bz = bField;
+ }
+
+ DocaTrackParameters(double bField, Particle p)
+ {
+ this(bField, p.getMomentum(), p.getOrigin(), p.getType().getCharge());
+ }
+
+ DocaTrackParameters(double bField, Hep3Vector momentum, Hep3Vector x, double q)
+ {
+ reset();
+ m_Bz = bField;
+ calculateDoca(momentum, x, q);
+ }
+
+ DocaTrackParameters(double bField, double[] momentum, double[] x, double q)
+ {
+ reset();
+ m_Bz = bField;
+ calculateDoca(momentum, x, q);
+ }
+
+ DocaTrackParameters(double bField, double[] momentum, double[] x, double q, double[][] errorMatrix)
+ {
+ this(bField, momentum, x, q);
+ fillErrorMatrix(errorMatrix);
+ }
+
+ DocaTrackParameters(double bField, double[] parameters)
+ {
+ m_Bz = bField;
+ m_parm = parameters;
+ }
+
+ DocaTrackParameters(double bField, double[] parameters, double[][] errorMatrix)
+ {
+ this(bField, parameters);
+ fillErrorMatrix(errorMatrix);
+ }
+
+ DocaTrackParameters(double bField, double[] parameters, double[][] errorMatrix, double chi2)
+ {
+ this(bField, parameters);
+ fillErrorMatrix(errorMatrix);
+ setChi2(chi2);
+ }
+
+ DocaTrackParameters(double bField, double[] parameters, double[][] errorMatrix, double chi2, int ndf)
+ {
+ this(bField, parameters);
+ fillErrorMatrix(errorMatrix);
+ setChi2(chi2);
+ setNDF(ndf);
+ }
+
+ public double getD0()
+ {
+ return m_parm[0];
+ }
+
+ /**
+ * Get the error matrix as a 2-D array
+ * @see #getTrackParameter
+ */
+ public double[][] getErrorMatrix()
+ {
+ return m_err;
+ }
+
+ /**
+ * Get an individual element of the error matrix
+ * @see #getTrackParameter
+ */
+ public double getErrorMatrixElement(int i, int j)
+ {
+ return (i < j) ? m_err[j][i] : m_err[i][j];
+ }
+
+ /**
+ * Get momentum at DOCA.
+ */
+ public double[] getMomentum()
+ {
+ return new double[] { getPX(), getPY(), getPZ() };
+ }
+
+ public double getOmega()
+ {
+ return m_parm[2];
+ }
+
+ public double getPX()
+ {
+ return getPt() * Math.cos(m_parm[1]);
+ }
+
+ public double getPY()
+ {
+ return getPt() * Math.sin(m_parm[1]);
+ }
+
+ public double getPZ()
+ {
+ return getPt() * m_parm[4];
+ }
+
+ public double getPhi0()
+ {
+ return m_parm[1];
+ }
+
+ /**
+ * Get total momentum at DOCA.
+ */
+ public double getPtot()
+ {
+ return Math.sqrt((getPX() * getPX()) + (getPY() * getPY()) + (getPZ() * getPZ()));
+ }
+
+ public double getTanL()
+ {
+ return m_parm[4];
+ }
+
+ /**
+ * Get an individual track parameter. <br>
+ *
+ * The track parameters for LCD are defined as follows
+ * <table>
+ * <tr><th>Index</th><th>Meaning</th></tr>
+ * <tr><td> 0 </td><td> d0 = XY impact parameter </td><tr>
+ * <tr><td> 1 </td><td> phi0 </td><tr> </td><tr>
+ * <tr><td> 2 </td><td> omega = 1/curv.radius (negative for negative tracks) </td><tr>
+ * <tr><td> 3 </td><td> z0 = z of track (z impact parameter) </td><tr>
+ * <tr><td> 4 </td><td> s = tan lambda </td><tr>
+ * </table>
+ * @param i The index of the track parameter
+ * @return The track parameter with the specified index
+ *
+ * All parameters are given at the DOCA.
+ */
+ public double getTrackParameter(int i)
+ {
+ return m_parm[i];
+ }
+
+ /**
+ * Get the track parameters as an array
+ * @see #getTrackParameter
+ */
+ public double[] getTrackParameters()
+ {
+ return m_parm;
+ }
+
+ /**
+ * Get the unit charge, ie +1, 0, -1.
+ */
+ public int getUnitCharge()
+ {
+ if (m_Bz != 0)
+ {
+ return (int) Math.signum(m_parm[2]);
+ }
+ else
+ {
+ return 0;
+ }
+ }
+
+ public double getZ0()
+ {
+ return m_parm[3];
+ }
+
+ /**
+ * Get cos(Theta) as calculated from the
+ * momentum vector at the DOCA. <br>
+ *
+ * Note: This is the same as getCosTheta()
+ */
+
+ public double magnitude()
+ {
+ return Math.sqrt(VecOp.dot(getDocaVec(), getDocaVec()));
+ }
+ public double magnitudeSquared()
+ {
+ return VecOp.dot(getDocaVec(), getDocaVec());
+ }
+
+ public double[] v()
+ {
+ return getDoca();
+ }
+
+ // IHep3Vector methods
+ public double x()
+ {
+ return getDocaX();
+ }
+
+ public double y()
+ {
+ return getDocaY();
+ }
+
+ public double z()
+ {
+ return getDocaZ();
+ }
+
+ /**
+ * Store the chi2.
+ */
+ void setChi2(double chi2)
+ {
+ m_chi2 = chi2;
+ }
+
+ /**
+ * Return the chi2 from smearing. <br>
+ *
+ * Note: The chi2 is to be calculated and
+ * stored by the smearing routine
+ * using setChi2().
+ */
+ double getChi2()
+ {
+ return m_chi2;
+ }
+
+ /**
+ * get cos(theta) at DOCA.
+ */
+ double getCosTheta()
+ {
+ return getPZ() / getPtot();
+ }
+
+ /**
+ * Get coordinates of DOCA.
+ */
+ double[] getDoca()
+ {
+ return new double[] { getDocaX(), getDocaY(), getDocaZ() };
+ }
+
+ double[] getDocaMomentum(double[] refPoint)
+ {
+ return getDocaMomentumVec(refPoint).v();
+ }
+
+ /**
+ * Calculate and get Doca momentum on track with respect to any space point.
+ */
+ Hep3Vector getDocaMomentumVec(Hep3Vector refPoint)
+ {
+ if ((refPoint.x() != 0.) && (refPoint.y() != 0.))
+ {
+ checkCalcDoca(refPoint);
+
+ return m_pdoca_ref;
+ }
+ else
+ {
+ return getMomentumVec();
+ }
+ }
+
+ Hep3Vector getDocaMomentumVec(double[] refPoint)
+ {
+ return getDocaMomentumVec(new BasicHep3Vector(refPoint[0], refPoint[1], refPoint[2]));
+ }
+
+ double[] getDocaPosition(double[] refPoint)
+ {
+ return getDocaPositionVec(refPoint).v();
+ }
+
+ //====================================================
+ //
+ // methods
+ //
+ //====================================================
+
+ /**
+ * Calculate and get Doca position on track with respect to
+ * any space point.
+ */
+ Hep3Vector getDocaPositionVec(Hep3Vector refPoint)
+ {
+ if ((refPoint.x() != 0.) && (refPoint.y() != 0.))
+ {
+ checkCalcDoca(refPoint);
+
+ return m_xdoca_ref;
+ }
+ else
+ {
+ return getDocaVec();
+ }
+ }
+
+ Hep3Vector getDocaPositionVec(double[] refPoint)
+ {
+ return getDocaPositionVec(new BasicHep3Vector(refPoint[0], refPoint[1], refPoint[2]));
+ }
+
+ /**
+ * Calculate and get path length on track for a doca to any space point
+ * in respect to the track defining doca (with respect to the origin).
+ * The length l is given in the transverse plane. <br>
+ * Use L = l*tan(lambda) to convert.
+ */
+ double getDocaTransversePathLength(Hep3Vector refPoint)
+ {
+ if ((refPoint.x() != 0.) && (refPoint.y() != 0))
+ {
+ checkCalcDoca(refPoint);
+
+ return m_l_ref;
+ }
+ else
+ {
+ return 0.;
+ }
+ }
+
+ double getDocaTransversePathLength(double[] refPoint)
+ {
+ return getDocaTransversePathLength(new BasicHep3Vector(refPoint[0], refPoint[1], refPoint[2]));
+ }
+
+ /**
+ * Get coordinates of DOCA.
+ */
+ Hep3Vector getDocaVec()
+ {
+ return new BasicHep3Vector(getDocaX(), getDocaY(), getDocaZ());
+ }
+
+ double getDocaX()
+ {
+ return (-m_parm[0] * Math.sin(m_parm[1]));
+ }
+
+ double getDocaY()
+ {
+ return (m_parm[0] * Math.cos(m_parm[1]));
+ }
+
+ double getDocaZ()
+ {
+ return (m_parm[3]);
+ }
+
+ /**
+ * Set the (transverse) path length l0 to original track vertex.
+ */
+ void setL0(double l0)
+ {
+ m_l0 = l0;
+ }
+
+ /**
+ * Get the (transverse) path length l0 to original track vertex.
+ */
+ double getL0()
+ {
+ return m_l0;
+ }
+
+ double[] getMomentum(double l)
+ {
+ return getMomentumVec(l).v();
+ }
+
+ /**
+ * Calculate and get momentum on track with respect to
+ * any path length l on track (l in xy plane).
+ */
+ Hep3Vector getMomentumVec(double l)
+ {
+ double phi0 = m_parm[1];
+ double omega = m_parm[2];
+ double tanl = m_parm[4];
+
+ int iq = getUnitCharge();
+
+ double phi = phi0 + (omega * l);
+ double pt = Constants.fieldConversion * iq * m_Bz / omega;
+
+ double px = pt * Math.cos(phi);
+ double py = pt * Math.sin(phi);
+ double pz = pt * tanl;
+
+ // System.out.println("l: "+l+" p: ("+px+", "+py+", "+pz+")");
+ return new BasicHep3Vector(px, py, pz);
+ }
+
+ Hep3Vector getMomentumVec()
+ {
+ return new BasicHep3Vector(getPX(), getPY(), getPZ());
+ }
+
+ /**
+ * Change the number degrees of freedom.
+ */
+ void setNDF(int ndf)
+ {
+ m_ndf = ndf;
+ }
+
+ /**
+ * Get the number degrees of freedom.
+ *
+ * Default is 5 unless changed with setNDF().
+ */
+ int getNDF()
+ {
+ return m_ndf;
+ }
+
+ double[] getPosition(double l)
+ {
+ return getPositionVec(l).v();
+ }
+
+ /**
+ * Calculate and get position on track with respect to
+ * any path length l on track (l in xy plane).
+ */
+ Hep3Vector getPositionVec(double l)
+ {
+ double d0 = m_parm[0];
+ double phi0 = m_parm[1];
+ double omega = m_parm[2];
+ double z0 = m_parm[3];
+ double tanl = m_parm[4];
+
+ double phi = phi0 + omega;
+ double rho = 1 / omega;
+
+ double x = (rho * Math.sin(phi)) - ((rho + d0) * Math.sin(phi0));
+ double y = (-rho * Math.cos(phi)) + ((rho + d0) * Math.cos(phi0));
+ double z = z0 + (l * tanl);
+
+ return new BasicHep3Vector(x, y, z);
+ }
+
+ /**
+ * Get transverse momentum at DOCA.
+ */
+ double getPt()
+ {
+ if (m_parm[2] != 0.)
+ {
+ return Math.abs(Constants.fieldConversion * m_Bz / m_parm[2]);
+ }
+ else
+ {
+ return 0.;
+ }
+ }
+
+ /**
+ * Get theta angle at DOCA.
+ */
+ double getTheta()
+ {
+ return Math.atan2(getPt(), getPZ());
+ }
+
+ /**
+ * Calculate the error matrix for the momentum
+ * for a point on the track specified by l.
+ * Result is given as a 3x3 array for the matrix.
+ */
+ double[][] calcMomentumErrorMatrix(double l)
+ {
+ double rho = 1. / getOmega();
+ double phi = getPhi0() + (getOmega() * l);
+ double tanl = getTanL();
+ double c = Constants.fieldConversion * Math.abs(m_Bz);
+ double sphi = Math.sin(phi);
+ double cphi = Math.cos(phi);
+
+ Matrix tMatrix = new Matrix(5, 3, 0.);
+ tMatrix.set(1, 0, -c * rho * sphi);
+ tMatrix.set(1, 1, c * rho * cphi);
+ tMatrix.set(2, 0, (-c * rho * rho * cphi) - (c * rho * l * sphi));
+ tMatrix.set(2, 1, (-c * rho * rho * sphi) + (c * rho * l * cphi));
+ tMatrix.set(2, 2, -c * rho * rho * tanl);
+ tMatrix.set(4, 2, c * rho);
+
+ Matrix errorMatrix = new Matrix(getErrorMatrix());
+ Matrix pErrorMatrix = tMatrix.transpose().times(errorMatrix.times(tMatrix));
+
+ return pErrorMatrix.getArrayCopy();
+ }
+
+ /**
+ * Calculate the error matrix for the position coordinates
+ * for a point on the track specified by l.
+ * Result is given as a 3x3 array for the matrix.
+ */
+ double[][] calcPositionErrorMatrix(double l)
+ {
+ double d0 = getD0();
+ double rho = 1. / getOmega();
+ double phi = getPhi0() + (getOmega() * l);
+ double sphi0 = Math.sin(getPhi0());
+ double cphi0 = Math.cos(getPhi0());
+ double sphi = Math.sin(phi);
+ double cphi = Math.cos(phi);
+
+ Matrix tMatrix = new Matrix(5, 3, 0.);
+ tMatrix.set(0, 0, -sphi0);
+ tMatrix.set(0, 1, cphi0);
+ tMatrix.set(1, 0, (rho * cphi) - ((rho + d0) * cphi0));
+ tMatrix.set(1, 1, (rho * sphi) - ((rho + d0) * sphi0));
+ tMatrix.set(2, 0, (rho * l * cphi) - (rho * rho * (sphi - sphi0)));
+ tMatrix.set(2, 1, (rho * l * sphi) + (rho * rho * (cphi - cphi0)));
+ tMatrix.set(3, 2, 1.);
+ tMatrix.set(4, 2, l);
+
+ Matrix errorMatrix = new Matrix(getErrorMatrix());
+
+ // MyContext.println(MyContext.getHeader());
+ // MyContext.printMatrix("Error matrix:",errorMatrix,10,15);
+ // MyContext.printMatrix("Transf matrix:",tMatrix,10,15);
+ Matrix xErrorMatrix = tMatrix.transpose().times(errorMatrix.times(tMatrix));
+
+ return xErrorMatrix.getArrayCopy();
+ }
+
+ void fillErrorMatrix(double[][] errorMatrix)
+ {
+ m_err = errorMatrix;
+
+ // fill upper triangle from lower triangle
+ for (int i = 0; i < 5; i++)
+ for (int j = 0; j < i; j++)
+ m_err[j][i] = m_err[i][j];
+ }
+
+ //====================================================
+ //
+ // private methods
+ //
+ //====================================================
+
+ /*
+ * Calculate the DOCA for a set of parameters with respect to the origin.
+ */
+ private void calculateDoca(double[] momentum, double[] trackPoint, double q)
+ {
+ Hep3Vector p = new BasicHep3Vector(momentum[0], momentum[1], momentum[2]);
+ Hep3Vector x = new BasicHep3Vector(trackPoint[0], trackPoint[1], trackPoint[2]);
+ calculateDoca(p, x, q);
+ }
+
+ /*
+ * Calculate the DOCA for a set of parameters in vectors with respect to the origin.
+ */
+ private void calculateDoca(Hep3Vector momentum, Hep3Vector trackPoint, double charge)
+ {
+ reset();
+
+ Hep3Vector xOrigin = new BasicHep3Vector(0., 0., 0.);
+
+ Hep3Vector[] result = calculateDoca(momentum, trackPoint, charge, xOrigin);
+ Hep3Vector xdoca = result[0];
+ Hep3Vector pdoca = result[1];
+ Hep3Vector dphdl = result[2];
+
+ int iq = (int) (charge / Math.abs(charge));
+ double pt = Math.sqrt((pdoca.x() * pdoca.x()) + (pdoca.y() * pdoca.y()));
+
+ // now calculate track parameters
+ double d0 = Math.sqrt((xdoca.x() * xdoca.x()) + (xdoca.y() * xdoca.y()));
+ if (VecOp.cross(xdoca, pdoca).z() > 0)
+ {
+ d0 = -d0;
+ }
+
+ double phi0 = Math.atan2(pdoca.y(), pdoca.x());
+ if (phi0 > Math.PI)
+ {
+ phi0 -= (2 * Math.PI);
+ }
+ if (phi0 < -Math.PI)
+ {
+ phi0 += (2 * Math.PI);
+ }
+
+ double omega = Constants.fieldConversion * iq * m_Bz / pt;
+ double tanl = pdoca.z() / pt;
+ double z0 = xdoca.z();
+
+ // now fill trackparameters
+ m_parm[0] = d0;
+ m_parm[1] = phi0;
+ m_parm[2] = omega;
+ m_parm[3] = z0;
+ m_parm[4] = tanl;
+
+ // save the distance to orignial track vertex
+ m_l0 = -dphdl.y();
+
+ // System.out.println("DocaTrackParameters: xdoca = ("+
+ // xdoca.x()+", "+xdoca.y()+", "+xdoca.z()+")");
+ // System.out.println("DocaTrackParameters: pdoca = ("+
+ // pdoca.x()+", "+pdoca.y()+", "+pdoca.z()+")");
+ // System.out.println("DocaTrackParameters: d0: "+m_parm[0]+
+ // " phi0: "+m_parm[1]+" omega: "+m_parm[2]+
+ // " z0: "+m_parm[3]+" tanl: "+m_parm[4]);
+ // System.out.println("DocaTrackParameters: m_l0 = "+m_l0);
+ }
+
+ /*
+ * Calculate DOCA position and momentum vectors with respect to any
+ * space point.
+ */
+ private Hep3Vector[] calculateDoca(Hep3Vector momentum, Hep3Vector trackPoint, double charge, Hep3Vector refPoint)
+ {
+ // subtract refPoint
+ Hep3Vector xp = VecOp.sub(trackPoint, refPoint);
+
+ int iq = (int) (charge / Math.abs(charge));
+ double pt = Math.sqrt((momentum.x() * momentum.x()) + (momentum.y() * momentum.y()));
+ double tanl = momentum.z() / pt;
+ double rho = pt / (m_Bz * Constants.fieldConversion);
+
+ BasicHep3Vector xdoca;
+ Hep3Vector pdoca;
+ Hep3Vector dphdl;
+
+ // System.out.println("calculateDoca: m_flip: "+m_flip+" iq: "+iq);
+ if (xp.magnitude() > 0.) // no need for calculation if xp = (0,0,0) !
+ {
+ // calculate position and momentum at doca
+ Hep3Vector nzv = new BasicHep3Vector(0., 0., iq);
+ Hep3Vector xxc = new BasicHep3Vector(VecOp.cross(momentum, nzv).x(), VecOp.cross(momentum, nzv).y(), 0.);
+ Hep3Vector nxxc = VecOp.unit(xxc);
+ BasicHep3Vector xc = (BasicHep3Vector) VecOp.add(xp, VecOp.mult(rho, nxxc));
+ xc.setV(xc.x(), xc.y(), 0.);
+
+ Hep3Vector nxc = VecOp.unit(xc);
+
+ BasicHep3Vector catMC = (BasicHep3Vector) VecOp.cross(nzv, nxc);
+ catMC.setV(catMC.x(), catMC.y(), 0.);
+
+ Hep3Vector ncatMC = VecOp.unit(catMC);
+
+ pdoca = new BasicHep3Vector(pt * ncatMC.x(), pt * ncatMC.y(), momentum.z());
+
+ double dphi = Math.asin(VecOp.cross(nxxc, nxc).z());
+ double dl = -dphi * rho * iq;
+
+ xdoca = (BasicHep3Vector) VecOp.add(xc, VecOp.mult(-rho, nxc));
+ xdoca.setV(xdoca.x(), xdoca.y(), xp.z() + (dl * tanl));
+
+ // save dphi and dl
+ dphdl = new BasicHep3Vector(dphi, dl, 0.);
+ }
+ else
+ {
+ xdoca = (BasicHep3Vector) xp;
+ pdoca = momentum;
+ dphdl = new BasicHep3Vector();
+ }
+
+ // add refPoint back in again
+ xdoca = (BasicHep3Vector) VecOp.add(xdoca, refPoint);
+
+ return new Hep3Vector[] { xdoca, pdoca, dphdl };
+ }
+
+ /*
+ * Calculate Doca for this track with respect to any space point,
+ * if this has not been done before.
+ */
+ private void checkCalcDoca(Hep3Vector refPoint)
+ {
+ // hassle with calculation only, if not done before yet!
+ if ((m_xref == null) || (refPoint.x() != m_xref.x()) || (refPoint.y() != m_xref.y()) || (refPoint.z() != m_xref.z()))
+ {
+ m_xref = refPoint;
+
+ // find doca vectors
+ Hep3Vector xdoca = getDocaVec();
+ Hep3Vector pdoca = getMomentumVec();
+ int iq = getUnitCharge();
+
+ // calculate doca to refPoint
+ Hep3Vector[] result = calculateDoca(pdoca, xdoca, (double) iq, refPoint);
+ m_xdoca_ref = result[0];
+ m_pdoca_ref = result[1];
+
+ // distance along track to original point
+ m_l_ref = result[2].y();
+ }
+ }
+
+ private void reset()
+ {
+ for (int i = 0; i < 5; i++)
+ {
+ m_parm[i] = 0;
+ for (int j = 0; j < 5; j++)
+ {
+ m_err[i][j] = 0;
+ }
+ }
+ m_l0 = 0.;
+ }
+}
lcsim/src/org/lcsim/contrib/garfield/tester
diff -N ExampleDriver.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ExampleDriver.java 23 Feb 2006 19:10:10 -0000 1.1
@@ -0,0 +1,38 @@
+package org.lcsim.contrib.garfield.tester;
+
+import java.util.*;
+import org.lcsim.event.*;
+import org.lcsim.recon.cluster.nn.NearestNeighborClusterDriver;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+import org.lcsim.contrib.garfield.*;
+
+/**
+ * An example driver showing how to use calorimeter-assisted track finder and
+ * the tester package.
+ *
+ * @author D. Onoprienko
+ * @version $Id: ExampleDriver.java,v 1.1 2006/02/23 19:10:10 onoprien Exp $
+ * InternalRelease 2.02.01 Feb 23, 2006
+ */
+public class ExampleDriver extends Driver {
+
+ private AIDA aida = AIDA.defaultInstance();
+
+ public ExampleDriver() {
+
+ add(new NearestNeighborClusterDriver(2, 2, 1, 2, 0.));
+ add(new EmcalMipStubs("GarfieldMipStubs",0));
+ add(new GarfieldTrackFinder(0));
+
+ TrackingTest test = new TrackingTest("Example test");
+ test.set(TrackingTest.Parameter.TRACK_COLLECTION_NAME, "GarfieldTracks");
+ add(test);
+
+ TrackingTest testKs = new TrackingTestKs("Example KS0 test");
+ testKs.set(TrackingTest.Parameter.TRACK_COLLECTION_NAME, "GarfieldTracks");
+ testKs.setPtCut(0.5*Const.GeV);
+ add(testKs);
+ }
+
+}
\ No newline at end of file
lcsim/src/org/lcsim/contrib/garfield/tester
diff -N MCTrack.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MCTrack.java 23 Feb 2006 19:10:10 -0000 1.1
@@ -0,0 +1,185 @@
+package org.lcsim.contrib.garfield.tester;
+
+import java.util.*;
+import org.lcsim.event.*;
+import org.lcsim.contrib.garfield.*;
+
+
+/**
+ * This class represents a Monte Carlo track - Monte Carlo particle that produced
+ * hits in a tracker. Provides access to the underlying {@link MCParticle} object
+ * and to associated tracker hits and reconstructed tracks.
+ * <p>
+ * An object of this class can be assigned any number of integer status codes through
+ * the call to {@link #setStatus(int,int) setStatus(int category, int status)} method.
+ * This allows users to mark Monte Carlo tracks that belong to a certain group, for example.
+ * Assigned statuses are retrieved by {@link #getStatus(int) getStatus(int category)}.
+ *
+ * @author D. Onoprienko
+ * @version $Id: MCTrack.java,v 1.1 2006/02/23 19:10:10 onoprien Exp $
+ * InternalRelease 2.02.01 Feb 23, 2006
+ */
+public class MCTrack implements Track {
+
+// -- Private fields : --------------------------------------------------------
+
+ private MCParticle _mcParticle;
+ private ArrayList<TrackerHit> _trackerHits;
+ private ArrayList<Track> _tracks;
+ private boolean _hitsSorted;
+
+ HashMap<Integer,Integer> _status;
+
+ static private double[] _cacheTrackParameters;
+ static private MCTrack _trackInCache;
+
+// -- Constructors : ----------------------------------------------------------
+
+ public MCTrack(MCParticle mcParticle) {
+ _mcParticle = mcParticle;
+ _trackerHits = new ArrayList<TrackerHit>(10);
+ _tracks = new ArrayList<Track>(1);
+ _hitsSorted = false;
+ _status = new HashMap<Integer,Integer>(2);
+ }
+
+// -- Getters : ---------------------------------------------------------------
+
+ /** Returns MCParticle that produced this track. */
+ public MCParticle getMCParticle() {return _mcParticle;}
+
+ /** Returns a list of associated reconstructed tracks. */
+ public List<Track> getTracks() {return _tracks;}
+
+ /** Returns a list of associated tracker hits. */
+ public List<TrackerHit> getTrackerHits() {
+ if (!_hitsSorted) {
+ sortHits();
+ _hitsSorted = true;
+ }
+ return _trackerHits;
+ }
+
+ /**
+ * Returns a status flag that has been previously assigned to this track through
+ * the call to {@link #setStatus(int,int) setStatus(int category, int status)}.
+ */
+ public int getStatus(int category) {
+ return _status.get(category);
+ }
+
+ /** Returns transverse momentum of the track. */
+ public double getPt() {
+ double[] p = getMomentum();
+ return Math.sqrt(p[0]*p[0]+p[1]*p[1]);
+ }
+
+ /** Returns momentum of the track. */
+ public double getP() {
+ double[] p = getMomentum();
+ return Math.sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
+ }
+
+
+// -- Setters : ---------------------------------------------------------------
+
+ /** Add track to the list of reconstructed tracks associated with this Monte Carlo track. */
+ public void addTrack(Track track) {_tracks.add(track);}
+
+ /** Add hit to the list of tracker hits associated with this Monte Carlo track. */
+ public void addHit(TrackerHit hit) {_trackerHits.add(hit);}
+
+ /** Add hits to the list of tracker hits associated with this Monte Carlo track. */
+ public void addHits(List<TrackerHit> hitList) {
+ for (TrackerHit hit : hitList) if (!_trackerHits.contains(hit)) _trackerHits.add(hit);
+ }
+
+ /**
+ * Set an integer status flag for this track.
+ *
+ * @param category Identifier of a status category. A track can be assigned
+ * statuses in any number of categories.
+ * @param status Status flag that can be later retrieved through the call to
+ * {@link #getStatus(int) getStatus(int category)}
+ */
+ public void setStatus(int category, int status) {_status.put(category, status);}
+
+// -- Implement org.lcsim.event.Track : ----------------------------------------
+
+ public double getTrackParameter(int param) {
+ return getTrackParameters()[param];
+ }
+
+ public double[] getTrackParameters() {
+ if (_trackInCache != this) {
+ _cacheTrackParameters = (new DocaTrackParameters(Const.bField, _mcParticle)).getTrackParameters();
+ _trackInCache = this;
+ }
+ return _cacheTrackParameters;
+ }
+
+ public boolean isReferencePointPCA() {return false;}
+ public double getReferencePointX() {return _mcParticle.getOriginX();}
+ public double getReferencePointY() {return _mcParticle.getOriginY();}
+ public double getReferencePointZ() {return _mcParticle.getOriginZ();}
+ public double[] getReferencePoint() {return _mcParticle.getOrigin().v();}
+
+ public double getdEdx() {return 0.;}
+ public double getdEdxError() {return 0.;}
+ public int getType() {return 0;}
+ public boolean fitSuccess() {return true;}
+ public double getChi2() {return 0.;}
+ public int getNDF() {return 0;}
+ public double getErrorMatrixElement(int param, int param1) {return 0.;}
+ public double[][] getErrorMatrix() {return new double[1][1];}
+
+ public int[] getSubdetectorHitNumbers() {
+ return new int[]{0,0,0,0};
+ }
+
+ public int getCharge() {
+ double charge = _mcParticle.getCharge();
+ int iCharge;
+ if (charge < -0.7) iCharge = -1;
+ else if (charge > 0.7) iCharge = 1;
+ else iCharge = 0;
+ return iCharge;
+ }
+
+ public double[] getMomentum() {return _mcParticle.getMomentum().v();}
+ public double getPX() {return _mcParticle.getPX();}
+ public double getPY() {return _mcParticle.getPY();}
+ public double getPZ() {return _mcParticle.getPZ();}
+
+ public double getRadiusOfInnermostHit() {
+ double rad = 0.;
+ if (!_trackerHits.isEmpty()) {
+ if (!_hitsSorted) {
+ sortHits();
+ _hitsSorted = true;
+ }
+ double[] pos = _trackerHits.get(0).getPosition();
+ rad = Math.hypot(pos[0], pos[1]);
+ }
+ return rad;
+ }
+
+// -- Private helper methods : ------------------------------------------------
+
+ /**
+ * Sort hits in the order of increasing distance from the detector center.
+ */
+ private void sortHits() {
+ Collections.sort(_trackerHits, new Comparator<TrackerHit>() {
+ public int compare(TrackerHit hit1, TrackerHit hit2) {
+ double[] pos = hit1.getPosition();
+ double r1 = pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2];
+ pos = hit2.getPosition();
+ double r2 = pos[0]*pos[0]+pos[1]*pos[1]+pos[2]*pos[2];
+ return (int)Math.signum(r1-r2);
+ }
+ });
+
+ }
+
+}
lcsim/src/org/lcsim/contrib/garfield/tester
diff -N RatedTrack.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ RatedTrack.java 23 Feb 2006 19:10:10 -0000 1.1
@@ -0,0 +1,113 @@
+package org.lcsim.contrib.garfield.tester;
+
+import java.util.*;
+import org.lcsim.event.*;
+import org.lcsim.contrib.garfield.*;
+
+/**
+ * This class wraps a {@link Track}, providing access to the underlying reconstructed
+ * track, associated Monte Carlo track, and allowing the user to associate any number
+ * of integer status flags with the track.
+ *
+ * @author D. Onoprienko
+ * @version $Id: RatedTrack.java,v 1.1 2006/02/23 19:10:10 onoprien Exp $
+ * InternalRelease 2.02.01 Feb 23, 2006
+ */
+public class RatedTrack implements Track {
+
+// -- Constructors : ----------------------------------------------------------
+
+ /** Create RatedTrack by wrapping a reconstructed track. */
+ public RatedTrack(Track track) {
+ _track = track;
+ _status = new HashMap<Integer,Integer>(2);
+ }
+
+// -- Getters : ---------------------------------------------------------------
+
+ /** Returns underlying reconstructed Track object. */
+ public Track getTrack() {return _track;}
+
+ /** Returns associated Monte Carlo track. */
+ public MCTrack getMCTrack() {return _mcTrack;}
+
+ /**
+ * Returns a status flag that has been previously assigned to this track through
+ * the call to {@link #setStatus(int,int) setStatus(int category, int status)}.
+ */
+ public int getStatus(int category) {return _status.get(category);}
+
+ /** Returns momentum of the track. */
+ public double getP() {
+ double mom;
+ try {
+ mom = ((GarfieldTrack)_track).getP(Const.bField);
+ } catch (Exception e) {
+ double p[] = _track.getMomentum();
+ mom = Math.sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
+ }
+ return mom;
+ }
+
+ /** Returns transverse momentum of the track. */
+ public double getPt() {
+ double pt;
+ try {
+ pt = ((GarfieldTrack)_track).getPt(Const.bField);
+ } catch (Exception e) {
+ double p[] = _track.getMomentum();
+ pt = Math.sqrt(p[0]*p[0]+p[1]*p[1]);
+ }
+ return pt;
+ }
+
+// -- Setters : ---------------------------------------------------------------
+
+ /** Associate this track with a Monte Carlo track. */
+ public void setMCTrack(MCTrack mcTrack) {_mcTrack = mcTrack;}
+
+ /**
+ * Set an integer status flag for this track.
+ *
+ * @param category Identifier of a status category. A track can be assigned
+ * statuses in any number of categories.
+ * @param status Status flag that can be later retrieved through the call to
+ * {@link #getStatus(int) getStatus(int category)}
+ */
+ public void setStatus(int category, int status) {_status.put(category, status);}
+
+// -- Implementing Track: forward calls to associated reconstructed track : ---
+
+ public double getTrackParameter(int param) {return _track.getTrackParameter(param);}
+ public boolean isReferencePointPCA() {return _track.isReferencePointPCA();}
+ public double getdEdxError() {return _track.getdEdxError();}
+ public double getdEdx() {return _track.getdEdx();}
+ public int getType() {return _track.getType();}
+ public List<Track> getTracks() {return _track.getTracks();}
+ public List<TrackerHit> getTrackerHits() {return _track.getTrackerHits();}
+ public double[] getTrackParameters() {return _track.getTrackParameters();}
+ public int[] getSubdetectorHitNumbers() {return _track.getSubdetectorHitNumbers();}
+ public boolean fitSuccess() {return _track.fitSuccess();}
+ public int getCharge() {return _track.getCharge();}
+ public double getChi2() {return _track.getChi2();}
+ public double[][] getErrorMatrix() {return _track.getErrorMatrix();}
+ public double getErrorMatrixElement(int param, int param1) {return _track.getErrorMatrixElement(param, param1);}
+ public double[] getMomentum() {return _track.getMomentum();}
+ public int getNDF() {return _track.getNDF();}
+ public double getPX() {return _track.getPX();}
+ public double getPY() {return _track.getPY();}
+ public double getPZ() {return _track.getPZ();}
+ public double getRadiusOfInnermostHit() {return _track.getRadiusOfInnermostHit();}
+ public double[] getReferencePoint() {return _track.getReferencePoint();}
+ public double getReferencePointX() {return _track.getReferencePointX();}
+ public double getReferencePointY() {return _track.getReferencePointY();}
+ public double getReferencePointZ() {return _track.getReferencePointZ();}
+
+// -- Private data : ----------------------------------------------------------
+
+ Track _track;
+
+ MCTrack _mcTrack;
+ HashMap<Integer,Integer> _status;
+
+}
lcsim/src/org/lcsim/contrib/garfield/tester
diff -N TrackingTest.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TrackingTest.java 23 Feb 2006 19:10:10 -0000 1.1
@@ -0,0 +1,122 @@
+package org.lcsim.contrib.garfield.tester;
+
+import java.util.*;
+import java.lang.reflect.*;
+import org.lcsim.event.*;
+import org.lcsim.event.MCParticle.SimulatorStatus;
+import org.lcsim.event.EventHeader.LCMetaData;
+import org.lcsim.util.Driver;
+import hep.physics.particle.Particle;
+import org.lcsim.util.aida.AIDA;
+import hep.aida.ITree;
+
+import javax.swing.*;
+
+
+/**
+ * A general example of track finding performance test.
+ *
+ * @author D. Onoprienko
+ * @version $Id: TrackingTest.java,v 1.1 2006/02/23 19:10:10 onoprien Exp $
+ * InternalRelease 2.02.01 Feb 23, 2006
+ */
+public class TrackingTest extends TrackingTestBase {
+
+// -- Indices into the statistics (per event and cumulative) arrays : ---------
+
+ protected final int _I_nMCTracks;
+ protected final int _I_eMCTracks;
+ protected final int _I_nRecoTracks;
+ protected final int _I_nRecoTracksFake;
+ protected final int _I_nMCTracksReconstructed;
+ protected final int _I_eMCTracksReconstructed;
+
+// -- Constructors : ----------------------------------------------------------
+
+ public TrackingTest(String name) {
+ this();
+ set(Parameter.TEST_NAME, name);
+ }
+
+ public TrackingTest() {
+ _I_nMCTracks = createIndex("Number of MC Tracks in class");
+ _I_eMCTracks = createIndex("Total momentum of MC Tracks in class", IndexType.NOPRINT);
+ _I_nMCTracksReconstructed = createIndex("Number of MC Tracks in class that have been reconstructed");
+ _I_eMCTracksReconstructed = createIndex("Total momentum of MC Tracks in class that have been reconstructed", IndexType.NOPRINT);
+ _I_nRecoTracks = createIndex("Number of reconstructed tracks", IndexType.NOPRINT);
+ _I_nRecoTracksFake = createIndex("Total number of fakes", IndexType.PRINTEVENT);
+
+ _ptCut = 0.;
+ }
+
+// -- Setters : ---------------------------------------------------------------
+
+ public void setPtCut(double ptCut) { _ptCut = ptCut;}
+
+// -- Processing event : ------------------------------------------------------
+
+ protected void analyzeEvent(EventHeader event) {
+
+ // Number and momentum of tracks in class and those of them that have been reconstructed
+
+ for (MCTrack mcTrack : _mcTrackList) {
+ if (mcTrack.getStatus(IN_CLASS) == YES) {
+ _statE[_I_nMCTracks] += 1.;
+ _statE[_I_eMCTracks] += mcTrack.getP();
+ if (mcTrack.getTracks().size() > 0) {
+ _statE[_I_nMCTracksReconstructed] += 1.;
+ _statE[_I_eMCTracksReconstructed] += mcTrack.getP();
+ }
+ }
+ }
+
+ // Number of reconstracted tracks and fakes
+
+ _statE[_I_nRecoTracks] = _ratedTrackList.size();
+ for (RatedTrack ratedTrack : _ratedTrackList) {
+ if (ratedTrack.getMCTrack() == null) _statE[_I_nRecoTracksFake] += 1.;
+ }
+
+ // Fill histograms for reconstructed tracks in class :
+
+ for (RatedTrack rTrack : _ratedTrackList) {
+ MCTrack mcTrack = rTrack.getMCTrack();
+ if (mcTrack != null && mcTrack.getStatus(IN_CLASS) == YES) {
+ double p = mcTrack.getP();
+ double dp = rTrack.getP()-p;
+ _aida.cloud1D(_testName + " deltaP over P").fill(dp/p);
+ p = mcTrack.getPt();
+ dp = rTrack.getPt()-p;
+ _aida.cloud1D(_testName + " deltaPt over Pt").fill(dp/p);
+ }
+ }
+
+ // Print statistics
+
+ printStatistics(event);
+ }
+
+
+ protected boolean mcTrackFilter(MCTrack mcTrack) {
+ MCParticle mcPart = mcTrack.getMCParticle();
+ if (mcPart.getGeneratorStatus() != Particle.FINAL_STATE) return false;
+ if (Math.abs(mcPart.getCharge()) < .8) return false;
+ if (mcTrack.getPt() < _ptCut) return false;
+ return true;
+ }
+
+ protected void suspend() {
+ String[] messages = new String[] {
+ "Reconstruction efficiency for tracks in class: " + _statC[_I_nMCTracksReconstructed]/_statC[_I_nMCTracks],
+ "Fraction of reconstructed momentum for tracks in class: " + _statC[_I_eMCTracksReconstructed]/_statC[_I_eMCTracks],
+ "Purity (all tracks): " + (1.-_statC[_I_nRecoTracksFake]/_statC[_I_nRecoTracks])
+ };
+ printCumulativeStatistics(messages);
+ }
+
+ protected void endOfData() {suspend();}
+
+// -- Private fields : --------------------------------------------------------
+
+ double _ptCut;
+}
lcsim/src/org/lcsim/contrib/garfield/tester
diff -N TrackingTestBase.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TrackingTestBase.java 23 Feb 2006 19:10:10 -0000 1.1
@@ -0,0 +1,448 @@
+package org.lcsim.contrib.garfield.tester;
+
+import java.util.*;
+import java.lang.reflect.*;
+import javax.swing.*;
+import org.lcsim.event.*;
+import org.lcsim.event.MCParticle.SimulatorStatus;
+import org.lcsim.event.EventHeader.LCMetaData;
+import org.lcsim.util.Driver;
+import hep.physics.particle.Particle;
+import org.lcsim.util.aida.AIDA;
+import hep.aida.ITree;
+
+
+/**
+ * Abstract class that provides infrastructure and utility methods for testing track finders.
+ * See {@link TrackingTest} and {@link TrackingTestKs} for examples of use.
+ * <p>
+ * Classes implementing particular tests should extend this class, overwriting
+ * {@link #analyzeEvent(EventHeader)} and, if necessary,
+ * {@link #eventFilter(EventHeader)},
+ * {@link #matchTracks(Track track, MCTrack mcTrack)},
+ * {@link #mcTrackFilter(MCTrack mcTrack)},
+ * {@link #suspend()}, and {@link #endOfData()} methods.
+ * <p>
+ * Subclusses extending this driver should be added to the event processing chain
+ * after track findind stage. When the driver's {@link #process(EventHeader event)}
+ * method is called by the framework, it creates a list of Monte Carlo tracks,
+ * associates hits and reconstructed tracks with them, and then calls
+ * {@link #analyzeEvent(EventHeader)} method that should be implemented by subclasses.
+ * <p>
+ * Only events for which {@link #eventFilter(EventHeader)} returns <code>true</code>
+ * are processed. Reconstructed track is associated with Monte Carlo track if
+ * {@link #matchTracks(Track track, MCTrack mcTrack)} returns <code>true</code>.
+ * Monte Carlo tracks for which {@link #mcTrackFilter(MCTrack mcTrack)} returns
+ * <code>true</code> are marked as being in the requested class (subsequent calls to
+ * {@link MCTrack#getStatus(int) mcTrack.getStatus(TrackingTestBase.IN_CLASS)}
+ * will return <code>TrackingTestBase.YES</code>. Default implementations are provided
+ * for these methods, but subclasses will need to overwrite some or all of them to
+ * perform a meaningful test.
+ * <p>
+ * Any number of variables for accumulating statistics can be created with
+ * {@link #createIndex(String name, IndexType... args)}. Data put into
+ * these variables is automatically accumulated at the end of event processing,
+ * and can be printed out with {@link #printStatistics(EventHeader event)} and
+ * {@link #printCumulativeStatistics()}.
+ * <p>
+ * A few utility methods that might be useful in testing tracking performance are also provided.
+ *
+ * @author D. Onoprienko
+ * @version $Id: TrackingTestBase.java,v 1.1 2006/02/23 19:10:10 onoprien Exp $
+ * InternalRelease 2.02.01 Feb 23, 2006
+ */
+abstract public class TrackingTestBase extends Driver {
+
+// -- Constants defining status categories and values: ------------------------
+
+ public static final int IN_CLASS = -1;
+ public static final int FAKE = -2;
+
+ public static final int YES = -1;
+ public static final int NO = -2;
+
+// -- Constructors : ----------------------------------------------------------
+
+ public TrackingTestBase() {
+ _testName = "";
+ _aida = AIDA.defaultInstance();
+ _tree = _aida.tree();
+ _trackCollectionName = "";
+ _trackerHitCollectionName = "";
+ _nEvents = 0;
+ }
+
+ // -- Setters : -------------------------------------------------------------
+
+ public enum Parameter {
+ TEST_NAME,
+ TRACK_COLLECTION_NAME,
+ HIT_COLLECTION_NAME
+ }
+
+ /**
+ * Set test parameters.
+ * @param par Indicates what parameter is being changed.
+ * @param value Value to be assigned to this parameter.
+ */
+ public void set(Parameter par, String value) {
+ switch (par) {
+ case TEST_NAME:
+ _testName = value;
+ break;
+ case TRACK_COLLECTION_NAME:
+ _trackCollectionName = value;
+ break;
+ case HIT_COLLECTION_NAME:
+ _trackerHitCollectionName = value;
+ break;
+ default: throw new Error("No such parameter: " + par);
+ }
+ }
+
+// -- Processing event : ------------------------------------------------------
+
+ /** Process event. */
+ public final void process(EventHeader event) {
+
+ // Skip undesired events
+
+ if (!eventFilter(event)) return;
+
+ // Reset "per event" statistics
+
+ if (_statE != null) Arrays.fill(_statE, 0.);
+
+ // Create list of Monte Carlo Tracks
+
+ List<List<TrackerHit>> hitCollections;
+ List<TrackerHit> hitList;
+ if (_trackerHitCollectionName == "") {
+ hitCollections = event.get(TrackerHit.class);
+ hitList = new ArrayList<TrackerHit>(100);
+ for (List<TrackerHit> list : hitCollections) hitList.addAll(list);
+ } else {
+ hitList = event.get(TrackerHit.class, _trackerHitCollectionName);
+ }
+
+ _mcTrackList = new ArrayList<MCTrack>(100);
+ List<MCParticle> mcParticleList = event.getMCParticles();
+ for (MCParticle mcPart : mcParticleList) {
+ List<TrackerHit> relatedHits = findRelatedHits(mcPart, hitList);
+ if (relatedHits.size() > 0) {
+ MCTrack mcTrack = new MCTrack(mcPart);
+ mcTrack.addHits(relatedHits);
+ mcTrack.setStatus(IN_CLASS, mcTrackFilter(mcTrack) ? YES : NO);
+ _mcTrackList.add(mcTrack);
+ }
+ }
+ hitList = null;
+ _mcTrackList.trimToSize();
+
+ // Associate reconstructed tracks with Monte Carlo Tracks
+
+ _recoTrackList = event.get(Track.class, _trackCollectionName);
+ _ratedTrackList = new ArrayList<RatedTrack>(_recoTrackList.size());
+
+ for (Track track : _recoTrackList) {
+ _ratedTrackList.add(createRatedTrack(track));
+ }
+
+ // User's analysis
+
+ analyzeEvent(event);
+
+ // Update cumulative statistics
+
+ if (_statE != null) {
+ for (int i=0; i<_statE.length; i++) _statC[i] += _statE[i];
+ }
+ _nEvents++;
+
+ // Free memory
+
+ _mcTrackList = null;
+ _recoTrackList = null;
+ _ratedTrackList = null;
+ }
+
+ /**
+ * Returns <code>true</code> for events that should be included in analysis.
+ * Default implementation always returns <code>true</code>.
+ * Subclasses should overwrite this method if some events need to be excluded.
+ */
+ protected boolean eventFilter(EventHeader event) {
+ return true;
+ }
+
+ /**
+ * Returns <code>true</code> if reconstructed track and Monte Carlo track are deemed to describe the same track.
+ * Default implementation returns <code>true</code> if all hits associated with the reconstructed track are
+ * also associated with the MC track. Subclasses should overwrite this method if different algorithm is desired.
+ */
+ protected boolean matchTracks(Track track, MCTrack mcTrack) {
+ List<TrackerHit> recoTrackHitList = track.getTrackerHits();
+ List<TrackerHit> mcTrackHitList = mcTrack.getTrackerHits();
+ return mcTrackHitList.containsAll(recoTrackHitList);
+ }
+
+ /**
+ * Returns <code>true</code> if Monte Carlo track belongs to the class that needs to be included in analysis.
+ * Default implementation always returns <code>true</code>.
+ * Subclasses should overwrite this method if more narrow class is desired (for example, requesting only
+ * tracks produced by final state charged particles is usually a good idea at the moment).
+ */
+ protected boolean mcTrackFilter(MCTrack mcTrack) {
+ return true;
+ }
+
+ /** User analysis of the event. */
+ abstract protected void analyzeEvent(EventHeader event);
+
+ /**
+ * Called by the framework when the event processing is suspended.
+ * Unless overwritten in the subclass, prints out standard cumulative statistics.
+ */
+ protected void suspend() {printCumulativeStatistics();}
+
+ /**
+ * Called when all data processing is finished.
+ * Unless overwritten in the subclass, prints out standard cumulative statistics.
+ */
+ protected void endOfData() {printCumulativeStatistics();}
+
+ /**
+ * Called before the first event is processed, or after a rewind.
+ * Unless overwritten in the subclass, resets standard cumulative statistics.
+ */
+ protected void startOfData() {
+ _nEvents = 0;
+ if (_statE != null) Arrays.fill(_statC, 0.);
+ }
+
+// -- Public helper methods : -------------------------------------------------
+
+ /**
+ * Returns a list of hits from the list supplied as an argument that were produced by the given MCParticle.
+ */
+ public static List<TrackerHit> findRelatedHits(MCParticle mcParticle, List<TrackerHit> hitList) {
+ ArrayList<TrackerHit> outputList = new ArrayList<TrackerHit>(1);
+ for (TrackerHit hit : hitList) {
+ List<MCParticle> mcList = extractMCParticles(hit);
+ if (mcList == null || mcList.isEmpty()) {
+ mcList = new ArrayList<MCParticle>(1);
+ List rawHitList = hit.getRawHits();
+ for (Object rawHit : rawHitList) {
+ List<MCParticle> list = extractMCParticles(rawHit);
+ if (list != null) mcList.addAll(list);
+ }
+ }
+ if (mcList.contains(mcParticle)) outputList.add(hit);
+ }
+ return outputList;
+ }
+
+ /**
+ * Returns a list of MCParticles associated with the given object.
+ * If the object has either <code>getMCParticle()</code> or <code>getMCParticles()</code>
+ * methods, their output is returned unless it is <code>null</code> (in which case empty list is
+ * returned). If the object has none of these methods, <code>null</code> is returned.
+ */
+ public static List<MCParticle> extractMCParticles(Object object) {
+ List<MCParticle> outputList;
+ Class objectClass = object.getClass();
+ try {
+ Method method = objectClass.getMethod("getMCParticle");
+ MCParticle mcPart = (MCParticle) method.invoke(object);
+ outputList = new ArrayList<MCParticle>(1);
+ if (mcPart != null) {
+ outputList.add(mcPart);
+ }
+ } catch (Exception ex) {
+ try {
+ Method method = objectClass.getMethod("getMCParticles");
+ outputList = (List<MCParticle>) method.invoke(object);
+ if (outputList == null) outputList = new ArrayList<MCParticle>(1);
+ } catch (Exception exx) {
+ outputList = null;
+ }
+ }
+ return outputList;
+ }
+
+// -- Protected helper methods: -----------------------------------------------
+
+ /**
+ * Creates RatedTrack from a reconstructed track and associates a Monte Carlo track with it.
+ * This method calls <code>matchTracks(Track track, MCTrack mcTrack)</code>
+ * method that can be overwritten by subclasses to reflect the
+ * desired algorithm for associating Monte Carlo tracks with reconstructed tracks.
+ */
+ protected RatedTrack createRatedTrack(Track track) {
+
+ RatedTrack ratedTrack = new RatedTrack(track);
+
+ MCTrack matchingMCTrack = null;
+ for (MCTrack mcTrack : _mcTrackList) {
+ if (matchTracks(track, mcTrack)) {
+ matchingMCTrack = mcTrack;
+ break;
+ }
+ }
+ if (matchingMCTrack != null) {
+ ratedTrack.setMCTrack(matchingMCTrack);
+ matchingMCTrack.addTrack(ratedTrack);
+ }
+
+ return ratedTrack;
+ }
+
+ protected enum IndexType {NOSUM, SUM, SUM2, NOPRINT, PRINTEVENT, PRINTSUM, PRINTALL}
+
+ /**
+ * Creates a new statistical variable and returns an index that can be used to access it.
+ * "Per event" data can be put into <code>_statE[index]</code>, these variables are reset
+ * to zero at the beginning of each event. At the end of the event, statistics is
+ * automatically accumulated into <code>_statC[index]</code>. The name of the variable
+ * is accessable through <code>_statName[index]</code>.
+ *
+ * @param name Name to be assigned to the variable.
+ * @param args Optional flags defining accumulation type and printing preferences. Possible flags and their meanings:<p>
+ * NOSUM - do not accumulate this variable.<br>
+ * SUM - <code>_statE[index]</code> added to <code>_statC[index]</code> at the end of event.<br>
+ * SUM2 - square of <code>_statE[index]</code> added to <code>_statC[index]</code> at the end of event.<br>
+ * If no flag is given, SUM is assumed.<p>
+ * NOPRINT - variable is not printed by {@link #printStatistics(EventHeader event)} and {@link #printCumulativeStatistics()}<br>
+ * PRINTEVENT - variable is printed by {@link #printStatistics(EventHeader event)} but not by {@link #printCumulativeStatistics()}<br>
+ * PRINTSUM - variable is printed by {@link #printCumulativeStatistics()} but not by {@link #printStatistics(EventHeader event)}<br>
+ * PRINTALL - variable is printed by both {@link #printCumulativeStatistics()} and {@link #printStatistics(EventHeader event)}<br>
+ * If no flag is given, PRINTALL is assumed.<p>
+ * @return Index into {@link #_statE}, {@link #_statC}, and {@link #_statName} arrays that can be used to access the variable being created.
+ */
+ protected int createIndex(String name, IndexType... args) {
+
+ int size;
+ if (_statE == null) {
+ _statE = new double[] {0.};
+ _statC = new double[] {0.};
+ _statName = new String[] {name};
+ _statAccumulateFlag = new IndexType[] {IndexType.SUM};
+ _statPrintFlag = new IndexType[] {IndexType.PRINTALL};
+ size = 0;
+ } else {
+
+ size = _statE.length;
+
+ double[] tempD = new double[size+1];
+ System.arraycopy(_statE, 0, tempD, 0, size);
+ _statE = tempD;
+ _statE[size] = 0.;
+
+ tempD = new double[size+1];
+ System.arraycopy(_statC, 0, tempD, 0, size);
+ _statC = tempD;
+ _statC[size] = 0.;
+
+ String[] tempS = new String[size+1];
+ System.arraycopy(_statName, 0, tempS, 0, size);
+ _statName = tempS;
+ _statName[size] = name;
+
+ IndexType[] tempI = new IndexType[size+1];
+ System.arraycopy(_statAccumulateFlag, 0, tempI, 0, size);
+ _statAccumulateFlag = tempI;
+ _statAccumulateFlag[size] = IndexType.SUM;
+
+ tempI = new IndexType[size+1];
+ System.arraycopy(_statPrintFlag, 0, tempI, 0, size);
+ _statPrintFlag = tempI;
+ _statPrintFlag[size] = IndexType.PRINTALL;
+ }
+
+ if (args != null) {
+ for (IndexType type : args) {
+ switch (type) {
+ case NOSUM:
+ case SUM:
+ case SUM2:
+ _statAccumulateFlag[size] = type; break;
+ case NOPRINT:
+ case PRINTEVENT:
+ case PRINTSUM:
+ case PRINTALL:
+ _statPrintFlag[size] = type; break;
+ }
+ }
+ }
+
+ return size;
+ }
+
+ /** Print statistics at the end of an event. */
+ protected void printStatistics(EventHeader event) {
+ System.out.println("");
+ System.out.println("Tracking test: " + _testName + ", statistics for event " + event.getEventNumber());
+ for (int i=0; i<_statE.length; i++) {
+ if (_statPrintFlag[i] == IndexType.PRINTEVENT || _statPrintFlag[i] == IndexType.PRINTALL) {
+ System.out.println(" " + _statName[i] + " " + _statE[i]);
+ }
+ }
+ System.out.println("");
+ }
+
+ /**
+ * Print cumulative statistics.
+ * @param messages Array of strings that should be printed before standard statistics.
+ * @param printStandard Flag to indicate whether to print standard statistics.
+ */
+ protected void printCumulativeStatistics(String[] messages, boolean printStandard) {
+ ArrayList<String> msgList = new ArrayList<String>(20);
+ msgList.add("");
+ if (messages != null) for (String msg : messages) msgList.add(" " + msg);
+ if (printStandard && _nEvents > 0) {
+ for (int i=0; i<_statE.length; i++) {
+ if (_statPrintFlag[i] == IndexType.PRINTSUM || _statPrintFlag[i] == IndexType.PRINTALL) {
+ msgList.add(" Average " + _statName[i] + " " + _statC[i]/_nEvents);
+ }
+ }
+ }
+ msgList.add("");
+ Object[] msgArray = msgList.toArray();
+ JOptionPane.showMessageDialog(null, msgArray,
+ "Tracking test: " + _testName + " cumulative statistics (" + _nEvents + " events)",
+ JOptionPane.INFORMATION_MESSAGE);
+ }
+
+ /**
+ * Print standard cumulative statistics.
+ */
+ protected void printCumulativeStatistics() {printCumulativeStatistics(null, true);}
+
+ /**
+ * Print cumulative statistics.
+ * @param messages Array of strings that should be printed before standard statistics.
+ */
+ protected void printCumulativeStatistics(String[] messages) {printCumulativeStatistics(messages, true);}
+
+// -- Private fields : --------------------------------------------------------
+
+ protected String _testName;
+
+ protected ArrayList<MCTrack> _mcTrackList;
+ protected List<Track> _recoTrackList;
+ protected ArrayList<RatedTrack> _ratedTrackList;
+
+ protected String _trackCollectionName;
+ protected String _trackerHitCollectionName;
+
+ protected int _nEvents; // Number of events analyzed
+ protected double[] _statE; // Array to hold per event statistics
+ protected double[] _statC; // Array to hold cumulative statistics
+ protected String[] _statName; // Array to hold names of statical variables
+ protected IndexType[] _statAccumulateFlag; // Array to hold flags indicating accumulation methods
+ protected IndexType[] _statPrintFlag ; // Array to hold flags indicating accumulation printing options
+
+ protected AIDA _aida;
+ protected ITree _tree;
+
+}
lcsim/src/org/lcsim/contrib/garfield/tester
diff -N TrackingTestKs.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ TrackingTestKs.java 23 Feb 2006 19:10:10 -0000 1.1
@@ -0,0 +1,69 @@
+package org.lcsim.contrib.garfield.tester;
+
+import java.util.*;
+import java.lang.reflect.*;
+import org.lcsim.event.*;
+import org.lcsim.event.MCParticle.SimulatorStatus;
+import org.lcsim.event.EventHeader.LCMetaData;
+import org.lcsim.util.Driver;
+import hep.physics.particle.Particle;
+import org.lcsim.util.aida.AIDA;
+import hep.aida.ITree;
+
+
+/**
+ * An example of track finding performance test that looks at reconstruction of
+ * charged pions produced in K short decays.
+ *
+ * @author D. Onoprienko
+ * @version $Id: TrackingTestKs.java,v 1.1 2006/02/23 19:10:10 onoprien Exp $
+ * InternalRelease 2.02.01 Feb 23, 2006
+ */
+public class TrackingTestKs extends TrackingTest {
+
+// -- Indices into the statistics (per event and cumulative) arrays : ---------
+
+ protected final int _I_eKshortsPFraction;
+
+// -- Constructors : ----------------------------------------------------------
+
+ public TrackingTestKs(String name) {
+ this();
+ set(Parameter.TEST_NAME, name);
+ }
+
+ public TrackingTestKs() {
+ _I_eKshortsPFraction = createIndex("Fraction of energy carried by Kshorts");
+ }
+
+// -- Processing event : ------------------------------------------------------
+
+ protected void analyzeEvent(EventHeader event) {
+
+ // Calculate fraction of charged track momenta carried by Kshorts
+
+ List<MCParticle> mcParticleList = event.getMCParticles();
+ double eTotal = 0.;
+ double eKs = 0.;
+ for (MCParticle mcPart : mcParticleList) {
+ if (mcPart.getGeneratorStatus() == Particle.FINAL_STATE) eTotal += mcPart.getEnergy();
+ if (mcPart.getPDGID() == 310) eKs += mcPart.getEnergy();
+ }
+ _statE[_I_eKshortsPFraction] = eKs / eTotal;
+
+ super.analyzeEvent(event);
+
+ }
+
+ protected boolean mcTrackFilter(MCTrack mcTrack) {
+
+ if (!super.mcTrackFilter(mcTrack)) return false;
+
+ MCParticle mcPart = mcTrack.getMCParticle();
+ List<MCParticle> parentList = mcPart.getParents();
+ return (((mcPart.getPDGID() == 211) || (mcPart.getPDGID() == -211))
+ && (parentList != null) && (parentList.size() == 1)
+ && (parentList.get(0).getPDGID() == 310));
+ }
+
+}
lcsim/src/org/lcsim/contrib/garfield/tester
diff -N package-info.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ package-info.java 23 Feb 2006 19:10:10 -0000 1.1
@@ -0,0 +1,5 @@
+/**
+ * Contains classes necessary for running simple performance tests on track finders.
+ */
+package org.lcsim.contrib.garfield.tester;
+
CVSspam 0.2.8
