Commit in hps-java/src/main/java/org/lcsim/hps/users/meeg on MAIN | |||
KinkAnalysisDriver.java | +177 | added 1.1 |
simple analysis driver for track kink angles
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; + } +}
Use REPLY-ALL to reply to list
To unsubscribe from the LCD-CVS list, click the following link:
https://listserv.slac.stanford.edu/cgi-bin/wa?SUBED1=LCD-CVS&A=1