PlaneGeo.h

Go to the documentation of this file.

#ifndef LARCOREALG_GEOMETRY_PLANEGEO_H
#define LARCOREALG_GEOMETRY_PLANEGEO_H

// LArSoft libraries
#include "larcorealg/CoreUtils/DereferenceIterator.h"
#include "larcorealg/Geometry/GeoObjectSorter.h"
#include "larcorealg/Geometry/SimpleGeo.h"
#include "larcorealg/Geometry/LocalTransformationGeo.h"
#include "larcorealg/Geometry/BoxBoundedGeo.h"
#include "larcorealg/Geometry/Decomposer.h"
#include "larcorealg/Geometry/WireGeo.h"
#include "larcorealg/Geometry/geo_vectors_utils.h" // geo::vect
#include "larcoreobj/SimpleTypesAndConstants/geo_types.h"

// ROOT libraries
#include "Math/GenVector/Cartesian2D.h"
#include "Math/GenVector/PositionVector2D.h"
#include "Math/GenVector/DisplacementVector2D.h"

// ROOT libraries
#include "TGeoMatrix.h" // TGeoHMatrix

// C/C++ standard libraries
#include <cmath> // std::atan2()
#include <vector>
#include <string>


class TGeoNode;
class TVector3;

namespace geo {

  namespace details {
    struct ActiveAreaCalculator;
  } // namespace details
  
  //......................................................................
  
  // Note: SignalType() and SetSignalType() have been removed.
  //       Use `geo::GeometryCore::SignalType` instead.
  //       (see LArSoft issue #14365 at https://cdcvs.fnal.gov/redmine/issues/14365 )
  class PlaneGeo {
    
    using DefaultVector_t = TVector3; // ... not for long
    using DefaultPoint_t = TVector3; // ... not for long
    
  public:
    
    using GeoNodePath_t = std::vector<TGeoNode const*>;
    

    struct PlaneGeoCoordinatesTag {};
    
    using LocalPoint_t = geo::Point3DBase_t<PlaneGeoCoordinatesTag>;
    
    using LocalVector_t = geo::Vector3DBase_t<PlaneGeoCoordinatesTag>;
    
    
    
    struct WireCoordinateReferenceTag {};
    
    using WireCoordProjection_t = ROOT::Math::DisplacementVector2D
      <ROOT::Math::Cartesian2D<double>, WireCoordinateReferenceTag>;
    
    using WireDecomposer_t = geo::Decomposer
      <geo::Vector_t, geo::Point_t, WireCoordProjection_t>;
    
    using WireDecomposedVector_t = WireDecomposer_t::DecomposedVector_t;
    
    
    
    
    struct WidthDepthReferenceTag {};
    
    using WidthDepthProjection_t = ROOT::Math::DisplacementVector2D
      <ROOT::Math::Cartesian2D<double>, WidthDepthReferenceTag>;
    
    using WidthDepthDisplacement_t = ROOT::Math::DisplacementVector2D
      <ROOT::Math::Cartesian2D<double>, WidthDepthReferenceTag>;
    
    using WidthDepthDecomposer_t = geo::Decomposer
      <geo::Vector_t, geo::Point_t, WidthDepthProjection_t>;
    
    using WDDecomposedVector_t = WidthDepthDecomposer_t::DecomposedVector_t;
    
    
    
    using Rect = lar::util::simple_geo::Rectangle<double>;
    
    
    PlaneGeo(GeoNodePath_t& path, size_t depth);
    
    
    
    View_t View()                                             const { return fView;          }
    
    Orient_t Orientation()                                    const { return fOrientation;   }

    double ThetaZ()                                           const;
    
    double PhiZ()                                             const
      { return std::atan2(fSinPhiZ, fCosPhiZ); }
    
    double SinPhiZ()                                          const { return fSinPhiZ; }
    
    double CosPhiZ()                                          const { return fCosPhiZ; }
    
    geo::PlaneID const& ID() const { return fID; }
    
    
    
    template <typename Vector = DefaultVector_t>
    Vector WidthDir() const { return geo::vect::convertTo<Vector>(fDecompFrame.MainDir()); }
    
    template <typename Vector = DefaultVector_t>
    Vector DepthDir() const { return geo::vect::convertTo<Vector>(fDecompFrame.SecondaryDir()); }
    
    double Width() const { return fFrameSize.Width(); }
    
