WireGeo.h

Go to the documentation of this file.

#ifndef LARCOREALG_GEOMETRY_WIREGEO_H
#define LARCOREALG_GEOMETRY_WIREGEO_H

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

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

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


// forward declarations
class TGeoNode;


namespace geo {
  
  struct WireID; // forward declaration
  
  
  class WireGeo {
    
    using DefaultVector_t = TVector3; // ... not for long
    using DefaultPoint_t = TVector3; // ... not for long
    
  public:
    
    using GeoNodePath_t = std::vector<TGeoNode const*>;
    

    struct WireGeoCoordinatesTag {};
    
    using LocalPoint_t = geo::Point3DBase_t<WireGeoCoordinatesTag>;
    
    using LocalVector_t = geo::Vector3DBase_t<WireGeoCoordinatesTag>;
    
    
    WireGeo(GeoNodePath_t const& path, size_t depth);
    
    
    
    double RMax() const;
    
    double HalfL() const { return fHalfL; }
    
    double RMin() const;
    
    void GetCenter(double* xyz, double localz=0.0) const;
    
    void GetStart(double* xyz) const { GetCenter(xyz, -fHalfL); }
    
    void GetEnd(double* xyz) const { GetCenter(xyz, +fHalfL); }

    template <typename Point = DefaultPoint_t>
    Point GetPositionFromCenter(double localz) const
      { return GetPositionFromCenterUnbounded<Point>(capLength(localz)); }
    
    template <typename Point = DefaultPoint_t>
    Point GetPositionFromCenterUnbounded(double localz) const;
    
    template <typename Point = DefaultPoint_t>
    Point GetCenter() const { return geo::vect::convertTo<Point>(fCenter); }
    
    template <typename Point = DefaultPoint_t>
    Point GetStart() const
      { return GetPositionFromCenterUnbounded<Point>(-HalfL()); }
    
    template <typename Point = DefaultPoint_t>
    Point GetEnd() const
      { return GetPositionFromCenterUnbounded<Point>(+HalfL()); }
    
    double Length() const { return 2. * HalfL(); }
    
    
    
    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(geo::WireGeo const& wire) const
      {
        return // parallel if the dot product of the directions is about +/- 1
          std::abs(std::abs(Direction<geo::Vector_t>().Dot(wire.Direction<geo::Vector_t>())) - 1.) < 1e-5;
      }
    
    template <typename Vector = DefaultVector_t>
    Vector Direction() const;
    
    template <typename Stream>
    void PrintWireInfo
      (Stream&& out, std::string indent = "", unsigned int verbosity = 1) const;
    
    static constexpr unsigned int MaxVerbosity = 4;
    
    
    
    

    void LocalToWorld(const double* wire, double* world) const
      { fTrans.LocalToWorld(wire, world); }
      
    geo::Point_t toWorldCoords(LocalPoint_t const& local) const
      { return fTrans.toWorldCoords(local); }
    
    void LocalToWorldVect(const double* wire, double* world) const
      { fTrans.LocalToWorldVect(wire, world); }
    
    geo::Vector_t toWorldCoords(LocalVector_t const& local) const
      { return fTrans.toWorldCoords(local); }
    
    void WorldToLocal(const double* world, double* wire) const
      { fTrans.WorldToLocal(world, wire); }
    
    LocalPoint_t toLocalCoords(geo::Point_t const& world) const
      { return fTrans.toLocalCoords(world); }
    
    void WorldToLocalVect(const double* world, double* wire) const
      { fTrans.WorldToLocalVect(world, wire); }
    
    LocalVector_t toLocalCoords(geo::Vector_t const& world) const
      { return fTrans.toLocalCoords(world); }
    
    
    
    const TGeoNode*     Node() const { return fWireNode; }
    
    double ComputeZatY0() const
      { return fCenter.Z() - fCenter.Y() / TanThetaZ(); }
    
    double DistanceFrom(geo::WireGeo const& wire) const;
    
    
    void UpdateAfterSorting(geo::WireID const&, bool flip);
    
    static double WirePitch(geo::WireGeo const& w1, geo::WireGeo const& w2)
      { return std::abs(w2.DistanceFrom(w1)); }
    
  private:
    using LocalTransformation_t
      = geo::LocalTransformationGeo<TGeoHMatrix, LocalPoint_t, LocalVector_t>;
    
    const TGeoNode*    fWireNode;  
    double             fThetaZ;    
    double             fHalfL;     
    geo::Point_t       fCenter;    
    LocalTransformation_t fTrans;  
    bool               flipped;    
    
    bool isFlipped() const { return flipped; }
    
    double relLength(double local) const { return isFlipped()? -local: local; }
    
    double capLength(double local) const
      { return std::min(+HalfL(), std::max(-HalfL(), local)); }
    
    double capRelLength(double local) const
      { return capLength(relLength(local)); }
    
    void Flip();
    
    
    static double gausSum(double a, double b) { return std::sqrt(a*a + b*b); }
    
  };
}


//------------------------------------------------------------------------------
//--- template implementation
//---
//------------------------------------------------------------------------------
template <typename Point>
Point geo::WireGeo::GetPositionFromCenterUnbounded(double localz) const {
  return geo::vect::convertTo<Point>
    (toWorldCoords(LocalPoint_t{ 0.0, 0.0, relLength(localz) }));
} // geo::WireGeo::GetPositionFromCenterImpl()


//------------------------------------------------------------------------------
template <typename Vector /* = DefaultVector_t */>
Vector geo::WireGeo::Direction() const {
  // maybe (GetCenter() - GetStart()) / HalfL() would be faster;
  // strangely, TVector3 does not implement operator/ (double).
  return geo::vect::convertTo<Vector>(GetEnd<geo::Point_t>() - GetStart<geo::Point_t>()) * (1.0 / Length());
} // geo::WireGeo::Direction()


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

  //----------------------------------------------------------------------------
  out << "wire from " << GetStart<geo::Point_t>()
    << " to " << GetEnd<geo::Point_t>();
  
  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<geo::Point_t>() << " cm";
    
  if (verbosity-- <= 0) return; // 3
  
  //----------------------------------------------------------------------------
  out << ", direction: " << Direction<geo::Vector_t>();
  if (isHorizontal()) out << " (horizontal)";
  if (isVertical()) out << " (vertical)";
    
//  if (verbosity-- <= 0) return; // 4
  
  //----------------------------------------------------------------------------
} // geo::WireGeo::PrintWireInfo()


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

#endif // LARCOREALG_GEOMETRY_WIREGEO_H