DoubleHitSpacePointBuilder.ipp

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file DoubleHitSpacePointBuilder.ipp

// This file is part of the Acts project.
//
// Copyright (C) 2018 CERN for the benefit of the Acts project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.

#include <cmath>
#include <limits>
#include "Acts/Utilities/Helpers.hpp"

namespace Acts {
namespace detail {
struct SpacePointParameters {
  Vector3D q;
  Vector3D r;
  Vector3D s;
  Vector3D t;
  Vector3D qs;
  Vector3D rt;
  double qmag = 0.;
  double m = 0.;
  double n = 0.;
  double limit = 1.;
  double limitExtended = 0.;
};

double differenceOfClustersChecked(const Vector3D& pos1, const Vector3D& pos2,
                                   const Vector3D& posVertex,
                                   const double maxDistance,
                                   const double maxAngleTheta2,
                                   const double maxAnglePhi2) {
  // Check if measurements are close enough to each other
  if ((pos1 - pos2).norm() > maxDistance) {
    return -1.;
  }

  // Calculate the angles of the vectors
  double phi1, theta1, phi2, theta2;
  phi1 = VectorHelpers::phi(pos1 - posVertex);
  theta1 = VectorHelpers::theta(pos1 - posVertex);
  phi2 = VectorHelpers::phi(pos2 - posVertex);
  theta2 = VectorHelpers::theta(pos2 - posVertex);

  // Calculate the squared difference between the theta angles
  double diffTheta2 = (theta1 - theta2) * (theta1 - theta2);
  if (diffTheta2 > maxAngleTheta2) {
    return -1.;
  }

  // Calculate the squared difference between the phi angles
  double diffPhi2 = (phi1 - phi2) * (phi1 - phi2);
  if (diffPhi2 > maxAnglePhi2) {
    return -1.;
  }

  // Return the squared distance between both vector
  return diffTheta2 + diffPhi2;
}

std::pair<Vector2D, Vector2D> findLocalTopAndBottomEnd(
    const Vector2D& local, const CartesianSegmentation* segment) {
  auto& binData = segment->binUtility().binningData();
  auto& boundariesX = binData[0].boundaries();
  auto& boundariesY = binData[1].boundaries();

  // Search the x-/y-bin of the cluster
  size_t binX = binData[0].searchLocal(local);
  size_t binY = binData[1].searchLocal(local);

  // Storage of the local top (first) and bottom (second) end
  std::pair<Vector2D, Vector2D> topBottomLocal;

  if (boundariesX[binX + 1] - boundariesX[binX] <
      boundariesY[binY + 1] - boundariesY[binY]) {
    // Set the top and bottom end of the strip in local coordinates
    topBottomLocal.first = {(boundariesX[binX] + boundariesX[binX + 1]) / 2,
                            boundariesY[binY + 1]};
    topBottomLocal.second = {(boundariesX[binX] + boundariesX[binX + 1]) / 2,
                             boundariesY[binY]};
  } else {
    // Set the top and bottom end of the strip in local coordinates
    topBottomLocal.first = {boundariesX[binX],
                            (boundariesY[binY] + boundariesY[binY + 1]) / 2};
    topBottomLocal.second = {boundariesX[binX + 1],
                             (boundariesY[binY] + boundariesY[binY + 1]) / 2};
  }
  return topBottomLocal;
}

double calcPerpendicularProjection(const Vector3D& a, const Vector3D& c,
                                   const Vector3D& q, const Vector3D& r) {
  Vector3D ac = c - a;
  double qr = q.dot(r);
  double denom = q.dot(q) - qr * qr;

  // Check for numerical stability
  if (fabs(denom) > 10e-7) {
    // Return lambda0
    return (ac.dot(r) * qr - ac.dot(q) * r.dot(r)) / denom;
  }
  // lambda0 is in the interval [-1,0]. This return serves as error check.
  return 1.;
}

bool recoverSpacePoint(SpacePointParameters& spaPoPa,
                       double stripLengthGapTolerance) {
  // Check if the limits are allowed to be increased
  if (stripLengthGapTolerance <= 0.) {
    return false;
  }
  spaPoPa.qmag = spaPoPa.q.norm();
  // Increase the limits. This allows a check if the point is just slightly
  // outside the SDE
  spaPoPa.limitExtended =
      spaPoPa.limit + stripLengthGapTolerance / spaPoPa.qmag;
  // Check if m is just slightly outside
  if (fabs(spaPoPa.m) > spaPoPa.limitExtended) {
    return false;
  }
  // Calculate n if not performed previously
  if (spaPoPa.n == 0.) {
    spaPoPa.n = -spaPoPa.t.dot(spaPoPa.qs) / spaPoPa.r.dot(spaPoPa.qs);
  }
  // Check if n is just slightly outside
  if (fabs(spaPoPa.n) > spaPoPa.limitExtended) {
    return false;
  }


  // Calculate the scaling factor to project lengths of the second SDE on the
  // first SDE
  double secOnFirstScale =
      spaPoPa.q.dot(spaPoPa.r) / (spaPoPa.qmag * spaPoPa.qmag);
  // Check if both overshoots are in the same direction
  if (spaPoPa.m > 1. && spaPoPa.n > 1.) {
    // Calculate the overshoots
    double mOvershoot = spaPoPa.m - 1.;
    double nOvershoot =
        (spaPoPa.n - 1.) * secOnFirstScale;  // Perform projection
    // Resolve worse overshoot
    double biggerOvershoot = std::max(mOvershoot, nOvershoot);
    // Move m and n towards 0
    spaPoPa.m -= biggerOvershoot;
    spaPoPa.n -= (biggerOvershoot / secOnFirstScale);
    // Check if this recovered the space point
    return fabs(spaPoPa.m) < spaPoPa.limit && fabs(spaPoPa.n) < spaPoPa.limit;
  }
  // Check if both overshoots are in the same direction
  if (spaPoPa.m < -1. && spaPoPa.n < -1.) {
    // Calculate the overshoots
    double mOvershoot = -(spaPoPa.m + 1.);
    double nOvershoot =
        -(spaPoPa.n + 1.) * secOnFirstScale;  // Perform projection
    // Resolve worse overshoot
    double biggerOvershoot = std::max(mOvershoot, nOvershoot);
    // Move m and n towards 0
    spaPoPa.m += biggerOvershoot;
    spaPoPa.n += (biggerOvershoot / secOnFirstScale);
    // Check if this recovered the space point
    return fabs(spaPoPa.m) < spaPoPa.limit && fabs(spaPoPa.n) < spaPoPa.limit;
  }
  // No solution could be found
  return false;
}

bool calculateSpacePoint(const std::pair<Vector3D, Vector3D>& stripEnds1,
                         const std::pair<Vector3D, Vector3D>& stripEnds2,
                         const Vector3D& posVertex,
                         SpacePointParameters& spaPoPa,
                         const double stripLengthTolerance) {

  spaPoPa.s = stripEnds1.first + stripEnds1.second - 2 * posVertex;
  spaPoPa.t = stripEnds2.first + stripEnds2.second - 2 * posVertex;
  spaPoPa.qs = spaPoPa.q.cross(spaPoPa.s);
  spaPoPa.rt = spaPoPa.r.cross(spaPoPa.t);
  spaPoPa.m = -spaPoPa.s.dot(spaPoPa.rt) / spaPoPa.q.dot(spaPoPa.rt);

  // Set the limit for the parameter
  if (spaPoPa.limit == 1. && stripLengthTolerance != 0.) {
    spaPoPa.limit = 1. + stripLengthTolerance;
  }

  // Check if m and n can be resolved in the interval (-1, 1)
  return (fabs(spaPoPa.m) <= spaPoPa.limit &&
          fabs(spaPoPa.n = -spaPoPa.t.dot(spaPoPa.qs) /
                           spaPoPa.r.dot(spaPoPa.qs)) <= spaPoPa.limit);
}
}  // namespace detail
}  // namespace Acts

