lcsim/src/org/lcsim/contrib/uiowa/structural/likelihood
diff -u -r1.3 -r1.4
--- MiscUtilities.java 18 Oct 2005 17:42:14 -0000 1.3
+++ MiscUtilities.java 18 Oct 2005 22:32:57 -0000 1.4
@@ -116,23 +116,76 @@
double pointR = Math.sqrt(point.x()*point.x() + point.y()*point.y());
double pointZ = point.z();
if (pointR > innerR && pointR < outerR && pointZ > minZ && pointZ < maxZ) {
- // Yes.
- if (subdet.isCalorimeter()) {
- //System.out.println("DEBUG: Point in calorimeter: "+point+" is inside system '"+name+"' which is of class '"+subdet.getClass().getName()+"'");
- if (subdet.isLayered()) {
- Layering layers = subdet.getLayering();
- }
- return subdet;
- } else {
- //System.out.println("DEBUG: Point not in calorimeter: "+point+" is inside system '"+name+"' which is of class '"+subdet.getClass().getName()+"'");
- return subdet;
- }
+ return subdet;
}
- } else {
- //System.out.println("DEBUG: Subdetector '"+name+"' of class '"+subdet.getClass().getName()+"' is not cylindrical.");
}
}
- return null; // dummy
+ return null; // Couldn't find a cylindrical detector containing this point
+ }
+
+ static protected org.lcsim.geometry.layer.Layer findLayer(Hep3Vector point, org.lcsim.geometry.Detector det)
+ {
+ org.lcsim.geometry.Subdetector subdet = MiscUtilities.findComponent(point, det);
+ if (subdet == null) {
+ // Point doesn't fall inside a cylindrical subdetector
+ return null;
+ }
+ if ( ! (subdet instanceof CylindricalSubdetector) ) {
+ throw new AssertionError("Asked to find layer info for subdetector "+subdet+" of class "+subdet.getClass().getName()+", which is not a cylinder");
+ }
+ if (!subdet.isLayered()) {
+ throw new AssertionError("Asked to find layer info for subdetector "+subdet+" of class "+subdet.getClass().getName()+", which is not layered");
+ }
+ CylindricalSubdetector cylinder = (CylindricalSubdetector) (subdet);
+ Layering layers = subdet.getLayering();
+ // Barrel: Spans from -DeltaZ/2 to +DeltaZ/2
+ // Endcap: Spans from +z to +z+DeltaZ or from -z to -Z-Deltaz
+ boolean isEndcap = (cylinder.getZMin() * cylinder.getZMax() > 0.0);
+ // Now, if this is a barrel-type detector, the layers are at constant r and go outwards in steps of z.
+ // If this is an endcap-type detector, the layers are at constant z and go outwards in steps of r.
+ double longitudinalCoord = 0.0;
+ double sign = 1.0;
+ if (isEndcap) {
+ longitudinalCoord = Math.sqrt(point.x()*point.x() + point.y()*point.y());
+ } else {
+ longitudinalCoord = point.z();
+ if (longitudinalCoord<0.0) {
+ sign = -1.0;
+ }
+ }
+ // Which layer are we in?
+ // Note that we have to be careful about the sign in the case where z<0...
+ int numLayers = layers.getLayerCount();
+ org.lcsim.geometry.layer.Layer matchingLayer = null;
+ for (int iLayer=0; iLayer<numLayers; iLayer++) {
+ double layerAbsMin = sign*layers.getDistanceToLayer(iLayer);
+ double layerAbsMax = layerAbsMin + layers.getLayer(iLayer).getThickness();
+ if (sign*longitudinalCoord>=layerAbsMin && sign*longitudinalCoord<layerAbsMax) {
+ matchingLayer = layers.getLayer(iLayer);
+ break;
+ }
+ }
+ // OK, we've found the layer (or not). Do some sanity checks:
+ // Total thickness must match sum of thicknesses, and layer distances
+ // must be monotonic increasing:
+ for (int iLayer=0; iLayer<numLayers-1; iLayer++) {
+ double thickness = layers.getLayer(iLayer).getThickness();
+ double thisLayerAbsMin = sign*layers.getDistanceToLayer(iLayer);
+ double thisLayerAbsMax = thisLayerAbsMin + thickness;
+ double nextLayerAbsMin = sign*layers.getDistanceToLayer(iLayer+1);
+ if (thisLayerAbsMin > thisLayerAbsMax) {
+ throw new AssertionError("Layer "+iLayer+" min = "+thisLayerAbsMin+" and layer "+iLayer+" max = "+thisLayerAbsMax+" -- not allowed.");
+ }
+ if (thisLayerAbsMax > nextLayerAbsMin) {
+ throw new AssertionError("Layer "+(iLayer+1)+" min = "+nextLayerAbsMin+" and layer "+iLayer+" max = "+thisLayerAbsMax+" -- not allowed.");
+ }
+ // Check that there are no holes:
+ double gap = (nextLayerAbsMin - thisLayerAbsMax);
+ if ( Math.abs(gap) > 0.001) {
+ throw new AssertionError("Layer "+iLayer+" max = "+thisLayerAbsMax+" and layer "+iLayer+" min = "+nextLayerAbsMin+" so gap = "+gap+" -- not allowed.");
+ }
+ }
+ return matchingLayer;
}
}
lcsim/src/org/lcsim/contrib/uiowa/structural/likelihood
diff -u -r1.1 -r1.2
--- TrackToTrackIntermediateHits.java 29 Sep 2005 21:04:50 -0000 1.1
+++ TrackToTrackIntermediateHits.java 18 Oct 2005 22:32:57 -0000 1.2
@@ -1,10 +1,15 @@
package structural.likelihood; // package org.lcsim.recon.cluster.structural.likelihood;
+import java.util.Set;
+import java.util.HashSet;
+
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.event.Cluster;
import org.lcsim.util.swim.Trajectory;
import org.lcsim.util.swim.Line;
-import org.lcsim.event.Cluster;
import hep.physics.vec.BasicHep3Vector;
import hep.physics.vec.Hep3Vector;
+import hep.physics.vec.VecOp;
abstract public class TrackToTrackIntermediateHits implements StructuralLikelihoodQuantity
{
@@ -15,9 +20,65 @@
transient protected IntermediateHitFinder.InfoLookForHits m_infoParallel = null;
// Figure out if tracks are parallel
- protected boolean areTracksParallel(Cluster clus1, Cluster clus2)
+ protected boolean areTracksParallel(Cluster track1, Cluster track2)
{
- throw new AssertionError("Not implemented");
+ // First step: Are these guys nearly the same?
+ // They are if:
+ // 1) Their directions are nearly the same:
+ // dot product of directions > 0.95
+ // 2) They are close to each other
+ // DOCA < 1.0
+ // 3) They are ~ disjoint in layer
+ // "Overlap" = hits from both MIPs in the same layer in the same detector component
+ // Require < 2 overlaps.
+
+ // Obtain basic quantities:
+ Line line1 = MiscUtilities.makeLine(track1);
+ Line line2 = MiscUtilities.makeLine(track2);
+ double[] distancesAlongLinesToPOCAs = Line.getPOCAOfLines(line1, line2);
+ Hep3Vector poca1 = line1.getPointAtDistance(distancesAlongLinesToPOCAs[0]);
+ Hep3Vector poca2 = line2.getPointAtDistance(distancesAlongLinesToPOCAs[1]);
+ double doca = VecOp.sub(poca1,poca2).magnitude();
+ Hep3Vector dir1 = line1.getDirection();
+ Hep3Vector dir2 = line2.getDirection();
+ double dotProduct = VecOp.dot(dir1, dir2);
+
+ // Check if direction/proximity/disjoint criteria are met:
+ boolean tracksAreParallel = (dotProduct > 0.95);
+ boolean tracksMeet = (doca < 1.0);
+ int numOverlaps = countOverlaps(track1, track2);
+ boolean disjoint = (numOverlaps < 2);
+
+ return (tracksAreParallel && tracksMeet && disjoint);
+ }
+
+ /**
+ * Count the number of distinct layers in which both
+ * clusters have a hit.
+ */
+ protected int countOverlaps(Cluster clus1, Cluster clus2)
+ {
+ Set<org.lcsim.geometry.layer.Layer> overlapLayers = new HashSet<org.lcsim.geometry.layer.Layer>();
+ // For each pair of hits, check if they come from the same subdetector and have the same
+ // layer index:
+ for (CalorimeterHit hit1 : clus1.getCalorimeterHits()) {
+ org.lcsim.geometry.Subdetector subdet1 = hit1.getSubdetector();
+ hit1.getIDDecoder().setID(hit1.getCellID());
+ int layer1 = hit1.getIDDecoder().getLayer();
+ for (CalorimeterHit hit2 : clus2.getCalorimeterHits()) {
+ org.lcsim.geometry.Subdetector subdet2 = hit2.getSubdetector();
+ hit2.getIDDecoder().setID(hit2.getCellID());
+ int layer2 = hit2.getIDDecoder().getLayer();
+ if (layer1==layer2 && subdet1==subdet2) {
+ // Overlap
+ org.lcsim.geometry.layer.Layer layerObject = subdet1.getLayering().getLayer(layer1);
+ overlapLayers.add(layerObject);
+ }
+ }
+ }
+
+ // How many distinct layers have overlaps?
+ return overlapLayers.size();
}
// Make info for divergent (non-parallel) tracks
