lcsim/sandbox/RobKutschke/TKNHits
diff -u -r1.5 -r1.6
--- ReadTKNRawHitsDriverV1.java 6 Nov 2007 00:30:56 -0000 1.5
+++ ReadTKNRawHitsDriverV1.java 16 Nov 2007 16:57:03 -0000 1.6
@@ -1,11 +1,10 @@
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.TreeMap;
import org.lcsim.util.Driver;
import org.lcsim.event.EventHeader;
@@ -15,10 +14,10 @@
import org.lcsim.event.base.BaseRawTrackerHit;
import org.lcsim.geometry.Detector;
-import org.lcsim.geometry.compact.Subdetector;
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;
@@ -30,53 +29,60 @@
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.detector.driver.SimTrackerHitIdentifierReadoutDriver;
import org.lcsim.contrib.RobKutschke.TKNHits.TKNRawHitsDriverV1;
-import org.lcsim.contrib.RobKutschke.TKNHits.TrackerIdentifierIndexCache;
/**
*
* Driver to illustrate use of the RawTrackerHits created by Tim's code.
* In this version, all of the hits are returned as a single unsorted vector
- * per system.
+ * 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.
*
- * This version also contains tests to compared 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: ReadTKNRawHitsDriverV1.java,v 1.5 2007/11/06 00:30:56 kutschke Exp $
+ *@version $Id: ReadTKNRawHitsDriverV1.java,v 1.6 2007/11/16 16:57:03 kutschke Exp $
*
- * Date $Date: 2007/11/06 00:30:56 $
+ * Date $Date: 2007/11/16 16:57:03 $
*
* Notes:
*
- * 0) This example is known to work with the file:
+ * 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 .
*
- * 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.
- *
- * 2) I am not aware of anything in the event model that informs the user that correct
+ * 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().
- *
- * 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.
+ * VecOp( local, uaxis) simplifies to local.x(). I believe that taking the dot product
+ * will be robust for all strip/pixel orientations.
*
- * 5) I have only tested the code with tracker barrel hits. I am told that it should work
+ * 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 code does not yet work with pixel devices.
+ * 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.
*
*/
@@ -84,10 +90,7 @@
{
private AIDA aida = AIDA.defaultInstance();
- // Cache of index values that are expensive to look up. See note 1.
- private TrackerIdentifierIndexCache index = null;
-
- // Names of collections on which to work. See note 5.
+ // Names of collections on which to work. See note 4.
private String[] sim = {"SiTrackerBarrel_RO"};
private String[] raw = {"RawBarrelTrackerHits"};
@@ -114,9 +117,10 @@
*/
public void detectorChanged( Detector detector) {
- // Initialize the index cache.
+ // Get the helper.
Subdetector tracker = findSubdetector( detector, "SiTrackerBarrel");
- index = new TrackerIdentifierIndexCache( tracker, true );
+ SiTrackerIdentifierHelper helper1 =
+ (SiTrackerIdentifierHelper)tracker.getDetectorElement().getIdentifierHelper();
}
/**
@@ -132,27 +136,31 @@
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 System Layer Module Sensor Side Strip Time ADC NSim\n");
+ 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>> sensormap = new HashMap<SiSensor,Map<Integer,RawTrackerHit>>();
+
// Loop over the raw hits.
for ( RawTrackerHit raw : rawhits ){
- // 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 = raw.getExpandedIdentifier();
-
- // Hit strip.
- int strip = id.getValue(index.getStrip());
-
// Information about the wafer on which this hit lives.
IDetectorElement de = raw.getDetectorElement();
SiSensor sensor = (SiSensor) de;
IGeometryInfo geom = de.getGeometry();
+ SiTrackerIdentifierHelper helper = (SiTrackerIdentifierHelper)de.getIdentifierHelper();
+
+ // Hit strip.
+ int strip = helper.getStripValue(raw.getIdentifier());
+
+ // Get the charge on the strip in units of electrons.
+ double charge = getChargeFromKPix( raw );
// The pattern of strips or pixels is encoded within this class. See note 2.
SiSensorElectrodes electrode = sensor.getReadoutElectrodes(ChargeCarrier.HOLE);
@@ -171,120 +179,119 @@
// 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> badsims = 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 ){
-
- // True crossing of track with local z=0 plane.
- Hep3Vector local_gen = getTrueZeroCrossing( sims, sensor);
+ // Links back to the SimTrackerHits that contributed to this hit.
+ List<SimTrackerHit> sims = raw.getSimTrackerHit();
- // The measured quantity associated with this point.
- double meas_gen = VecOp.dot( local_gen, uaxis);
+ // True crossing of track with local z=0 plane.
+ Hep3Vector local_gen = getTrueZeroCrossing( sims, sensor );
- // Difference between measured and generated quantities.
- double du = meas_reco-meas_gen;
+ // The true measured quantity associated with this point.
+ double meas_gen = VecOp.dot( local_gen, uaxis);
- // 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());
+ // 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.
+ // Some printout that exercises the identifier helper.
if ( nprint < 6 ){
- System.out.printf ( "%5d %6d %5d %6d %6d %4d %5d %4d %5d %4d\n",
+
+ System.out.printf ( "%5d %6d %6d %5d %6d %6d %4d %5d %4d %10.1f %4d \n",
header.getEventNumber(),
- id.getValue(index.getSystem()),
- id.getValue(index.getLayer()),
- id.getValue(index.getModule()),
- id.getValue(index.getSensor()),
- id.getValue(index.getSide()),
- id.getValue(index.getStrip()),
+ 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(),
- adc,
- badsims.size()
+ charge,
+ sims.size()
);
-
- // Uncomment this block to inspect the corrupted list of SimTrackerHits.
- /*
- for ( SimTrackerHit badsim : badsims ){
- double[] pos = badsim.getPoint();
- System.out.printf ( "SimTrackerHit: %9.3f %9.3f %9.3f | %20.10f | ",
- pos[0], pos[1], pos[2], badsim.getTime() );
- System.out.println (badsim);
- }
- */
-
-
}
-
+
// Some more histograms.
- aida.cloud1D("Raw ADC spectrum").fill(adc);
+ 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());
- }
- } // process(EventHeader)
+ // The time field is not yet filled in; it is always zero.
+ aida.cloud1D("Time").fill(raw.getTime());
- /**
- * 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>();
+ // See note 5.
+ if ( genmcps.size() == 0 ){
+ aida.histogram1D("Muon Only/Number of parent SimTrackerHits",30,0.,30.).fill(sims.size());
+ }
- // List of all SimTrackerHits in the event.
- List<SimTrackerHit> eventHits = header.get(SimTrackerHit.class, sim[0]);
+ // Add this hit to the per sensor map. See note 6.
+ Map<Integer,RawTrackerHit> rawmap = sensormap.get(sensor);
+ if ( rawmap == null ){
+ rawmap = new TreeMap<Integer,RawTrackerHit>();
+ sensormap.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()
+ );
+ }
- // 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();
+ } // loop over rawhits
- for ( SimTrackerHit sim : eventHits ){
+ // Inspect the organized collection.
+ for ( SiSensor sensor: sensormap.keySet() ){
- // The position of this hit in global coordinates.
- Hep3Vector pos = new BasicHep3Vector( sim.getPoint() );
+ // Histogram of the number of hits per sensor.
+ aida.cloud1D("Hits per sensor").fill(sensormap.get(sensor).size());
- // Position of this hit in sensor-local coordinates.
- Hep3Vector local = geom.transformGlobalToLocal(pos);
+ if ( nprint < 6 ){
+
+ // Hits on this sensor.
+ Map<Integer,RawTrackerHit> rawmap = sensormap.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());
- // 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);
+ System.out.println ("Sensor: " + module + " " + i );
+ }
}
+
}
- return list;
- }
+ } // process(EventHeader)
+
/**
*
- * Find true coordinates of a hit.
+ * 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 track of interest. The list must be non-null.
+ * and are associated with the RawTrackerHit of interest.
*
* The second argument is the sensor of interest.
*
@@ -299,7 +306,7 @@
// and compute the local coord of that point.
SimTrackerHit closest_sim = null;
Hep3Vector closest_local = null;
- double mindz = 1.e10;
+ double mindz = Double.MAX_VALUE;
for ( SimTrackerHit sim : list ){
// The position of this hit in global coordinates.
@@ -308,6 +315,7 @@
// 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;
@@ -317,32 +325,31 @@
}
- // 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.
+ // 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 tp = rot.rotated(t);
+ 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()/tp.z();
+ 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, tp));
-
+ Hep3Vector intersection = VecOp.add( closest_local, VecOp.mult( s, tlocal));
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();
@@ -356,4 +363,28 @@
return tracker;
}
+ /**
+ *
+ * Use the Kpix code to compute the charge on a strip, specified as
+ * a RawTrackerHit. Tim says that this complexity will eventualy be wrapped
+ * up so that users do not need to see it.
+ *
+ * For Kpix chips, the short[] return of getADCValues() is not actually an array
+ * of adc values; it contains chip state information that includes the adc value
+ * plus additional information.
+ *
+ * @return Charge on the strip in units of electrons.
+ *
+ */
+ private double getChargeFromKPix( RawTrackerHit raw){
+
+ short adc = raw.getADCValues()[1];
+ Kpix.ControlRegisters control_registers =
+ Kpix.ControlRegisters.decoded(raw.getADCValues()[0]);
+ Kpix.KpixChannel.ReadoutRegisters readout_registers =
+ Kpix.KpixChannel.ReadoutRegisters.decoded(raw.getADCValues()[1]);
+ double charge = Kpix.computeCharge(readout_registers, control_registers);
+ return charge;
+ }
+
}
