PlaneGeo.h

Go to the documentation of this file.

#ifndef LARCOREALG_GEOMETRY_PLANEGEO_H
#define LARCOREALG_GEOMETRY_PLANEGEO_H

// LArSoft libraries
#include "larcorealg/Geometry/BoxBoundedGeo.h"
#include "larcorealg/Geometry/Decomposer.h"
#include "larcorealg/Geometry/GeoObjectSorter.h"
#include "larcorealg/Geometry/LocalTransformationGeo.h"
#include "larcorealg/Geometry/SimpleGeo.h"
#include "larcorealg/Geometry/TransformationMatrix.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/DisplacementVector2D.h"
#include "Math/GenVector/PositionVector2D.h"

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

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

class TGeoNode;

namespace geo {

  class WireGeo;
  using WirePtr = WireGeo const*; // this will become an object in the future

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

    using ElementIteratorBox = WireCollection_t const&;


    struct PlaneGeoCoordinatesTag {};

    using LocalPoint_t = Point3DBase_t<PlaneGeoCoordinatesTag>;

    using LocalVector_t = Vector3DBase_t<PlaneGeoCoordinatesTag>;



    struct WireCoordinateReferenceTag {};

    using WireCoordProjection_t =
      ROOT::Math::DisplacementVector2D<ROOT::Math::Cartesian2D<double>, WireCoordinateReferenceTag>;

    using WireDecomposer_t = Decomposer<Vector_t, 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 = Decomposer<Vector_t, Point_t, WidthDepthProjection_t>;

    using WDDecomposedVector_t = WidthDepthDecomposer_t::DecomposedVector_t;


    using Rect = lar::util::simple_geo::Rectangle<double>;

    PlaneGeo(TGeoNode const* node, TransformationMatrix&& trans, WireCollection_t&& wires);

