3 removed files
lcsim/src/org/lcsim/contrib/uiowa
diff -N LocalHelixExtrapolationTrackClusterMatcher.java
--- LocalHelixExtrapolationTrackClusterMatcher.java 26 Mar 2008 19:44:04 -0000 1.8
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,206 +0,0 @@
-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;
-
-/**
- * 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. We use truth information
- * to find the tracker hits, so we need to look up truth information for the Track.
- *
- * 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();
- 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. */
- 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;
- m_extrap.process(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) {
- 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 = 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()); }
- // Look through clusters, starting with the nearest ones:
- List<Cluster> nearestClusters = findNearestClusters(point, clusters);
- for (Cluster nearbyCluster : nearestClusters) {
- double separation = proximity(point, nearbyCluster);
- if (m_debug) { System.out.println("DEBUG: "+this.getClass().getName()+": comparing track with p="+(new BasicHep3Vector(tr.getMomentum())).magnitude()+" to a cluster of "+nearbyCluster.getCalorimeterHits().size()+" hits at a separation of "+separation); }
- if (separation > m_cutSeparation) {
- // This cluster (and therefore all subsequent ones) are too far away to pass
- if (m_debug) { System.out.println("DEBUG: "+this.getClass().getName()+": for track with p="+(new BasicHep3Vector(tr.getMomentum())).magnitude()+", remaining clusters are too far away (cut="+m_cutSeparation+") => no match."); }
- 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) {
- if (m_debug) { System.out.println("DEBUG: "+this.getClass().getName()+": comparing track with p="+(new BasicHep3Vector(tr.getMomentum())).magnitude()+" to a cluster of "+nearbyCluster.getCalorimeterHits().size()+" hits at a separation of "+separation+": First hit in ECAL is in layer "+getLayer(firstHitInECAL)+" -- OK!"); }
- // First hit layer OK.
- if (m_extraCut == null || m_extraCut.valid(tr,nearbyCluster) ) {
- // Extra cut not specified or passed
- // All cuts passed.
- if (m_debug) { System.out.println("DEBUG: "+this.getClass().getName()+": comparing track with p="+(new BasicHep3Vector(tr.getMomentum())).magnitude()+" to a cluster of "+nearbyCluster.getCalorimeterHits().size()+" hits at a separation of "+separation+": All cuts passed => accept"); }
- return nearbyCluster;
- } else {
- if (m_debug) { System.out.println("DEBUG: "+this.getClass().getName()+": comparing track with p="+(new BasicHep3Vector(tr.getMomentum())).magnitude()+" to a cluster of "+nearbyCluster.getCalorimeterHits().size()+" hits at a separation of "+separation+": Failed extra cut ("+m_extraCut.getClass().getName()+")"); }
- }
- } else {
- if (m_debug) { System.out.println("DEBUG: "+this.getClass().getName()+": comparing track with p="+(new BasicHep3Vector(tr.getMomentum())).magnitude()+" to a cluster of "+nearbyCluster.getCalorimeterHits().size()+" hits at a separation of "+separation+": Failed to find hit in ECAL layer < "+m_cutFirstLayer); }
- }
- }
- }
- } else {
- if (m_debug) { System.out.println("DEBUG: "+this.getClass().getName()+": failed to extrapolate track => no match"); }
- }
- // No valid match
- return null;
- }
-
- 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 LocalHelixExtrapolator m_extrap = null;
-
- // Utility routines
- // ----------------
-
- 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 boolean m_debug = false;
- public void setDebug(boolean debug) { m_debug = debug ; }
-}
-
lcsim/src/org/lcsim/contrib/uiowa
diff -N LocalHelixExtrapolationTrackMIPClusterMatcher.java
--- LocalHelixExtrapolationTrackMIPClusterMatcher.java 25 Mar 2008 22:46:01 -0000 1.3
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,81 +0,0 @@
-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;
-import org.lcsim.event.SimTrackerHit;
-import org.lcsim.event.SimCalorimeterHit;
-
-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 = m_extrap.performExtrapolation(tr);
- if (point != null) {
- // Extrapolated to ECAL OK.
- // Now, what is the tangent at that point?
- Hep3Vector tangent = m_extrap.getTangent();
- // 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());
- if (m_debug) {
- System.out.println("DEBUG: LocalHelixExtrapolationTrackMIPClusterMatcher: Computed dot product as clusterDir=("+clusterDir.x()+", "+clusterDir.y()+", "+clusterDir.z()+"), tangent=("+tangent.x()+", "+tangent.y()+", "+tangent.z()+") using "+copy.getCalorimeterHits().size()+" hits.");
- }
- return unitDotProduct;
- }
-
- protected double m_cutDotProduct = 0.85;
-}
lcsim/src/org/lcsim/contrib/uiowa
diff -N LocalHelixExtrapolator.java
--- LocalHelixExtrapolator.java 25 Jul 2008 23:02:08 -0000 1.3
+++ /dev/null 1 Jan 1970 00:00:00 -0000
@@ -1,458 +0,0 @@
-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;
- boolean m_track_dphi_negative = false;
- 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;
- m_track_dphi_negative = (dphi < 0);
- // 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;
- if (m_track_dphi_negative) { 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;
- }
-
- public Hep3Vector getTangent(Hep3Vector v, Track tr){
-
- double x0 = v.x();
- double y0 = v.y();
- double phi0 = Math.atan2(y0-m_trackParam_xc, x0-m_trackParam_yc); // in the range -pi through +pi
-
- boolean _debug =false;
- if(_debug){
- Hep3Vector P = new BasicHep3Vector(tr.getMomentum());
- Double pT = Math.sqrt(P.x()*P.x()+P.y()*P.y());
- Hep3Vector ip = new BasicHep3Vector();
- double zField = m_event.getDetector().getFieldMap().getField(P).z();
- double r = 1000*pT/(0.3*zField);
-
- System.out.println("R from track = " + m_trackParam_radius);
- System.out.println("R from p/0.3B= " + r + " magnetic field= " + zField);
- System.out.println("Charge: " + tr.getCharge());
- if(m_track_dphi_negative) System.out.println("T");
- else System.out.println("F");
- }
-
- double dphi = 0.01;
- if (m_track_dphi_negative) { dphi = -0.01; }
- double dx = m_trackParam_radius * ( Math.cos(phi0+dphi) - Math.cos(phi0) );
- double dy = m_trackParam_radius * ( Math.sin(phi0+dphi) - Math.sin(phi0) );
- 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;
- }
- }
-}
CVSspam 0.2.8
