WireGeo.h

Go to the documentation of this file.

#ifndef LARCOREALG_GEOMETRY_WIREGEO_H
#define LARCOREALG_GEOMETRY_WIREGEO_H

// LArSoft
#include "larcorealg/Geometry/LineClosestPoint.h"
#include "larcorealg/Geometry/LocalTransformationGeo.h"
#include "larcorealg/Geometry/TransformationMatrix.h"
#include "larcorealg/Geometry/geo_vectors_utils.h" // geo::vect

// ROOT
#include "Math/GenVector/Transform3D.h"
#include "TVector3.h"

// C/C++ libraries
#include <cmath> // std::sin(), ...
#include <string>
#include <type_traits> // std::is_nothrow_move_constructible<>
#include <vector>

class TGeoNode;

namespace geo {

  class WireGeo;
  struct WireID;

  Point_t WiresIntersection(WireGeo const& wireA, WireGeo const& wireB);

  IntersectionPointAndOffsets<Point_t> WiresIntersectionAndOffsets(WireGeo const& wireA,
                                                                   WireGeo const& wireB);

  class WireGeo {
  public:
    using ID_t = WireID;
    using GeoNodePath_t = std::vector<TGeoNode const*>;

    // -- BEGIN -- Types for geometry-local reference vectors ------------------

    struct WireGeoCoordinatesTag {};

    using LocalPoint_t = Point3DBase_t<WireGeoCoordinatesTag>;

    using LocalVector_t = Vector3DBase_t<WireGeoCoordinatesTag>;

    // -- END ---- Types for geometry-local reference vectors ------------------

    WireGeo(TGeoNode const* node, TransformationMatrix&& trans);

    // -- BEGIN -- Size and coordinates ----------------------------------------

    double RMax() const;

    double HalfL() const { return fHalfL; }

    double RMin() const;


    Point_t GetPositionFromCenter(double localz) const
    {
      return GetPositionFromCenterUnbounded(capLength(localz));
    }


    Point_t GetPositionFromCenterUnbounded(double localz) const
    {
      return toWorldCoords(LocalPoint_t{0.0, 0.0, relLength(localz)});
    }

    Point_t const& GetCenter() const { return fCenter; }

    Point_t GetStart() const { return GetPositionFromCenterUnbounded(-HalfL()); }

    Point_t GetEnd() const { return GetPositionFromCenterUnbounded(+HalfL()); }

    double Length() const { return 2. * HalfL(); }

    // -- END ---- Size and coordinates  ---------------------------------------

    // -- BEGIN -- Orientation and angles --------------------------------------

    double ThetaZ() const { return fThetaZ; }


    double ThetaZ(bool degrees) const;

    double CosThetaZ() const { return std::cos(ThetaZ()); }
    double SinThetaZ() const { return std::sin(ThetaZ()); }
    double TanThetaZ() const { return std::tan(ThetaZ()); }

    bool isHorizontal() const { return std::abs(SinThetaZ()) < 1e-5; }

    bool isVertical() const { return std::abs(CosThetaZ()) < 1e-5; }

    bool isParallelTo(WireGeo const& wire) const
    {
      return // parallel if the dot product of the directions is about +/- 1
        std::abs(std::abs(Direction().Dot(wire.Direction())) - 1.) < 1e-5;
    }

    Vector_t Direction() const
    {
      // maybe (GetCenter() - GetStart()) / HalfL() would be faster;
      return (GetEnd() - GetStart()) / Length();
    }

    // -- END ---- Orientation and angles --------------------------------------

    // -- BEGIN -- Printing ----------------------------------------------------

    template <typename Stream>
    void PrintWireInfo(Stream&& out, std::string indent = "", unsigned int verbosity = 1) const;

    std::string WireInfo(std::string indent = "", unsigned int verbosity = 1) const;

    static constexpr unsigned int MaxVerbosity = 4;

    // -- END ---- Printing ----------------------------------------------------

    // -- BEGIN -- Coordinate transformation -----------------------------------

