lcsim/sandbox/RobKutschke/TKNHits
diff -N ReadTKNClusterDriverV1.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ReadTKNClusterDriverV1.java 21 Nov 2007 22:56:48 -0000 1.1
@@ -0,0 +1,370 @@
+import org.lcsim.util.aida.AIDA;
+
+import java.util.List;
+import java.util.ArrayList;
+import java.util.Set;
+import java.util.HashSet;
+import java.util.Map;
+import java.util.HashMap;
+import java.util.Iterator;
+
+import org.lcsim.util.Driver;
+
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.event.RawTrackerHit;
+import org.lcsim.event.TrackerHit;
+import org.lcsim.event.base.BaseTrackerHit;
+import org.lcsim.event.base.BaseTrackerHitMC;
+
+import org.lcsim.geometry.Detector;
+import org.lcsim.detector.IDetectorElement;
+import org.lcsim.detector.IGeometryInfo;
+import org.lcsim.detector.IRotation3D;
+import org.lcsim.geometry.compact.Subdetector;
+import org.lcsim.detector.identifier.IExpandedIdentifier;
+import org.lcsim.detector.solids.ISolid;
+import org.lcsim.detector.solids.Box;
+
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.BasicHep3Vector;
+import hep.physics.vec.BasicHep3Matrix;
+import hep.physics.vec.VecOp;
+import hep.physics.matrix.Matrix;
+import hep.physics.matrix.SymmetricMatrix;
+import hep.physics.matrix.MatrixOp;
+
+import org.lcsim.detector.tracker.silicon.SiSensor;
+
+import org.lcsim.detector.driver.SimTrackerHitIdentifierReadoutDriver;
+import org.lcsim.contrib.RobKutschke.TKNHits.TKNRawHitsDriverV1;
+import org.lcsim.contrib.RobKutschke.TKNHits.TKNClusterMakerDriverV1;
+
+/**
+ *
+ * Driver to illustrate use of the clusters created by Tim's code.
+ * The clusters are stored as SimTrackerHits.
+ *
+ *@author $Author: kutschke $
+ *@version $Id: ReadTKNClusterDriverV1.java,v 1.1 2007/11/21 22:56:48 kutschke Exp $
+ *
+ * Date $Date: 2007/11/21 22:56:48 $
+ *
+ * Notes:
+ *
+ * 1) This example is known to work with the file:
+ * http://www.lcsim.org/test/lcio/mu-_10GeV_SiTrackerBarrelTest00.slcio
+ * It has 100 events of single muons, all with an azimuth in a small cone around
+ * -90 degrees. This file uses the detector description: SiTrackerBarrelTest00 .
+ *
+ * 2) I have only tested the code with tracker barrel hits. I am told that it should work
+ * with tracker endcap hits if waferized versions of the tracker endcap are are present
+ * in the geometry file. The hit simulation code does not yet work with pixel devices.
+ * In this code, there is one container per type of hit (sim/raw/tracker) per subdetector,
+ * where subdetector = ( tracker/vertex/forward x barrel/endcap).
+ *
+ * 3) I have only checked this code with axial strips. I don't know yet how the code will
+ * deal with stereo strips or with strips that are not parallel to a side of a wafer.
+ * In those cases, does local x refer to the measurement direction or to some feature
+ * of the wafer, independent of the stip direction. If the later, the knowledge of the
+ * stereo angle is encoded in xy correlation coefficient in the local covariance matrix
+ * of the hit position.
+ *
+ * 4) In some distribitions there are outliers that come from confusion caused by the presence
+ * of two particles in the event. In most cases I have remade the plot only for events
+ * with a single muon. This cleans up the outliers.
+ *
+ */
+
+public class ReadTKNClusterDriverV1 extends Driver
+{
+ private AIDA aida = AIDA.defaultInstance();
+
+ // Names of collections on which to work. See note 2.
+ private String[] sim = {"SiTrackerBarrel_RO"};
+ private String[] raw = {"RawBarrelTrackerHits"};
+ private String[] clus = {"BarrelTrackerHits"};
+
+ // Printout limit control.
+ private int nprint = 0;
+
+ /**
+ * Constructor:
+ * Add the necessary sub-drivers.
+ *
+ */
+ public ReadTKNClusterDriverV1(){
+
+ // Attach SimTrackerHits to their detector elements.
+ super.add( new SimTrackerHitIdentifierReadoutDriver( sim ) );
+
+ // Create Raw Trackerhits and put them in the event.
+ super.add( new TKNRawHitsDriverV1( sim, raw) );
+
+ // Create Raw Trackerhits and put them in the event.
+ super.add( new TKNClusterMakerDriverV1( raw, clus) );
+ }
+
+ /**
+ * Loop over the RawTrackerHits; make printout and histograms.
+ *
+ */
+ public void process(EventHeader header)
+ {
+ // Create the list of TrackerHits and add it to the event.
+ super.process(header);
+
+ // Extract list of TrackerHits from the event.
+ List<TrackerHit> trkhits = header.get(TrackerHit.class, clus[0] );
+ aida.cloud1D("Number of clusters per event",-1).fill(trkhits.size());
+
+ // List of generated MCParticles. Used in diagnostics only.
+ List<MCParticle> genmcps = header.get(MCParticle.class,header.MC_PARTICLES);
+
+ // Per sensor hit collections, built by this code.
+ Map<SiSensor,List<TrackerHit>> sensorhits = new HashMap<SiSensor,List<TrackerHit>>();
+
+ // Loop over the TrackerHits.
+ for ( TrackerHit thit : trkhits ){
+
+ // Information about this cluster.
+ int cluster_size = thit.getRawHits().size();
+ double dedx = thit.getdEdx();
+
+ // So far these are always zero.
+ double time = thit.getTime();
+ int type = thit.getType();
+
+ // Information about the sensor on which this cluster lies.
+ SiSensor sensor = extractSiSensor( thit );
+ IGeometryInfo geom = sensor.getGeometry();
+
+ // Hit position, in global and local coordinate systems. See note 3.
+ Hep3Vector global = new BasicHep3Vector( thit.getPosition() );
+ Hep3Vector local = geom.transformGlobalToLocal(global);
+
+ // Covariance matrix of hit position in global and local coordinates.
+ Matrix gcov = new SymmetricMatrix ( 3, thit.getCovMatrix(), true);
+ Matrix lcov = covToLocal( geom, gcov );
+
+ // Errors in local x and local y. See note 3.
+ double sx = Math.sqrt(lcov.e(0,0));
+ double sy = Math.sqrt(lcov.e(1,1));
+
+ // Local XY correlation coefficient. See note 3.
+ double rhoxy = lcov.e(0,1);
+ if ( sx >0 && sy >0 ){
+ rhoxy = rhoxy/sx/sy;
+ }
+
+ // Add this hit to the per sensor hit collections.
+ List<TrackerHit> thits = sensorhits.get(sensor);
+ if ( thits == null ){
+ thits = new ArrayList<TrackerHit>();
+ sensorhits.put(sensor,thits);
+ }
+ thits.add(thit);
+
+ // Diagnostics.
+ double r = Math.sqrt(global.x()*global.x() + global.y()*global.y() );
+ aida.cloud1D("Cluster size").fill( cluster_size );
+ aida.cloud2D("Cluster size vs Hit radius (mm)",-1).fill(r,cluster_size);
+ aida.cloud1D("dEdx per TrackerHit (GeV)").fill(dedx);
+ aida.cloud1D("Time (ps)").fill(time);
+ aida.cloud1D("Type").fill(type);
+ aida.histogram1D("Sigma x (mm)",100,0., 0.010).fill(sx);
+ aida.histogram1D("Sigma y (mm)",100,0.,50. ).fill(sy);
+ aida.cloud1D("XY Correlation Coefficient").fill(rhoxy);
+ if( genmcps.size() == 1 ){
+ aida.cloud1D("Cluster size, single muon only").fill( cluster_size );
+ aida.cloud2D("Cluster size vs Hit radius (mm), single muon only",-1).fill(r,cluster_size);
+ aida.cloud1D("dEdx per TrackerHit (GeV), single muon only").fill(dedx);
+ }
+
+ ISolid solid = geom.getLogicalVolume().getSolid();
+ Box box = (Box)solid;
+ double xhalf = box.getXHalfLength();
+ double yhalf = box.getYHalfLength();
+ double zhalf = box.getZHalfLength();
+
+ // Make sure the cluster is actually on the sensor.
+ aida.histogram1D("Local x (units of xhalf)",100,-1.,1.).fill(local.x()/xhalf);
+ aida.histogram1D("Local y (units of yhalf)",100,-1.,1.).fill(local.y()/yhalf);
+ aida.histogram1D("Local z (units of zhalf)",100,-1.,1.).fill(local.z()/zhalf);
+
+ // Sigma**2 in z; should be zero within roundoff errors, so it can be negative.
+ double s2z = lcov.e(2,2);
+ aida.histogram1D("Sigma**2 z (micron^2)",100,-1.,1.).fill(s2z*1.e6);
+
+ // How many different MCParticles contributed to this hit?
+ BaseTrackerHitMC mc_hit = (BaseTrackerHitMC)thit;
+ aida.cloud1D("Number of MCParticles contributing to cluster",-1).fill( mc_hit.mcParticles().size() );
+
+ // Check the resolution.
+ List<SimTrackerHit> sims = getSimTrackerHits(thit);
+ Hep3Vector genpos = getTrueZeroCrossing( sims, sensor );
+ double dx = local.x()-genpos.x();
+ aida.cloud1D("Number of SimTrackerHits per TrackerHit").fill( sims.size() );
+ aida.cloud1D("Measured - Generated position").fill(dx);
+ if( genmcps.size() == 1 ){
+ aida.cloud1D("Measured - Generated position, single muon only").fill(dx);
+ if ( sx > 0 ){
+ aida.cloud1D("Pull in Measured position, single muon only").fill(dx/sx);
+ }
+ }
+
+ } // loop over TrackerHits
+
+
+ // Diagnostics from the per Sensor hit collections.
+ aida.cloud1D("Number of sensors with clusters").fill( sensorhits.keySet().size());
+ for ( SiSensor sensor: sensorhits.keySet() ){
+ aida.cloud1D("Clusters per sensor").fill(sensorhits.get(sensor).size());
+ if ( genmcps.size() == 1 ){
+ aida.cloud1D("Clusters per sensor, single muon only").fill(sensorhits.get(sensor).size());
+ }
+ }
+
+ } // process(EventHeader)
+
+
+ /**
+ *
+ * Naive method for finding the true coordinates of a hit.
+ *
+ * The true value is defined to be the point at which the track passes through
+ * the local z=0 plane of the sensor. The algorithm can be confused if more
+ * than one true track contributes energy to one raw trackerhit. So the result
+ * should only be trusted 100% for events that contain only a single track.
+ *
+ * @return the point, in sensor-local coordinates, at which the generated track
+ * passed through the local z=0 plane of the specified sensor.
+ *
+ * The first argument is a list of SimTrackerHits that are on the sensor
+ * and are associated with the RawTrackerHit or TrackerHit of interest.
+ *
+ * The second argument is the sensor of interest.
+ *
+ */
+ private Hep3Vector getTrueZeroCrossing( List<SimTrackerHit> list, SiSensor sensor ){
+
+ // Geometry information about this sensor.
+ IGeometryInfo geom = sensor.getGeometry();
+ IRotation3D rot = geom.getGlobalToLocal().getRotation();
+
+ // Find the SimTrackerHit that is closest to the local z=0 plane
+ // and compute the local coord of that point.
+ SimTrackerHit closest_sim = null;
+ Hep3Vector closest_local = null;
+ double mindz = Double.MAX_VALUE;
+ for ( SimTrackerHit sim : list ){
+
+ // The position of this hit in global coordinates.
+ Hep3Vector pos = new BasicHep3Vector( sim.getPoint() );
+
+ // Position of this hit in sensor-local coordinates.
+ Hep3Vector local = geom.transformGlobalToLocal(pos);
+
+ // Find hit that is closest to the z=0 plane of the sensor.
+ double dz = Math.abs(local.z());
+ if ( dz < mindz ){
+ mindz = dz;
+ closest_sim = sim;
+ closest_local = local;
+ }
+ }
+
+ // Extrapolate the track, as a straight line, to the point at which it crosses the
+ // local z=0 plane. Start at the SimTrackerHit that is closest to the z=0 plane.
+ // Typically this is step length less than 100 microns but it can be much longer for
+ // a track with a direction that is almost in the plane of the sensor.
+
+ // Unit vector along the track direction (in global coord).
+ Hep3Vector t = VecOp.unit(new BasicHep3Vector( closest_sim.getMomentum() ));
+
+ // Unit vector along track direction (in local coord).
+ Hep3Vector tlocal = rot.rotated(t);
+
+ // Step length along the track required to move the track to the zlocal=0 plane.
+ double s = -closest_local.z()/tlocal.z();
+
+ // Local coord of intersection between the track and the plane.
+ Hep3Vector intersection = VecOp.add( closest_local, VecOp.mult( s, tlocal));
+ return intersection;
+
+ } // getTrueZeroCrossing();
+
+ /**
+ *
+ * Utility routine to construct a list of SimTrackerHits that contribute to this TrackerHit.
+ * A SimTrackerHit may contribute to more than one RawTrackerHit; if this occurs, this code
+ * will return a single reference to any duplicated SimTrackerHits.
+ *
+ */
+ private List<SimTrackerHit> getSimTrackerHits( TrackerHit hit ){
+
+ // The set of SimTrackerHits that belong to this TrackerHit.
+ Set<SimTrackerHit> simhits = new HashSet<SimTrackerHit>();
+
+ // For each RawTrackerHit, add the SimTrackHits underneath it.
+ for ( Iterator ihit=hit.getRawHits().iterator(); ihit.hasNext(); ){
+ RawTrackerHit raw = (RawTrackerHit)ihit.next();
+ simhits.addAll( raw.getSimTrackerHit() );
+ }
+
+ return new ArrayList<SimTrackerHit>(simhits);
+
+ }
+
+ /**
+ *
+ * Utility routine to extract the SiSensor from the TrackerHit.
+ * Also does a consistency check.
+ *
+ */
+ private SiSensor extractSiSensor( TrackerHit hit ){
+
+ // List of individual hit strip/pixel.
+ List raws = hit.getRawHits();
+
+ // Extract sensor from first hit strip/pixel.
+ RawTrackerHit raw0 = (RawTrackerHit)raws.get(0);
+ IDetectorElement de0 = raw0.getDetectorElement();
+
+ // Consistency check.
+ for ( Iterator ihit=raws.iterator(); ihit.hasNext(); ){
+ RawTrackerHit raw = (RawTrackerHit)ihit.next();
+ IDetectorElement de = raw.getDetectorElement();
+
+ if ( de != de0 ){
+ System.out.println ("Confused ReadTKNClusterDriverV1. TrackerHit spread across sensors.");
+ }
+ }
+
+ return (SiSensor) de0;
+ }
+
+ /**
+ *
+ * Utility routine to transform a covariance matrix from the global to the
+ * local coordinate systems.
+ *
+ */
+ private Matrix covToLocal( IGeometryInfo geom, Matrix gcov ){
+
+ // Global to local rotation matrix, and its transpose (which is also its inverse).
+ BasicHep3Matrix R = (BasicHep3Matrix)geom.getGlobalToLocal().getRotation().getRotationMatrix();
+ BasicHep3Matrix RT = (BasicHep3Matrix)geom.getLocalToGlobal().getRotation().getRotationMatrix();
+
+ // The last line is equivalent to the following:
+ // No idea if there is a speed difference between these two ways of accesing the transpose.
+ //BasicHep3Matrix RT = new BasicHep3Matrix(R);
+ //RT.transpose();
+
+ // Rotate covariance matrix from global to local.
+ return MatrixOp.mult(R,MatrixOp.mult(gcov,RT));
+
+ }
+
+}