    double Depth() const { return fFrameSize.Depth(); }

    
    geo::BoxBoundedGeo BoundingBox() const;
    
    
    
    
    unsigned int Nwires()                                     const { return fWire.size();   }
    unsigned int NElements()                                  const { return Nwires();       }
    

    bool HasWire(unsigned int iwire) const { return iwire < Nwires(); }
    bool HasElement(unsigned int iwire) const { return HasWire(iwire); }
    

    bool HasWire(geo::WireID const& wireid) const
      { return HasWire(wireid.Wire); }
    bool HasElement(geo::WireID const& wireid) const
      { return HasWire(wireid); }
    
    WireGeo const& Wire(unsigned int iwire) const;
    

    WireGeo const& Wire(WireID const& wireid) const
      { return Wire(wireid.Wire); }
    WireGeo const& GetElement(WireID const& wireid) const
      { return Wire(wireid); }
    
    WireGeo const* WirePtr(unsigned int iwire) const
      { return HasWire(iwire)? &(fWire[iwire]): nullptr; }
    

    WireGeo const* WirePtr(WireID const& wireid) const
      { return WirePtr(wireid.Wire); }
    WireGeo const* GetElementPtr(WireID const& wireid) const
      { return WirePtr(wireid); }
    
    
    const WireGeo& FirstWire()                                const { return Wire(0);        }
    
    const WireGeo& MiddleWire()                               const { return Wire(Nwires()/2); }
    
    const WireGeo& LastWire()                                 const { return Wire(Nwires()-1); }
    
    auto IterateWires() const
      { return fWire; }
    
    
    
    
    double WirePitch() const { return fWirePitch; }
    
    bool WireIDincreasesWithZ() const;
    
    template <typename Vector = DefaultVector_t>
    Vector GetNormalDirection() const { return geo::vect::convertTo<Vector>(fNormal); }

    template <typename Vector = DefaultVector_t>
    Vector GetIncreasingWireDirection() const
      { return geo::vect::convertTo<Vector>(fDecompWire.SecondaryDir()); }
    
    
    template <typename Point = DefaultPoint_t>
    Point GetCenter() const
      { return geo::vect::convertTo<Point>(fCenter); }
    
    template <typename Point = DefaultPoint_t>
    Point GetBoxCenter() const
      { return geo::vect::convertTo<Point>(toWorldCoords(LocalPoint_t{ 0.0, 0.0, 0.0 })); }
    
    template <typename Vector = DefaultVector_t>
    Vector GetWireDirection() const { return geo::vect::convertTo<Vector>(fDecompWire.MainDir()); }
    
    

    geo::WireID NearestWireID(geo::Point_t const& pos) const;
    geo::WireID NearestWireID(TVector3 const& pos) const
      { return NearestWireID(geo::vect::toPoint(pos)); }
    
    
    geo::WireGeo const& NearestWire(geo::Point_t const& pos) const;
    
    
    geo::WireID ClosestWireID(geo::WireID::WireID_t wireNo) const;
    
    
    geo::WireID ClosestWireID(geo::WireID const& wireid) const; // inline
    
    

    double DistanceFromPlane(geo::Point_t const& point) const
      { return fDecompWire.PointNormalComponent(point); }
    double DistanceFromPlane(TVector3 const& point) const
      { return DistanceFromPlane(geo::vect::toPoint(point)); }
    
    

    void DriftPoint(geo::Point_t& position, double distance) const
      { position -= distance * GetNormalDirection<geo::Vector_t>(); }
    void DriftPoint(TVector3& position, double distance) const
      { position -= distance * GetNormalDirection(); }
    

    void DriftPoint(geo::Point_t& position) const
      { DriftPoint(position, DistanceFromPlane(position)); }
    void DriftPoint(TVector3& position) const
      { DriftPoint(position, DistanceFromPlane(position)); }
    
    
    Rect const& ActiveArea() const { return fActiveArea; }
    
    lar::util::simple_geo::Volume<> Coverage() const;
    
    template <typename Stream>
    void PrintPlaneInfo
      (Stream&& out, std::string indent = "", unsigned int verbosity = 1) const;
    
    static constexpr unsigned int MaxVerbosity = 6;
    
    
    
    