    Point_t toWorldCoords(LocalPoint_t const& local) const { return fTrans.toWorldCoords(local); }

    Vector_t toWorldCoords(LocalVector_t const& local) const { return fTrans.toWorldCoords(local); }

    LocalPoint_t toLocalCoords(Point_t const& world) const { return fTrans.toLocalCoords(world); }

    LocalVector_t toLocalCoords(Vector_t const& world) const { return fTrans.toLocalCoords(world); }

    // -- END ---- Coordinate transformation -----------------------------------

    const TGeoNode* Node() const { return fWireNode; }

    // -- BEGIN -- Geometric properties and algorithms -------------------------

    double ComputeZatY0() const { return fCenter.Z() - fCenter.Y() / TanThetaZ(); }

    double DistanceFrom(WireGeo const& wire) const;

    Point_t IntersectionWith(WireGeo const& other) const { return WiresIntersection(*this, other); }

    IntersectionPointAndOffsets<Point_t> IntersectionAndOffsetsWith(WireGeo const& other) const
    {
      auto const [point, ofsA, ofsB] = WiresIntersectionAndOffsets(*this, other);
      return {point, ofsA, ofsB};
    }

    // -- END ---- Geometric properties and algorithms -------------------------

    void UpdateAfterSorting(WireID const&, bool flip);

    static double WirePitch(WireGeo const& w1, WireGeo const& w2)
    {
      return std::abs(w2.DistanceFrom(w1));
    }

  private:
    using LocalTransformation_t =
      LocalTransformationGeo<ROOT::Math::Transform3D, LocalPoint_t, LocalVector_t>;

    const TGeoNode* fWireNode;    
    double fThetaZ;               
    double fHalfL;                
    Point_t fCenter;              
    LocalTransformation_t fTrans; 
    bool flipped;                 

    double relLength(double local) const { return flipped ? -local : local; }

    double capLength(double local) const { return std::min(+HalfL(), std::max(-HalfL(), local)); }

    void Flip();

  }; // class WireGeo

  static_assert(std::is_move_assignable_v<WireGeo>);
  static_assert(std::is_move_constructible_v<WireGeo>);

} // namespace geo

//------------------------------------------------------------------------------
//--- template implementation
//---
//------------------------------------------------------------------------------
template <typename Stream>
void geo::WireGeo::PrintWireInfo(Stream&& out,
                                 std::string indent /* = "" */,
                                 unsigned int verbosity /* = 1 */
) const
{
  //----------------------------------------------------------------------------
  out << "wire from " << GetStart() << " to " << GetEnd();

  if (verbosity-- <= 0) return; // 0

  //----------------------------------------------------------------------------
  out << " (" << Length() << " cm long)";

  if (verbosity-- <= 0) return; // 1

  //----------------------------------------------------------------------------
  out << ", theta(z)=" << ThetaZ() << " rad";

  if (verbosity-- <= 0) return; // 2

  //----------------------------------------------------------------------------
  out << "\n" << indent << "  center at " << GetCenter() << " cm";

  if (verbosity-- <= 0) return; // 3

  //----------------------------------------------------------------------------
  out << ", direction: " << Direction();
  if (isHorizontal()) out << " (horizontal)";
  if (isVertical()) out << " (vertical)";
}

//------------------------------------------------------------------------------
inline geo::IntersectionPointAndOffsets<geo::Point_t> geo::WiresIntersectionAndOffsets(
  WireGeo const& wireA,
  WireGeo const& wireB)
{

  return LineClosestPointAndOffsetsWithUnitVectors(
    wireA.GetCenter(), wireA.Direction(), wireB.GetCenter(), wireB.Direction());
}

//------------------------------------------------------------------------------
inline geo::Point_t geo::WiresIntersection(WireGeo const& wireA, WireGeo const& wireB)
{
  return LineClosestPointWithUnitVectors(
    wireA.GetCenter(), wireA.Direction(), wireB.GetCenter(), wireB.Direction());
}

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

#endif // LARCOREALG_GEOMETRY_WIREGEO_H