BoxBoundedGeo.cxx

Go to the documentation of this file.

#include "larcorealg/Geometry/BoxBoundedGeo.h"

// LArSoft libraries
#include "larcorealg/Geometry/geo_vectors_utils.h"

// ROOT libraries
#include "Math/GenVector/Cartesian3D.h"
#include "Math/GenVector/DisplacementVector3D.h"
#include "Math/GenVector/PositionVector3D.h"

#include <array>
#include <type_traits>
#include <vector>

namespace geo {
  //----------------------------------------------------------------------------
  bool BoxBoundedGeo::ContainsPosition(TVector3 const& point, double wiggle /* = 1.0 */) const
  {
    return ContainsPosition(geo::vect::toPoint(point), wiggle);
  }

  //----------------------------------------------------------------------------
  bool BoxBoundedGeo::ContainsPosition(double const* point, double wiggle /* = 1.0 */) const
  {
    return ContainsPosition(geo::vect::makePointFromCoords(point), wiggle);
  }

  //----------------------------------------------------------------------------
  std::vector<geo::Point_t> BoxBoundedGeo::GetIntersections(
    geo::Point_t const& TrajectoryStart,
    geo::Vector_t const& TrajectoryDirect) const
  {

    std::vector<geo::Point_t> IntersectionPoints;
    std::vector<double> LineParameters;

    // Generate normal vectors and offsets for every plane of the box
    // All normal vectors are headed outwards
    // BUG the double brace syntax is required to work around clang bug 21629
    // (https://bugs.llvm.org/show_bug.cgi?id=21629)
    static std::array<geo::Vector_t, 6U> const NormalVectors = {{
      -geo::Xaxis<geo::Vector_t>(),
      geo::Xaxis<geo::Vector_t>(), // anode, cathode,
      -geo::Yaxis<geo::Vector_t>(),
      geo::Yaxis<geo::Vector_t>(), // bottom, top,
      -geo::Zaxis<geo::Vector_t>(),
      geo::Zaxis<geo::Vector_t>() // upstream, downstream
    }};
    // BUG the double brace syntax is required to work around clang bug 21629
    // (https://bugs.llvm.org/show_bug.cgi?id=21629)
    std::array<geo::Point_t, 6U> const NormalVectorOffsets = {{
      geo::Point_t{Min().X(), Min().Y(), Min().Z()}, // Anode offset
      geo::Point_t{Max().X(), Min().Y(), Min().Z()}, // Cathode
      geo::Point_t{Min().X(), Min().Y(), Min().Z()}, // Bottom
      geo::Point_t{Min().X(), Max().Y(), Min().Z()}, // Top
      geo::Point_t{Min().X(), Min().Y(), Min().Z()}, // upstream
      geo::Point_t{Min().X(), Min().Y(), Max().Z()}  // downstream
    }};

    // Loop over all surfaces of the box
    for (unsigned int face_no = 0; face_no < NormalVectors.size(); face_no++) {
      // Check if trajectory and surface are not parallel
      if (NormalVectors[face_no].Dot(TrajectoryDirect)) {
        // Calculate the line parameter for the intersection points
        LineParameters.push_back(
          NormalVectors[face_no].Dot(NormalVectorOffsets.at(face_no) - TrajectoryStart) /
          NormalVectors[face_no].Dot(TrajectoryDirect));
      }
      else
        continue;

      // Calculate intersection point using the line parameter
      IntersectionPoints.push_back(TrajectoryStart + LineParameters.back() * TrajectoryDirect);

      // Coordinate which should be ignored when checking for limits added by Christoph Rudolf von Rohr 05/21/2016
      unsigned int NoCheckCoord;

      // Calculate NoCheckCoord out of the face_no
      if (face_no % 2) {
        // Convert odd face number to coordinate
        NoCheckCoord = (face_no - 1) / 2;
      }
      else {
        // Convert even face number to coordinate
        NoCheckCoord = face_no / 2;
      }

      // Loop over all three space coordinates
      unsigned int coord = 0;
      for (auto extractCoord : geo::vect::coordReaders<geo::Point_t>()) {
        auto const lastPointCoord = geo::vect::bindCoord(IntersectionPoints.back(), extractCoord);
        auto const minCoord = geo::vect::bindCoord(c_min, extractCoord);
        auto const maxCoord = geo::vect::bindCoord(c_max, extractCoord);

        // Changed by Christoph Rudolf von Rohr 05/21/2016
        // Then check if point is not within the surface limits at this coordinate, without looking
        // at the plane normal vector coordinate. We can assume, that our algorithm already found this coordinate correctily.
        // In rare cases, looking for boundaries in this coordinate causes unexpected behavior due to floating point inaccuracies.
        if (coord++ != NoCheckCoord &&
            ((lastPointCoord() > maxCoord()) || (lastPointCoord() < minCoord))) {
          // if off limits, get rid of the useless data and break the coordinate loop
          LineParameters.pop_back();
          IntersectionPoints.pop_back();
          break;
        }
      } // coordinate loop
    }   // Surcaces loop

    // sort points according to their parameter value (first is entry, second is exit)
    if (LineParameters.size() == 2 && LineParameters.front() > LineParameters.back()) {
      std::swap(IntersectionPoints.front(), IntersectionPoints.back());
    }

    return IntersectionPoints;
  } // GetIntersections()

  //----------------------------------------------------------------------------
  std::vector<TVector3> BoxBoundedGeo::GetIntersections(TVector3 const& TrajectoryStart,
                                                        TVector3 const& TrajectoryDirect) const
  {
    std::vector<TVector3> intersections;
    for (auto const& point :
         GetIntersections(geo::Point_t(TrajectoryStart), geo::Vector_t(TrajectoryDirect)))
      intersections.emplace_back(point.X(), point.Y(), point.Z());
    return intersections;
  } // GetIntersections(TVector3)

  //----------------------------------------------------------------------------
  void BoxBoundedGeo::SortCoordinates()
  {

    for (auto coordMan : geo::vect::coordManagers<geo::Point_t>()) {
      auto min = geo::vect::bindCoord(c_min, coordMan);
      auto max = geo::vect::bindCoord(c_max, coordMan);
      if (min() > max()) {
        auto temp = min();
        min = max();
        max = temp;
      }
    } // for

  } // BoxBoundedGeo::SortCoordinates()

  //----------------------------------------------------------------------------

} // namespace geo