1 added + 4 modified, total 5 files
lcsim/src/org/lcsim/contrib/uiowa
diff -N LocalHelixExtrapolator.java
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ LocalHelixExtrapolator.java 25 Mar 2008 22:46:01 -0000 1.1
@@ -0,0 +1,424 @@
+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;
+import org.lcsim.event.base.*;
+import org.lcsim.geometry.IDDecoder;
+
+/**
+ * Support class for local helix extrapolation, used to estimate
+ * intercept point of simulated track on an arbitrary ECAL layer.
+ * This will at some point be replaced by full tracking so it's
+ * not a good idea to depend heavily on it.
+ *
+ */
+
+public class LocalHelixExtrapolator extends Driver
+{
+
+ 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 Vector<Double> m_ECAL_barrel_layering_r;
+ protected Vector<Double> m_ECAL_endcap_layering_z;
+ protected double m_ECAL_endcap_rmin;
+ protected double m_ECAL_endcap_rmax;
+ protected double[] m_fieldStrength;
+ 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;
+ double m_track_z0 = 0.0;
+ double m_track_phi0 = 0.0;
+ double m_track_phi1 = 0.0;
+ EventHeader m_event = null;
+
+ /** Process this event, storing geometry info if needed. */
+ public void process(EventHeader event) {
+ m_event = event;
+ initGeometry(event);
+ }
+
+ /** 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_ECAL_barrel_layering_r = new Vector<Double>();
+ for (int iLayer=0; iLayer<emb.getLayering().getLayerCount(); iLayer++) {
+ double r = emb.getLayering().getDistanceToLayerSensorMid(iLayer);
+ m_ECAL_barrel_layering_r.add(new Double(r));
+ }
+ m_ECAL_endcap_layering_z = new Vector<Double>();
+ for (int iLayer=0; iLayer<eme.getLayering().getLayerCount(); iLayer++) {
+ double z = eme.getLayering().getDistanceToLayerSensorMid(iLayer);
+ m_ECAL_endcap_layering_z.add(new Double(z));
+ }
+ m_init = true;
+ }
+
+ protected double m_cutSeparation = 30.0; // 3cm
+ public void setCutSeparation(double cutSeparation) { m_cutSeparation = cutSeparation; }
+
+ // Utility routines
+ // ----------------
+
+ // ONLY "performExtrapolation" is allowed to modify this!
+ private BasicHep3Vector m_lastInterceptPoint = null;
+
+ public BasicHep3Vector getLastInterceptPoint() {
+ if (m_lastInterceptPoint == null) {
+ return null;
+ } else {
+ return new BasicHep3Vector(m_lastInterceptPoint.x(), m_lastInterceptPoint.y(), m_lastInterceptPoint.z());
+ }
+ }
+
+ public 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.
+ m_lastInterceptPoint = null;
+ return m_lastInterceptPoint;
+ }
+ // 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;
+ m_track_phi0 = phi0;
+ m_track_phi1 = phi1;
+ m_track_z0 = z0;
+ if (m_endcapValid) {
+ m_trackPoint_z = found_endcap_z;
+ m_trackPoint_phi = found_endcap_phi;
+ m_lastInterceptPoint = new BasicHep3Vector(found_endcap_x, found_endcap_y, found_endcap_z);
+ return m_lastInterceptPoint;
+ }
+ if (m_barrelValid) {
+ m_trackPoint_z = found_barrel_z;
+ m_trackPoint_phi = found_barrel_phi;
+ m_lastInterceptPoint = new BasicHep3Vector(found_barrel_x, found_barrel_y, found_barrel_z);
+ return m_lastInterceptPoint;
+ }
+
+ // No solution
+ m_lastInterceptPoint = null;
+ return m_lastInterceptPoint;
+ }
+
+ private Collection<SimTrackerHit> findHits(Track tr) {
+ // Find truth particle of track:
+ MCParticle truth = null;
+ if (tr instanceof org.lcsim.mc.fast.tracking.ReconTrack) {
+ org.lcsim.mc.fast.tracking.ReconTrack reconTr = (org.lcsim.mc.fast.tracking.ReconTrack) (tr);
+ truth = (MCParticle)(reconTr.getMCParticle());
+ } else if (tr instanceof BaseTrackMC) {
+ truth = ((BaseTrackMC)(tr)).getMCParticle();
+ }
+ // 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() == truth) {
+ hitsMatched.add(hit);
+ }
+ }
+ return hitsMatched;
+ }
+
+ public Long extendToECALLayerAndFindCell(int layer) {
+ Hep3Vector point = extendToECALLayer(layer);
+ IDDecoder id = null;
+ if (m_barrelValid) {
+ id = m_event.getDetector().getDecoder("EcalBarrHits");
+ if (id == null) { throw new AssertionError("Failed to find barrel ID"); }
+ } else if (m_endcapValid) {
+ id = m_event.getDetector().getDecoder("EcalEndcapHits");
+ if (id == null) { throw new AssertionError("Failed to find endcap ID"); }
+ }
+ if (id != null) {
+ long cell = id.findCellContainingXYZ(point);
+ id.setID(cell);
+ Hep3Vector cellCenter = id.getPositionVector();
+ return new Long(cell);
+ } else {
+ return null;
+ }
+ }
+
+ /** Assumes extrapolation has already been done. */
+ public Hep3Vector getTangent() {
+ 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));
+ return tangent;
+ }
+
+ /** Assumes extrapolation has already been done. */
+ public Hep3Vector extendToECALLayer(int layer) {
+ double dphi = (m_track_phi0 - m_track_phi1);
+ while (dphi < -Math.PI) { dphi += 2.0*Math.PI; }
+ while (dphi > Math.PI) { dphi -= 2.0*Math.PI; }
+ double x_c = m_trackParam_xc;
+ double y_c = m_trackParam_yc;
+ double radius = m_trackParam_radius;
+ double dz_by_dphi = m_trackParam_dz_by_dphi;
+
+ if (m_barrelValid) {
+ double a = radius;
+ double b = Math.sqrt(x_c*x_c + y_c*y_c);
+ double r = m_ECAL_barrel_layering_r.get(layer);
+ 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;
+ if (cos_phi_minus_theta < -1.0) {
+ // No solution
+ return null;
+ } else if (cos_phi_minus_theta == -1.0) {
+ // Unique solution: phi = theta + pi
+ dphi_to_barrel = theta + Math.PI - m_track_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; }
+ } else if (cos_phi_minus_theta == 1.0) {
+ // Unique solution: phi = theta
+ dphi_to_barrel = theta - m_track_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; }
+ } else if (cos_phi_minus_theta > 1.0) {
+ // No solution
+ return null;
+ } 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 - m_track_phi0;
+ double dphi_to_barrel_secondSolution = phi_second_solution - m_track_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; }
+ }
+ double found_barrel_z = m_track_z0 + dz_by_dphi * dphi_to_barrel;
+ double found_barrel_phi = m_track_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);
+ boolean validSolution =(found_barrel_z >= m_ECAL_barrel_zmin-m_cutSeparation && found_barrel_z <= m_ECAL_barrel_zmax+m_cutSeparation);
+ if (validSolution) {
+ return new BasicHep3Vector(found_barrel_x, found_barrel_y, found_barrel_z);
+ } else {
+ return null;
+ }
+ } else if (m_endcapValid) {
+ double previous_z = m_trackPoint_z;
+ double new_z = Math.abs(m_ECAL_endcap_layering_z.get(layer));
+ if (previous_z < 0) { new_z = -new_z; }
+ double dz = new_z - previous_z;
+ double deltaPhi = dz / m_trackParam_dz_by_dphi;
+ double phi = m_trackPoint_phi + deltaPhi;
+ double found_endcap_x = x_c + radius * Math.cos(phi);
+ double found_endcap_y = y_c + radius * Math.sin(phi);
+ double found_endcap_z = new_z;
+ double found_endcap_polar_r = Math.sqrt(found_endcap_x*found_endcap_x + found_endcap_y*found_endcap_y);
+ boolean validSolution = (found_endcap_polar_r >= m_ECAL_endcap_rmin-m_cutSeparation && found_endcap_polar_r <= m_ECAL_endcap_rmax+m_cutSeparation);
+ if (validSolution) {
+ return new BasicHep3Vector(found_endcap_x, found_endcap_y, found_endcap_z);
+ } else {
+ return null;
+ }
+ } else {
+ // No solution
+ return null;
+ }
+ }
+}
lcsim/src/org/lcsim/contrib/uiowa
diff -u -r1.6 -r1.7
--- LocalHelixExtrapolationTrackClusterMatcher.java 8 Mar 2008 20:51:30 -0000 1.6
+++ LocalHelixExtrapolationTrackClusterMatcher.java 25 Mar 2008 22:46:01 -0000 1.7
@@ -39,12 +39,15 @@
/** Simple constructor. */
public LocalHelixExtrapolationTrackClusterMatcher() {
super();
+ m_extrap= new LocalHelixExtrapolator();
}
/** Constructor, specifying additional cut to apply. */
public LocalHelixExtrapolationTrackClusterMatcher(DecisionMakerPair<Track,Cluster> extraCut) {
super();
m_extraCut = extraCut;
+ m_extrap= new LocalHelixExtrapolator();
+ add(m_extrap);
}
/** Another way to specify the additional cut. */
@@ -55,7 +58,8 @@
/** Process this event, storing geometry info if needed. */
public void process(EventHeader event) {
m_event = event;
- initGeometry(event);
+ m_extrap.process(event);
+
}
/**
@@ -65,7 +69,7 @@
*/
public Cluster matchTrackToCluster(Track tr, List<Cluster> clusters) {
if (m_debug) { System.out.println("DEBUG: "+this.getClass().getName()+" trying to match track with p="+(new BasicHep3Vector(tr.getMomentum())).magnitude()+" to a list of "+clusters.size()+" clusters."); }
- Hep3Vector point = performExtrapolation(tr);
+ Hep3Vector point = m_extrap.performExtrapolation(tr);
if (point != null) {
// Extrapolated to ECAL OK.
if (m_debug) { System.out.println("DEBUG: "+this.getClass().getName()+": extrapolated track with p="+(new BasicHep3Vector(tr.getMomentum())).magnitude()+" to r="+(Math.sqrt(point.x()*point.x()+point.y()*point.y()))+", z="+point.z()); }
@@ -104,270 +108,21 @@
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_ECAL_barrel_layering_r = new Vector<Double>();
- for (int iLayer=0; iLayer<emb.getLayering().getLayerCount(); iLayer++) {
- double r = emb.getLayering().getDistanceToLayerSensorMid(iLayer);
- m_ECAL_barrel_layering_r.add(new Double(r));
- }
- m_ECAL_endcap_layering_z = new Vector<Double>();
- for (int iLayer=0; iLayer<eme.getLayering().getLayerCount(); iLayer++) {
- double z = eme.getLayering().getDistanceToLayerSensorMid(iLayer);
- m_ECAL_endcap_layering_z.add(new Double(z));
- }
- m_init = true;
- }
-
- /**
- * Perform track matching.
- *
- * Return value is the associated cluster.
- *
- * The argumets <code>tr</code> and <code>clusters</code> are inputs.
- *
- * The argument <code>interceptPoint</code> will be changed to the (x,y,z)
- * co-ordinates of the estimated intercept point of the track on the calorimeter
- * surface. If no intercept point is found, the intercept point will be set
- * to (NaN, NaN, NaN). If a null <code>interceptPoint</code> is supplied by
- * the caller, this will cause a null pointer exception.
- */
- public Cluster matchTrackToClusterQuotingInterceptPoint(Track tr, List<Cluster> clusters, BasicHep3Vector interceptPoint) {
- Cluster clus = matchTrackToCluster(tr, clusters);
- BasicHep3Vector foundInterceptPoint = getLastInterceptPoint();
- if (foundInterceptPoint != null) {
- interceptPoint.setV(foundInterceptPoint.x(), foundInterceptPoint.y(), foundInterceptPoint.z());
- } else {
- interceptPoint.setV(Double.NaN, Double.NaN, Double.NaN);
- }
- return clus;
- }
-
- public void setCutSeparation(double cutSeparation) { m_cutSeparation = cutSeparation; }
+ public void setCutSeparation(double cutSeparation) { m_cutSeparation = cutSeparation; m_extrap.setCutSeparation(cutSeparation); }
public void setCutFirstLayer(int cutFirstLayer) { m_cutFirstLayer = cutFirstLayer; }
+ public LocalHelixExtrapolator getExtrapolator() { return m_extrap; }
+ public void setExtrapolator(LocalHelixExtrapolator extrap) { m_extrap = extrap; }
+
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 Vector<Double> m_ECAL_barrel_layering_r;
- protected Vector<Double> m_ECAL_endcap_layering_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;
- double m_track_z0 = 0.0;
- double m_track_phi0 = 0.0;
- double m_track_phi1 = 0.0;
+ protected LocalHelixExtrapolator m_extrap = null;
// Utility routines
// ----------------
-
- // ONLY "performExtrapolation" is allowed to modify this!
- private BasicHep3Vector m_lastInterceptPoint = null;
- // Don't make this public!
- protected BasicHep3Vector getLastInterceptPoint() {
- if (m_lastInterceptPoint == null) {
- return null;
- } else {
- return new BasicHep3Vector(m_lastInterceptPoint.x(), m_lastInterceptPoint.y(), m_lastInterceptPoint.z());
- }
- }
-
- 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.
- m_lastInterceptPoint = null;
- return m_lastInterceptPoint;
- }
- // 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;
- m_track_phi0 = phi0;
- m_track_phi1 = phi1;
- m_track_z0 = z0;
- if (m_endcapValid) {
- m_trackPoint_z = found_endcap_z;
- m_trackPoint_phi = found_endcap_phi;
- m_lastInterceptPoint = new BasicHep3Vector(found_endcap_x, found_endcap_y, found_endcap_z);
- return m_lastInterceptPoint;
- }
- if (m_barrelValid) {
- m_trackPoint_z = found_barrel_z;
- m_trackPoint_phi = found_barrel_phi;
- m_lastInterceptPoint = new BasicHep3Vector(found_barrel_x, found_barrel_y, found_barrel_z);
- return m_lastInterceptPoint;
- }
-
- // No solution
- m_lastInterceptPoint = null;
- return m_lastInterceptPoint;
- }
-
protected List<Cluster> findNearestClusters(Hep3Vector point, List<Cluster> clusterList) {
Map<Cluster,Double> mapClusterToDistance = new HashMap<Cluster, Double>();
List<Cluster> sortedListOfClusters = new Vector<Cluster>();
@@ -446,142 +201,8 @@
Map<T,Double> m_map;
}
- protected Collection<SimTrackerHit> findHits(Track tr) {
- // Find truth particle of track:
- MCParticle truth = null;
- if (tr instanceof org.lcsim.mc.fast.tracking.ReconTrack) {
- org.lcsim.mc.fast.tracking.ReconTrack reconTr = (org.lcsim.mc.fast.tracking.ReconTrack) (tr);
- truth = (MCParticle)(reconTr.getMCParticle());
- } else if (tr instanceof BaseTrackMC) {
- truth = ((BaseTrackMC)(tr)).getMCParticle();
- }
- // 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() == truth) {
- hitsMatched.add(hit);
- }
- }
- return hitsMatched;
- }
-
protected boolean m_debug = false;
public void setDebug(boolean debug) { m_debug = debug ; }
- protected Long extendToLayerAndFindCell(int layer) {
- Hep3Vector point = extendToLayer(layer);
- IDDecoder id = null;
- if (m_barrelValid) {
- id = m_event.getDetector().getDecoder("EcalBarrHits");
- if (id == null) { throw new AssertionError("Failed to find barrel ID"); }
- } else if (m_endcapValid) {
- id = m_event.getDetector().getDecoder("EcalEndcapHits");
- if (id == null) { throw new AssertionError("Failed to find endcap ID"); }
- }
- if (id != null) {
- long cell = id.findCellContainingXYZ(point);
- id.setID(cell);
- Hep3Vector cellCenter = id.getPositionVector();
- return new Long(cell);
- } else {
- return null;
- }
- }
-
- protected Hep3Vector extendToLayer(int layer) {
- double dphi = (m_track_phi0 - m_track_phi1);
- while (dphi < -Math.PI) { dphi += 2.0*Math.PI; }
- while (dphi > Math.PI) { dphi -= 2.0*Math.PI; }
- double x_c = m_trackParam_xc;
- double y_c = m_trackParam_yc;
- double radius = m_trackParam_radius;
- double dz_by_dphi = m_trackParam_dz_by_dphi;
-
- if (m_barrelValid) {
- double a = radius;
- double b = Math.sqrt(x_c*x_c + y_c*y_c);
- double r = m_ECAL_barrel_layering_r.get(layer);
- 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;
- if (cos_phi_minus_theta < -1.0) {
- // No solution
- return null;
- } else if (cos_phi_minus_theta == -1.0) {
- // Unique solution: phi = theta + pi
- dphi_to_barrel = theta + Math.PI - m_track_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; }
- } else if (cos_phi_minus_theta == 1.0) {
- // Unique solution: phi = theta
- dphi_to_barrel = theta - m_track_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; }
- } else if (cos_phi_minus_theta > 1.0) {
- // No solution
- return null;
- } 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 - m_track_phi0;
- double dphi_to_barrel_secondSolution = phi_second_solution - m_track_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; }
- }
- double found_barrel_z = m_track_z0 + dz_by_dphi * dphi_to_barrel;
- double found_barrel_phi = m_track_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);
- boolean validSolution =(found_barrel_z >= m_ECAL_barrel_zmin-m_cutSeparation && found_barrel_z <= m_ECAL_barrel_zmax+m_cutSeparation);
- if (validSolution) {
- return new BasicHep3Vector(found_barrel_x, found_barrel_y, found_barrel_z);
- } else {
- return null;
- }
- } else if (m_endcapValid) {
- double previous_z = m_trackPoint_z;
- double new_z = Math.abs(m_ECAL_endcap_layering_z.get(layer));
- if (previous_z < 0) { new_z = -new_z; }
- double dz = new_z - previous_z;
- double deltaPhi = dz / m_trackParam_dz_by_dphi;
- double phi = m_trackPoint_phi + deltaPhi;
- double found_endcap_x = x_c + radius * Math.cos(phi);
- double found_endcap_y = y_c + radius * Math.sin(phi);
- double found_endcap_z = new_z;
- double found_endcap_polar_r = Math.sqrt(found_endcap_x*found_endcap_x + found_endcap_y*found_endcap_y);
- boolean validSolution = (found_endcap_polar_r >= m_ECAL_endcap_rmin-m_cutSeparation && found_endcap_polar_r <= m_ECAL_endcap_rmax+m_cutSeparation);
- if (validSolution) {
- return new BasicHep3Vector(found_endcap_x, found_endcap_y, found_endcap_z);
- } else {
- return null;
- }
- } else {
- // No solution
- return null;
- }
- }
}
lcsim/src/org/lcsim/contrib/uiowa
diff -u -r1.2 -r1.3
--- LocalHelixExtrapolationTrackMIPClusterMatcher.java 19 Nov 2007 21:04:42 -0000 1.2
+++ LocalHelixExtrapolationTrackMIPClusterMatcher.java 25 Mar 2008 22:46:01 -0000 1.3
@@ -26,83 +26,14 @@
}
public Cluster matchTrackToCluster(Track tr, List<Cluster> clusters) {
- Hep3Vector point = performExtrapolation(tr);
+ Hep3Vector point = m_extrap.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));
+ Hep3Vector tangent = m_extrap.getTangent();
// Look through clusters, starting with the nearest ones:
List<Cluster> nearestClusters = findNearestClusters(point, clusters);
for (Cluster nearbyCluster : nearestClusters) {
- if (m_debug) {
- double debug_separation = proximity(point, nearbyCluster);
- int debug_layer = -1;
- CalorimeterHit debug_firstHitInECAL = findInnermostHitInECAL(nearbyCluster);
- if (debug_firstHitInECAL != null) {
- debug_layer = getLayer(debug_firstHitInECAL);
- }
- double debug_unitDotProduct = findUnitDotProduct(tangent, nearbyCluster);
- System.out.println("DEBUG: LocalHelixExtrapolationTrackMIPClusterMatcher: trying to connect a track to a MIP ("+nearbyCluster.getCalorimeterHits().size()+" hits) with separation="+debug_separation+" and first ECAL layer="+debug_layer+" and unit dot product="+debug_unitDotProduct);
- System.out.println("DEBUG: LocalHelixExtrapolationTrackMIPClusterMatcher: Extrapolated point = ("+point.x()+", "+point.y()+", "+point.z()+")");
- System.out.println("DEBUG: LocalHelixExtrapolationTrackMIPClusterMatcher: Cross-check: point = ("+
- (m_trackParam_xc + m_trackParam_radius * Math.cos(m_trackPoint_phi))+", "+
- (m_trackParam_yc + m_trackParam_radius * Math.sin(m_trackPoint_phi))+", "+
- m_trackPoint_z+")");
- System.out.println("DEBUG: LocalHelixExtrapolationTrackMIPClusterMatcher: ... since xc="+m_trackParam_xc+" and yc="+m_trackParam_yc+" and radius="+m_trackParam_radius+" and phi="+m_trackPoint_phi);
- Collection<SimTrackerHit> trackerHits = findHits(tr);
- 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;
- }
- }
- List<SimTrackerHit> otherHits = new Vector<SimTrackerHit>();
- List<SimTrackerHit> tmpHits = new Vector<SimTrackerHit>();
- tmpHits.add(hit0);
- tmpHits.add(hit1);
- tmpHits.add(hit2);
- for (SimTrackerHit hit : trackerHits) {
- if (hit==hit0 || hit==hit1 || hit==hit2) {
- // ...
- } else {
- otherHits.add(hit);
- }
- }
- for (SimTrackerHit hit : tmpHits) {
- double x = hit.getPoint()[0];
- double y = hit.getPoint()[1];
- double z = hit.getPoint()[2];
- double r = Math.sqrt(x*x + y*y);
- System.out.println("DEBUG: LocalHelixExtrapolationTrackMIPClusterMatcher: SimHit at x="+x+", y="+y+", z="+z+", r="+r);
- }
- for (SimTrackerHit hit : otherHits) {
- double x = hit.getPoint()[0];
- double y = hit.getPoint()[1];
- double z = hit.getPoint()[2];
- double r = Math.sqrt(x*x + y*y);
- System.out.println("DEBUG: LocalHelixExtrapolationTrackMIPClusterMatcher: SimHit at x="+x+", y="+y+", z="+z+", r="+r+" [extra]");
- }
- if (new Double(debug_unitDotProduct).isNaN()) {
- System.out.println("DEBUG: Dot product is NaN => dumping hits...");
- for (CalorimeterHit hit : nearbyCluster.getCalorimeterHits()) {
- double[] pos = hit.getPosition();
- System.out.println("DEBUG: Layer "+getLayer(hit)+", x="+pos[0]+", y="+pos[1]+", z="+pos[2]);
- }
- }
- }
double separation = proximity(point, nearbyCluster);
if (separation > m_cutSeparation) {
// This cluster (and therefore all subsequent ones) are too far away to pass
lcsim/src/org/lcsim/contrib/uiowa
diff -u -r1.16 -r1.17
--- ReclusterDriver.java 8 Mar 2008 20:51:30 -0000 1.16
+++ ReclusterDriver.java 25 Mar 2008 22:46:01 -0000 1.17
@@ -32,7 +32,7 @@
*
* This version is PRELIMINARY.
*
- * @version $Id: ReclusterDriver.java,v 1.16 2008/03/08 20:51:30 mcharles Exp $
+ * @version $Id: ReclusterDriver.java,v 1.17 2008/03/25 22:46:01 mcharles Exp $
* @author Mat Charles
*/
@@ -2744,7 +2744,7 @@
double pt = Math.sqrt(px*px + py*py);
double cosTheta = pz/pmag;
String printme = new String("Track with p="+p.magnitude()+" (cosTheta="+cosTheta+")");
- LocalHelixExtrapolationTrackClusterMatcher extrap = new LocalHelixExtrapolationTrackClusterMatcher();
+ LocalHelixExtrapolator extrap = new LocalHelixExtrapolator();
extrap.process(m_event);
Hep3Vector point = extrap.performExtrapolation(tr);
if (point == null) {
@@ -2842,7 +2842,7 @@
}
printme += " (dominant particle: "+maxParticle.getType().getPDGID()+" with "+clusterTruth.get(maxParticle).size()+" hits)";
- LocalHelixExtrapolationTrackClusterMatcher extrap = new LocalHelixExtrapolationTrackClusterMatcher();
+ LocalHelixExtrapolator extrap = new LocalHelixExtrapolator();
extrap.process(m_event);
Hep3Vector point = extrap.performExtrapolation(tr);
printme += " -- extrap point is ";
lcsim/src/org/lcsim/contrib/uiowa
diff -u -r1.12 -r1.13
--- ReclusterDTreeDriver.java 20 Mar 2008 02:29:18 -0000 1.12
+++ ReclusterDTreeDriver.java 25 Mar 2008 22:46:01 -0000 1.13
@@ -259,7 +259,6 @@
clustersMatchedToTracks.put(clus, tracksOfThisCluster);
tracksMatchedToClusters.put(tr, clus);
if (m_photonDebug) {
- debugTestPhoton(tr, clus);
double radius50 = radiusToCoverFractionOfPhoton(clus, 0.5);
double radius67 = radiusToCoverFractionOfPhoton(clus, 0.67);
double radius90 = radiusToCoverFractionOfPhoton(clus, 0.9);
@@ -588,7 +587,7 @@
// In any case, it clearly pointed to the calorimeter so we shouldn't
// add on the track momentum (that would be double-counting)
} else {
- Hep3Vector interceptPoint = debugTrackMatch.extendToLayer(0);
+ Hep3Vector interceptPoint = debugTrackMatch.getExtrapolator().extendToECALLayer(0);
if (interceptPoint == null) {
// No valid extrap to calorimeter
unmatchedTracksThatDontReachCalorimeter.add(tr);
@@ -2396,105 +2395,6 @@
}
}
- private void debugTestPhoton(Track tr, Cluster matchedCluster) {
- List<CalorimeterHit> allHitsEcalBarrel = m_event.get(CalorimeterHit.class, "EcalBarrDigiHits");
- List<CalorimeterHit> allHitsEcalEndcap = m_event.get(CalorimeterHit.class, "EcalEndcapDigiHits");
- Set<Long> allHitsEcal = new HashSet<Long>();
- for (CalorimeterHit hit : allHitsEcalBarrel) { allHitsEcal.add(hit.getCellID()); }
- for (CalorimeterHit hit : allHitsEcalEndcap) { allHitsEcal.add(hit.getCellID()); }
- PhotonClusterEnergyCalculator ronPhotonCalib = new PhotonClusterEnergyCalculator();
- LocalHelixExtrapolationTrackClusterMatcher debugTrackMatch = new LocalHelixExtrapolationTrackClusterMatcher();
- debugTrackMatch.process(m_event);
- List<Cluster> tmpList = new Vector<Cluster>();
- tmpList.add(matchedCluster);
- Cluster tmpMatchedCluster = debugTrackMatch.matchTrackToCluster(tr, tmpList);
- Hep3Vector interceptPoint = debugTrackMatch.extendToLayer(0);
- BaseTrackMC mctr = (BaseTrackMC)(tr);
- System.out.println("DEBUG: Matched track (ID="+mctr.getMCParticle().getType()+") with p="+mctr.getMCParticle().getMomentum().magnitude()+" to photon cluster with "+matchedCluster.getCalorimeterHits().size()+" hits E = "+ronPhotonCalib.getEnergy(matchedCluster));
- Map<MCParticle, List<SimCalorimeterHit>> truthMap = truthFromCluster(matchedCluster);
- for (MCParticle truthPart : truthMap.keySet()) {
- List<SimCalorimeterHit> hitsOfTruthPart = truthMap.get(truthPart);
- System.out.println("DEBUG: * Contribution of "+hitsOfTruthPart.size()+" hits from truth particle (ID="+truthPart.getType()+") with p = "+truthPart.getMomentum().magnitude());
- }
- Cluster coreSubCluster = matchedCluster.getClusters().get(0);
- {
- BasicCluster copyOfCoreSubCluster = new BasicCluster();
- copyOfCoreSubCluster.addCluster(coreSubCluster);
- TensorClusterPropertyCalculator calc = new TensorClusterPropertyCalculator();
- copyOfCoreSubCluster.setPropertyCalculator(calc);
- copyOfCoreSubCluster.calculateProperties();
- double[][]axes = calc.getPrincipleAxis();
- Hep3Vector coreDirection = new BasicHep3Vector(axes[0][0], axes[0][1], axes[0][2]);
- Hep3Vector corePosition = new BasicHep3Vector(calc.getPosition());
- Line line = new Line(corePosition, coreDirection);
- if (interceptPoint != null) {
- double s = line.getDistanceToPoint(interceptPoint);
- Hep3Vector poca = line.getPointAtDistance(s);
- double doca = VecOp.sub(poca, interceptPoint).magnitude();
- System.out.println("DEBUG: Distance from core centroid to intercept point: "+doca);
- }
- Hep3Vector origin = new BasicHep3Vector(0,0,0);
- double sOrigin = line.getDistanceToPoint(origin);
- Hep3Vector pocaToOrigin = line.getPointAtDistance(sOrigin);
- double docaToOrigin = VecOp.sub(pocaToOrigin, origin).magnitude();
- System.out.println("DEBUG: Impact parameter of photon core to origin: "+docaToOrigin);
- }
- Map<MCParticle, List<SimCalorimeterHit>> truthMapCore = truthFromCluster(coreSubCluster);
- for (MCParticle truthPart : truthMapCore.keySet()) {
- List<SimCalorimeterHit> hitsOfTruthPart = truthMapCore.get(truthPart);
- System.out.println("DEBUG: * CORE: Contribution of "+hitsOfTruthPart.size()+" hits from truth particle (ID="+truthPart.getType()+") with p = "+truthPart.getMomentum().magnitude());
- }
- Clusterer oldMipfinder = new TrackClusterDriver();
- List<Cluster> mipClustersInsidePhoton = oldMipfinder.createClusters(matchedCluster.getCalorimeterHits());
- System.out.println("DEBUG: * Found "+mipClustersInsidePhoton.size()+" MIPs inside photon:");
- for (Cluster clus : mipClustersInsidePhoton) {
- double energySum = 0.0;
- double hitCount = 0.0;
- int outerLayer=-1;
- int innerLayer=-1;
- boolean firstHit = true;
- for (CalorimeterHit hit : clus.getCalorimeterHits()) {
- energySum += hit.getCorrectedEnergy();
- hitCount++;
- int layer = getLayer(hit);
- if (firstHit) {
- outerLayer = layer;
- innerLayer = layer;
- firstHit = false;
- } else {
- if (layer > outerLayer) {
- outerLayer = layer;
- }
- if (layer < innerLayer) {
- innerLayer = layer;
- }
- }
- }
- double meanEnergyPerHit = energySum / hitCount;
- System.out.println("DEBUG: * "+clus.getCalorimeterHits().size()+" hits ("+innerLayer+"-"+outerLayer+") with mean energy "+meanEnergyPerHit);
- }
- Set<Long> allClusterHits = new HashSet<Long>();
- Set<Long> allCoreHits = new HashSet<Long>();
- for (CalorimeterHit hit : coreSubCluster.getCalorimeterHits()) {
- allCoreHits.add(hit.getCellID());
- }
- for (CalorimeterHit hit : matchedCluster.getCalorimeterHits()) {
- allClusterHits.add(hit.getCellID());
- }
- for (int iLayer=0; iLayer<5; iLayer++) {
- Long cellID = debugTrackMatch.extendToLayerAndFindCell(iLayer);
- if (allCoreHits.contains(cellID)) {
- System.out.println("DEBUG: * Extrap hit from layer "+iLayer+" in core");
- } else if (allClusterHits.contains(cellID)) {
- System.out.println("DEBUG: * Extrap hit from layer "+iLayer+" in cluster but not core");
- } else if (allHitsEcal.contains(cellID)) {
- System.out.println("DEBUG: * Extrap hit from layer "+iLayer+" not in cluster (but in ECAL)");
- } else {
- System.out.println("DEBUG: * Extrap hit from layer "+iLayer+" missing");
- }
- }
- }
-
private double electronEnergyNormalizedResidual(Track tr, Cluster clus) {
double energyAssumingElectron = energy(clus, m_photonCalib);
double trackMomentum = (new BasicHep3Vector(tr.getMomentum())).magnitude();
@@ -2508,18 +2408,15 @@
}
private int countHitsInClusterInFirstLayers(Track tr, Cluster clus, int nLayers) {
- LocalHelixExtrapolationTrackClusterMatcher debugTrackMatch = new LocalHelixExtrapolationTrackClusterMatcher();
- debugTrackMatch.process(m_event);
- List<Cluster> tmpList = new Vector<Cluster>();
- tmpList.add(clus);
- Cluster tmpMatchedCluster = debugTrackMatch.matchTrackToCluster(tr, tmpList);
+ LocalHelixExtrapolator extrap = new LocalHelixExtrapolator();
+ extrap.process(m_event);
Set<Long> allClusterHits = new HashSet<Long>();
for (CalorimeterHit hit : clus.getCalorimeterHits()) {
allClusterHits.add(hit.getCellID());
}
int countMatches = 0;
for (int iLayer=0; iLayer<nLayers; iLayer++) {
- Long cellID = debugTrackMatch.extendToLayerAndFindCell(iLayer);
+ Long cellID = extrap.extendToECALLayerAndFindCell(iLayer);
if (allClusterHits.contains(cellID)) {
countMatches++;
}
@@ -2528,18 +2425,15 @@
}
private int countHitsInCoreInFirstLayers(Track tr, Cluster clus, int nLayers) {
- LocalHelixExtrapolationTrackClusterMatcher debugTrackMatch = new LocalHelixExtrapolationTrackClusterMatcher();
- debugTrackMatch.process(m_event);
- List<Cluster> tmpList = new Vector<Cluster>();
- tmpList.add(clus);
- Cluster tmpMatchedCluster = debugTrackMatch.matchTrackToCluster(tr, tmpList);
+ LocalHelixExtrapolator extrap = new LocalHelixExtrapolator();
+ extrap.process(m_event);
Set<Long> coreClusterHits = new HashSet<Long>();
for (CalorimeterHit hit : clus.getClusters().get(0).getCalorimeterHits()) {
coreClusterHits.add(hit.getCellID());
}
int countMatches = 0;
for (int iLayer=0; iLayer<nLayers; iLayer++) {
- Long cellID = debugTrackMatch.extendToLayerAndFindCell(iLayer);
+ Long cellID = extrap.extendToECALLayerAndFindCell(iLayer);
if (coreClusterHits.contains(cellID)) {
countMatches++;
}
@@ -2566,12 +2460,9 @@
}
private double distanceFromTrackToPhotonCore(Track tr, Cluster clus) {
- LocalHelixExtrapolationTrackClusterMatcher debugTrackMatch = new LocalHelixExtrapolationTrackClusterMatcher();
- debugTrackMatch.process(m_event);
- List<Cluster> tmpList = new Vector<Cluster>();
- tmpList.add(clus);
- Cluster tmpMatchedCluster = debugTrackMatch.matchTrackToCluster(tr, tmpList);
- Hep3Vector interceptPoint = debugTrackMatch.extendToLayer(0);
+ LocalHelixExtrapolator extrap = new LocalHelixExtrapolator();
+ extrap.process(m_event);
+ Hep3Vector interceptPoint = extrap.performExtrapolation(tr);
if (interceptPoint != null) {
Cluster coreSubCluster = clus.getClusters().get(0);
BasicCluster copyOfCoreSubCluster = new BasicCluster();
@@ -2913,9 +2804,9 @@
}
Cluster seedClusToUse = null;
Cluster tmpMatchedClusterLayer0 = debugTrackMatch.matchTrackToCluster(tr, tmpListLayer0);
+ Hep3Vector interceptPointLayer0 = debugTrackMatch.getExtrapolator().extendToECALLayer(0);
Cluster tmpMatchedClusterLayer1 = debugTrackMatch.matchTrackToCluster(tr, tmpListLayer1);
- Hep3Vector interceptPointLayer0 = debugTrackMatch.extendToLayer(0);
- Hep3Vector interceptPointLayer1 = debugTrackMatch.extendToLayer(1);
+ Hep3Vector interceptPointLayer1 = debugTrackMatch.getExtrapolator().extendToECALLayer(1);
if (tmpMatchedClusterLayer0 == null && tmpMatchedClusterLayer1 == null) {
// No seed found for this track -- maybe extrapolation is too imprecise.
// Handle it the usual way.
@@ -2943,7 +2834,7 @@
seedHitToUse = seedClusToUse.getCalorimeterHits().iterator().next();
// Require within a certain distance (sanity check)
Hep3Vector positionOfSeedHit = new BasicHep3Vector(seedHitToUse.getPosition());
- Hep3Vector trackExtrapPointInLayer = debugTrackMatch.extendToLayer(getLayer(seedHitToUse));
+ Hep3Vector trackExtrapPointInLayer = debugTrackMatch.getExtrapolator().extendToECALLayer(getLayer(seedHitToUse));
double distForCut = VecOp.sub(positionOfSeedHit, trackExtrapPointInLayer).magnitude();
if (distForCut < 10.0) {
// Within 1cm
CVSspam 0.2.8