    double PlaneCoordinateFrom
      (geo::Point_t const& point, geo::WireGeo const& refWire) const
      {
        return
          fDecompWire.VectorSecondaryComponent(point - geo::vect::toPoint(refWire.GetCenter()));
      }
    double PlaneCoordinateFrom
      (TVector3 const& point, geo::WireGeo const& refWire) const
      { return PlaneCoordinateFrom(geo::vect::toPoint(point), refWire); }
    
    

    double PlaneCoordinate(geo::Point_t const& point) const
      { return fDecompWire.PointSecondaryComponent(point); }
    double PlaneCoordinate(TVector3 const& point) const
      { return PlaneCoordinate(geo::vect::toPoint(point)); }
    
    template <typename Point = DefaultPoint_t>
    double WireCoordinate(Point const& point) const
      { return PlaneCoordinate(point) / WirePitch(); }
    
    

    WireDecomposedVector_t DecomposePoint(geo::Point_t const& point) const
      { return fDecompWire.DecomposePoint(point); }
    WireDecomposedVector_t DecomposePoint(TVector3 const& point) const
      { return DecomposePoint(geo::vect::toPoint(point)); }
      
    
    template <typename Point = DefaultPoint_t>
    Point ProjectionReferencePoint() const
      { return geo::vect::convertTo<Point>(fDecompWire.ReferencePoint()); }
    

    WireCoordProjection_t Projection(geo::Point_t const& point) const
      { return fDecompWire.ProjectPointOnPlane(point); }
    WireCoordProjection_t PointProjection(geo::Point_t const& point) const
      { return Projection(point); }
    WireCoordProjection_t PointProjection(TVector3 const& point) const
      { return PointProjection(geo::vect::toPoint(point)); }
    

    WireCoordProjection_t Projection(geo::Vector_t const& v) const
      { return fDecompWire.ProjectVectorOnPlane(v); }
    WireCoordProjection_t VectorProjection(geo::Vector_t const& v) const
      { return Projection(v); }
    WireCoordProjection_t VectorProjection(TVector3 const& v) const
      { return VectorProjection(geo::vect::toVector(v)); }
    
    template <typename Point = DefaultPoint_t>
    Point ComposePoint(WireDecomposedVector_t const& decomp) const
      { return geo::vect::convertTo<Point>(fDecompWire.ComposePoint(decomp)); }
    
    template <typename Point = DefaultPoint_t>
    Point ComposePoint
      (double distance, WireCoordProjection_t const& proj) const
      { return geo::vect::convertTo<Point>(fDecompWire.ComposePoint(distance, proj)); }
    
    
    
    
    

    WDDecomposedVector_t DecomposePointWidthDepth(geo::Point_t const& point) const
      { return fDecompFrame.DecomposePoint(point); }
    WDDecomposedVector_t DecomposePointWidthDepth(TVector3 const& point) const
      { return DecomposePointWidthDepth(geo::vect::toPoint(point)); }
    

    WidthDepthProjection_t PointWidthDepthProjection
      (geo::Point_t const& point) const
      { return fDecompFrame.ProjectPointOnPlane(point); }
    WidthDepthProjection_t PointWidthDepthProjection
      (TVector3 const& point) const
      { return PointWidthDepthProjection(geo::vect::toPoint(point)); }
    

    WidthDepthProjection_t VectorWidthDepthProjection
      (geo::Vector_t const& v) const
      { return fDecompFrame.ProjectVectorOnPlane(v); }
    WidthDepthProjection_t VectorWidthDepthProjection
      (TVector3 const& v) const
      { return VectorWidthDepthProjection(geo::vect::toVector(v)); }
    

    bool isProjectionOnPlane(geo::Point_t const& point) const;
    bool isProjectionOnPlane(TVector3 const& point) const
      { return isProjectionOnPlane(geo::vect::toPoint(point)); }
    
    WidthDepthProjection_t DeltaFromPlane
      (WidthDepthProjection_t const& proj, double wMargin, double dMargin)
       const;
    
    WidthDepthProjection_t DeltaFromPlane
      (WidthDepthProjection_t const& proj, double margin = 0.0) const
      { return DeltaFromPlane(proj, margin, margin); }
    
    WidthDepthProjection_t DeltaFromActivePlane
      (WidthDepthProjection_t const& proj, double wMargin, double dMargin)
      const;
    
    WidthDepthProjection_t DeltaFromActivePlane
      (WidthDepthProjection_t const& proj, double margin = 0.0) const
      { return DeltaFromActivePlane(proj, margin, margin); }
    
