RecoHGCal/TICL: optimize calculation for comparison in HGCDoublet::ar…

…eAligned By squaring cosTheta and multiplying by the squares of the magnitudes an equivalent comparison is made without the division and square root which are costly FP ops.
cms-sw · Jul 17, 2020 · a9840a9 · a9840a9
1 parent 92997c6
commit a9840a9
Showing 1 changed file with 15 additions and 13 deletions.
diff --git a/RecoHGCal/TICL/plugins/HGCDoublet.cc b/RecoHGCal/TICL/plugins/HGCDoublet.cc
@@ -82,13 +82,13 @@ int HGCDoublet::areAligned(double xi,
   // inner product
   auto dot = dx1 * dx2 + dy1 * dy2 + dz1 * dz2;
   // magnitudes
-  auto mag1 = std::sqrt(dx1 * dx1 + dy1 * dy1 + dz1 * dz1);
-  auto mag2 = std::sqrt(dx2 * dx2 + dy2 * dy2 + dz2 * dz2);
-  // angle between the vectors
-  auto cosTheta = dot / (mag1 * mag2);
+  auto mag1sq = dx1 * dx1 + dy1 * dy1 + dz1 * dz1;
+  auto mag2sq = dx2 * dx2 + dy2 * dy2 + dz2 * dz2;
   if (debug) {
-    LogDebug("HGCDoublet") << "-- Are Aligned -- dot: " << dot << " mag1: " << mag1 << " mag2: " << mag2
-                           << " cosTheta: " << cosTheta << " isWithinLimits: " << (cosTheta > minCosTheta) << std::endl;
+    LogDebug("HGCDoublet") << "-- Are Aligned -- dotsq: " << dot * dot << " mag1sq: " << mag1sq << " mag2sq: " << mag2sq
+                           << "minCosTheta_sq:" << minCosTheta * minCosTheta
+                           << " isWithinLimits: " << (dot * dot > minCosTheta * minCosTheta * mag1sq * mag2sq)
+                           << std::endl;
   }
 
   // Now check the compatibility with the pointing origin.
@@ -102,15 +102,17 @@ int HGCDoublet::areAligned(double xi,
   const GlobalVector pointingDir = (seedIndex_ == -1) ? GlobalVector(innerX(), innerY(), innerZ()) : refDir;
 
   auto dot_pointing = pointingDir.dot(firstDoublet);
-  auto mag_pointing = sqrt(pointingDir.mag2());
-  auto cosTheta_pointing = dot_pointing / (mag2 * mag_pointing);
+  auto mag_pointing_sq = pointingDir.mag2();
   if (debug) {
-    LogDebug("HGCDoublet") << "-- Are Aligned -- dot_pointing: " << dot_pointing << " mag_pointing: " << mag_pointing
-                           << " mag2: " << mag2 << " cosTheta_pointing: " << cosTheta_pointing
-                           << " isWithinLimits: " << (cosTheta_pointing > minCosPointing) << std::endl;
+    LogDebug("HGCDoublet") << "-- Are Aligned -- dot_pointing_sq: " << dot_pointing * dot_pointing
+                           << " mag_pointing_sq: " << mag_pointing_sq << " mag2sq: " << mag2sq << " isWithinLimits: "
+                           << (dot_pointing * dot_pointing > minCosPointing * minCosPointing * mag_pointing_sq * mag2sq)
+                           << std::endl;
   }
-
-  return (cosTheta > minCosTheta) && (cosTheta_pointing > minCosPointing);
+  // by squaring cosTheta and multiplying by the squares of the magnitudes
+  // an equivalent comparison is made without the division and square root which are costly FP ops.
+  return (dot * dot > minCosTheta * minCosTheta * mag1sq * mag2sq) &&
+         (dot_pointing * dot_pointing > minCosPointing * minCosPointing * mag_pointing_sq * mag2sq);
 }
 
 void HGCDoublet::findNtuplets(std::vector<HGCDoublet> &allDoublets,