lcsim/src/org/lcsim/contrib/uiowa
diff -N MIPGeometryHandler.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ MIPGeometryHandler.java 15 Aug 2008 17:32:33 -0000 1.1
@@ -0,0 +1,119 @@
+package org.lcsim.contrib.uiowa;
+
+import java.util.*;
+import hep.physics.vec.*;
+import org.lcsim.event.Track;
+
+/**
+ * A simple abstract class designed to make it easy
+ * to get the showering point and track tangent at
+ * that point. Implementation is delegated to
+ * subclasses.
+ *
+ * @version $Id: MIPGeometryHandler.java,v 1.1 2008/08/15 17:32:33 mcharles Exp $
+ */
+
+public abstract class MIPGeometryHandler {
+
+ protected Map<Track, Hep3Vector> m_cachePoint;
+ protected Map<Track, Hep3Vector> m_cacheTangent;
+
+ /** Simple constructor. */
+ protected MIPGeometryHandler() {
+ m_cachePoint = new HashMap<Track, Hep3Vector>();
+ m_cacheTangent = new HashMap<Track, Hep3Vector>();
+ }
+
+ /**
+ * Quote the showering point of the supplied track.
+ * Returns null if not known or if track does not
+ * reach calorimeter.
+ */
+ public Hep3Vector getShowerPoint(Track tr) {
+ try {
+ BasicHep3Vector outputPoint = new BasicHep3Vector();
+ BasicHep3Vector outputTangentUnit = new BasicHep3Vector();
+ findPointAndTangent(tr, outputPoint, outputTangentUnit);
+ return outputPoint;
+ } catch (ExtrapolationFailureException x) {
+ // Failed to make extrapolation
+ return null;
+ }
+ }
+
+ /**
+ * Quote the tangent unit vector of the supplied track at
+ * its showering point.
+ * Returns null if not known or if track does not
+ * reach calorimeter.
+ */
+ public Hep3Vector getTangentUnit(Track tr) {
+ try {
+ BasicHep3Vector outputPoint = new BasicHep3Vector();
+ BasicHep3Vector outputTangentUnit = new BasicHep3Vector();
+ findPointAndTangent(tr, outputPoint, outputTangentUnit);
+ return outputTangentUnit;
+ } catch (ExtrapolationFailureException x) {
+ // Failed to make extrapolation
+ return null;
+ }
+ }
+
+ /**
+ * Internal helper routine. This handles caching and sanity checks, but
+ * delegates the actual calculation to the subclass. The outputs (3-vectors
+ * for the position and direction at the showering point) are passed back
+ * by over-writing the two supplied Hep3Vectors. If the calculation fails
+ * for some reason, an exception is thrown and the values of the two
+ * Hep3Vectors are undefined.
+ *
+ * @param tr The track whose showering point and tangent are to be calculated.
+ * @param outputPoint A non-null vector whose contents will be over-written with the showering point.
+ * @param outputTangentUnit A non-null vector whose contents will be over-written with the tangent unit vector at the showering point.
+ */
+ protected void findPointAndTangent(Track tr, BasicHep3Vector outputPoint, BasicHep3Vector outputTangentUnit) throws ExtrapolationFailureException {
+ Hep3Vector point = m_cachePoint.get(tr);
+ Hep3Vector tangent = m_cacheTangent.get(tr);
+ if (point == null && tangent == null) {
+ BasicHep3Vector tmpPoint = new BasicHep3Vector();
+ BasicHep3Vector tmpTangent = new BasicHep3Vector();
+ try {
+ findPointAndTangentNoCache(tr, tmpPoint, tmpTangent);
+ m_cachePoint.put(tr, tmpPoint);
+ m_cacheTangent.put(tr, tmpTangent);
+ point = tmpPoint;
+ tangent = tmpTangent;
+ } catch (ExtrapolationFailureException x) {
+ // Exception due to failure to extrapolate track to ECAL
+ throw x;
+ }
+ } else if (point!=null && tangent!=null) {
+ // retrieved from cache OK
+ } else {
+ throw new AssertionError("Cache error");
+ }
+
+ outputPoint.setV(point.x(), point.y(), point.z());
+ outputTangentUnit.setV(tangent.x(), tangent.y(), tangent.z());
+ return;
+ }
+
+ /**
+ * Calculate the showering point and unit tangent at that point given a track.
+ * Implementing class must over-ride this. The outputs (3-vectors
+ * for the position and direction at the showering point) are passed back
+ * by over-writing the two supplied Hep3Vectors. If the calculation fails
+ * for some reason, an exception is thrown and the values of the two
+ * Hep3Vectors are undefined.
+ *
+ * @param tr The track whose showering point and tangent are to be calculated.
+ * @param outputPoint A non-null vector whose contents will be over-written with the showering point.
+ * @param outputTangentUnit A non-null vector whose contents will be over-written with the tangent unit vector at the showering point.
+ */
+ protected abstract void findPointAndTangentNoCache(Track tr, BasicHep3Vector outputPoint, BasicHep3Vector outputTangentUnit) throws ExtrapolationFailureException;
+
+ protected class ExtrapolationFailureException extends java.lang.Exception {
+ public ExtrapolationFailureException(String m) { super(m); }
+ }
+
+}
lcsim/src/org/lcsim/contrib/uiowa
diff -N LayerBasedMIPGeometryHandler.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ LayerBasedMIPGeometryHandler.java 15 Aug 2008 17:32:33 -0000 1.1
@@ -0,0 +1,228 @@
+package org.lcsim.contrib.uiowa;
+
+import java.util.*;
+import org.lcsim.util.*;
+import org.lcsim.event.util.*;
+import org.lcsim.event.*;
+import hep.physics.vec.*;
+import org.lcsim.recon.pfa.identifier.LocalHelixExtrapolator;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import org.lcsim.geometry.*;
+
+/**
+ * This class provides a convenient wrapper to get the
+ * showering point and the tangent at that point. It
+ * takes a map from tracks to the pre-shower MIP cluster
+ * as input. The end user should call the
+ * getShowerPoint() and getTangentUnit() routines defined
+ * in the parent class to access these.
+ *
+ * The calculation is based on looking through the hits in
+ * the MIP cluster, splitting them up by subdetector,
+ * finding which subdetector is furthest along the track
+ * from the IP, then finding the hits in the outermost layer
+ * of that subdetector.
+ *
+ * @version $Id: LayerBasedMIPGeometryHandler.java,v 1.1 2008/08/15 17:32:33 mcharles Exp $
+ */
+
+public class LayerBasedMIPGeometryHandler extends MIPGeometryHandler {
+
+ private Map<Track, BasicCluster> m_newMapMIP;
+ private LocalHelixExtrapolator m_extrap;
+
+ public LayerBasedMIPGeometryHandler(Map<Track, BasicCluster> mapTrackToMIP, LocalHelixExtrapolator extrap) {
+ super();
+ m_extrap = extrap;
+ m_newMapMIP = mapTrackToMIP;
+ }
+
+ protected void findPointAndTangentNoCache(Track tr, BasicHep3Vector outputPoint, BasicHep3Vector outputTangentUnit) throws ExtrapolationFailureException {
+ // Make sure the extrapolation is done
+ m_extrap.performExtrapolation(tr);
+
+ // Find the MIP trace for the track
+ BasicCluster mip = m_newMapMIP.get(tr);
+ if (mip == null) {
+ throw new AssertionError("Track of class "+tr.getClass().getName()+" with p="+(new BasicHep3Vector(tr.getMomentum())).magnitude()+" has null mip!");
+ } else if (mip.getCalorimeterHits().size()==0) {
+ Hep3Vector interceptPoint = m_extrap.getLastInterceptPoint();
+ if (interceptPoint == null) {
+ throw new ExtrapolationFailureException("Failure to extrapolate");
+ } else {
+ outputPoint.setV(interceptPoint.x(), interceptPoint.y(), interceptPoint.z());
+ }
+ Hep3Vector tangent = m_extrap.getTangent();
+ outputTangentUnit.setV(tangent.x(), tangent.y(), tangent.z());
+ return;
+ }
+
+ // Find the outermost hit(s)
+
+ // 1) Split hits up by subdetector
+ Map<Subdetector, List<CalorimeterHit>> hitsInEachSubdetector = splitHitsBySubdetector(mip);
+ // 2) In each subdetector, find outermost layer hit & hits in it
+ Map<Subdetector, List<CalorimeterHit>> outermostHitsInEachSubdetector = findOutermostHitsInEachSubdetector(hitsInEachSubdetector);
+ // 3) Find which subdetector has the furthest-out hits (checked with |z|)
+ Subdetector outermostSubdet = findSubdetectorWithOutermostHits(outermostHitsInEachSubdetector);
+ // 4) For the hits in the furthest-out layer in that subdetector, find the mean position:
+ Hep3Vector outermostPosition = findMeanPositionOfHits(outermostHitsInEachSubdetector.get(outermostSubdet));
+ outputPoint.setV(outermostPosition.x(), outermostPosition.y(), outermostPosition.z());
+
+ // Check layer index:
+ int lastLayer = getLayer(outermostHitsInEachSubdetector.get(outermostSubdet).get(0));
+ for (CalorimeterHit hit : outermostHitsInEachSubdetector.get(outermostSubdet)) {
+ if (getLayer(hit) != lastLayer) { throw new AssertionError("Book-keeping error: layer mis-match"); }
+ if (hit.getSubdetector() != outermostSubdet) { throw new AssertionError("Book-keeping error"); }
+ }
+
+ // OK, now we have the position & layer. Next, we need the tangent at/near that point.
+ try {
+ Hep3Vector tangentUnit = getTangentUnit(lastLayer, outermostSubdet);
+ outputTangentUnit.setV(tangentUnit.x(), tangentUnit.y(), tangentUnit.z());
+ } catch (ExtrapolationFailureException x) {
+ // Failure...
+ Hep3Vector tangentUnit = m_extrap.getTangent();
+ if (tangentUnit != null) {
+ // Use tangent at ECAL front surface instead (iffy...)
+ outputTangentUnit.setV(tangentUnit.x(), tangentUnit.y(), tangentUnit.z());
+ } else {
+ // Complete failure
+ throw new ExtrapolationFailureException("Failure to compute tangent");
+ }
+ }
+
+ // All done
+ return;
+ }
+
+ protected Hep3Vector getTangentUnit(int lastLayer, Subdetector outermostSubdet) throws ExtrapolationFailureException {
+ // We need the tangent in the layer of interest.
+ // To do this, we extrapolate the track to the corresponding layer, and the previous one.
+ // (Or, if this is the first layer, use the next one).
+ Hep3Vector trackPointInLayer_N = null;
+ Hep3Vector trackPointInLayer_NminusOne = null;
+ int layerN = lastLayer;
+ if (lastLayer == 0) {
+ layerN = lastLayer + 1;
+ }
+
+ // To do the extrapolation, we need to know if we're looking at a barrel or an endcap subdetector.
+ String calName = outermostSubdet.getName();
+ if (calName.compareTo("HADBarrel")==0) {
+ trackPointInLayer_N = m_extrap.extendToHCALBarrelLayer(layerN);
+ trackPointInLayer_NminusOne = m_extrap.extendToHCALBarrelLayer(layerN-1);
+ } else if (calName.compareTo("HADEndcap")==0) {
+ trackPointInLayer_N = m_extrap.extendToHCALEndcapLayer(layerN);
+ trackPointInLayer_NminusOne = m_extrap.extendToHCALEndcapLayer(layerN-1);
+ } else if (calName.compareTo("EMBarrel")==0) {
+ trackPointInLayer_N = m_extrap.extendToECALBarrelLayer(layerN);
+ trackPointInLayer_NminusOne = m_extrap.extendToECALBarrelLayer(layerN-1);
+ } else if (calName.compareTo("EMEndcap")==0) {
+ trackPointInLayer_N = m_extrap.extendToECALEndcapLayer(layerN);
+ trackPointInLayer_NminusOne = m_extrap.extendToECALEndcapLayer(layerN-1);
+ } else {
+ throw new AssertionError("Calorimeter component "+calName+" not recognized!");
+ }
+
+ if (trackPointInLayer_N == null || trackPointInLayer_NminusOne==null) {
+ System.out.println("ERROR: Extrapolation failure when computing tangent for lastLayer="+lastLayer+" and subdet="+outermostSubdet.getClass().getName());
+ System.out.println("Tried to extrapolate to layer "+layerN+" and layer "+(layerN-1)+" of "+calName);
+ if (trackPointInLayer_N == null) {
+ System.out.println(" -- point in layer "+layerN+" was NULL!");
+ } else {
+ System.out.println(" -- point in layer "+layerN+" found OK: ("+trackPointInLayer_N.x()+", "+trackPointInLayer_N.y()+", "+trackPointInLayer_N.z()+")");
+ }
+ if (trackPointInLayer_NminusOne == null) {
+ System.out.println(" -- point in layer "+(layerN-1)+" was NULL!");
+ } else {
+ System.out.println(" -- point in layer "+(layerN-1)+" found OK: ("+trackPointInLayer_NminusOne.x()+", "+trackPointInLayer_NminusOne.y()+", "+trackPointInLayer_NminusOne.z()+")");
+ }
+ throw new ExtrapolationFailureException("Unable to compute tangent");
+ }
+
+ Hep3Vector tangentUnit = VecOp.unit(VecOp.sub(trackPointInLayer_N, trackPointInLayer_NminusOne));
+ return tangentUnit;
+ }
+
+ // Split hits up by subdetector
+ protected Map<Subdetector, List<CalorimeterHit>> splitHitsBySubdetector(Cluster mip) {
+ if (mip.getCalorimeterHits().size()==0) { throw new AssertionError("MIP has no hits!"); }
+ Map<Subdetector, List<CalorimeterHit>> hitsInEachSubdetector = new HashMap<Subdetector, List<CalorimeterHit>>();
+ for (CalorimeterHit hit : mip.getCalorimeterHits()) {
+ Subdetector det = hit.getSubdetector();
+ List<CalorimeterHit> hitList = hitsInEachSubdetector.get(det);
+ if (hitList == null) {
+ hitList = new Vector<CalorimeterHit>();
+ hitsInEachSubdetector.put(det, hitList);
+ }
+ hitList.add(hit);
+ }
+ return hitsInEachSubdetector;
+ }
+
+ // In each subdetector, find outermost layer hit & hits in it
+ Map<Subdetector, List<CalorimeterHit>> findOutermostHitsInEachSubdetector(Map<Subdetector, List<CalorimeterHit>> hitsInEachSubdetector) {
+ Map<Subdetector, List<CalorimeterHit>> outermostHitsInEachSubdetector = new HashMap<Subdetector, List<CalorimeterHit>>();
+ if (hitsInEachSubdetector.keySet().size()==0) { throw new AssertionError("No subdetectors supplied!"); }
+ for (Subdetector det : hitsInEachSubdetector.keySet()) {
+ List<CalorimeterHit> hits = hitsInEachSubdetector.get(det);
+ if (hits.size()==0) { throw new AssertionError("Book-keeping error"); }
+ int layerMax = -1;
+ List<CalorimeterHit> hitsInMaxLayer = new Vector<CalorimeterHit>();
+ for (CalorimeterHit hit : hits) {
+ int layer = getLayer(hit);
+ if (layer > layerMax) {
+ layerMax = layer;
+ hitsInMaxLayer.clear();
+ }
+ if (layer == layerMax) {
+ hitsInMaxLayer.add(hit);
+ }
+ }
+ if (layerMax<0) { throw new AssertionError("Book-keeping error"); }
+ if (hitsInMaxLayer.size()==0) { throw new AssertionError("Book-keeping error"); }
+ outermostHitsInEachSubdetector.put(det, hitsInMaxLayer);
+ }
+ return outermostHitsInEachSubdetector;
+ }
+
+ // Find which subdetector has the furthest-out hits (checked with |z|)
+ Subdetector findSubdetectorWithOutermostHits(Map<Subdetector, List<CalorimeterHit>> outermostHitsInEachSubdetector) {
+ if (outermostHitsInEachSubdetector.keySet().size()==0) { throw new AssertionError("No subdetectors supplied!"); }
+ Subdetector subdetectorWithMaxZ = null;
+ double maxZ = 0;
+ for (Subdetector det : outermostHitsInEachSubdetector.keySet()) {
+ List<CalorimeterHit> hits = outermostHitsInEachSubdetector.get(det);
+ if (hits.size()==0) { throw new AssertionError("No hits supplied for this subdetector!"); }
+ for (CalorimeterHit hit : hits) {
+ double modZ = Math.abs(hit.getPosition()[2]);
+ if (modZ > maxZ) {
+ maxZ = modZ;
+ subdetectorWithMaxZ = det;
+ }
+ }
+ }
+ if (subdetectorWithMaxZ == null) { throw new AssertionError("Book-keeping error when scanning over "+outermostHitsInEachSubdetector.keySet().size()+" subdetectors (maxZ="+maxZ+")"); }
+ return subdetectorWithMaxZ;
+ }
+
+ // For the hits in the furthest-out layer in that subdetector, find the mean position:
+ Hep3Vector findMeanPositionOfHits(Collection<CalorimeterHit> hitsInOutermostLayer) {
+ Hep3Vector meanPosition = new BasicHep3Vector(0.0, 0.0, 0.0);
+ for (CalorimeterHit hit : hitsInOutermostLayer) {
+ Hep3Vector position = new BasicHep3Vector(hit.getPosition());
+ double scaleFactor = 1.0 / ((double)(hitsInOutermostLayer.size()));
+ Hep3Vector contribution = VecOp.mult(scaleFactor, position);
+ meanPosition = VecOp.add(meanPosition, contribution);
+ }
+ return meanPosition;
+ }
+
+ protected int getLayer(CalorimeterHit hit) {
+ org.lcsim.geometry.IDDecoder id = hit.getIDDecoder();
+ id.setID(hit.getCellID());
+ int layer = id.getLayer();
+ return layer;
+ }
+}
lcsim/src/org/lcsim/contrib/uiowa
diff -N HelixTangentMIPGeometryHandler.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ HelixTangentMIPGeometryHandler.java 15 Aug 2008 17:32:33 -0000 1.1
@@ -0,0 +1,92 @@
+package org.lcsim.contrib.uiowa;
+
+import java.util.*;
+import org.lcsim.util.*;
+import org.lcsim.event.util.*;
+import org.lcsim.event.*;
+import hep.physics.vec.*;
+import org.lcsim.recon.pfa.identifier.LocalHelixExtrapolator;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import org.lcsim.geometry.*;
+
+/**
+ * This class provides a convenient wrapper to get the
+ * showering point and the tangent at that point. It
+ * takes a map from tracks to the pre-shower MIP cluster
+ * as input. The end user should call the
+ * getShowerPoint() and getTangentUnit() routines defined
+ * in the parent class to access these.
+ *
+ * The calculation is based on identifying the outermost hit
+ * and then checking the track helix near that point.
+ *
+ * @version $Id: HelixTangentMIPGeometryHandler.java,v 1.1 2008/08/15 17:32:33 mcharles Exp $
+ */
+
+public class HelixTangentMIPGeometryHandler extends MIPGeometryHandler {
+
+ private Map<Track, BasicCluster> m_newMapMIP;
+ private LocalHelixExtrapolator m_extrap;
+
+ public HelixTangentMIPGeometryHandler(Map<Track, BasicCluster> trkmipmap, LocalHelixExtrapolator extrap) {
+ super();
+ m_extrap = extrap;
+ m_newMapMIP = trkmipmap;
+ }
+
+ protected void findPointAndTangentNoCache(Track tr, BasicHep3Vector outputPoint, BasicHep3Vector outputTangentUnit) throws ExtrapolationFailureException {
+ m_extrap.performExtrapolation(tr); // Make sure that the extrapolation is valid.
+ BasicCluster newmip = m_newMapMIP.get(tr);
+ List<Hep3Vector> lastTwoPositions = new ArrayList<Hep3Vector>();
+ Hep3Vector last0pos = null;
+ Hep3Vector tangent = null;
+ if(newmip.getCalorimeterHits().size() < 3){ //These is no mip close to track. There should at least two hits.
+ last0pos = m_extrap.getLastInterceptPoint();
+ if (last0pos != null) {
+ tangent = m_extrap.getTangent(last0pos,tr); //tangent obtained using extrapolated track.
+ } else {
+ // Track didn't reach calorimeter
+ if (newmip.getCalorimeterHits().size()>=2) {
+ // 2+ hits => use last two hits to get position & direction
+ Hep3Vector last1pos = new BasicHep3Vector(newmip.getCalorimeterHits().get(newmip.getCalorimeterHits().size()-2).getPosition());
+ last0pos = new BasicHep3Vector(newmip.getCalorimeterHits().get(newmip.getCalorimeterHits().size()-1).getPosition());
+ //tangent = VecOp.sub(last0pos, last1pos);
+ tangent = m_extrap.getTangent(last0pos,tr); //tangent obtained using extrapolated track.
+ } else if (newmip.getCalorimeterHits().size()==1) {
+ Hep3Vector last1pos = new BasicHep3Vector(0,0,0);
+ last0pos = new BasicHep3Vector(newmip.getCalorimeterHits().get(newmip.getCalorimeterHits().size()-1).getPosition());
+ //tangent = VecOp.sub(last0pos, last1pos);
+ tangent = m_extrap.getTangent(last0pos,tr); //tangent obtained using extrapolated track.
+ } else {
+ // We don't have any information on the track position or intercept at all
+ throw new ExtrapolationFailureException("No extrapolation");
+ }
+ }
+ }else{ //There is a mip attached to track.
+ for(int i=0; i<2 ; i++){
+ CalorimeterHit hit = newmip.getCalorimeterHits().get(newmip.getCalorimeterHits().size()-1-i);
+ IDDecoder id = hit.getIDDecoder();
+ id.setID(hit.getCellID());
+ double x = id.getPosition()[0];
+ double y = id.getPosition()[1];
+ double z = id.getPosition()[2];
+ Hep3Vector v = new BasicHep3Vector(x,y,z);
+ lastTwoPositions.add(v);
+ }
+ last0pos = lastTwoPositions.get(0);
+ Hep3Vector last1pos = lastTwoPositions.get(1);
+ //Chose one
+ boolean usingExtrap = true;
+ if(usingExtrap) tangent = m_extrap.getTangent(last0pos, tr); //tangent obtained using extrapoltated track
+ else tangent = VecOp.sub(last0pos, last1pos);
+ }
+
+ Hep3Vector tangentUnit = VecOp.unit(tangent);
+ Hep3Vector showerPoint = last0pos;
+
+ outputPoint.setV(last0pos.x(), last0pos.y(), last0pos.z());
+ outputTangentUnit.setV(tangentUnit.x(), tangentUnit.y(), tangentUnit.z());
+ return;
+ }
+
+}
