diff -N KinkAnalysisDriver.java
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ KinkAnalysisDriver.java	14 Aug 2013 23:23:02 -0000	1.1
@@ -0,0 +1,177 @@

+package org.lcsim.hps.users.meeg;
+
+import hep.aida.IHistogram1D;
+import hep.aida.IHistogram2D;
+import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
+import java.util.ArrayList;
+import java.util.Collections;
+import java.util.HashMap;
+import java.util.List;
+import java.util.Map;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.MCParticle;
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.geometry.Detector;
+import org.lcsim.util.Driver;
+import org.lcsim.util.aida.AIDA;
+
+/*
+ * @author Sho Uemura <[log in to unmask]>
+ * @version $Id: KinkAnalysisDriver.java,v 1.1 2013/08/14 23:23:02 meeg Exp $
+ */
+public class KinkAnalysisDriver extends Driver {
+
+    private boolean debug = false;
+//    private double bfield = 0.5;
+    private AIDA aida = AIDA.defaultInstance();
+//    IHistogram2D sinLayer;
+    IHistogram2D dThetaLayer;
+//    IHistogram1D maxSin;
+    IHistogram1D maxDTheta;
+//    IHistogram2D sinLayerScatter;
+    IHistogram2D dThetaLayerLargest;
+    IHistogram2D dThetaLayerScatter;
+    IHistogram2D dThetaLayerNoScatter;
+    IHistogram1D maxEventDTheta;
+    IHistogram1D maxDThetaLayer;
+    int neventsWithScatter = 0;
+    int hardScatterIsLargestScatter = 0;
+
+    @Override
+    protected void detectorChanged(Detector detector) {
+//        sinLayer = aida.histogram2D("deflection vs. layer", 12, 0.5, 12.5, 1000, 0.0, 0.1);
+        dThetaLayer = aida.histogram2D("deflection in Y vs. layer", 12, 0.5, 12.5, 1000, -0.02, 0.02);
+//        maxSin = aida.histogram1D("largest deflection", 1000, 0.0, 0.1);
+        maxDTheta = aida.histogram1D("largest deflection in Y in track", 1000, -0.02, 0.02);
+//        sinLayerScatter = aida.histogram2D("deflection vs. layer with hard scatter", 12, 0.5, 12.5, 1000, 0.0, 0.1);
+        dThetaLayerScatter = aida.histogram2D("deflection in Y vs. layer, with hard scatter", 12, 0.5, 12.5, 1000, -0.02, 0.02);
+        dThetaLayerNoScatter = aida.histogram2D("deflection in Y vs. layer, no hard scatter", 12, 0.5, 12.5, 1000, -0.02, 0.02);
+        dThetaLayerLargest = aida.histogram2D("deflection in Y vs. layer, largest scatter in event", 12, 0.5, 12.5, 1000, -0.02, 0.02);
+        maxEventDTheta = aida.histogram1D("largest deflection in Y in event", 1000, -0.02, 0.02);
+        maxDThetaLayer = aida.histogram1D("layer of largest deflection in Y in event", 12, 0.5, 12.5);
+    }
+
+    @Override
+    public void process(EventHeader event) {
+        List<MCParticle> MCParticles = event.getMCParticles();
+        List<Hep3Vector> hardScatters = new ArrayList<Hep3Vector>();
+
+        for (MCParticle particle : MCParticles) {
+            if (particle.getOrigin().magnitude() > 10.0) {
+                hardScatters.add(VecOp.neg(particle.getOrigin()));
+            }
+        }
+
+
+        List<SimTrackerHit> trackerHits = event.getSimTrackerHits("TrackerHits");
+
+        Map<MCParticle, List<SimTrackerHit>> hitMap = new HashMap<MCParticle, List<SimTrackerHit>>();
+        for (SimTrackerHit hit : trackerHits) {
+            List hitList = hitMap.get(hit.getMCParticle());
+            if (hitList == null) {
+                hitList = new ArrayList<SimTrackerHit>();
+                hitMap.put(hit.getMCParticle(), hitList);
+            }
+            hitList.add(hit);
+        }
+        double maxEventDT = 0.0;
+        int maxEventDTLayer = 0;
+        boolean maxEventDTHardScatter = false;
+
+        trackloop:
+        for (MCParticle particle : hitMap.keySet()) {
+            List<SimTrackerHit> hits = hitMap.get(particle);
+            Map<Integer, SimTrackerHit> layerMap = new HashMap<Integer, SimTrackerHit>();
+            for (SimTrackerHit hit : hits) {
+                int layer = hit.getIdentifierFieldValue("layer");
+                if (layerMap.containsKey(layer)) {
+//                    System.out.format("Double hit in layer %d\n", layer);
+                    if (layerMap.get(layer).getPathLength() < hit.getPathLength()) {
+                        continue;
+                    }
+                }
+                layerMap.put(layer, hit);
+            }
+            List<Integer> layers = new ArrayList<Integer>(layerMap.keySet());
+            if (layers.size() < 8) {
+                continue;
+            }
+            Collections.sort(layers);
+            double maxDT = 0.0;
+            int maxDTLayer = 0;
+            boolean maxDTHardScatter = false;
+
+            for (int i = 0; i < layers.size() - 1; i++) {
+                SimTrackerHit hit = layerMap.get(layers.get(i));
+                
+                boolean nearHardScatter = false;
+                for (Hep3Vector scatter : hardScatters) {
+                    if (VecOp.add(hit.getPositionVec(), scatter).magnitude() < 5.0) {
+                        nearHardScatter = true;
+                    }
+                }
+
+                double y = hit.getPosition()[1];
+                double z = hit.getPosition()[2];
+
+                double lastY = 0.0;
+                double lastZ = 0.0;
+                if (i > 0) {
+                    SimTrackerHit lastHit = layerMap.get(layers.get(i - 1));
+                    lastY = lastHit.getPosition()[1];
+                    lastZ = lastHit.getPosition()[2];
+                }
+
+                SimTrackerHit nextHit = layerMap.get(layers.get(i + 1));
+                double nextY = nextHit.getPosition()[1];
+                double nextZ = nextHit.getPosition()[2];
+
+                double inAngle = Math.atan2((y - lastY), (z - lastZ));
+                double outAngle = Math.atan2((nextY - y), (nextZ - z));
+
+                double deflectionY = (outAngle - inAngle) * Math.signum(inAngle);
+//                    System.out.format("in: %f, out: %f, delta: %f\n", inAngle,outAngle,deflectionY);
+                dThetaLayer.fill(layers.get(i), deflectionY);
+                if (nearHardScatter) {
+                    dThetaLayerScatter.fill(layers.get(i), deflectionY);
+                } else {
+                    dThetaLayerNoScatter.fill(layers.get(i), deflectionY);
+                }
+                if (Math.abs(outAngle - inAngle) > Math.abs(maxDT)) {
+                    maxDT = deflectionY;
+                    maxDTLayer = layers.get(i);
+                    maxDTHardScatter = nearHardScatter;
+                }
+            }
+            if (maxDT != 0.0) {
+                maxDTheta.fill(maxDT);
+            }
+
+            if (Math.abs(maxDT) > Math.abs(maxEventDT)) {
+                maxEventDT = maxDT;
+                maxEventDTLayer = maxDTLayer;
+                maxEventDTHardScatter = maxDTHardScatter;
+            }
+        }
+
+        if (maxEventDT != 0.0) {
+            maxEventDTheta.fill(maxEventDT);
+            maxDThetaLayer.fill(maxEventDTLayer);
+            dThetaLayerLargest.fill(maxEventDTLayer, maxEventDT);
+            neventsWithScatter++;
+            if (maxEventDTHardScatter) {
+                hardScatterIsLargestScatter++;
+            }
+        }
+    }
+
+    @Override
+    protected void endOfData() {
+        System.out.format("%d events, %d had hard scatter as largest scatter\n", neventsWithScatter, hardScatterIsLargestScatter);
+    }
+
+    public void setDebug(boolean debug) {
+        this.debug = debug;
+    }
+}