Author: [log in to unmask]
Date: Sun Dec 6 13:05:52 2015
New Revision: 4011
Log:
updating cluster/track matching
Modified:
java/branches/nathan-dev/recon/src/main/java/org/hps/recon/particle/ReconParticleDriver.java
java/branches/nathan-dev/recon/src/main/java/org/hps/recon/utils/TrackClusterMatcher.java
Modified: java/branches/nathan-dev/recon/src/main/java/org/hps/recon/particle/ReconParticleDriver.java
=============================================================================
--- java/branches/nathan-dev/recon/src/main/java/org/hps/recon/particle/ReconParticleDriver.java (original)
+++ java/branches/nathan-dev/recon/src/main/java/org/hps/recon/particle/ReconParticleDriver.java Sun Dec 6 13:05:52 2015
@@ -267,30 +267,41 @@
((BaseReconstructedParticle) particle).setParticleIdUsed(new SimpleParticleID(11, 0, 0, 0));
}
+ // the cluster matched to the track:
Cluster matchedCluster = null;
-
- // Track the best matching cluster for the track.
- // TODO: This should find the best match not just the first match.
- clusterLoop:
- for (Cluster cluster : clusters) {
- // Check if the cluster and track are a valid match.
- if (matcher.isMatch(cluster, track)) {
-
- // Store the matched cluster to matched track.
- clusterToTrack.put(cluster, track);
-
- // Store the matched cluster.
- matchedCluster = cluster;
-
- // Since a match has been found, the loop can be
- // terminated.
- break clusterLoop;
- }
+
+ // normalized cluster-track distance required for qualifying as a match:
+ final double maximumNSigma=5;
+
+ // only match clusters to GBL tracks:
+ if (track.getType() >= 32) {
+
+ // normalized distance of the closest match:
+ double smallestNSigma=Double.MAX_VALUE;
+
+ for (Cluster cluster : clusters) {
+
+ // normalized distance between this cluster and track:
+ double thisNSigma=matcher.getNSigmaFromMeanGBL(cluster, track);
+
+ // ignore if distance doesn't make the cut:
+ if (thisNSigma > maximumNSigma) continue;
+
+ // keep the cluster with the smallest normalized distance to the track:
+ if (thisNSigma < smallestNSigma) {
+ smallestNSigma = thisNSigma;
+ matchedCluster = cluster;
+ }
+ }
}
-
+
// If a cluster was found that matches the track...
if (matchedCluster != null) {
- particle.addCluster(matchedCluster);
+
+ // store the matching for later cluster corrections:
+ clusterToTrack.put(matchedCluster, track);
+
+ particle.addCluster(matchedCluster);
int pid = particle.getParticleIDUsed().getPDG();
if (Math.abs(pid) == 11) {
@@ -645,7 +656,7 @@
// Beam size variables.
// The beamsize array is in the tracking frame
/* TODO mg-May 14, 2014: the the beam size from the conditions db...also beam position! */
- protected double[] beamSize = {0.001, 0.130, 0.050}; //rough estimate from harp scans during engineering run production running
+ protected double[] beamSize = {0.001, 0.2, 0.02};
protected double bField;
// flipSign is a kludge...
Modified: java/branches/nathan-dev/recon/src/main/java/org/hps/recon/utils/TrackClusterMatcher.java
=============================================================================
--- java/branches/nathan-dev/recon/src/main/java/org/hps/recon/utils/TrackClusterMatcher.java (original)
+++ java/branches/nathan-dev/recon/src/main/java/org/hps/recon/utils/TrackClusterMatcher.java Sun Dec 6 13:05:52 2015
@@ -64,6 +64,104 @@
private double bottomClusterTrackMatchDeltaYLow = -28; // mm
private double bottomClusterTrackMatchDeltaYHigh = 24; // mm
+
+ /**
+ * Rafo's parameterization of cluster-seed x/y position residuals as function of energy.
+ * Derived using GBL tracks, non-analytic extrapolation, and uncorrected cluster positions.
+ */
+ private static final double dxMeanTopPosiGBL[] = { 6.67414,-9.57296, 5.70647, 27.4523,-28.1103,-9.11424 };
+ private static final double dxSigmTopPosiGBL[] = { 52.6437,-478.805, 1896.73,-3761.48, 3676.77,-1408.31 };
+ private static final double dxMeanBotPosiGBL[] = { 4.13802, 15.8887,-74.2844,-9.78944, 308.541,-287.668 };
+ private static final double dxSigmBotPosiGBL[] = { 37.6513,-294.851, 1002.15,-1639.08, 1228.02,-308.754 };
+ private static final double dxMeanTopElecGBL[] = {-1.6473, 5.58701, 25.3977,-17.1523,-121.025, 145.969 };
+ private static final double dxSigmTopElecGBL[] = { 48.7018,-423.599, 1592.66,-2959.99, 2668.97,-919.876 };
+ private static final double dxMeanBotElecGBL[] = {-6.63558, 83.7763,-460.451, 1275.63,-1702.83, 873.913 };
+ private static final double dxSigmBotElecGBL[] = { 47.0029,-411.784, 1586.52,-3083.37, 2985.58,-1145.53 };
+ private static final double dyMeanTopPosiGBL[] = { 0.71245, 5.57585,-6.50267,-8.21688, 39.8607,-43.9661 };
+ private static final double dySigmTopPosiGBL[] = { 33.0213,-275.174, 1168.77,-2642.34, 3045.52,-1406.21 };
+ private static final double dyMeanBotPosiGBL[] = {-5.532, 74.9738,-383.972, 977.849,-1250.28, 637.75 };
+ private static final double dySigmBotPosiGBL[] = { 19.019, -83.9253, 133.813, 119.883,-546.951, 405.207 };
+ private static final double dyMeanTopElecGBL[] = { 2.88498,-20.4101, 62.9689, 25.6386,-259.957, 207.145 };
+ private static final double dySigmTopElecGBL[] = { 8.65583, 120.676,-1166.43, 3811.72,-5383.19, 2787.42 };
+ private static final double dyMeanBotElecGBL[] = {-9.02276, 112.329,-489.761, 953.037,-829.96, 260.772 };
+ private static final double dySigmBotElecGBL[] = { 23.4856,-108.19, 158.7, 189.261,-682.034, 459.15 };
+
+ public double getNSigmaFromMeanGBL(Cluster cluster,Track track) {
+
+ if (this.useAnalyticExtrapolator)
+ throw new RuntimeException("This is to be used with non-analytic extrapolator only.");
+
+ if (track.getType() < 32)
+ throw new RuntimeException("This is to be used with GBL tracks only.");
+
+ // Get the cluster position:
+ Hep3Vector cPos = new BasicHep3Vector(cluster.getPosition());
+
+ // Get the extrapolated track position at the calorimeter:
+ TrackState trackStateAtEcal = TrackUtils.getTrackStateAtECal(track);
+ Hep3Vector tPos = new BasicHep3Vector(trackStateAtEcal.getReferencePoint());
+ tPos = CoordinateTransformations.transformVectorToDetector(tPos);
+
+ // whether track is in top half of detector:
+ final boolean isTopTrack = track.getTrackStates().get(0).getTanLambda() > 0;
+
+ // track and cluster in different halves:
+ if (isTopTrack != cPos.y()>0) return Double.MAX_VALUE;
+
+ // choose which parameterization of mean and sigma to use:
+ double dxMean[],dyMean[],dxSigm[],dySigm[];
+ if (track.getCharge()==1) {
+ if (isTopTrack) {
+ dxMean=dxMeanTopPosiGBL;
+ dxSigm=dxSigmTopPosiGBL;
+ dyMean=dyMeanTopPosiGBL;
+ dySigm=dySigmTopPosiGBL;
+ }
+ else {
+ dxMean=dxMeanBotPosiGBL;
+ dxSigm=dxSigmBotPosiGBL;
+ dyMean=dyMeanBotPosiGBL;
+ dySigm=dySigmBotPosiGBL;
+ }
+ }
+ else if (track.getCharge()==-1) {
+ if (isTopTrack) {
+ dxMean=dxMeanTopElecGBL;
+ dxSigm=dxSigmTopElecGBL;
+ dyMean=dyMeanTopElecGBL;
+ dySigm=dySigmTopElecGBL;
+ }
+ else {
+ dxMean=dxMeanBotElecGBL;
+ dxSigm=dxSigmBotElecGBL;
+ dyMean=dyMeanBotElecGBL;
+ dySigm=dySigmBotElecGBL;
+ }
+ }
+ else {
+ return Double.MAX_VALUE;
+ }
+
+ // get particle energy:
+ final double pp[] = track.getMomentum();
+ double ee = Math.sqrt( pp[0]*pp[0] + pp[1]*pp[1] + pp[2]*pp[2] );
+
+ // force the parameterizations to be constant above a certain energy:
+ if (ee > 0.6) ee=0.6;
+
+ // calculate measured mean and sigma of deltaX and deltaY for this energy:
+ double aDxMean=0,aDxSigm=0,aDyMean=0,aDySigm=0;
+ for (int ii=dxMean.length-1; ii>=0; ii--) aDxMean = dxMean[ii] + ee*aDxMean;
+ for (int ii=dxSigm.length-1; ii>=0; ii--) aDxSigm = dxSigm[ii] + ee*aDxSigm;
+ for (int ii=dyMean.length-1; ii>=0; ii--) aDyMean = dyMean[ii] + ee*aDyMean;
+ for (int ii=dySigm.length-1; ii>=0; ii--) aDySigm = dySigm[ii] + ee*aDySigm;
+
+ // calculate nSigma between track and cluster:
+ final double nSigmaX = Math.abs(cPos.x() - tPos.x() - aDxMean) / aDxSigm;
+ final double nSigmaY = Math.abs(cPos.y() - tPos.y() - aDyMean) / aDySigm;
+ return Math.sqrt(nSigmaX*nSigmaX + nSigmaY*nSigmaY);
+ }
+
/**
* Z position to start extrapolation from
*/
@@ -194,7 +292,7 @@
tPos = new BasicHep3Vector(trackStateAtEcal.getReferencePoint());
tPos = CoordinateTransformations.transformVectorToDetector(tPos);
}
-
+
return Math.sqrt(Math.pow(cPos.x()-tPos.x(),2)+Math.pow(cPos.y()-tPos.y(),2));
}
