Author: [log in to unmask]
Date: Thu Jul 23 12:15:16 2015
New Revision: 3275
Log:
add version of extrapolator that uses field edge from geometry
Modified:
java/trunk/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java
Modified: java/trunk/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java
=============================================================================
--- java/trunk/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java (original)
+++ java/trunk/tracking/src/main/java/org/hps/recon/tracking/TrackUtils.java Thu Jul 23 12:15:16 2015
@@ -36,6 +36,7 @@
import org.lcsim.fit.helicaltrack.HelixUtils;
import org.lcsim.fit.helicaltrack.HitUtils;
import org.lcsim.fit.helicaltrack.MultipleScatter;
+import org.lcsim.geometry.Detector;
import org.lcsim.geometry.subdetector.BarrelEndcapFlag;
import org.lcsim.recon.tracking.seedtracker.SeedCandidate;
import org.lcsim.recon.tracking.seedtracker.SeedTrack;
@@ -43,8 +44,7 @@
/**
* Assorted helper functions for the track and helix objects in lcsim. Re-use as
- * much of HelixUtils
- * as possible.
+ * much of HelixUtils as possible.
*
* @author Omar Moreno <[log in to unmask]>
*/
@@ -60,8 +60,7 @@
/**
* Extrapolate track to a position along the x-axis. Turn the track into a
- * helix object in
- * order to use HelixUtils.
+ * helix object in order to use HelixUtils.
*
* @param track
* @param x
@@ -130,9 +129,8 @@
// ==========================================================================
/**
* Calculate the point of interception between the helix and a plane in
- * space. Uses an iterative procedure.
- * This function makes assumptions on the sign and convecntion of the
- * B-field. Be careful.
+ * space. Uses an iterative procedure. This function makes assumptions on
+ * the sign and convecntion of the B-field. Be careful.
*
* @param helfit - helix
* @param unit_vec_normal_to_plane - unit vector normal to the plane
@@ -146,18 +144,19 @@
//WTrack wtrack = new WTrack(helfit, -1.0*bfield); //
Hep3Vector B = new BasicHep3Vector(0, 0, 1);
WTrack wtrack = new WTrack(helfit, bfield); //
- if (debug)
+ if (debug) {
System.out.printf("getHelixPlaneIntercept:find intercept between plane defined by point on plane %s, unit vec %s, bfield %.3f, h=%s and WTrack \n%s \n", point_on_plane.toString(), unit_vec_normal_to_plane.toString(), bfield, B.toString(), wtrack.toString());
+ }
Hep3Vector intercept_point = wtrack.getHelixAndPlaneIntercept(point_on_plane, unit_vec_normal_to_plane, B);
- if (debug)
+ if (debug) {
System.out.printf("getHelixPlaneIntercept: found intercept point at %s\n", intercept_point.toString());
+ }
return intercept_point;
}
/**
* Calculate the point of interception between the helix and a plane in
- * space. Uses an
- * iterative procedure.
+ * space. Uses an iterative procedure.
*
* @param helfit - helix
* @param strip - strip cluster that will define the plane
@@ -220,8 +219,8 @@
*/
public static Hep3Vector extrapolateTrack(Track track, double z) {
- Hep3Vector trackPosition = null;
- double dz = 0;
+ Hep3Vector trackPosition;
+ double dz;
if (z >= BeamlineConstants.DIPOLE_EDGE_TESTRUN) {
trackPosition = extrapolateHelixToXPlane(track, BeamlineConstants.DIPOLE_EDGE_TESTRUN);
dz = z - BeamlineConstants.DIPOLE_EDGE_TESTRUN;
@@ -233,6 +232,50 @@
// System.out.printf("detVec %s\n", detVecTracking.toString());
return new BasicHep3Vector(detVecTracking.y(), detVecTracking.z(), detVecTracking.x());
}
+
+ // Get the track azimuthal angle
+ double phi = getPhi(track, trackPosition);
+
+ // Find the distance to the point of interest
+ double r = dz / (getSinTheta(track) * Math.cos(phi));
+ double dx = r * getSinTheta(track) * Math.sin(phi);
+ double dy = r * getCosTheta(track);
+
+ // Find the track position at the point of interest
+ double x = trackPosition.y() + dx;
+ double y = trackPosition.z() + dy;
+
+ return new BasicHep3Vector(x, y, z);
+ }
+
+ /**
+ * Extrapolate track to given position, using dipole position from geometry.
+ *
+ * @param helix - to be extrapolated
+ * @param track - position along the x-axis of the helix in lcsim
+ * coordiantes
+ * @return
+ */
+ public static Hep3Vector extrapolateTrack(Track track, double z, Detector detector) {
+
+ Hep3Vector trackPosition;
+// <constant name="dipoleMagnetPositionZ" value="45.72*cm"/>
+// <constant name="dipoleMagnetLength" value="108*cm"/>
+
+ double magnetLength = detector.getConstants().get("dipoleMagnetLength").getValue();
+ double magnetZ = detector.getConstants().get("dipoleMagnetPositionZ").getValue();
+ double magnetDownstreamEdge = magnetZ + magnetLength / 2;
+ double magnetUpstreamEdge = magnetZ - magnetLength / 2;
+ if (z >= magnetDownstreamEdge) {
+ trackPosition = extrapolateHelixToXPlane(track, magnetDownstreamEdge);
+ } else if (z <= magnetUpstreamEdge) {
+ trackPosition = extrapolateHelixToXPlane(track, magnetUpstreamEdge);
+ } else {
+ Hep3Vector detVecTracking = extrapolateHelixToXPlane(track, z);
+ // System.out.printf("detVec %s\n", detVecTracking.toString());
+ return new BasicHep3Vector(detVecTracking.y(), detVecTracking.z(), detVecTracking.x());
+ }
+ double dz = z - trackPosition.x();
// Get the track azimuthal angle
double phi = getPhi(track, trackPosition);
@@ -294,8 +337,9 @@
// Check if we are on the plane
d = VecOp.dot(VecOp.sub(pos, origin), normal);
dx += -1.0 * d / 2.0;
- if (debug)
+ if (debug) {
System.out.printf("%d d %.10f pos [%.10f %.10f %.10f] dx %.10f\n", nIter, d, pos.x(), pos.y(), pos.z(), dx);
+ }
nIter += 1;
}
return pos;
@@ -322,10 +366,12 @@
double sinPhi = (xc - x) / R;
double phi_at_x = Math.asin(sinPhi);
double dphi_at_x = phi_at_x - phi0;
- if (dphi_at_x > Math.PI)
+ if (dphi_at_x > Math.PI) {
dphi_at_x -= 2.0 * Math.PI;
- if (dphi_at_x < -Math.PI)
+ }
+ if (dphi_at_x < -Math.PI) {
dphi_at_x += 2.0 * Math.PI;
+ }
double s_at_x = -1.0 * dphi_at_x * R;
double y = dca * Math.cos(phi0) - R * Math.cos(phi0) + R * Math.cos(phi_at_x);
double z = z0 + s_at_x * slope;
@@ -333,8 +379,9 @@
// System.out.printf("pos %s xc %f phi_at_x %f dphi_at_x %f s_at_x %f\n",
// pos.toString(),xc,phi_at_x,dphi_at_x,s_at_x);
Hep3Vector posXCheck = TrackUtils.extrapolateHelixToXPlane(helix, x);
- if (VecOp.sub(pos, posXCheck).magnitude() > 0.0000001)
+ if (VecOp.sub(pos, posXCheck).magnitude() > 0.0000001) {
throw new RuntimeException(String.format("ERROR the helix propagation equations do not agree? (%f,%f,%f) vs (%f,%f,%f) in HelixUtils", pos.x(), pos.y(), pos.z(), posXCheck.x(), posXCheck.y(), posXCheck.z()));
+ }
return pos;
}
@@ -354,56 +401,56 @@
Box sensorSolid = (Box) sensor.getGeometry().getLogicalVolume().getSolid();
Polygon3D sensorFace = sensorSolid.getFacesNormalTo(new BasicHep3Vector(0, 0, 1)).get(0);
- if (debug)
+ if (debug) {
System.out.println("sensorContainsTrack: Track Position: " + trackPosition.toString());
+ }
List<Point3D> vertices = new ArrayList<Point3D>();
- for (int index = 0; index < 4; index++)
+ for (int index = 0; index < 4; index++) {
vertices.add(new Point3D());
- for (Point3D vertex : sensorFace.getVertices())
+ }
+ for (Point3D vertex : sensorFace.getVertices()) {
if (vertex.y() < 0 && vertex.x() > 0) {
localToGlobal.transform(vertex);
// vertices.set(0, new Point3D(vertex.y() + sensorPos.x(), vertex.x() +
// sensorPos.y(), vertex.z() + sensorPos.z()));
vertices.set(0, new Point3D(vertex.x(), vertex.y(), vertex.z()));
- if (debug)
+ if (debug) {
System.out.println("sensorContainsTrack: Vertex 1 Position: " + vertices.get(0).toString()); // System.out.println("sensorContainsTrack: Transformed Vertex 1 Position: " +
- // localToGlobal.transformed(vertex).toString());
+ } // localToGlobal.transformed(vertex).toString());
} else if (vertex.y() > 0 && vertex.x() > 0) {
localToGlobal.transform(vertex);
// vertices.set(1, new Point3D(vertex.y() + sensorPos.x(), vertex.x() +
// sensorPos.y(), vertex.z() + sensorPos.z()));
vertices.set(1, new Point3D(vertex.x(), vertex.y(), vertex.z()));
- if (debug)
+ if (debug) {
System.out.println("sensorContainsTrack: Vertex 2 Position: " + vertices.get(1).toString()); // System.out.println("sensorContainsTrack: Transformed Vertex 2 Position: " +
- // localToGlobal.transformed(vertex).toString());
+ } // localToGlobal.transformed(vertex).toString());
} else if (vertex.y() > 0 && vertex.x() < 0) {
localToGlobal.transform(vertex);
// vertices.set(2, new Point3D(vertex.y() + sensorPos.x(), vertex.x() +
// sensorPos.y(), vertex.z() + sensorPos.z()));
vertices.set(2, new Point3D(vertex.x(), vertex.y(), vertex.z()));
- if (debug)
+ if (debug) {
System.out.println("sensorContainsTrack: Vertex 3 Position: " + vertices.get(2).toString()); // System.out.println("sensorContainsTrack: Transformed Vertex 3 Position: " +
- // localToGlobal.transformed(vertex).toString());
+ } // localToGlobal.transformed(vertex).toString());
} else if (vertex.y() < 0 && vertex.x() < 0) {
localToGlobal.transform(vertex);
// vertices.set(3, new Point3D(vertex.y() + sensorPos.x(), vertex.x() +
// sensorPos.y(), vertex.z() + sensorPos.z()));
vertices.set(3, new Point3D(vertex.x(), vertex.y(), vertex.z()));
- if (debug)
+ if (debug) {
System.out.println("sensorContainsTrack: Vertex 4 Position: " + vertices.get(3).toString()); // System.out.println("sensorContainsTrack: Transformed Vertex 4 Position: " +
- // localToGlobal.transformed(vertex).toString());
- }
+ } // localToGlobal.transformed(vertex).toString());
+ }
+ }
double area1 = TrackUtils.findTriangleArea(vertices.get(0).x(), vertices.get(0).y(), vertices.get(1).x(), vertices.get(1).y(), trackPosition.y(), trackPosition.z());
double area2 = TrackUtils.findTriangleArea(vertices.get(1).x(), vertices.get(1).y(), vertices.get(2).x(), vertices.get(2).y(), trackPosition.y(), trackPosition.z());
double area3 = TrackUtils.findTriangleArea(vertices.get(2).x(), vertices.get(2).y(), vertices.get(3).x(), vertices.get(3).y(), trackPosition.y(), trackPosition.z());
double area4 = TrackUtils.findTriangleArea(vertices.get(3).x(), vertices.get(3).y(), vertices.get(0).x(), vertices.get(0).y(), trackPosition.y(), trackPosition.z());
- if ((area1 > 0 && area2 > 0 && area3 > 0 && area4 > 0) || (area1 < 0 && area2 < 0 && area3 < 0 && area4 < 0))
- return true;
-
- return false;
+ return (area1 > 0 && area2 > 0 && area3 > 0 && area4 > 0) || (area1 < 0 && area2 < 0 && area3 < 0 && area4 < 0);
}
public static Map<String, Double> calculateTrackHitResidual(HelicalTrackHit hth, HelicalTrackFit track, boolean includeMS) {
@@ -475,7 +522,7 @@
HelicalTrackStripGbl stripGbl = new HelicalTrackStripGbl(strip, true);
return calculateLocalTrackHitResiduals(track, hth, stripGbl, bFieldInZ);
}
-
+
public static Map<String, Double> calculateLocalTrackHitResiduals(Track track, HelicalTrackHit hth, HelicalTrackStripGbl strip, double bFieldInZ) {
SeedTrack st = (SeedTrack) track;
@@ -489,21 +536,19 @@
public static Map<String, Double> calculateLocalTrackHitResiduals(HelicalTrackFit _trk, HelicalTrackStrip strip, double bFieldInZ) {
HelicalTrackStripGbl stripGbl = new HelicalTrackStripGbl(strip, true);
- return calculateLocalTrackHitResiduals( _trk, stripGbl, 0.0, 0.0, bFieldInZ);
- }
-
-
-
+ return calculateLocalTrackHitResiduals(_trk, stripGbl, 0.0, 0.0, bFieldInZ);
+ }
+
public static Map<String, Double> calculateLocalTrackHitResiduals(HelicalTrackFit _trk, HelicalTrackStripGbl strip, double msdrdphi, double msdz, double bFieldInZ) {
boolean debug = false;
boolean includeMS = true;
if (debug) {
- System.out.printf("calculateLocalTrackHitResiduals: for strip on sensor %s \n",
- ((RawTrackerHit)strip.getStrip().rawhits().get(0)).getDetectorElement().getName());
- }
-
+ System.out.printf("calculateLocalTrackHitResiduals: for strip on sensor %s \n",
+ ((RawTrackerHit) strip.getStrip().rawhits().get(0)).getDetectorElement().getName());
+ }
+
Hep3Vector u = strip.u();
Hep3Vector corigin = strip.origin();
@@ -511,16 +556,16 @@
Hep3Vector trkpos = TrackUtils.getHelixPlaneIntercept(_trk, strip, Math.abs(bFieldInZ));
if (debug) {
- System.out.printf("calculateLocalTrackHitResiduals: strip u %s origin %s \n", u.toString(),corigin.toString());
+ System.out.printf("calculateLocalTrackHitResiduals: strip u %s origin %s \n", u.toString(), corigin.toString());
System.out.printf("calculateLocalTrackHitResiduals: found interception point with sensor at %s \n", trkpos.toString());
}
-
+
if (Double.isNaN(trkpos.x()) || Double.isNaN(trkpos.y()) || Double.isNaN(trkpos.z())) {
System.out.printf("calculateLocalTrackHitResiduals: failed to get interception point (%s) \n", trkpos.toString());
System.out.printf("calculateLocalTrackHitResiduals: track params\n%s\n", _trk.toString());
System.out.printf("calculateLocalTrackHitResiduals: track pT=%.3f chi2=[%.3f][%.3f] \n", _trk.pT(bFieldInZ), _trk.chisq()[0], _trk.chisq()[1]);
- trkpos = TrackUtils.getHelixPlaneIntercept(_trk, strip, bFieldInZ);
- System.exit(1);
+// trkpos = TrackUtils.getHelixPlaneIntercept(_trk, strip, bFieldInZ);
+ throw new RuntimeException();
}
double xint = trkpos.x();
@@ -547,12 +592,11 @@
double wmeas = 0;
double wError = 10.0 / Math.sqrt(12); // 0.001;
-
- if (debug) {
- System.out.printf("calculateLocalTrackHitResiduals: vdiffTrk %s vdiff %s umc %f umeas %f du %f\n",
- vdiffTrk.toString(),vdiff.toString(),umc, umeas, umeas-umc);
- }
-
+ if (debug) {
+ System.out.printf("calculateLocalTrackHitResiduals: vdiffTrk %s vdiff %s umc %f umeas %f du %f\n",
+ vdiffTrk.toString(), vdiff.toString(), umc, umeas, umeas - umc);
+ }
+
Map<String, Double> res = new HashMap<String, Double>();
res.put("ures", umeas - umc);
res.put("ureserr", includeMS ? Math.sqrt(uError * uError + msuError * msuError) : uError);
@@ -573,38 +617,37 @@
HelicalTrackHit hth = (HelicalTrackHit) hit;
//===> if (SvtUtils.getInstance().isTopLayer((SiSensor) ((RawTrackerHit) hth.getRawHits().get(0)).getDetectorElement())) {
HpsSiSensor sensor = ((HpsSiSensor) ((RawTrackerHit) hth.getRawHits().get(0)).getDetectorElement());
- if (sensor.isTopLayer())
+ if (sensor.isTopLayer()) {
n[0] = n[0] + 1;
- else
+ } else {
n[1] = n[1] + 1;
+ }
}
return n;
}
public static boolean isTopTrack(Track track, int minhits) {
- return isTopOrBottomTrack(track, minhits) == 1 ? true : false;
+ return isTopOrBottomTrack(track, minhits) == 1;
}
public static boolean isBottomTrack(Track track, int minhits) {
- return isTopOrBottomTrack(track, minhits) == 0 ? true : false;
+ return isTopOrBottomTrack(track, minhits) == 0;
}
public static int isTopOrBottomTrack(Track track, int minhits) {
int nhits[] = getHitsInTopBottom(track);
- if (nhits[0] >= minhits && nhits[1] == 0)
+ if (nhits[0] >= minhits && nhits[1] == 0) {
return 1;
- else if (nhits[1] >= minhits && nhits[0] == 0)
+ } else if (nhits[1] >= minhits && nhits[0] == 0) {
return 0;
- else
+ } else {
return -1;
+ }
}
public static boolean hasTopBotHit(Track track) {
int nhits[] = getHitsInTopBottom(track);
- if (nhits[0] > 0 && nhits[1] > 0)
- return true;
- else
- return false;
+ return nhits[0] > 0 && nhits[1] > 0;
}
public static boolean isSharedHit(TrackerHit hit, List<Track> othertracks) {
@@ -613,9 +656,10 @@
List<TrackerHit> hitsOnTrack = track.getTrackerHits();
for (TrackerHit loop_hit : hitsOnTrack) {
HelicalTrackHit loop_hth = (HelicalTrackHit) loop_hit;
- if (hth.equals(loop_hth))
- // System.out.printf("share hit at layer %d at %s (%s) with track w/ chi2=%f\n",hth.Layer(),hth.getCorrectedPosition().toString(),loop_hth.getCorrectedPosition().toString(),track.getChi2());
+ if (hth.equals(loop_hth)) // System.out.printf("share hit at layer %d at %s (%s) with track w/ chi2=%f\n",hth.Layer(),hth.getCorrectedPosition().toString(),loop_hth.getCorrectedPosition().toString(),track.getChi2());
+ {
return true;
+ }
}
}
return false;
@@ -627,22 +671,25 @@
// System.out.printf("look for another track with chi2=%f and px=%f \n",track.getChi2(),track.getTrackStates().get(0).getMomentum()[0]);
for (Track t : tracklist) {
// System.out.printf("add track with chi2=%f and px=%f ?\n",t.getChi2(),t.getTrackStates().get(0).getMomentum()[0]);
- if (t.equals(track))
- // System.out.printf("NOPE\n");
+ if (t.equals(track)) // System.out.printf("NOPE\n");
+ {
continue;
+ }
// System.out.printf("YEPP\n");
tracks.add(t);
}
List<TrackerHit> hitsOnTrack = track.getTrackerHits();
int n_shared = 0;
- for (TrackerHit hit : hitsOnTrack)
- if (isSharedHit(hit, tracks))
+ for (TrackerHit hit : hitsOnTrack) {
+ if (isSharedHit(hit, tracks)) {
++n_shared;
+ }
+ }
return n_shared;
}
public static boolean hasSharedHits(Track track, List<Track> tracklist) {
- return numberOfSharedHits(track, tracklist) == 0 ? false : true;
+ return numberOfSharedHits(track, tracklist) != 0;
}
public static void cut(int cuts[], EventQuality.Cut bit) {
@@ -651,22 +698,27 @@
public static boolean isGoodTrack(Track track, List<Track> tracklist, EventQuality.Quality trk_quality) {
int cuts = passTrackSelections(track, tracklist, trk_quality);
- return cuts == 0 ? true : false;
+ return cuts == 0;
}
public static int passTrackSelections(Track track, List<Track> tracklist, EventQuality.Quality trk_quality) {
int cuts[] = {0};
if (trk_quality.compareTo(Quality.NONE) != 0) {
- if (track.getTrackStates().get(0).getMomentum()[0] < EventQuality.instance().getCutValue(EventQuality.Cut.PZ, trk_quality))
+ if (track.getTrackStates().get(0).getMomentum()[0] < EventQuality.instance().getCutValue(EventQuality.Cut.PZ, trk_quality)) {
cut(cuts, EventQuality.Cut.PZ);
- if (track.getChi2() >= EventQuality.instance().getCutValue(EventQuality.Cut.CHI2, trk_quality))
+ }
+ if (track.getChi2() >= EventQuality.instance().getCutValue(EventQuality.Cut.CHI2, trk_quality)) {
cut(cuts, EventQuality.Cut.CHI2);
- if (numberOfSharedHits(track, tracklist) > ((int) Math.round(EventQuality.instance().getCutValue(EventQuality.Cut.SHAREDHIT, trk_quality))))
+ }
+ if (numberOfSharedHits(track, tracklist) > ((int) Math.round(EventQuality.instance().getCutValue(EventQuality.Cut.SHAREDHIT, trk_quality)))) {
cut(cuts, EventQuality.Cut.SHAREDHIT);
- if (hasTopBotHit(track))
+ }
+ if (hasTopBotHit(track)) {
cut(cuts, EventQuality.Cut.TOPBOTHIT);
- if (track.getTrackerHits().size() < ((int) Math.round(EventQuality.instance().getCutValue(EventQuality.Cut.NHITS, trk_quality))))
+ }
+ if (track.getTrackerHits().size() < ((int) Math.round(EventQuality.instance().getCutValue(EventQuality.Cut.NHITS, trk_quality)))) {
cut(cuts, EventQuality.Cut.NHITS);
+ }
}
return cuts[0];
}
@@ -681,8 +733,7 @@
/**
* Transform MCParticle into a Helix object. Note that it produces the helix
- * parameters at
- * nominal x=0 and assumes that there is no field at x<0
+ * parameters at nominal x=0 and assumes that there is no field at x<0
*
* @param mcp MC particle to be transformed
* @return helix object based on the MC particle
@@ -696,8 +747,9 @@
Hep3Vector org = CoordinateTransformations.transformVectorToTracking(mcp.getOrigin());
Hep3Vector p = CoordinateTransformations.transformVectorToTracking(mcp.getMomentum());
- if (debug)
+ if (debug) {
System.out.printf("mcp org %s mc p %s (trans)\n", org.toString(), p.toString());
+ }
// Move to x=0 if needed
double targetX = BeamlineConstants.DIPOLE_EDGELOW_TESTRUN;
@@ -708,15 +760,17 @@
double y = delta_x * dydx + org.y();
double z = delta_x * dzdx + org.z();
double x = org.x() + delta_x;
- if (Math.abs(x - targetX) > 1e-8)
+ if (Math.abs(x - targetX) > 1e-8) {
throw new RuntimeException("Error: origin is not zero!");
+ }
org = new BasicHep3Vector(x, y, z);
// System.out.printf("org %s p %s -> org %s\n",
// old.toString(),p.toString(),org.toString());
}
- if (debug)
+ if (debug) {
System.out.printf("mcp org %s mc p %s (trans2)\n", org.toString(), p.toString());
+ }
HelixParamCalculator helixParamCalculator = new HelixParamCalculator(p, org, -1 * ((int) mcp.getCharge()), Bz);
double par[] = new double[5];
@@ -732,17 +786,19 @@
}
public static HelicalTrackFit getHTF(Track track) {
- if (track.getClass().isInstance(SeedTrack.class))
+ if (track.getClass().isInstance(SeedTrack.class)) {
return ((SeedTrack) track).getSeedCandidate().getHelix();
- else
+ } else {
return getHTF(track.getTrackStates().get(0).getParameters());
+ }
}
public static HelicalTrackFit getHTF(double par[]) {
// need to have matrix that makes sense? Really?
SymmetricMatrix cov = new SymmetricMatrix(5);
- for (int i = 0; i < cov.getNRows(); ++i)
+ for (int i = 0; i < cov.getNRows(); ++i) {
cov.setElement(i, i, 1.);
+ }
HelicalTrackFit htf = new HelicalTrackFit(par, cov, new double[2], new int[2], null, null);
return htf;
}
@@ -755,8 +811,9 @@
if (useFringe) {
// broken because you need ot provide the Field Map to get this...
// pos1 = hpstrk1.getPositionAtZMap(100.0, zVal, 5.0)[0];
- } else
+ } else {
pos1 = TrackUtils.extrapolateTrack(trk1, zVal);
+ }
// System.out.printf("%s: Position1 at edge of fringe %s\n",this.getClass().getSimpleName(),pos1.toString());
Helix traj = (Helix) hpstrk1.getTrajectory();
if (traj == null) {
@@ -774,19 +831,22 @@
Map<MCParticle, Integer> particlesOnTrack = new HashMap<MCParticle, Integer>();
- if (debug)
+ if (debug) {
System.out.printf("getMatchedTruthParticle: getmatched mc particle from %d tracker hits on the track \n", track.getTrackerHits().size());
+ }
for (TrackerHit hit : track.getTrackerHits()) {
List<MCParticle> mcps = ((HelicalTrackHit) hit).getMCParticles();
- if (mcps == null)
+ if (mcps == null) {
System.out.printf("getMatchedTruthParticle: warning, this hit (layer %d pos=%s) has no mc particles.\n", ((HelicalTrackHit) hit).Layer(), ((HelicalTrackHit) hit).getCorrectedPosition().toString());
- else {
- if (debug)
+ } else {
+ if (debug) {
System.out.printf("getMatchedTruthParticle: this hit (layer %d pos=%s) has %d mc particles.\n", ((HelicalTrackHit) hit).Layer(), ((HelicalTrackHit) hit).getCorrectedPosition().toString(), mcps.size());
+ }
for (MCParticle mcp : mcps) {
- if (!particlesOnTrack.containsKey(mcp))
+ if (!particlesOnTrack.containsKey(mcp)) {
particlesOnTrack.put(mcp, 0);
+ }
int c = particlesOnTrack.get(mcp);
particlesOnTrack.put(mcp, c + 1);
}
@@ -795,23 +855,26 @@
if (debug) {
System.out.printf("Track p=[ %f, %f, %f] \n", track.getTrackStates().get(0).getMomentum()[0], track.getTrackStates().get(0).getMomentum()[1], track.getTrackStates().get(0).getMomentum()[1]);
System.out.printf("Found %d particles\n", particlesOnTrack.size());
- for (Map.Entry<MCParticle, Integer> entry : particlesOnTrack.entrySet())
+ for (Map.Entry<MCParticle, Integer> entry : particlesOnTrack.entrySet()) {
System.out.printf("%d hits assigned to %d p=%s \n", entry.getValue(), entry.getKey().getPDGID(), entry.getKey().getMomentum().toString());
+ }
}
Map.Entry<MCParticle, Integer> maxEntry = null;
- for (Map.Entry<MCParticle, Integer> entry : particlesOnTrack.entrySet())
- if (maxEntry == null || entry.getValue().compareTo(maxEntry.getValue()) > 0)
+ for (Map.Entry<MCParticle, Integer> entry : particlesOnTrack.entrySet()) {
+ if (maxEntry == null || entry.getValue().compareTo(maxEntry.getValue()) > 0) {
maxEntry = entry; //if ( maxEntry != null ) {
- // if(entry.getValue().compareTo(maxEntry.getValue()) < 0) continue;
- //}
+ } // if(entry.getValue().compareTo(maxEntry.getValue()) < 0) continue;
+ } //}
//maxEntry = entry;
- if (debug)
- if (maxEntry != null)
+ if (debug) {
+ if (maxEntry != null) {
System.out.printf("Matched particle with pdgId=%d and mom %s to track with charge %d and momentum [%f %f %f]\n",
maxEntry.getKey().getPDGID(), maxEntry.getKey().getMomentum().toString(),
track.getCharge(), track.getTrackStates().get(0).getMomentum()[0], track.getTrackStates().get(0).getMomentum()[1], track.getTrackStates().get(0).getMomentum()[2]);
- else
+ } else {
System.out.printf("No truth particle found on this track\n");
+ }
+ }
return maxEntry == null ? null : maxEntry.getKey();
}
@@ -835,8 +898,9 @@
// now get the hits and make them helicaltrackhits
List<TrackerHit> rth = track.getTrackerHits();
List<HelicalTrackHit> hth = new ArrayList<>();
- for (TrackerHit hit : rth)
+ for (TrackerHit hit : rth) {
hth.add(makeHelicalTrackHitFromTrackerHit(hit));
+ }
// SeedCandidate(List<HelicalTrackHit> , SeedStrategy strategy, HelicalTrackFit helix, double bfield) ;
SeedCandidate scand = new SeedCandidate(hth, null, htf, 0.24);
SeedTrack st = new SeedTrack();
@@ -845,10 +909,10 @@
}
/*
- cast TrackerHit as HTH...this is pretty dumb; just a rearrangment of information
- in TrackerHit. The important information that's in HTH but not
- in Tracker hit is the HelicalTrackCrosses (and from there the individual strip clusters)
- is lost; some work to get them back.
+ cast TrackerHit as HTH...this is pretty dumb; just a rearrangment of information
+ in TrackerHit. The important information that's in HTH but not
+ in Tracker hit is the HelicalTrackCrosses (and from there the individual strip clusters)
+ is lost; some work to get them back.
*/
public static HelicalTrackHit makeHelicalTrackHitFromTrackerHit(TrackerHit hit) {
Hep3Vector pos = new BasicHep3Vector(hit.getPosition());
@@ -862,7 +926,7 @@
BarrelEndcapFlag beflag = BarrelEndcapFlag.BARREL;
return new HelicalTrackHit(pos, hitcov, dedx, time, type, rhits, detname, layer, beflag);
}
-
+
public static RelationalTable getHitToStripsTable(EventHeader event) {
RelationalTable hitToStrips = new BaseRelationalTable(RelationalTable.Mode.MANY_TO_MANY, RelationalTable.Weighting.UNWEIGHTED);
List<LCRelation> hitrelations = event.get(LCRelation.class, "HelicalTrackHitRelations");
@@ -900,7 +964,7 @@
meanTime /= nStrips;
return meanTime;
}
-
+
public static boolean hasSharedStrips(Track track1, Track track2, RelationalTable hitToStrips, RelationalTable hitToRotated) {
Set<TrackerHit> track1hits = new HashSet<TrackerHit>();
for (TrackerHit hit : track1.getTrackerHits()) {
@@ -919,7 +983,7 @@
public static int getLayer(TrackerHit strip) {
return ((RawTrackerHit) strip.getRawHits().get(0)).getLayerNumber();
}
-
+
/**
* Compute strip isolation, defined as the minimum distance to another strip
* in the same sensor. Strips are only checked if they formed valid crosses
@@ -936,13 +1000,32 @@
for (TrackerHit cross : (Set<TrackerHit>) hitToStrips.allTo(otherStrip)) {
for (TrackerHit crossStrip : (Set<TrackerHit>) hitToStrips.allFrom(cross)) {
if (crossStrip != strip && crossStrip != otherStrip) {
+ int stripMin = Integer.MAX_VALUE;
+ int stripMax = Integer.MIN_VALUE;
+ int crossMin = Integer.MAX_VALUE;
+ int crossMax = Integer.MIN_VALUE;
+ for (Object rth : strip.getRawHits()) {
+ int hitStrip = ((RawTrackerHit) rth).getIdentifierFieldValue("strip");
+ stripMin = Math.min(stripMin, hitStrip);
+ stripMax = Math.max(stripMax, hitStrip);
+ }
+ for (Object rth : crossStrip.getRawHits()) {
+ int hitStrip = ((RawTrackerHit) rth).getIdentifierFieldValue("strip");
+ crossMin = Math.min(crossMin, hitStrip);
+ crossMax = Math.max(crossMax, hitStrip);
+ }
+ if (stripMin - crossMax <= 1 && crossMin - stripMax <= 1) {
+ continue; //adjacent strips don't count
+ }
Hep3Vector stripPosition = new BasicHep3Vector(strip.getPosition());
Hep3Vector crossStripPosition = new BasicHep3Vector(crossStrip.getPosition());
double distance = VecOp.sub(stripPosition, crossStripPosition).magnitude();
-// System.out.format("%s, %s, %s, %f\n", stripPosition, crossStripPosition, VecOp.sub(stripPosition, crossStripPosition), distance);
- if (distance < nearestDistance) {
- nearestDistance = distance;
- }
+// if (distance<=0.0601) continue; //hack to avoid counting adjacent strips that didn't get clustered together
+// if (distance<0.1) System.out.format("%d, %d, %f\n",strip.getRawHits().size(),crossStrip.getRawHits().size(), distance);
+// if (distance < 0.1) {
+// System.out.format("%s, %s, %s, %f\n", stripPosition, crossStripPosition, VecOp.sub(stripPosition, crossStripPosition), distance);
+// }
+ nearestDistance = Math.min(nearestDistance, distance);
}
}
}
|