template <typename Cluster>
Acts::SpacePointBuilder<Acts::SpacePoint<Cluster>>::SpacePointBuilder(
    DoubleHitSpacePointConfig cfg)
    : m_cfg(std::move(cfg)) {}

template <typename Cluster>
Acts::Vector2D Acts::SpacePointBuilder<Acts::SpacePoint<Cluster>>::localCoords(
    const Cluster& cluster) const {
  // Local position information
  auto par = cluster.parameters();
  Acts::Vector2D local(par[Acts::ParDef::eLOC_0], par[Acts::ParDef::eLOC_1]);
  return local;
}

template <typename Cluster>
Acts::Vector3D Acts::SpacePointBuilder<Acts::SpacePoint<Cluster>>::globalCoords(
    const GeometryContext& gctx, const Cluster& cluster) const {
  // Receive corresponding surface
  auto& clusterSurface = cluster.referenceSurface();

  // Transform local into global position information
  Acts::Vector3D pos, mom;
  clusterSurface.localToGlobal(gctx, localCoords(cluster), mom, pos);

  return pos;
}

template <typename Cluster>
void Acts::SpacePointBuilder<Acts::SpacePoint<Cluster>>::makeClusterPairs(
    const GeometryContext& gctx,
    const std::vector<const Cluster*>& clustersFront,
    const std::vector<const Cluster*>& clustersBack,
    std::vector<std::pair<const Cluster*, const Cluster*>>& clusterPairs)
    const {
  // Return if no clusters are given in a vector
  if (clustersFront.empty() || clustersBack.empty()) {
    return;
  }

  // Declare helper variables
  double currentDiff;
  double diffMin;
  unsigned int clusterMinDist;

  // Walk through all clusters on both surfaces
  for (unsigned int iClustersFront = 0; iClustersFront < clustersFront.size();
       iClustersFront++) {
    // Set the closest distance to the maximum of double
    diffMin = std::numeric_limits<double>::max();
    // Set the corresponding index to an element not in the list of clusters
    clusterMinDist = clustersBack.size();
    for (unsigned int iClustersBack = 0; iClustersBack < clustersBack.size();
         iClustersBack++) {
      // Calculate the distances between the hits
      currentDiff = detail::differenceOfClustersChecked(
          globalCoords(gctx, *(clustersFront[iClustersFront])),
          globalCoords(gctx, *(clustersBack[iClustersBack])), m_cfg.vertex,
          m_cfg.diffDist, m_cfg.diffPhi2, m_cfg.diffTheta2);
      // Store the closest clusters (distance and index) calculated so far
      if (currentDiff < diffMin && currentDiff >= 0.) {
        diffMin = currentDiff;
        clusterMinDist = iClustersBack;
      }
    }

    // Store the best (=closest) result
    if (clusterMinDist < clustersBack.size()) {
      std::pair<const Cluster*, const Cluster*> clusterPair;
      clusterPair = std::make_pair(clustersFront[iClustersFront],
                                   clustersBack[clusterMinDist]);
      clusterPairs.push_back(clusterPair);
    }
  }
}

