lcsim/sandbox/RobKutschke/TKNHits
diff -N ReadTKNRawHitsDriverV2.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ ReadTKNRawHitsDriverV2.java 20 Nov 2007 20:21:55 -0000 1.1
@@ -0,0 +1,498 @@
+import org.lcsim.util.aida.AIDA;
+
+import java.util.List;
+import java.util.Set;
+import java.util.Map;
+import java.util.HashMap;
+import java.util.TreeMap;
+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.base.BaseRawTrackerHit;
+
+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.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.tracker.silicon.SiTrackerIdentifierHelper;
+
+import org.lcsim.detector.DetectorIdentifierHelper;
+
+import org.lcsim.contrib.SiStripSim.Kpix;
+import org.lcsim.contrib.SiStripSim.ReadoutChip;
+
+
+import org.lcsim.detector.driver.SimTrackerHitIdentifierReadoutDriver;
+import org.lcsim.contrib.RobKutschke.TKNHits.TKNRawHitsDriverV1;
+
+/**
+ *
+ * Driver to illustrate use of the RawTrackerHits created by Tim's code.
+ *
+ *
+ * This version does the same work as V1 but with some additional diagnostics.
+ * The model to be used by general usesrs is V1.
+ *
+ * The diagnostics show that Tim's hits pass a sanity check:
+ * - they are all on the correct sensor;
+ * - the outliers in the measured-generated plots are due to confusion
+ * caused by 2 tracks contributing to the same RawTrackerHit.
+ * - A given SimTrackerHit may contribute to multiple RawTrackerHits
+ * - A given RawTrackerHit may have contributions from multiple SimTrackerHits.
+ *
+ * In this version, all of the hits are returned as a single unsorted vector
+ * per system, where system means one of Tracker Barrel, Tracker Endcap,
+ * Vertex Barrel, Vertex Endcap, Forward Tracker.
+ *
+ * This version contains tests to compare measured and generated positions.
+ *
+ * It also creates a map which organizes RawTrackerHits by sensor; this is just
+ * one of many possible organizations.
+ *
+ *@author $Author: kutschke $
+ *@version $Id: ReadTKNRawHitsDriverV2.java,v 1.1 2007/11/20 20:21:55 kutschke Exp $
+ *
+ * Date $Date: 2007/11/20 20:21:55 $
+ *
+ * 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 am not aware of anything in the event model that informs the user that the 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(). I believe that taking the dot product
+ * will be robust for all strip/pixel orientations.
+ *
+ * 4) 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.
+ *
+ * 5) In some events in the standard input file the muons decay or interact in the detector,
+ * which means that the code to compute the mc true value of the measurement may sometimes
+ * give the wrong answer. This gives rise to outliers in the (measured-generated) plots.
+ * To verify that the code is working correctly, I have also made some plots only for events
+ * that contain a single muon and nothing else; in these plots the tails and outliers are gone.
+ *
+ * 6) The collection of hits on each sensor are organized as a TreeMap. Why?
+ * - Choose a Map so that one can address hits by strip number
+ * - Choose a TreeMap so that the keyset is sorted in order of ascending strip number.
+ *
+ * 7) There are two ways to get the helper:
+ * - From the subdetector
+ * - From the RawTrackerHit
+ * I do not know yet if the same helper will work for all of Tracker/Vertex/Forward
+ * and endcap/barrel or if a separate helper is needed for each of these. If all the same then
+ * we should get it in the detectorChanged() method. If different we should get it from the
+ * RawTrackerHit. I coded the second option for now.
+ */
+
+public class ReadTKNRawHitsDriverV2 extends Driver
+{
+ private AIDA aida = AIDA.defaultInstance();
+
+ // Names of collections on which to work. See note 4.
+ private String[] sim = {"SiTrackerBarrel_RO"};
+ private String[] raw = {"RawBarrelTrackerHits"};
+
+ // Printout limit control.
+ private int nprint = 0;
+
+ // Object that can interpret raw data to return charge, noise etc.
+ // Someday this will be (sub)detector dependent.
+ ReadoutChip readout = new Kpix();
+
+ /**
+ * Constructor:
+ * Add the necessary sub-drivers using super.
+ *
+ */
+ public ReadTKNRawHitsDriverV2(){
+
+ // 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) );
+ }
+
+ /**
+ * Update any geometry related information when the detector changes.
+ *
+ */
+ public void detectorChanged( Detector detector) {
+
+ // An alternate way to get the helper. See note 7.
+ Subdetector tracker = findSubdetector( detector, "SiTrackerBarrel");
+ SiTrackerIdentifierHelper helper1 =
+ (SiTrackerIdentifierHelper)tracker.getDetectorElement().getIdentifierHelper();
+ }
+
+ /**
+ * Loop over the RawTrackerHits; make printout and histograms.
+ *
+ */
+ public void process(EventHeader header)
+ {
+ // 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());
+
+ // List of generated MCParticles. Used in diagnostics only.
+ List<MCParticle> genmcps = header.get(MCParticle.class,header.MC_PARTICLES);
+
+ // Header for the tabluar printout.
+ if ( ++nprint == 1 ){
+ System.out.printf ( "Event Barrel System Layer Module Sensor Side Strip Time Charge NSim\n");
+ }
+
+ // Container to organize RawTrackerHits on each sensor.
+ Map<SiSensor,Map<Integer,RawTrackerHit>> sensorhits = new HashMap<SiSensor,Map<Integer,RawTrackerHit>>();
+
+ // For diagnostics: Sum of SimTrackerHits associated with each RawTrackerhits.
+ int sumsim = 0;
+
+ // For diagnostics: number of RawTrackerHits associated with each SimTrackerHit.
+ Map<SimTrackerHit,Integer> simToRaw = new HashMap<SimTrackerHit,Integer>();
+ for (SimTrackerHit s : header.get(SimTrackerHit.class, sim[0]) ){
+ simToRaw.put( s, 0);
+ }
+
+ // Loop over the raw hits.
+ for ( RawTrackerHit raw : rawhits ){
+
+ // Information about the wafer on which this hit lives.
+ IDetectorElement de = raw.getDetectorElement();
+ SiSensor sensor = (SiSensor) de;
+ IGeometryInfo geom = de.getGeometry();
+
+ // See note 7.
+ SiTrackerIdentifierHelper helper = (SiTrackerIdentifierHelper)de.getIdentifierHelper();
+
+ // The pattern of strips or pixels is encoded within this class. See note 2.
+ SiSensorElectrodes electrode = sensor.getReadoutElectrodes(ChargeCarrier.HOLE);
+
+ // Hit strip.
+ int strip = helper.getStripValue(raw.getIdentifier());
+
+ // Get the charge on the strip in units of electrons.
+ double charge = readout.decodeCharge(raw);
+
+ // The expected noise level on this strip, in units of electrons.
+ double noise = readout.getChannel(strip).computeNoise(electrode.getCapacitance(strip));
+
+ // 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.
+ List<SimTrackerHit> sims = raw.getSimTrackerHit();
+ sumsim += sims.size();
+
+ // True crossing of track with local z=0 plane.
+ Hep3Vector local_gen = getTrueZeroCrossing( sims, sensor, genmcps );
+
+ // The true 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;
+
+ // A diagnostic plot.
+ aida.histogram1D("Measured-Generated coordinate (mm)",100,-0.15,0.15).fill(du);
+
+ // Some additional diagnostic plots for events with only a single generated track. See note 5.
+ if ( genmcps.size() ==1 ){
+ aida.histogram1D("Muon Only/Measured-Generated coordinate (mm)",100,-0.15,0.15).fill(du);
+ aida.cloud2D("Muon Only/Measured y vs x positions (mm)",-1).fill(global.x(),global.y());
+ aida.cloud2D("Muon Only/Global y vs (Measured-Generated coordinate) (mm)",-1).fill(du,global.y());
+ }
+
+ // Some printout that exercises the identifier helper.
+ if ( nprint < 6 ){
+
+ System.out.printf ( "%5d %6d %6d %5d %6d %6d %4d %5d %4d %10.1f %4d \n",
+ header.getEventNumber(),
+ helper.getSystem(raw.getIdentifier()),
+ helper.getBarrel(raw.getIdentifier()),
+ helper.getLayer(raw.getIdentifier()),
+ helper.getModuleValue(raw.getIdentifier()),
+ helper.getSensorValue(raw.getIdentifier()),
+ helper.getSideValue(raw.getIdentifier()),
+ helper.getStripValue(raw.getIdentifier()),
+ raw.getTime(),
+ charge,
+ sims.size()
+ );
+ }
+
+ // Some more histograms.
+ aida.cloud1D("Raw charge (electrons)").fill(charge);
+ aida.histogram1D("Raw charge, detail (electrons)",50,0.,40000.).fill(charge);
+ aida.cloud1D("Strip number").fill(strip);
+ aida.cloud1D("Reconstructed local position (mm)").fill(meas_reco);
+ aida.histogram1D("Number of parent SimTrackerHits",30,0.,30.).fill(sims.size());
+ aida.cloud1D("Expected Noise",-1).fill(noise);
+ if ( noise > 0. ){
+ aida.cloud1D("Signal to Noise Ratio",-1).fill(charge/noise);
+ }
+
+ // The time field is not yet filled in; it is always zero.
+ aida.cloud1D("Time").fill(raw.getTime());
+
+ // See note 5.
+ if ( genmcps.size() == 0 ){
+ aida.histogram1D("Muon Only/Number of parent SimTrackerHits",30,0.,30.).fill(sims.size());
+ }
+
+ // Add this hit to the per sensor map. See note 6.
+ Map<Integer,RawTrackerHit> rawmap = sensorhits.get(sensor);
+ if ( rawmap == null ){
+ rawmap = new TreeMap<Integer,RawTrackerHit>();
+ sensorhits.put(sensor,rawmap);
+ }
+ if ( rawmap.get(strip) == null ){
+ rawmap.put(strip,raw);
+ }else{
+ System.out.println ("Duplicate strip: " + strip
+ + " on module: " + helper.getModuleValue(raw.getIdentifier())
+ + " in event: " + header.getEventNumber()
+ );
+ }
+
+ // Identify RawTracker hits that have SimTrackerHits from more than 1 MCParticle.
+ FindDifferentMCPs( sims, genmcps );
+
+ // Increment counters for each SimTrackerHit used by this RawTrackerHit.
+ for ( SimTrackerHit s: sims ){
+ Integer n = simToRaw.get(s);
+ simToRaw.put(s,++n);
+ }
+
+ } // loop over rawhits
+
+ // Inspect the organized collection.
+ for ( SiSensor sensor: sensorhits.keySet() ){
+
+ // Histogram of the number of hits per sensor.
+ aida.cloud1D("Hits per sensor").fill(sensorhits.get(sensor).size());
+
+ if ( nprint < 6 ){
+
+ // Hits on this sensor.
+ Map<Integer,RawTrackerHit> rawmap = sensorhits.get(sensor);
+
+ // This traverses strips in ascending order or strip number.
+ for ( Integer i : rawmap.keySet() ){
+ RawTrackerHit raw = rawmap.get(i);
+ IDetectorElement de = raw.getDetectorElement();
+ SiTrackerIdentifierHelper helper = (SiTrackerIdentifierHelper)de.getIdentifierHelper();
+ int module = helper.getModuleValue(raw.getIdentifier());
+
+ System.out.println ("Sensor: " + module + " " + i );
+ }
+ }
+
+ }
+
+ // Plot the info in simToRaw.
+ for ( SimTrackerHit s : simToRaw.keySet() ){
+ Integer n = simToRaw.get(s);
+ aida.cloud1D("Number of RawTrackerHits per SimTrackerHit").fill(n.intValue());
+ }
+
+ } // 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 of interest.
+ *
+ * The second argument is the sensor of interest.
+ *
+ */
+ private Hep3Vector getTrueZeroCrossing( List<SimTrackerHit> list,
+ SiSensor sensor,
+ List<MCParticle> genmcps
+ ){
+
+ // Geometry information about this sensor.
+ IGeometryInfo geom = sensor.getGeometry();
+ IRotation3D rot = geom.getGlobalToLocal().getRotation();
+ ISolid solid = geom.getLogicalVolume().getSolid();
+ Box box = (Box)solid;
+ double xhalf = box.getXHalfLength();
+ double yhalf = box.getYHalfLength();
+ double zhalf = box.getZHalfLength();
+
+ // 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);
+
+ // Check that all SimTrackerHits are indeed inside the sensor.
+ aida.histogram1D( "True/Local X (units of xhalf)",100,-1.,1.).fill( local.x()/xhalf );
+ aida.histogram1D( "True/Local Y (units of yhalf)",100,-1.,1.).fill( local.y()/yhalf );
+ aida.histogram1D( "True/Local Z (units of zhalf)",100,-1.,1.).fill( local.z()/zhalf );
+
+ // 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));
+
+ // Additional diagnostics.
+ aida.histogram1D("True/z of closest SimTrackerHit (units of zhalf)",100,-1,1.).fill(closest_local.z()/zhalf);
+ aida.histogram1D("True/Step length to z=0 plane (mm)",100,-1,1.).fill(s);
+ if ( genmcps.size() == 1 ){
+ aida.histogram1D("True/Step length to z=0 plane (mm),muon only",100,-1,1.).fill(s);
+ }
+ aida.histogram1D("True/Local slope, all",100,0.1,1.0).fill(tlocal.z());
+ aida.cloud2D("True/Number of simtracker hits vs tz").fill( tlocal.z(), list.size());
+ aida.cloud2D("True/s vs tz").fill( tlocal.z(), s);
+ aida.cloud2D("True/s vs size").fill( list.size(), s);
+ if ( Math.abs(s) > 0.125 ){
+ aida.cloud1D("True/Local slope, large s").fill(tlocal.z());
+ aida.histogram1D("True/Number of parent SimTrackerHits, long s",20,0.,20.).fill(list.size());
+ if ( genmcps.size() == 1 ){
+ aida.cloud1D("True/Local slope, large s, mu only").fill(tlocal.z());
+ aida.histogram1D("True/Number of parent SimTrackerHits, long s,mu only",20,0.,20.).fill(list.size());
+ System.out.println("Long s, mu only: " + list.size() + " " + tlocal.z() + " " + s);
+
+ }
+ } else{
+ aida.histogram1D("True/Number of parent SimTrackerHits, short s",20,0.,20.).fill(list.size());
+ }
+
+ return intersection;
+
+ } // getTrueZeroCrossing();
+
+ /**
+ * Find a named subdector object within the detector.
+ *
+ * @return Subdetector object that is specified by the second argument.
+ *
+ */
+ private Subdetector findSubdetector( Detector detector, String sub ){
+ Map<String, Subdetector> subDetMap = detector.getSubdetectors();
+ Subdetector tracker = subDetMap.get(sub);
+ if ( tracker == null ){
+ System.out.println ("Cannot find subdetector: " + sub
+ + " in detector " + detector.getName() );
+ System.out.println ("Aborting now.");
+ System.exit(-1);
+ }
+ return tracker;
+ }
+
+
+ /**
+ *
+ * Identify RawTrackerHits that come from more than one MCParticle.
+ * Make plots and print some diagnostic info.
+ *
+ */
+ private void FindDifferentMCPs ( List<SimTrackerHit> list, List<MCParticle> genmcps ){
+
+ // Count number of times each MCParticle contributes to this RawTrackerHit.
+ Map<MCParticle,Integer> mcpToRaw = new HashMap<MCParticle,Integer>();
+ for ( SimTrackerHit sim : list ){
+ MCParticle mcp = sim.getMCParticle();
+ Integer count = mcpToRaw.get(mcp);
+ if ( count == null ){
+ count = new Integer(0);
+ mcpToRaw.put(mcp,count);
+ }
+ mcpToRaw.put(mcp,++count);
+ }
+ // Plot number of MCParticles that contribute.
+ aida.histogram1D("True/Number of contributing MCParticles",10,0,10).fill(mcpToRaw.size());
+
+ // If more than 1 MCParticle contributes, print info.
+ if ( mcpToRaw.size() > 1 ) {
+ System.out.printf ("Contributions from different MCParticles: " );
+ for ( MCParticle mcp : genmcps ){
+ Integer count = mcpToRaw.get(mcp);
+ if ( count != null ){
+ System.out.printf ( "( %3d : %3d )", genmcps.indexOf(mcp), count.intValue() );
+ }
+ }
+ System.out.printf (" | Size of input list: %3d\n", list.size());
+ }
+ }
+
+}