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 <string>
#include <vector>

class TGeoNode;

namespace geo {

  namespace details {
    struct ActiveAreaCalculator;
  } // namespace details

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

  // 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 {
  public:
    using ID_t = PlaneID;
    using WireCollection_t = std::vector<geo::WireGeo>;
    using GeoNodePath_t = std::vector<TGeoNode const*>;

    using ElementIteratorBox = WireCollection_t 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(TGeoNode const& node, geo::TransformationMatrix&& trans, WireCollection_t&& wires);


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




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

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

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


    geo::WirePtr WirePtr(WireID const& wireid) const { return WirePtr(wireid.Wire); }
    geo::WirePtr 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); }

    // @{
    ElementIteratorBox IterateElements() const { return fWire; }
    ElementIteratorBox IterateWires() const { return IterateElements(); }
    // @}



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


    geo::WireID NearestWireID(geo::Point_t const& pos) const;

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

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

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


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


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


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


    double InterWireProjectedDistance(WireCoordProjection_t const& projDir) const;


    double InterWireDistance(geo::Vector_t const& dir) const;


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

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

    // lar::util::simple_geo::Volume<Point_t> Coverage() const;

    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(geo::Point_t const& point, geo::WireGeo const& refWire) const
    {
      return fDecompWire.VectorSecondaryComponent(point - refWire.GetCenter());
    }


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

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


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


    Point_t ProjectionReferencePoint() const { return 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 Projection(geo::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(geo::Point_t const& point) const
    {
      return fDecompFrame.DecomposePoint(point);
    }


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


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


    bool isProjectionOnPlane(geo::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;


    geo::Point_t MovePointOverPlane(geo::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);
    }



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

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

    LocalPoint_t toLocalCoords(geo::Point_t const& world) const
    {
      return fTrans.toLocalCoords(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 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 LocalTransformation_t =
      geo::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

    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;

  }; // class PlaneGeo

} // namespace geo

//------------------------------------------------------------------------------
//--- inline implementation
//---
inline geo::WireID geo::PlaneGeo::ClosestWireID(geo::WireID::WireID_t wireNo) const
{
  return {ID(), std::min(wireNo, 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() << " 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();

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

  //----------------------------------------------------------------------------
} // geo::PlaneGeo::PrintPlaneInfo()

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

#endif // LARCOREALG_GEOMETRY_PLANEGEO_H