Commit in hps-java/src/main/java/org/lcsim/hps on MAIN | |||
examples/StarterAnalysisDriver.java | +8 | -3 | 1.3 -> 1.4 |
users/meeg/LCIOTrackAnalysis.java | +2 | -2 | 1.2 -> 1.3 |
/KinkAnalysisDriver.java | +176 | -41 | 1.2 -> 1.3 |
+186 | -46 |
minor changes to example analysis driver, final version of kink analysis driver
diff -u -r1.3 -r1.4 --- StarterAnalysisDriver.java 8 Feb 2013 22:00:14 -0000 1.3 +++ StarterAnalysisDriver.java 24 Oct 2013 18:11:43 -0000 1.4 @@ -13,6 +13,7 @@
import org.lcsim.event.base.BaseRelationalTable; import org.lcsim.event.base.BaseTrackState; import org.lcsim.geometry.Detector;
+import org.lcsim.geometry.compact.Field;
import org.lcsim.hps.users.meeg.LCIOTrackAnalysis; import org.lcsim.util.Driver; import org.lcsim.util.aida.AIDA;
@@ -20,7 +21,7 @@
/* * Example analysis driver. * @author Sho Uemura <[log in to unmask]>
- * @version $Id: StarterAnalysisDriver.java,v 1.3 2013/02/08 22:00:14 meeg Exp $
+ * @version $Id: StarterAnalysisDriver.java,v 1.4 2013/10/24 18:11:43 meeg Exp $
*/ public class StarterAnalysisDriver extends Driver {
@@ -34,8 +35,12 @@
@Override protected void detectorChanged(Detector detector) {
- trkPz = aida.histogram1D("Track momentum (Pz)", 50, 0, 10.0); - trkMC = aida.histogram2D("Track momentum vs. MCParticle momentum (Pz)", 50, 0, 10.0, 50, 0, 10.0);
+ for (String str : detector.getFields().keySet()) { + System.out.println(str); + ((Field) detector.getFields().get(str)).getName(); + } + trkPz = aida.histogram1D("Track momentum (Pz)", 200, 0, 10.0); + trkMC = aida.histogram2D("Track momentum vs. MCParticle momentum (Pz)", 200, 0, 10.0, 50, 0, 10.0);
} @Override
diff -u -r1.2 -r1.3 --- LCIOTrackAnalysis.java 8 Feb 2013 22:00:14 -0000 1.2 +++ LCIOTrackAnalysis.java 24 Oct 2013 18:11:43 -0000 1.3 @@ -18,7 +18,7 @@
/** * * @author Sho Uemura <[log in to unmask]>
- * @version $Id: LCIOTrackAnalysis.java,v 1.2 2013/02/08 22:00:14 meeg Exp $
+ * @version $Id: LCIOTrackAnalysis.java,v 1.3 2013/10/24 18:11:43 meeg Exp $
*/ public class LCIOTrackAnalysis {
@@ -96,7 +96,7 @@
_hitLocationPerLayer.put(layer, new BasicHep3Vector(cl.getPosition())); _nhitsNew++;
- boolean isAxial = SvtUtils.getInstance().isAxial(SvtUtils.getInstance().getSensor(module, layer));
+ boolean isAxial = SvtUtils.getInstance().isAxial(SvtUtils.getInstance().getSensor(module, layer - 1));
if (isAxial) { _nAxialhits++; } else {
diff -u -r1.2 -r1.3 --- KinkAnalysisDriver.java 14 Aug 2013 23:53:39 -0000 1.2 +++ KinkAnalysisDriver.java 24 Oct 2013 18:11:43 -0000 1.3 @@ -2,6 +2,7 @@
import hep.aida.IHistogram1D; import hep.aida.IHistogram2D;
+import hep.physics.vec.BasicHep3Vector;
import hep.physics.vec.Hep3Vector; import hep.physics.vec.VecOp; import java.util.ArrayList;
@@ -9,20 +10,23 @@
import java.util.HashMap; import java.util.List; import java.util.Map;
+import java.util.Set;
import org.lcsim.event.EventHeader; import org.lcsim.event.MCParticle;
+import org.lcsim.event.ReconstructedParticle;
import org.lcsim.event.SimTrackerHit;
+import org.lcsim.event.Vertex;
import org.lcsim.geometry.Detector;
+import org.lcsim.units.clhep.PhysicalConstants;
import org.lcsim.util.Driver; import org.lcsim.util.aida.AIDA; /* * @author Sho Uemura <[log in to unmask]>
- * @version $Id: KinkAnalysisDriver.java,v 1.2 2013/08/14 23:53:39 meeg Exp $
+ * @version $Id: KinkAnalysisDriver.java,v 1.3 2013/10/24 18:11:43 meeg Exp $
*/ public class KinkAnalysisDriver extends Driver {
- private boolean debug = false;
private AIDA aida = AIDA.defaultInstance(); IHistogram2D dThetaLayer; IHistogram2D dThetaLayerLargestT;
@@ -33,6 +37,12 @@
IHistogram1D dThetaLargestE; IHistogram1D layerLargestT; IHistogram1D layerLargestE;
+ IHistogram2D frontDPDT; + IHistogram1D frontDP; + IHistogram1D frontDT; + IHistogram1D bsZ; + IHistogram2D twoTrackFrontDT; + IHistogram1D frontDTSum;
int neventsWithScatter = 0; int hardScatterIsLargestScatter = 0;
@@ -47,10 +57,43 @@
dThetaLargestE = aida.histogram1D("largest deflection in Y in event", 1000, -0.02, 0.02); layerLargestT = aida.histogram1D("layer of largest deflection in Y in track", 12, 0.5, 12.5); layerLargestE = aida.histogram1D("layer of largest deflection in Y in event", 12, 0.5, 12.5);
+ frontDPDT = aida.histogram2D("deflection in Y vs. energy loss, front layers", 200, 0.0, 3.0, 1000, -0.02, 0.02); + frontDP = aida.histogram1D("energy loss, front layers", 200, 0.0, 3.0); + frontDT = aida.histogram1D("deflection in Y, front layers", 1000, -0.02, 0.02); + bsZ = aida.histogram1D("vertex Z, BS constrained", 200, -10.0, 20.0); + twoTrackFrontDT = aida.histogram2D("deflection in Y, two tracks", 200, -0.02, 0.02, 200, -0.02, 0.02); + frontDTSum = aida.histogram1D("sum of deflections in Y, two tracks", 500, -0.02, 0.02);
} @Override public void process(EventHeader event) {
+ +// for (List<ReconstructedParticle> list : event.get(ReconstructedParticle.class)) { +// for (ReconstructedParticle particle : list) { +// Vertex vertex = particle.getStartVertex(); +// if (vertex != null) { +// System.out.format("%s: z %f, chisq %f\n", event.getMetaData(list).getName(), vertex.getPosition().x(), vertex.getChi2()); +// } +// } +// } + + if (event.hasCollection(ReconstructedParticle.class, "AprimeBeamspotConstrained")) { + List<ReconstructedParticle> particles = event.get(ReconstructedParticle.class, "AprimeBeamspotConstrained"); + int nvertices = 0; + Vertex theVertex = null; + for (ReconstructedParticle particle : particles) { + Vertex vertex = particle.getStartVertex(); + if (vertex != null) { + nvertices++; + theVertex = vertex; +// System.out.println(vertex.getPosition().x()); + } + } + if (nvertices == 1) { + bsZ.fill(theVertex.getPosition().x()); + } + } +
List<MCParticle> MCParticles = event.getMCParticles(); List<Hep3Vector> hardScatters = new ArrayList<Hep3Vector>();
@@ -63,42 +106,107 @@
List<SimTrackerHit> trackerHits = event.getSimTrackerHits("TrackerHits");
- Map<MCParticle, List<SimTrackerHit>> hitMap = new HashMap<MCParticle, List<SimTrackerHit>>();
+// Map<MCParticle, List<SimTrackerHit>> hitMap = new HashMap<MCParticle, List<SimTrackerHit>>();
+ Map<MCParticle, Map<Integer, SimTrackerHit>> trackMap = new HashMap<MCParticle, Map<Integer, SimTrackerHit>>();
+
for (SimTrackerHit hit : trackerHits) {
- List hitList = hitMap.get(hit.getMCParticle()); - if (hitList == null) { - hitList = new ArrayList<SimTrackerHit>(); - hitMap.put(hit.getMCParticle(), hitList);
+// List hitList = hitMap.get(hit.getMCParticle()); +// if (hitList == null) { +// hitList = new ArrayList<SimTrackerHit>(); +// hitMap.put(hit.getMCParticle(), hitList); +// } +// hitList.add(hit); + + Map<Integer, SimTrackerHit> layerMap = trackMap.get(hit.getMCParticle()); + if (layerMap == null) { + layerMap = new HashMap<Integer, SimTrackerHit>(); + trackMap.put(hit.getMCParticle(), layerMap); + } + int layer = hit.getIdentifierFieldValue("layer"); + if (layerMap.containsKey(layer)) { + boolean nearHardScatter = false; + for (Hep3Vector scatter : hardScatters) { + if (VecOp.add(hit.getPositionVec(), scatter).magnitude() < 5.0) { + nearHardScatter = true; + } + } + if (!nearHardScatter) { + hardScatters.add(VecOp.neg(hit.getPositionVec())); + } +// System.out.format("Double hit in layer %d, %s\n", layer, nearHardScatter ? "near hard scatter" : "not near hard scatter"); + if (layerMap.get(layer).getPathLength() < hit.getPathLength()) { + continue; + } + } + layerMap.put(layer, hit); + } + + List<MCParticle> particlesWithoutTracks = new ArrayList<MCParticle>(); + for (MCParticle particle : trackMap.keySet()) { + Set<Integer> layers = trackMap.get(particle).keySet(); + int pairCount = 0; + for (Integer layer : layers) { + if (layer % 2 == 0 && layers.contains(layer - 1)) { + pairCount++; + } + } + if (pairCount < 4) { + particlesWithoutTracks.add(particle); + } + } + for (MCParticle particle : particlesWithoutTracks) { + trackMap.remove(particle); + } + + MCParticle electron = null, positron = null; + for (MCParticle particle : trackMap.keySet()) { + if (particle.getCharge() > 0) { + if (positron != null) { + return; + } + positron = particle; + } + if (particle.getCharge() < 0) { + if (electron != null) { + return; + } + electron = particle;
}
- hitList.add(hit);
}
+ if (electron == null || positron == null) { + return; + } +
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); - }
+ double deflection12_ele = deflection(trackMap.get(electron), 0, 4); + double deflection12_pos = deflection(trackMap.get(positron), 0, 4); + twoTrackFrontDT.fill(deflection12_ele, deflection12_pos); + frontDTSum.fill(deflection12_ele + deflection12_pos); + + for (MCParticle particle : trackMap.keySet()) { + Map<Integer, SimTrackerHit> layerMap = trackMap.get(particle);
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;
+ double angle1 = angle(layers, layerMap, 0, 1); + Hep3Vector p1 = particle.getMomentum(); + + double angle2 = angle(layers, layerMap, 4, 5); + Hep3Vector p2 = new BasicHep3Vector(layerMap.get(layers.get(1)).getMomentum()); + + double deflection12 = (angle2 - angle1) * Math.signum(angle1); + + frontDPDT.fill(p1.magnitude() - p2.magnitude(), deflection12); + frontDP.fill(p1.magnitude() - p2.magnitude()); + frontDT.fill(deflection12); + +
for (int i = 0; i < layers.size() - 1; i++) { SimTrackerHit hit = layerMap.get(layers.get(i));
@@ -109,23 +217,16 @@
} }
- double y = hit.getPosition()[1]; - double z = hit.getPosition()[2]; - - double lastY = 0.0; - double lastZ = 0.0;
+ SimTrackerHit lastHit = null;
if (i > 0) {
- SimTrackerHit lastHit = layerMap.get(layers.get(i - 1)); - lastY = lastHit.getPosition()[1]; - lastZ = lastHit.getPosition()[2];
+ lastHit = layerMap.get(layers.get(i - 1)); +// double theta = Constants.fieldConversion*0.5*PhysicalConstants.c_light;
} 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 inAngle = angle(lastHit, hit); + double outAngle = angle(hit, nextHit);
double deflectionY = (outAngle - inAngle) * Math.signum(inAngle); // System.out.format("in: %f, out: %f, delta: %f\n", inAngle,outAngle,deflectionY);
@@ -165,12 +266,46 @@
} }
+ private double angle(SimTrackerHit hit1, SimTrackerHit hit2) { + double y1 = hit2.getMCParticle().getOriginY(); +// double z1 = hit2.getMCParticle().getOriginZ(); + double s1 = hit2.getMCParticle().getProductionTime() * PhysicalConstants.c_light; + if (hit1 != null) { + y1 = hit1.getPosition()[1]; +// z1 = hit1.getPosition()[2]; + s1 = hit1.getTime() * PhysicalConstants.c_light; + } + + double y2 = hit2.getPosition()[1]; +// double z2 = hit2.getPosition()[2]; + double s2 = hit2.getTime() * PhysicalConstants.c_light; +// double outAngle = Math.atan2((nextY - y), (nextZ - z)); + + return Math.asin((y2 - y1) / (s2 - s1)); + } + + private double angle(List<Integer> layers, Map<Integer, SimTrackerHit> layerMap, int layer1, int layer2) { + SimTrackerHit hit1 = null; + if (layer1 > 0) { + hit1 = layerMap.get(layers.get(layer1 - 1)); + } + SimTrackerHit hit2 = layerMap.get(layers.get(layer2 - 1)); + + return angle(hit1, hit2); + } + + private double deflection(Map<Integer, SimTrackerHit> layerMap, int layer1, int layer2) { + List<Integer> layers = new ArrayList<Integer>(layerMap.keySet()); + Collections.sort(layers); + + double angle1 = angle(layers, layerMap, layer1, layer1 + 1); + double angle2 = angle(layers, layerMap, layer2, layer2 + 1); + + return (angle2 - angle1) * Math.signum(angle1); + } +
@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