    TGeoVolume const* TPCVolume() const { return fNode->GetMotherVolume(); }


    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; }

    PlaneID const& ID() const { return fID; }




    Vector_t const& WidthDir() const { return fDecompFrame.MainDir(); }


    Vector_t const& DepthDir() const { return fDecompFrame.SecondaryDir(); }

    double Width() const { return fFrameSize.Width(); }

    double Depth() const { return fFrameSize.Depth(); }

    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(WireID const& wireid) const { return HasWire(wireid.Wire); }
    bool HasElement(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); }

    auto WirePtr(unsigned int iwire) const { return HasWire(iwire) ? &fWire[iwire] : nullptr; }


    auto WirePtr(WireID const& wireid) const { return WirePtr(wireid.Wire); }
    auto GetElementPtr(WireID const& wireid) const { return WirePtr(wireid); }

    WireGeo const& FirstWire() const { return Wire(0); }

    WireGeo const& MiddleWire() const { return Wire(Nwires() / 2); }

    WireGeo const& LastWire() const { return Wire(Nwires() - 1); }



    double WirePitch() const { return fWirePitch; }

    bool WireIDincreasesWithZ() const;


    Vector_t const& GetNormalDirection() const { return fNormal; }


    Vector_t const& GetIncreasingWireDirection() const { return fDecompWire.SecondaryDir(); }


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


    Point_t GetBoxCenter() const { return toWorldCoords(LocalPoint_t{0.0, 0.0, 0.0}); }


    Vector_t const& GetWireDirection() const { return fDecompWire.MainDir(); }


    WireID NearestWireID(Point_t const& pos) const;

    WireGeo const& NearestWire(Point_t const& pos) const;

    WireID ClosestWireID(WireID::WireID_t wireNo) const; // inline

    WireID ClosestWireID(WireID const& wireid) const; // inline


    double DistanceFromPlane(Point_t const& point) const
    {
      return fDecompWire.PointNormalComponent(point);
    }


    void DriftPoint(Point_t& position, double distance) const
    {
      position -= distance * GetNormalDirection();
    }


    void DriftPoint(Point_t& position) const { DriftPoint(position, DistanceFromPlane(position)); }


    double InterWireProjectedDistance(WireCoordProjection_t const& projDir) const;


    double InterWireDistance(Vector_t const& dir) const;


    double InterWireProjectedDistance(Vector_t const& dir) const
    {
      return InterWireProjectedDistance(Projection(dir));
    }

    Rect const& ActiveArea() const { return fActiveArea; }

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

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

    static constexpr unsigned int MaxVerbosity = 6;




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


    double PlaneCoordinate(Point_t const& point) const
    {
      return fDecompWire.PointSecondaryComponent(point);
    }

    double WireCoordinate(Point_t const& point) const
    {
      return PlaneCoordinate(point) / WirePitch();
    }


    WireDecomposedVector_t DecomposePoint(Point_t const& point) const
    {
      return fDecompWire.DecomposePoint(point);
    }


    Point_t ProjectionReferencePoint() const { return fDecompWire.ReferencePoint(); }


    WireCoordProjection_t Projection(Point_t const& point) const
    {
      return fDecompWire.ProjectPointOnPlane(point);
    }


    WireCoordProjection_t Projection(Vector_t const& v) const
    {
      return fDecompWire.ProjectVectorOnPlane(v);
    }


    Vector_t ComposeVector(WireDecomposedVector_t const& decomp) const
    {
      return fDecompWire.ComposeVector(decomp);
    } //@}


    Vector_t ComposeVector(double distance, WireCoordProjection_t const& proj) const
    {
      return fDecompWire.ComposeVector(distance, proj);
    }


    Point_t ComposePoint(WireDecomposedVector_t const& decomp) const
    {
      return fDecompWire.ComposePoint(decomp);
    }


    Point_t ComposePoint(double distance, WireCoordProjection_t const& proj) const
    {
      return fDecompWire.ComposePoint(distance, proj);
    }




    WDDecomposedVector_t DecomposePointWidthDepth(Point_t const& point) const
    {
      return fDecompFrame.DecomposePoint(point);
    }


    WidthDepthProjection_t PointWidthDepthProjection(Point_t const& point) const
    {
      return fDecompFrame.ProjectPointOnPlane(point);
    }


    WidthDepthProjection_t VectorWidthDepthProjection(Vector_t const& v) const
    {
      return fDecompFrame.ProjectVectorOnPlane(v);
    }


    bool isProjectionOnPlane(Point_t const& point) const;

    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;


    Point_t MovePointOverPlane(Point_t const& point) const;


    Point_t ComposePoint(WDDecomposedVector_t const& decomp) const
    {
      return fDecompFrame.ComposePoint(decomp);
    }


    Point_t ComposePoint(double distance, WidthDepthProjection_t const& proj) const
    {
      return fDecompFrame.ComposePoint(distance, proj);
    }



    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); }



    void SetView(View_t view) { fView = view; }


    void SortWires(Compare<WireGeo> cmp);

    void UpdateAfterSorting(PlaneID planeid, BoxBoundedGeo const& TPCbox);

    static std::string ViewName(View_t view);

    static std::string OrientationName(Orient_t orientation);

  private:
    TGeoVolume const& Volume() const { return *fNode->GetVolume(); }

    void DetectGeometryDirections();

    Vector_t GetNormalAxis() const;

    void UpdatePlaneNormal(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(WireGeo const& wire) const;

    using LocalTransformation_t =
      LocalTransformationGeo<ROOT::Math::Transform3D, 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

    TGeoNode const* fNode;            
    LocalTransformation_t fTrans;     
    View_t fView{kUnknown};           
    Orient_t fOrientation{kVertical}; 
    WireCollection_t fWire;           
    double fWirePitch{0.};            
    double fSinPhiZ{0.};              
    double fCosPhiZ{0.};              

    Vector_t fNormal; 
    WireDecomposer_t fDecompWire{};
    WidthDepthDecomposer_t fDecompFrame{};
    RectSpecs fFrameSize; 
    Rect fActiveArea;
    Point_t fCenter{};

    PlaneID fID; 

    // friend struct details::ActiveAreaCalculator;

  }; // class PlaneGeo

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

} // namespace geo

//------------------------------------------------------------------------------
//--- inline implementation
//---
inline geo::WireID geo::PlaneGeo::ClosestWireID(WireID::WireID_t wireNo) const
{
  return {ID(), std::min(wireNo, Nwires())};
}

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

//------------------------------------------------------------------------------
//--- 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() << " 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();
  auto const& incrZdir = GetIncreasingWireDirection();
  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();
  auto const& widthDir = WidthDir();
  auto const& depthDir = DepthDir();
  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();
}

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

#endif // LARCOREALG_GEOMETRY_PLANEGEO_H