Print

Print


Commit in lcsim/sandbox/RobKutschke/TKNHits on MAIN
ReadTKNRawHitsDriverV2.java+498added 1.1
First version.

lcsim/sandbox/RobKutschke/TKNHits
ReadTKNRawHitsDriverV2.java added at 1.1
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());
+	}
+    }
+
+} 
CVSspam 0.2.8