lcsim/sandbox/RobKutschke/TKNHits
diff -u -r1.1 -r1.2
--- ReadTKNRawHitsDriverV1.java 1 Nov 2007 22:55:19 -0000 1.1
+++ ReadTKNRawHitsDriverV1.java 5 Nov 2007 21:33:27 -0000 1.2
@@ -1,6 +1,11 @@
-
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 org.lcsim.util.Driver;
import org.lcsim.event.EventHeader;
@@ -12,10 +17,19 @@
import org.lcsim.geometry.Detector;
import org.lcsim.detector.IDetectorElement;
import org.lcsim.detector.IGeometryInfo;
+import org.lcsim.detector.IRotation3D;
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.VecOp;
+
+import org.lcsim.detector.tracker.silicon.SiSensor;
+import org.lcsim.detector.tracker.silicon.SiSensorElectrodes;
+import org.lcsim.detector.tracker.silicon.ChargeCarrier;
+
import org.lcsim.detector.driver.SimTrackerHitIdentifierReadoutDriver;
import org.lcsim.contrib.RobKutschke.TKNHits.TKNRawHitsDriverV1;
@@ -27,24 +41,55 @@
* In this version, all of the hits are returned as a single unsorted vector
* per system.
*
+ * This version also contains tests to compared measured and generated positions.
+ *
*@author $Author: kutschke $
- *@version $Id: ReadTKNRawHitsDriverV1.java,v 1.1 2007/11/01 22:55:19 kutschke Exp $
+ *@version $Id: ReadTKNRawHitsDriverV1.java,v 1.2 2007/11/05 21:33:27 kutschke Exp $
+ *
+ * Date $Date: 2007/11/05 21:33:27 $
+ *
+ * Notes:
+ *
+ * 0) 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 .
*
- * Date $Date: 2007/11/01 22:55:19 $
+ * 1) There are a few magic numbers that I need to look up in a dictionary and I only
+ * want to do this when the dictionary changes. Ideally all of the code should be
+ * within detectorChanged() but I don't yet know how to do that. So the cache is
+ * invalidated when the detector changes and recomputed on the first event following
+ * a detector change.
+ *
+ * 2) I am not aware of anything in the event model that informs the user that correct
+ * argument to the getElectrodes() method is HOLES, not ELECTRONS. This is just something
+ * that you need to know a priori.
+ *
+ * 3) In some cases one can simplify the dot product. For example, for axial strips,
+ * VecOp( local, uaxis) simplifies to local.x().
+ *
+ * 4) The RawTrackerHit.getSimTrackerHit() method is currently broken. I have
+ * created a hack work around that only works for events with exactly on MCParticle
+ * ( ie does not work for events with decays in flight ). The Hack returns null
+ * for events with more than 1 MCParticle. If there is more than 1 MCParticle per event,
+ * and if two MCParticles hit the same sensor, then I have no way of matching
+ * SimTrackerHits with RawTrackerHits, so I skip the comparison of reconstructed and
+ * generated quantities for these events.
+ *
+ * 5) 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 code does not yet work with pixel devices.
*
*/
-
public class ReadTKNRawHitsDriverV1 extends Driver
{
private AIDA aida = AIDA.defaultInstance();
- // Cache of index values that are expensive to look up.
+ // Cache of index values that are expensive to look up. See note 1.
private TrackerIdentifierIndexCache index = null;
- // Names of collections on which to work.
- // As of 11/01/07 the example only works with the tracker barrel but, in the future
- // it will work with all Tracker, Vertex Detector and forward tracker hits.
+ // Names of collections on which to work. See note 5.
private String[] sim = {"SiTrackerBarrel_RO"};
private String[] raw = {"RawBarrelTrackerHits"};
@@ -53,12 +98,12 @@
/**
* Constructor:
- * Add the necessary drivers using super.
+ * Add the necessary sub-drivers using super.
*
*/
public ReadTKNRawHitsDriverV1(){
- // Attach SimTrackerHits to the detector.
+ // Attach SimTrackerHits to their detector elements.
super.add( new SimTrackerHitIdentifierReadoutDriver( sim ) );
// Create Raw Trackerhits and put them in the event.
@@ -71,24 +116,26 @@
*/
public void detectorChanged( Detector detector) {
- // Invalidate the cached indices when the detector changes.
+ // Invalidate the cached indices. See note 1.
index = null;
+ System.out.println ("Detector: " + detector.getName() );
+
}
/**
- * This does the main work of the class.
+ * Loop over the RawTrackerHits; make printout and histograms.
*
*/
public void process(EventHeader header)
{
- // Create the list of RawTrackerHits in the event.
+ // Create the list of RawTrackerHits and add it to the event.
super.process(header);
-
+
// Get the list of RawTrackerHits from the event.
List<RawTrackerHit> rawhits = header.get(RawTrackerHit.class, raw[0] );
aida.cloud1D("Number of RawTrackerHits per event").fill(rawhits.size());
- // Make sure cached indices are valid.
+ // Make sure cached indices are valid. See note 1.
if ( index == null ){
if ( rawhits.size() > 0 ){
index = new TrackerIdentifierIndexCache( rawhits.get(0), true );
@@ -99,39 +146,67 @@
// Header for the tabluar printout.
if ( ++nprint == 1 ){
- System.out.printf ( "Event System Layer Module Sensor Side Strip Time ADC NSim");
- System.out.printf ( " ( Local Origin )\n");
+ System.out.printf ( "Event System Layer Module Sensor Side Strip Time ADC NSim\n");
}
// Loop over the raw hits.
for ( RawTrackerHit raw : rawhits ){
- // Need concrete type to access some methods.
- // At present, this will not survive persistency.
- BaseRawTrackerHit base = (BaseRawTrackerHit) raw;
-
// For hits of this type, the adc array will always be of length 1, exactly.
short adc = raw.getADCValues()[0];
// Position in detector heirarchy.
- IExpandedIdentifier id = base.getExpandedIdentifier();
+ IExpandedIdentifier id = raw.getExpandedIdentifier();
+
+ // Hit strip.
+ int strip = id.getValue(index.getStrip());
// Information about the wafer on which this hit lives.
- IDetectorElement de = base.getDetectorElement();
+ IDetectorElement de = raw.getDetectorElement();
+ SiSensor sensor = (SiSensor) de;
IGeometryInfo geom = de.getGeometry();
- // Coordinates of wafer center in global system.
- Hep3Vector center = geom.getPosition();
+ // The pattern of strips or pixels is encoded within this class. See note 2.
+ SiSensorElectrodes electrode = sensor.getReadoutElectrodes(ChargeCarrier.HOLE);
+
+ // Convert the strip number to positions in local and global coord.
+ Hep3Vector local = electrode.getCellPosition(strip);
+ Hep3Vector global = geom.transformLocalToGlobal(local);
+
+ // Unit vector in measurement direction, in local coord.
+ Hep3Vector uaxis = electrode.getMeasuredCoordinate(0);
+
+ // For pixels there would be a second measurement direction, in local coord.
+ Hep3Vector vaxis = ( electrode.getNAxes() == 2 ) ?
+ electrode.getMeasuredCoordinate(1): null;
+
+ // Reconstructed measurement. See note 3. Assume strips not pixels.
+ double meas_reco = VecOp.dot( local, uaxis);
+
+ // Links back to the SimTrackerHits that contributed to this hit. See note 4.
+ //List<SimTrackerHit> sims = raw.getSimTrackerHit(); // Broken.
+ List<SimTrackerHit> sims = hackGetSimTrackerHits(header,sensor);
+
+ // Skip events for which the hack on the previous line does not work. See note 4.
+ if ( sims != null ){
- // Links back to the SimTrackerHits that contributed to this hit.
- // As of 11/01/07 this returns a corrupt list.
- List<SimTrackerHit> sims = raw.getSimTrackerHit();
+ // True crossing of track with local z=0 plane.
+ Hep3Vector local_gen = getTrueZeroCrossing( sims, sensor);
- // Coming soon: convert the strip number to a local position.
+ // The measured quantity associated with this point.
+ double meas_gen = VecOp.dot( local_gen, uaxis);
+
+ // Difference between measured and generated quantities.
+ double du = meas_reco-meas_gen;
+
+ // Some plots.
+ aida.histogram1D("Measured-Generated coordinate (mm)",100,-0.1,0.1).fill(du);
+ aida.cloud2D("Measured y vs x positions (mm)",-1).fill(global.x(),global.y());
+ }
// Some printout.
if ( nprint < 6 ){
- System.out.printf ( "%5d %6d %5d %6d %6d %4d %5d %4d %5d %4d (%9.3f %9.3f %9.3f)\n",
+ System.out.printf ( "%5d %6d %5d %6d %6d %4d %5d %4d %5d %4d\n",
header.getEventNumber(),
id.getValue(index.getSystem()),
id.getValue(index.getLayer()),
@@ -141,10 +216,7 @@
id.getValue(index.getStrip()),
raw.getTime(),
adc,
- sims.size(),
- center.x(),
- center.y(),
- center.z()
+ sims.size()
);
// Uncomment this block to inspect the corrupted list of SimTrackerHits.
@@ -159,9 +231,122 @@
}
- // And some histograms.
+ // Some more histograms.
aida.cloud1D("Raw ADC spectrum").fill(adc);
+ aida.cloud1D("Strip number").fill(strip);
+ aida.cloud1D("Reconstructed local position (mm)").fill(meas_reco);
}
} // process(EventHeader)
+
+ /**
+ * Hack to replace broken functionality of RawTrackerHit.getSimTrackerHit().
+ *
+ * @return a list of all SimTrackerHits that can make a hit on the specified sensor.
+ * If there is more than 1 MCParticle in the event, return a null list.
+ *
+ */
+
+ private List<SimTrackerHit> hackGetSimTrackerHits( EventHeader header, SiSensor sensor ){
+
+ // This hack only works reliably if there is exactly one particle in the event.
+ List<MCParticle> particles = header.get(MCParticle.class,header.MC_PARTICLES);
+ if ( particles.size() > 1 ) return null;
+
+ // Intialize the output.
+ List<SimTrackerHit> list = new ArrayList<SimTrackerHit>();
+
+ // List of all SimTrackerHits in the event.
+ List<SimTrackerHit> eventHits = header.get(SimTrackerHit.class, sim[0]);
+
+ // Geometry information about this sensor.
+ IGeometryInfo geom = sensor.getGeometry();
+ ISolid solid = geom.getLogicalVolume().getSolid();
+ Box box = (Box)solid;
+ double xhalf = box.getXHalfLength();
+ double yhalf = box.getYHalfLength();
+ double zhalf = box.getZHalfLength();
+
+ for ( SimTrackerHit sim : eventHits ){
+
+ // 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);
+
+ // Check if hit is on this sensor.
+ double dx = Math.abs(local.x());
+ double dy = Math.abs(local.y());
+ double dz = Math.abs(local.z());
+ if ( dx<=xhalf & dy<yhalf & dz<zhalf ){
+ list.add(sim);
+ }
+ }
+
+ return list;
+ }
+
+ /**
+ *
+ * Find true coordinates of a hit.
+ *
+ * @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 track of interest. The list must be non-null.
+ *
+ * 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 = 1.e10;
+ 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);
+
+ double dz = Math.abs(local.z());
+ if ( dz < mindz ){
+ mindz = dz;
+ closest_sim = sim;
+ closest_local = local;
+ }
+
+ }
+
+ // Extrapolate the track, as a straight line, from the closest point to the
+ // point at which it crosses the local z=0 plane. Typically this is step length
+ // less than 100 microns.
+
+ // 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 tp = rot.rotated(t);
+
+ // Step length along the track required to move the track to the zlocal=0 plane.
+ double s = -closest_local.z()/tp.z();
+
+ // Local coord of intersection between the track and the plane.
+ Hep3Vector intersection = VecOp.add( closest_local, VecOp.mult( s, tp));
+
+ return intersection;
+
+ }
+
+
}
