lcsim/src/org/lcsim/contrib/uiowa
diff -N CheckEoverP.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ CheckEoverP.java 9 May 2007 00:25:25 -0000 1.1
@@ -0,0 +1,42 @@
+package org.lcsim.contrib.uiowa;
+
+import java.util.*;
+import org.lcsim.util.decision.*;
+import org.lcsim.event.*;
+import org.lcsim.recon.cluster.util.ClusterEnergyCalculator;
+
+public class CheckEoverP implements DecisionMakerPair<Track,Cluster>
+{
+ public CheckEoverP(ClusterEnergyCalculator calib, double allowedVariation) {
+ m_calib = calib;
+ m_allowedVariation = allowedVariation;
+ if (allowedVariation < 0.0) { throw new AssertionError("Limit must be >=0"); }
+ }
+
+ public boolean valid(Track tr, Cluster clus) {
+ // Check energy and uncertainty from calorimeter:
+ double estimatedClusterEnergy = m_calib.getEnergy(clus);
+ double estimatedClusterEnergyUncertainty = 0.7 * Math.sqrt(estimatedClusterEnergy); // 70%/sqrt(E) for hadrons
+
+ // Check energy from track
+ double[] trackMomentum = tr.getMomentum();
+ double trackMomentumMagSq = trackMomentum[0]*trackMomentum[0] + trackMomentum[1]*trackMomentum[1] + trackMomentum[2]*trackMomentum[2];
+ double trackMomentumMag = Math.sqrt(trackMomentumMagSq);
+ double massPion = 0.14;
+ double massProton = 0.94;
+ double energyReleaseIfPion = Math.sqrt(trackMomentumMagSq + massPion*massPion);
+ double energyReleaseIfProton = trackMomentumMag;
+ double energyReleaseIfAntiproton = trackMomentumMag + massProton + massProton;
+
+ // Compare them
+ boolean energyDiffOK_pion = Math.abs((energyReleaseIfPion - estimatedClusterEnergy)/estimatedClusterEnergyUncertainty) < m_allowedVariation;
+ boolean energyDiffOK_proton = Math.abs((energyReleaseIfProton - estimatedClusterEnergy)/estimatedClusterEnergyUncertainty) < m_allowedVariation;
+ boolean energyDiffOK_antiproton = Math.abs((energyReleaseIfAntiproton - estimatedClusterEnergy)/estimatedClusterEnergyUncertainty) < m_allowedVariation;
+ boolean energyDiffOK = energyDiffOK_pion || energyDiffOK_proton || energyDiffOK_antiproton;
+
+ return energyDiffOK;
+ }
+
+ protected double m_allowedVariation;
+ protected ClusterEnergyCalculator m_calib = null;
+}
lcsim/src/org/lcsim/contrib/uiowa
diff -N LocalHelixExtrapolationTrackClusterMatcher.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ LocalHelixExtrapolationTrackClusterMatcher.java 9 May 2007 00:25:25 -0000 1.1
@@ -0,0 +1,388 @@
+package org.lcsim.contrib.uiowa;
+
+import java.util.*;
+import hep.physics.vec.*;
+import org.lcsim.util.decision.*;
+
+import hep.physics.particle.Particle;
+import org.lcsim.util.swim.HelixSwimmer;
+import org.lcsim.event.Cluster;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.geometry.subdetector.CylindricalCalorimeter;
+import org.lcsim.geometry.Detector;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import org.lcsim.recon.cluster.util.TensorClusterPropertyCalculator;
+import org.lcsim.event.EventHeader;
+import org.lcsim.event.Track;
+import org.lcsim.util.Driver;
+import org.lcsim.event.MCParticle;
+import org.lcsim.mc.fast.tracking.ReconTrack;
+import org.lcsim.event.SimCalorimeterHit;
+import org.lcsim.recon.pfa.identifier.TrackClusterMatcher;
+import org.lcsim.event.SimTrackerHit;
+import org.lcsim.mc.fast.tracking.ReconTrack;
+
+/**
+ * Attempt to match a Track to a Cluster. The Track is extrapolated to the inner surface
+ * of the ECAL with a helix derived from nearby tracker hits.
+ *
+ * Certain cuts are hard-coded, but additional cuts can be specified via a
+ * DecisionMakerPair<Track,Cluster>.
+ */
+
+
+public class LocalHelixExtrapolationTrackClusterMatcher extends Driver implements TrackClusterMatcher
+{
+ /** Simple constructor. */
+ public LocalHelixExtrapolationTrackClusterMatcher() {
+ super();
+ }
+
+ /** Constructor, specifying additional cut to apply. */
+ public LocalHelixExtrapolationTrackClusterMatcher(DecisionMakerPair<Track,Cluster> extraCut) {
+ super();
+ m_extraCut = extraCut;
+ }
+
+ /** Another way to specify the additional cut. */
+ public void setExtraCheck(DecisionMakerPair<Track,Cluster> extraCut) {
+ m_extraCut = extraCut;
+ }
+
+ /** Process this event, storing geometry info if needed. */
+ public void process(EventHeader event) {
+ m_event = event;
+ initGeometry(event);
+ }
+
+ /**
+ * Match this track to a cluster from the list supplied.
+ * Returns the best matched cluster, or null if there is no
+ * acceptable match.
+ */
+ public Cluster matchTrackToCluster(Track tr, List<Cluster> clusters) {
+ Hep3Vector point = performExtrapolation(tr);
+ if (point != null) {
+ // Extrapolated to ECAL OK.
+ // Look through clusters, starting with the nearest ones:
+ List<Cluster> nearestClusters = findNearestClusters(point, clusters);
+ for (Cluster nearbyCluster : nearestClusters) {
+ double separation = proximity(point, nearbyCluster);
+ if (separation > m_cutSeparation) {
+ // This cluster (and therefore all subsequent ones) are too far away to pass
+ break;
+ } else {
+ // Separation OK. Next, check that the cluster has a hit in the first n layers of the ECAL:
+ CalorimeterHit firstHitInECAL = findInnermostHitInECAL(nearbyCluster);
+ if (firstHitInECAL != null && getLayer(firstHitInECAL) < m_cutFirstLayer) {
+ // First hit layer OK.
+ if (m_extraCut == null || m_extraCut.valid(tr,nearbyCluster) ) {
+ // Extra cut not specified or passed
+ // All cuts passed.
+ return nearbyCluster;
+ }
+ }
+ }
+ }
+ }
+ // No valid match
+ return null;
+ }
+
+ /** Initialize geometry. This is done automatically if the class is run as a driver. */
+ public void initGeometry(EventHeader event) {
+ Detector det = event.getDetector();
+ CylindricalCalorimeter emb = ((CylindricalCalorimeter) det.getSubdetectors().get("EMBarrel"));
+ CylindricalCalorimeter eme = ((CylindricalCalorimeter) det.getSubdetectors().get("EMEndcap"));
+ m_ECAL_barrel_zmin = emb.getZMin();
+ m_ECAL_barrel_zmax = emb.getZMax();
+ m_ECAL_barrel_r = emb.getLayering().getDistanceToLayerSensorMid(0);
+ m_ECAL_endcap_z = eme.getLayering().getDistanceToLayerSensorMid(0);
+ m_ECAL_endcap_rmin = eme.getInnerRadius();
+ m_ECAL_endcap_rmax = eme.getOuterRadius();
+ double[] zero = {0, 0, 0};
+ m_fieldStrength = det.getFieldMap().getField(zero);
+ m_init = true;
+ }
+
+ protected DecisionMakerPair<Track,Cluster> m_extraCut = null;
+ protected double m_cutSeparation = 30.0; // 3cm
+ protected int m_cutFirstLayer = 5; // Cluster must have a hit in in layer 0,1,2,3, or 4.
+ protected EventHeader m_event;
+ protected boolean m_init = false;
+ protected double m_ECAL_barrel_zmin;
+ protected double m_ECAL_barrel_zmax;
+ protected double m_ECAL_barrel_r;
+ protected double m_ECAL_endcap_z;
+ protected double m_ECAL_endcap_rmin;
+ protected double m_ECAL_endcap_rmax;
+ protected double[] m_fieldStrength;
+ // Cache some geometry information for subclasses
+ boolean m_barrelValid = false;
+ boolean m_endcapValid = false;
+ double m_trackParam_xc = 0.0;
+ double m_trackParam_yc = 0.0;
+ double m_trackParam_radius = 0.0;
+ double m_trackParam_dz_by_dphi = 0.0;
+ double m_trackPoint_z = 0.0;
+ double m_trackPoint_phi = 0.0;
+
+ // Utility routines
+ // ----------------
+
+ protected Hep3Vector performExtrapolation(Track tr) {
+ // Find hits. For now these are SimTrackerHits because
+ // of the inability to get hold of TrackerHits.
+ Collection<SimTrackerHit> trackerHits = findHits(tr);
+ if (trackerHits.size() < 3) {
+ // Need at least 3 hits to make a helix.
+ return null;
+ }
+ // Hit 0 is the closest to the calorimeter (confusing!)
+ SimTrackerHit hit0 = null;
+ SimTrackerHit hit1 = null;
+ SimTrackerHit hit2 = null;
+ for (SimTrackerHit hit : trackerHits) {
+ if (hit0==null || Math.abs(hit.getPoint()[2])>Math.abs(hit0.getPoint()[2])) {
+ hit2 = hit1;
+ hit1 = hit0;
+ hit0 = hit;
+ } else if (hit1==null || Math.abs(hit.getPoint()[2])>Math.abs(hit1.getPoint()[2])) {
+ hit2 = hit1;
+ hit1 = hit;
+ } else if (hit2==null || Math.abs(hit.getPoint()[2])>Math.abs(hit2.getPoint()[2])) {
+ hit2 = hit;
+ }
+ }
+ // First look at the xy (circle) projection
+ double x0 = hit0.getPoint()[0];
+ double x1 = hit1.getPoint()[0];
+ double x2 = hit2.getPoint()[0];
+ double y0 = hit0.getPoint()[1];
+ double y1 = hit1.getPoint()[1];
+ double y2 = hit2.getPoint()[1];
+ double a1 = 2.0 * (x1-x0);
+ double a2 = 2.0 * (x2-x0);
+ double b1 = 2.0 * (y1-y0);
+ double b2 = 2.0 * (y2-y0);
+ double c1 = -x0*x0 -y0*y0 +x1*x1 +y1*y1;
+ double c2 = -x0*x0 -y0*y0 +x2*x2 +y2*y2;
+ // Circle center is at (x_c, y_c)
+ double x_c = (c1*b2-c2*b1)/(a1*b2-a2*b1);
+ double y_c = (c1*a2-c2*a1)/(b1*a2-b2*a1);
+ // Circle radius
+ double radiusSquared = (x_c-x0)*(x_c-x0) + (y_c-y0)*(y_c-y0);
+ double radius = Math.sqrt(radiusSquared);
+ // Now look at z/phi projection
+ // Watch out, this is phi around the circle, not the azimuthal angle phi
+ double z0 = hit0.getPoint()[2];
+ double z1 = hit1.getPoint()[2];
+ double z2 = hit2.getPoint()[2];
+ double phi0 = Math.atan2(y0-y_c, x0-x_c); // in the range -pi through +pi
+ double phi1 = Math.atan2(y1-y_c, x1-x_c);
+ double phi2 = Math.atan2(y2-y_c, x2-x_c);
+ double dz = (z0-z1);
+ double dphi = (phi0-phi1);
+ while (dphi < -Math.PI) { dphi += 2.0*Math.PI; }
+ while (dphi > Math.PI) { dphi -= 2.0*Math.PI; }
+ double dz_by_dphi = dz/dphi;
+ // Now, try to project along to the endcaps (fairly straightforward)
+ double dz_to_endcap = Math.abs(m_ECAL_endcap_z) - z0;
+ if (z0 < 0) {
+ dz_to_endcap = -Math.abs(m_ECAL_endcap_z) - z0;
+ }
+ double dphi_to_endcap = dz_to_endcap / dz_by_dphi;
+ double found_endcap_z = z0 + dz_to_endcap;
+ double found_endcap_phi = phi0 + dphi_to_endcap;
+ double found_endcap_x = x_c + radius * Math.cos(found_endcap_phi);
+ double found_endcap_y = y_c + radius * Math.sin(found_endcap_phi);
+ double found_endcap_polar_r = Math.sqrt(found_endcap_x*found_endcap_x + found_endcap_y*found_endcap_y);
+ double found_endcap_polar_phi = Math.atan2(found_endcap_y, found_endcap_x);
+ m_endcapValid = (found_endcap_polar_r >= m_ECAL_endcap_rmin-m_cutSeparation && found_endcap_polar_r <= m_ECAL_endcap_rmax+m_cutSeparation);
+ // Now project along to the barrel (harder!)
+ // We have phi such that (x-x_c)=a*cos(phi), (y-y_c)=a*sin(phi)
+ // Define theta such that x_c = b*cos(theta), y_c = b*sin(theta)
+ double a = radius;
+ double b = Math.sqrt(x_c*x_c + y_c*y_c);
+ double r = m_ECAL_barrel_r; // barrel radius
+ double cos_phi_minus_theta = (r*r - a*a - b*b) / (2.0*a*b); // Obviously, this blows up if a or b is zero
+ double theta = Math.atan2(y_c, x_c); // in the range (-pi, +pi)
+ double dphi_to_barrel = 0.0;
+ m_barrelValid = false;
+ if (cos_phi_minus_theta < -1.0) {
+ // No solution
+ } else if (cos_phi_minus_theta == -1.0) {
+ // Unique solution: phi = theta + pi
+ dphi_to_barrel = theta + Math.PI - phi0;
+ while (dphi_to_barrel < -Math.PI) { dphi_to_barrel += 2.0*Math.PI; }
+ while (dphi_to_barrel > Math.PI) { dphi_to_barrel -= 2.0*Math.PI; }
+ m_barrelValid = true;
+ } else if (cos_phi_minus_theta == 1.0) {
+ // Unique solution: phi = theta
+ dphi_to_barrel = theta - phi0;
+ while (dphi_to_barrel < -Math.PI) { dphi_to_barrel += 2.0*Math.PI; }
+ while (dphi_to_barrel > Math.PI) { dphi_to_barrel -= 2.0*Math.PI; }
+ m_barrelValid = true;
+ } else if (cos_phi_minus_theta > 1.0) {
+ // No solution
+ } else {
+ // Two solutions
+ double phi_minus_theta_first_solution = Math.acos(cos_phi_minus_theta); // in the range 0 through pi
+ double phi_minus_theta_second_solution = -phi_minus_theta_first_solution; // in the range -pi through 0
+ double phi_first_solution = phi_minus_theta_first_solution + theta; // in the range (-pi, 2pi)
+ double phi_second_solution = phi_minus_theta_second_solution + theta; // in the range (-2pi, pi)
+ double dphi_to_barrel_firstSolution = phi_first_solution - phi0;
+ double dphi_to_barrel_secondSolution = phi_second_solution - phi0;
+ while (dphi_to_barrel_firstSolution < -Math.PI) { dphi_to_barrel_firstSolution += 2.0*Math.PI; }
+ while (dphi_to_barrel_secondSolution < -Math.PI) { dphi_to_barrel_secondSolution += 2.0*Math.PI; }
+ while (dphi_to_barrel_firstSolution > Math.PI) { dphi_to_barrel_firstSolution -= 2.0*Math.PI; }
+ while (dphi_to_barrel_secondSolution > Math.PI) { dphi_to_barrel_secondSolution -= 2.0*Math.PI; }
+ // OK, now which of the two solutions is better?
+ double test_dphi1 = dphi_to_barrel_firstSolution * dphi;
+ double test_dphi2 = dphi_to_barrel_secondSolution * dphi;
+ while (test_dphi1 < 0) { test_dphi1 += 2.0*Math.PI; }
+ while (test_dphi2 < 0) { test_dphi2 += 2.0*Math.PI; }
+ if (test_dphi1 < test_dphi2) {
+ dphi_to_barrel = dphi_to_barrel_firstSolution;
+ } else {
+ dphi_to_barrel = dphi_to_barrel_secondSolution;
+ }
+ while (dphi_to_barrel < -Math.PI) { dphi_to_barrel += 2.0*Math.PI; }
+ while (dphi_to_barrel > Math.PI) { dphi_to_barrel -= 2.0*Math.PI; }
+ m_barrelValid = true;
+ }
+ double found_barrel_z = z0 + dz_by_dphi * dphi_to_barrel;
+ double found_barrel_phi = phi0 + dphi_to_barrel;
+ double found_barrel_x = x_c + radius * Math.cos(found_barrel_phi);
+ double found_barrel_y = y_c + radius * Math.sin(found_barrel_phi);
+ double found_barrel_polar_r = Math.sqrt(found_barrel_x*found_barrel_x + found_barrel_y*found_barrel_y);
+ double found_barrel_polar_phi = Math.atan2(found_barrel_y, found_barrel_x);
+ m_barrelValid = m_barrelValid && (found_barrel_z >= m_ECAL_barrel_zmin-m_cutSeparation && found_barrel_z <= m_ECAL_barrel_zmax+m_cutSeparation);
+
+ if (m_barrelValid && m_endcapValid) {
+ // Weird case: two possible solutions
+ // Look at which is closer to last tracker hit.
+ double distanceToBarrelSq = (x0-found_barrel_x)*(x0-found_barrel_x) + (y0-found_barrel_y)*(y0-found_barrel_y) + (z0-found_barrel_z)*(z0-found_barrel_z);
+ double distanceToEndcapSq = (x0-found_endcap_x)*(x0-found_endcap_x) + (y0-found_endcap_y)*(y0-found_endcap_y) + (z0-found_endcap_z)*(z0-found_endcap_z);
+ if (distanceToBarrelSq<distanceToEndcapSq) {
+ m_endcapValid = false;
+ } else {
+ m_barrelValid = false;
+ }
+ }
+
+ m_trackParam_xc = x_c;
+ m_trackParam_yc = y_c;
+ m_trackParam_radius = radius;
+ m_trackParam_dz_by_dphi = dz_by_dphi;
+ if (m_endcapValid) {
+ m_trackPoint_z = found_barrel_z;
+ m_trackPoint_phi = found_barrel_phi;
+ return new BasicHep3Vector(found_endcap_x, found_endcap_y, found_endcap_z);
+ }
+ if (m_barrelValid) {
+ m_trackPoint_z = found_endcap_z;
+ m_trackPoint_phi = found_endcap_phi;
+ return new BasicHep3Vector(found_barrel_x, found_barrel_y, found_barrel_z);
+ }
+
+ // No solution
+ return null;
+ }
+
+ protected List<Cluster> findNearestClusters(Hep3Vector point, List<Cluster> clusterList) {
+ Map<Cluster,Double> mapClusterToDistance = new HashMap<Cluster, Double>();
+ List<Cluster> sortedListOfClusters = new Vector<Cluster>();
+ for (Cluster clus : clusterList) {
+ double dist = proximity(point, clus);
+ mapClusterToDistance.put(clus, new Double(dist));
+ sortedListOfClusters.add(clus);
+ }
+ Comparator<Cluster> comp = new CompareMapping<Cluster>(mapClusterToDistance);
+ Collections.sort(sortedListOfClusters, comp);
+ return sortedListOfClusters;
+ }
+ protected CalorimeterHit findInnermostHitInECAL(Cluster clus) {
+ CalorimeterHit innermostHit = null;
+ for (CalorimeterHit hit : clus.getCalorimeterHits()) {
+ int layer = getLayer(hit);
+ org.lcsim.geometry.Subdetector subdet = hit.getSubdetector();
+ if ( ! subdet.isCalorimeter() ) { throw new AssertionError("Cluster hit outside calorimeter"); }
+ String name = subdet.getName();
+ if (name.compareTo("EMBarrel") == 0 || name.compareTo("EMEndcap") == 0) {
+ // EM -- OK
+ if (innermostHit==null || getLayer(innermostHit)>layer) {
+ innermostHit = hit;
+ }
+ }
+ }
+ return innermostHit;
+ }
+ protected double proximity(Hep3Vector point, Cluster clus) {
+ CalorimeterHit nearestHit = findNearestHit(point, clus);
+ return proximity(point, nearestHit);
+ }
+ protected CalorimeterHit findNearestHit(Hep3Vector point, Cluster clus) {
+ CalorimeterHit nearest = null;
+ double minDist = 0;
+ for (CalorimeterHit hit : clus.getCalorimeterHits()) {
+ Hep3Vector hitPosition = new BasicHep3Vector(hit.getPosition());
+ double distance = VecOp.sub(hitPosition, point).magnitude();
+ if (distance<minDist || nearest==null) {
+ nearest = hit;
+ minDist = distance;
+ }
+ }
+ return nearest;
+ }
+ protected double proximity(Hep3Vector point, CalorimeterHit hit) {
+ Hep3Vector hitPosition = new BasicHep3Vector(hit.getPosition());
+ double distance = VecOp.sub(hitPosition, point).magnitude();
+ return distance;
+ }
+ protected int getLayer(CalorimeterHit hit) {
+ org.lcsim.geometry.IDDecoder id = hit.getIDDecoder();
+ id.setID(hit.getCellID());
+ int layer = id.getLayer();
+ return layer;
+ }
+
+ private class CompareMapping<T> implements Comparator<T> {
+ public CompareMapping(Map<T,Double> map) {
+ m_map = map;
+ }
+ public int compare(Object o1, Object o2) {
+ Cluster c1 = (Cluster) o1;
+ Cluster c2 = (Cluster) o2;
+ Double D1 = m_map.get(c1);
+ Double D2 = m_map.get(c2);
+ if (D1.equals(D2)) {
+ // Equal
+ return 0;
+ } else if (D1.doubleValue() < D2.doubleValue()) {
+ return -1;
+ } else {
+ return +1;
+ }
+ }
+ Map<T,Double> m_map;
+ }
+
+ protected Collection<SimTrackerHit> findHits(Track tr) {
+ // Find truth particle of track:
+ org.lcsim.mc.fast.tracking.ReconTrack reconTr = (org.lcsim.mc.fast.tracking.ReconTrack) (tr);
+ // Look up all hits in event:
+ List<SimTrackerHit> trackerHits = new Vector<SimTrackerHit>();
+ trackerHits.addAll(m_event.get(org.lcsim.event.SimTrackerHit.class, "TkrBarrHits"));
+ trackerHits.addAll(m_event.get(org.lcsim.event.SimTrackerHit.class, "TkrEndcapHits"));
+ // Find hits that match track:
+ List<SimTrackerHit> hitsMatched = new Vector<org.lcsim.event.SimTrackerHit>();
+ for (SimTrackerHit hit : trackerHits) {
+ if (hit.getMCParticle() == reconTr.getMCParticle()) {
+ hitsMatched.add(hit);
+ }
+ }
+ return hitsMatched;
+ }
+}
+
lcsim/src/org/lcsim/contrib/uiowa
diff -N LocalHelixExtrapolationTrackMIPClusterMatcher.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ LocalHelixExtrapolationTrackMIPClusterMatcher.java 9 May 2007 00:25:25 -0000 1.1
@@ -0,0 +1,80 @@
+package org.lcsim.contrib.uiowa;
+
+import java.util.*;
+import hep.physics.vec.*;
+import org.lcsim.util.decision.*;
+
+import org.lcsim.event.Cluster;
+import org.lcsim.event.CalorimeterHit;
+import org.lcsim.geometry.Detector;
+import org.lcsim.recon.cluster.util.BasicCluster;
+import org.lcsim.recon.cluster.util.TensorClusterPropertyCalculator;
+import org.lcsim.event.EventHeader;
+import org.lcsim.util.Driver;
+import org.lcsim.event.Track;
+
+public class LocalHelixExtrapolationTrackMIPClusterMatcher extends LocalHelixExtrapolationTrackClusterMatcher
+{
+ public LocalHelixExtrapolationTrackMIPClusterMatcher() {
+ super();
+ }
+
+ public LocalHelixExtrapolationTrackMIPClusterMatcher(DecisionMakerPair<Track,Cluster> extraCut) {
+ super(extraCut);
+ }
+
+ public Cluster matchTrackToCluster(Track tr, List<Cluster> clusters) {
+ Hep3Vector point = performExtrapolation(tr);
+ if (point != null) {
+ // Extrapolated to ECAL OK.
+ // Now, what is the tangent at that point?
+ double dphi = 0.01;
+ double dx = m_trackParam_radius * ( Math.cos(m_trackPoint_phi+dphi) - Math.cos(m_trackPoint_phi) );
+ double dy = m_trackParam_radius * ( Math.sin(m_trackPoint_phi+dphi) - Math.sin(m_trackPoint_phi) );
+ double dz = m_trackParam_dz_by_dphi * dphi;
+ Hep3Vector tangent = VecOp.unit(new BasicHep3Vector(dx,dy,dz));
+ // Look through clusters, starting with the nearest ones:
+ List<Cluster> nearestClusters = findNearestClusters(point, clusters);
+ for (Cluster nearbyCluster : nearestClusters) {
+ double separation = proximity(point, nearbyCluster);
+ if (separation > m_cutSeparation) {
+ // This cluster (and therefore all subsequent ones) are too far away to pass
+ break;
+ } else {
+ // Separation OK. Next, check that the cluster has a hit in the first n layers of the ECAL:
+ CalorimeterHit firstHitInECAL = findInnermostHitInECAL(nearbyCluster);
+ if (firstHitInECAL != null && getLayer(firstHitInECAL) < m_cutFirstLayer) {
+ // First hit layer OK. Next, check the dot-product of directions:
+ double unitDotProduct = findUnitDotProduct(tangent, nearbyCluster);
+ if (Math.abs(unitDotProduct) > m_cutDotProduct) {
+ // Dot product OK. Next, check extra cut if specified:
+ if (m_extraCut == null || m_extraCut.valid(tr,nearbyCluster) ) {
+ // Extra cut not specified or passed
+ // All cuts passed.
+ return nearbyCluster;
+ }
+ }
+ }
+ }
+ }
+ }
+ // No valid match
+ return null;
+ }
+
+ protected double findUnitDotProduct(Hep3Vector tangent, Cluster clus) {
+ // Find the cluster direction
+ BasicCluster copy = new BasicCluster();
+ copy.addCluster(clus);
+ TensorClusterPropertyCalculator calc = new TensorClusterPropertyCalculator();
+ copy.setPropertyCalculator(calc);
+ copy.calculateProperties();
+ double[][]axes = calc.getPrincipleAxis();
+ Hep3Vector clusterDir = new BasicHep3Vector(axes[0][0], axes[0][1], axes[0][2]);
+ // Get the dot product:
+ double unitDotProduct = VecOp.dot(tangent, clusterDir) / (tangent.magnitude() * clusterDir.magnitude());
+ return unitDotProduct;
+ }
+
+ protected double m_cutDotProduct = 0.85;
+}