template <typename Cluster>
std::pair<Acts::Vector3D, Acts::Vector3D>
Acts::SpacePointBuilder<Acts::SpacePoint<Cluster>>::endsOfStrip(
    const GeometryContext& gctx, const Cluster& cluster) const {
  // Calculate the local coordinates of the cluster
  const Acts::Vector2D local = localCoords(cluster);

  // Receive the binning
  auto segment = dynamic_cast<const Acts::CartesianSegmentation*>(
      &(cluster.digitizationModule()->segmentation()));

  std::pair<Vector2D, Vector2D> topBottomLocal =
      detail::findLocalTopAndBottomEnd(local, segment);

  // Calculate the global coordinates of the top and bottom end of the strip
  Acts::Vector3D topGlobal, bottomGlobal, mom;  // mom is a dummy variable
  const auto* sur = &cluster.referenceSurface();
  sur->localToGlobal(gctx, topBottomLocal.first, mom, topGlobal);
  sur->localToGlobal(gctx, topBottomLocal.second, mom, bottomGlobal);

  // Return the top and bottom end of the strip in global coordinates
  return std::make_pair(topGlobal, bottomGlobal);
}

template <typename Cluster>
void Acts::SpacePointBuilder<Acts::SpacePoint<Cluster>>::calculateSpacePoints(
    const GeometryContext& gctx,
    const std::vector<std::pair<const Cluster*, const Cluster*>>& clusterPairs,
    std::vector<Acts::SpacePoint<Cluster>>& spacePoints) const {

  detail::SpacePointParameters spaPoPa;

  // Walk over every found candidate pair
  for (const auto& cp : clusterPairs) {
    // Calculate the ends of the SDEs
    const auto& ends1 = endsOfStrip(gctx, *(cp.first));
    const auto& ends2 = endsOfStrip(gctx, *(cp.second));

    spaPoPa.q = ends1.first - ends1.second;
    spaPoPa.r = ends2.first - ends2.second;

    // Fast skipping if a perpendicular projection should be used
    double resultPerpProj;
    if (m_cfg.usePerpProj) {
      resultPerpProj = detail::calcPerpendicularProjection(
          ends1.first, ends2.first, spaPoPa.q, spaPoPa.r);
      if (resultPerpProj <= 0.) {
        Acts::SpacePoint<Cluster> sp;
        sp.clusterModule.push_back(cp.first);
        sp.clusterModule.push_back(cp.second);
        Vector3D pos = ends1.first + resultPerpProj * spaPoPa.q;
        sp.vector = pos;
        spacePoints.push_back(std::move(sp));
        continue;
      }
    }

    if (calculateSpacePoint(ends1, ends2, m_cfg.vertex, spaPoPa,
                            m_cfg.stripLengthTolerance)) {
      // Store the space point
      Acts::SpacePoint<Cluster> sp;
      sp.clusterModule.push_back(cp.first);
      sp.clusterModule.push_back(cp.second);
      Vector3D pos = 0.5 * (ends1.first + ends1.second + spaPoPa.m * spaPoPa.q);
      // TODO: Clusters should deliver timestamp
      sp.vector = pos;
      spacePoints.push_back(std::move(sp));
    } else {
      // Check if a recovery the point(s) and store them if successful
      if (detail::recoverSpacePoint(spaPoPa, m_cfg.stripLengthGapTolerance)) {
        Acts::SpacePoint<Cluster> sp;
        sp.clusterModule.push_back(cp.first);
        sp.clusterModule.push_back(cp.second);
        Vector3D pos =
            0.5 * (ends1.first + ends1.second + spaPoPa.m * spaPoPa.q);
        sp.vector = pos;
        spacePoints.push_back(std::move(sp));
      }
    }
  }
}