    WidthDepthProjection_t MoveProjectionToPlane
      (WidthDepthProjection_t const& proj) const;
    

    geo::Point_t MovePointOverPlane(geo::Point_t const& point) const;
    TVector3 MovePointOverPlane(TVector3 const& point) const;
    
    template <typename Point = DefaultPoint_t>
    Point ComposePoint(WDDecomposedVector_t const& decomp) const
      { return geo::vect::convertTo<Point>(fDecompFrame.ComposePoint(decomp)); }
    
    template <typename Point = DefaultPoint_t>
    Point ComposePoint
      (double distance, WidthDepthProjection_t const& proj) const
      { return geo::vect::convertTo<Point>(fDecompFrame.ComposePoint(distance, proj)); }
    
    
    
    

    void LocalToWorld(const double* plane, double* world) const
      { fTrans.LocalToWorld(plane, world); }
      
    TVector3 LocalToWorld(const TVector3& local) const
      { return fTrans.LocalToWorld<TVector3>(local); }
    
    geo::Point_t toWorldCoords(LocalPoint_t const& local) const
      { return fTrans.toWorldCoords(local); }
    
    void LocalToWorldVect(const double* plane, double* world) const
      { fTrans.LocalToWorldVect(plane, world); }
    
    TVector3 LocalToWorldVect(const TVector3& local) const
      { return fTrans.LocalToWorldVect<TVector3>(local); }
    
    geo::Vector_t toWorldCoords(LocalVector_t const& local) const
      { return fTrans.toWorldCoords(local); }
    
    void WorldToLocal(const double* world, double* plane) const
      { fTrans.WorldToLocal(world, plane); }
    
    TVector3 WorldToLocal(TVector3 const& world) const
      { return fTrans.WorldToLocal<TVector3>(world); }
    
    LocalPoint_t toLocalCoords(geo::Point_t const& world) const
      { return fTrans.toLocalCoords(world); }
    
    void WorldToLocalVect(const double* world, double* plane) const
      { fTrans.WorldToLocalVect(world, plane); }
    
    TVector3 WorldToLocalVect(TVector3 const& world) const
      { return fTrans.WorldToLocalVect<TVector3>(world); }
    
    LocalVector_t toLocalCoords(geo::Vector_t const& world) const
      { return fTrans.toLocalCoords(world); }
    
    
    
    
    
    void SetView(geo::View_t view)                                  { fView = view; }
    
    
    void SortWires(geo::GeoObjectSorter const& sorter);
    
    void UpdateAfterSorting
      (geo::PlaneID planeid, geo::BoxBoundedGeo const& TPCbox);
    
    static std::string ViewName(geo::View_t view);
    
    static std::string OrientationName(geo::Orient_t orientation);
    
    
  private:
    
    void FindWire(GeoNodePath_t& path, size_t depth);
    void MakeWire(GeoNodePath_t& path, size_t depth);
    
    void DetectGeometryDirections();
    
    geo::Vector_t GetNormalAxis() const;
    
    void UpdatePlaneNormal(geo::BoxBoundedGeo const& TPCbox);
    
    void UpdateWidthDepthDir();
    
    void UpdateIncreasingWireDir();
    
    void UpdateWireDir();
    
    void UpdateOrientation();
    
    void UpdateWirePitch();
    
    void UpdateWirePlaneCenter();
    
    void UpdatePhiZ();
      
    void UpdateView();
    
    void UpdateWirePitchSlow();
    
    void UpdateDecompWireOrigin();
    
    void UpdateActiveArea();
    
    bool shouldFlipWire(geo::WireGeo const& wire) const;
    
  private:
    using WireCollection_t = std::vector<geo::WireGeo>;
    
    using LocalTransformation_t
      = geo::LocalTransformationGeo<TGeoHMatrix, LocalPoint_t, LocalVector_t>;
    
    struct RectSpecs {
      double halfWidth;
      double halfDepth;
      
      double HalfWidth() const { return halfWidth; }
      double HalfDepth() const { return halfDepth; }
      double Width() const { return 2.0 * HalfWidth(); }
      double Depth() const { return 2.0 * HalfDepth(); }
    }; // RectSpecs
    
    LocalTransformation_t fTrans;       
    TGeoVolume const*     fVolume;      
    View_t                fView;        
    Orient_t              fOrientation; 
    WireCollection_t      fWire;        
    double                fWirePitch;   
    double                fSinPhiZ;     
    double                fCosPhiZ;     
    
    geo::Vector_t         fNormal;      
    WireDecomposer_t      fDecompWire;
    WidthDepthDecomposer_t fDecompFrame;
    RectSpecs             fFrameSize;   
    Rect                  fActiveArea;
    geo::Point_t          fCenter;

    geo::PlaneID          fID;          
    
    friend struct details::ActiveAreaCalculator;
    
    template <typename T>
    static T boundedValue(T v, T min, T max)
      { return std::min(max, std::max(min, v)); }
    
  }; // class PlaneGeo
  
} // namespace geo


//------------------------------------------------------------------------------
//--- inline implementation
//---
inline geo::WireID geo::PlaneGeo::ClosestWireID
  (geo::WireID::WireID_t wireNo) const
{
  return { ID(), boundedValue<geo::WireID::WireID_t>(wireNo, 0, Nwires()) };
}

inline geo::WireID geo::PlaneGeo::ClosestWireID(geo::WireID const& wireid) const
{
  if (wireid.asPlaneID() != ID()) {
    geo::WireID invalid{ wireid };
    invalid.markInvalid();
    return invalid;
  }
  return ClosestWireID(wireid.Wire);
} // geo::PlaneGeo::ClosestWireID()

//------------------------------------------------------------------------------
//--- template implementation
//---
template <typename Stream>
void geo::PlaneGeo::PrintPlaneInfo(
  Stream&& out,
  std::string indent /* = "" */,
  unsigned int verbosity /* = 1 */
) const {
  
  //----------------------------------------------------------------------------
  out << "plane " << std::string(ID());
  
  if (verbosity-- <= 0) return; // 0
  
  //----------------------------------------------------------------------------
  out
    << " at " << GetCenter<geo::Vector_t>() << " cm"
    << ", theta: " << ThetaZ() << " rad";
  
  if (verbosity-- <= 0) return; // 1
  
  //----------------------------------------------------------------------------
  unsigned int const nWires = Nwires();
  
  out << "\n" << indent
    << "normal to wire: " << PhiZ() << " rad"
      << ", with orientation " << OrientationName(Orientation())
      << ", has " << nWires << " wires measuring " << ViewName(View())
      << " with a wire pitch of " << WirePitch() << " cm"
    ;
  
  if (verbosity-- <= 0) return; // 2
  
  //----------------------------------------------------------------------------
  auto const& normal = GetNormalDirection<geo::Vector_t>();
  auto const& incrZdir = GetIncreasingWireDirection<geo::Vector_t>();
  auto const& wireNormalDir = fDecompWire.NormalDir();
  out << "\n" << indent
    << "normal to plane: " << normal
    << ", direction of increasing wire number: " << incrZdir
    << " [wire frame normal: " << wireNormalDir << "]"
    << " (" << (WireIDincreasesWithZ()? "increases": "decreases") << " with z)";
    
  if (verbosity-- <= 0) return; // 3
  
  //----------------------------------------------------------------------------
  
  auto const& wireDir = GetWireDirection<geo::Vector_t>();
  auto const& widthDir = WidthDir<geo::Vector_t>();
  auto const& depthDir = DepthDir<geo::Vector_t>();
  auto const& frameNormalDir = fDecompFrame.NormalDir();
  
  out << "\n" << indent
    << "wire direction: " << wireDir
    << "; width " << Width() << " cm in direction: " << widthDir
    << ", depth " << Depth() << " cm in direction: " << depthDir
    << " [normal: " << frameNormalDir << "]"
    ;
    
  if (verbosity-- <= 0) return; // 4
  
  //----------------------------------------------------------------------------
  // get the area spanned by the wires
  out << "\n" << indent << "wires cover width "
    << ActiveArea().width.lower << " to " << ActiveArea().width.upper
    << ", depth "
    << ActiveArea().depth.lower << " to " << ActiveArea().depth.upper
    << " cm";
  if (verbosity-- <= 0) return; // 5
  
  //----------------------------------------------------------------------------
  // print also the containing box
  auto const box = BoundingBox();
  out << "\n" << indent
    << "bounding box: " << box.Min() << " -- " << box.Max();
  
//  if (verbosity-- <= 0) return; // 6
  
  //----------------------------------------------------------------------------
} // geo::PlaneGeo::PrintPlaneInfo()


#endif // LARCOREALG_GEOMETRY_PLANEGEO_H