TPCGeo.h

Go to the documentation of this file.

#ifndef LARCOREALG_GEOMETRY_TPCGEO_H
#define LARCOREALG_GEOMETRY_TPCGEO_H

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

// ROOT libraries
#include "TGeoVolume.h"
#include "TGeoMatrix.h"
#include "TVector3.h"

// C/C++ standard library
#include <set>
#include <vector>


class TGeoNode;

namespace geo {

  //......................................................................
  class TPCGeo: public BoxBoundedGeo {
    
    using DefaultVector_t = TVector3; // ... not for long
    using DefaultPoint_t = TVector3; // ... not for long
    
  public:
    
    using GeoNodePath_t = geo::WireGeo::GeoNodePath_t;
    

    struct TPCGeoCoordinatesTag {};
    
    using LocalPoint_t = geo::Point3DBase_t<TPCGeoCoordinatesTag>;
    
    using LocalVector_t = geo::Vector3DBase_t<TPCGeoCoordinatesTag>;
    
    
    using Projection_t = geo::PlaneGeo::WidthDepthProjection_t;
    
    using DecomposedVector_t = geo::PlaneGeo::WDDecomposedVector_t;
    
    
    // Construct a representation of a single plane of the detector
    TPCGeo(GeoNodePath_t& path, size_t depth);
    
    
    
    double            ActiveHalfWidth()                         const { return fActiveHalfWidth;        }
    double            ActiveWidth()                             const { return 2.0 * ActiveHalfWidth(); }
    double            ActiveHalfHeight()                        const { return fActiveHalfHeight;       }
    double            ActiveHeight()                            const { return 2.0 * ActiveHalfHeight(); }
    double            ActiveLength()                            const { return fActiveLength;           }
    double            ActiveHalfLength()                        const { return fActiveLength / 2.0;     }
    double            HalfWidth()                               const { return fHalfWidth;              }
    double            Width()                                   const { return 2.0 * HalfWidth();       }
    double            HalfHeight()                              const { return fHalfHeight;             }
    double            Height()                                  const { return 2.0 * HalfHeight();      }
    double            Length()                                  const { return fLength;                 }
    double            HalfLength()                              const { return fLength / 2.0;           }
    double            ActiveMass()                              const { return fActiveVolume->Weight(); }
    const TGeoVolume* ActiveVolume()                            const { return fActiveVolume;           }
    const TGeoVolume* TotalVolume()                             const { return fTotalVolume;            }
    
    template <typename Vector = DefaultVector_t>
    decltype(auto)    WidthDir()                                const { return geo::vect::convertTo<Vector>(fWidthDir); }
    
    template <typename Vector = DefaultVector_t>
    decltype(auto)    HeightDir()                               const { return geo::vect::convertTo<Vector>(fHeightDir); }
    
    template <typename Vector = DefaultVector_t>
    decltype(auto)    LengthDir()                               const { return geo::vect::convertTo<Vector>(fLengthDir); }
    
    DriftDirection_t  DriftDirection()                          const { return fDriftDirection;         }
    
    template <typename Vector = DefaultVector_t>
    Vector DriftDir()                                           const;
    
    double            DriftDistance()                           const { return ComputeDriftDistance();  }
    
    
    
    
    unsigned int    Nplanes()                                   const { return fPlanes.size();   }
    unsigned int    NElements()                                 const { return fPlanes.size();   }
    

    bool HasPlane(unsigned int iplane) const { return iplane < Nplanes(); }
    bool HasElement(unsigned int iplane) const { return HasPlane(iplane); }
    

    bool HasPlane(geo::PlaneID const& planeid) const
      { return HasPlane(planeid.Plane); }
    bool HasElement(geo::PlaneID const& planeid) const
      { return HasPlane(planeid); }
    
    PlaneGeo const&     Plane(geo::View_t  view)                const;
    
    
    PlaneGeo const&     Plane(unsigned int iplane)              const;
    

    const PlaneGeo&     Plane(PlaneID const& planeid)           const
      { return Plane(planeid.Plane); }
    const PlaneGeo&     GetElement(PlaneID const& planeid)            const
      { return Plane(planeid); }
    
    PlaneGeo const*     PlanePtr(unsigned int iplane)                const
      { return HasPlane(iplane)? &(fPlanes[iplane]): nullptr; }
    

    PlaneGeo const*     PlanePtr(PlaneID const& planeid)                   const
      { return PlanePtr(planeid.Plane); }
    PlaneGeo const*     GetElementPtr(PlaneID const& planeid)            const
      { return PlanePtr(planeid); }
    
    geo::PlaneGeo const& SmallestPlane() const;
    
    geo::PlaneGeo const& FirstPlane() const { return fPlanes[0]; }
    
    geo::PlaneGeo const& LastPlane() const { return fPlanes[Nplanes() - 1]; }
    
    unsigned int MaxWires() const;
    
    std::set<geo::View_t> Views() const;
    
    
    
    template <typename Point = DefaultPoint_t>
    Point GetCenter() const;
    
    template <typename Point = DefaultPoint_t>
    Point GetActiveVolumeCenter() const
      { return geo::vect::convertTo<Point>(fActiveCenter); }
    
    template <typename Point = DefaultPoint_t>
    Point GetCathodeCenter() const
      { return geo::vect::convertTo<Point>(GetCathodeCenterImpl()); }
    
    template <typename Point = geo::Point_t>
    Point GetFrontFaceCenter() const
      { return geo::vect::convertTo<Point>(GetFrontFaceCenterImpl()); }
    
    geo::BoxBoundedGeo const& BoundingBox() const
      { return *this; }
    
    geo::BoxBoundedGeo const& ActiveBoundingBox() const
      { return fActiveBox; }
    
    
    const double*     PlaneLocation(unsigned int p)             const; 
    double            Plane0Pitch(unsigned int p)               const;
    double            PlanePitch(unsigned int p1=0,
                                 unsigned int p2=1)             const;
    double WirePitch(unsigned plane = 0) const;
    
    geo::TPCID const& ID() const { return fID; }
    
    
    
    
    geo::PlaneGeo const& ReferencePlane() const { return FirstPlane(); }
    
    geo::PlaneID const& ReferencePlaneID() const
      { return ReferencePlane().ID(); }
    
    template <typename Vector = DefaultPoint_t>
    Vector RefWidthDir() const { return ReferencePlane().WidthDir<Vector>(); }
    
    template <typename Vector = DefaultPoint_t>
    Vector RefDepthDir() const { return ReferencePlane().DepthDir<Vector>(); }
    
    

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

    void DriftPoint(geo::Point_t& position, double distance) const
      { ReferencePlane().DriftPoint(position, distance); }
    void DriftPoint(TVector3& position, double distance) const
      { ReferencePlane().DriftPoint(position, distance); }
    
    

    DecomposedVector_t DecomposePoint(geo::Point_t const& point) const
      { return ReferencePlane().DecomposePointWidthDepth(point); }
    DecomposedVector_t DecomposePoint(TVector3 const& point) const
      { return DecomposePoint(geo::vect::toPoint(point)); }
    
    template <typename Point = DefaultPoint_t>
    Point ProjectionReferencePoint() const
      { return ReferencePlane().GetCenter<Point>(); }
    

    Projection_t Projection(geo::Point_t const& point) const
      { return ReferencePlane().PointWidthDepthProjection(point); }
    Projection_t PointProjection(geo::Point_t const& point) const
      { return Projection(point); }
    Projection_t PointProjection(TVector3 const& point) const
      { return Projection(geo::vect::toPoint(point)); }
    

    Projection_t Projection(geo::Vector_t const& v) const
      { return ReferencePlane().VectorWidthDepthProjection(v); }
    Projection_t VectorProjection(geo::Vector_t const& v) const
      { return Projection(v); }
    Projection_t VectorProjection(TVector3 const& v) const
      { return Projection(geo::vect::toVector(v)); }
    
    template <typename Point = DefaultPoint_t>
    Point ComposePoint(DecomposedVector_t const& decomp) const
      { return ReferencePlane().ComposePoint<Point>(decomp); }
    
    template <typename Point = DefaultPoint_t>
    Point ComposePoint(double distance, Projection_t const& proj) const
      { return ReferencePlane().ComposePoint<Point>(distance, proj); }
    
    
    
    
    void LocalToWorld(const double* tpc, double* world) const
      { fTrans.LocalToWorld(tpc, 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* tpc, double* world) const
      { fTrans.LocalToWorldVect(tpc, 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* tpc) const
      { fTrans.WorldToLocal(world, tpc); }
    
    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* tpc) const
      { fTrans.WorldToLocalVect(world, tpc); }
    
    TVector3 WorldToLocalVect(TVector3 const& world) const
      { return fTrans.WorldToLocalVect<TVector3>(world); }
    
    LocalVector_t toLocalCoords(geo::Vector_t const& world) const
      { return fTrans.toLocalCoords(world); }
    
    
    
    short int DetectDriftDirection() const;
    
    
    
    void              SortSubVolumes(geo::GeoObjectSorter const& sorter);
    
    void UpdateAfterSorting(geo::TPCID tpcid);
    
    
    template <typename Stream>
    void PrintTPCInfo
      (Stream&& out, std::string indent = "", unsigned int verbosity = 1) const;
    
    static constexpr unsigned int MaxVerbosity = 6;
    
    
    static bool CoordinateContained
      (double c, double min, double max, double wiggle = 1.)
      {
        return (c >= (min > 0? min / wiggle: min * wiggle))
          && (c <= (max < 0? max / wiggle: max * wiggle));
      } // CoordinateContained()
    
    static bool CoordinateContained
      (double c, double const* range, double wiggle = 1.)
      { return CoordinateContained(c, range[0], range[1], wiggle); }
    
    
  private:
    
    void FindPlane(GeoNodePath_t& path, size_t depth);
    void MakePlane(GeoNodePath_t& path, size_t depth);
    
  private:
    using LocalTransformation_t
      = geo::LocalTransformationGeo<TGeoHMatrix, LocalPoint_t, LocalVector_t>;
    
    LocalTransformation_t              fTrans;          
    
    std::vector<PlaneGeo>              fPlanes;         
    TGeoVolume*                        fActiveVolume;   
    TGeoVolume*                        fTotalVolume;    
    DriftDirection_t                   fDriftDirection; 
    std::vector<double>                fPlane0Pitch;    
    std::vector< std::vector<double> > fPlaneLocation;  
    geo::Point_t                       fActiveCenter;   

    double                             fActiveHalfWidth;  
    double                             fActiveHalfHeight; 
    double                             fActiveLength;     
    double                             fHalfWidth;        
    double                             fHalfHeight;       
    double                             fLength;           

    geo::Vector_t fWidthDir; 
    geo::Vector_t fHeightDir; 
    geo::Vector_t fLengthDir; 
    geo::Vector_t fDriftDir; 
    
    geo::BoxBoundedGeo                 fActiveBox;      
    
    geo::TPCID                         fID;             
    
    std::vector<geo::PlaneID::PlaneID_t> fViewToPlaneNumber;
  
    void ResetDriftDirection();
    
    double ComputeDriftDistance() const;
    
    void UpdatePlaneViewCache();
    
    void InitTPCBoundaries();
    
    void SortPlanes(std::vector<geo::PlaneGeo>&) const;

    geo::Point_t GetFrontFaceCenterImpl() const;
    geo::Point_t GetCathodeCenterImpl() const;
  
  };
}


//------------------------------------------------------------------------------
//--- template implementation
//---
//------------------------------------------------------------------------------
template <typename Vector /* = DefaultVector_t */>
Vector geo::TPCGeo::DriftDir() const
  { return geo::vect::convertTo<Vector>(fDriftDir); }


//------------------------------------------------------------------------------
template <typename Point /* = DefaultPoint_t */>
Point geo::TPCGeo::GetCenter() const {
  
  // convert the origin (default constructed TVector)
  return geo::vect::convertTo<Point>(toWorldCoords(LocalPoint_t{ 0.0, 0.0, 0.0 }));
  
} // geo::TPCGeo::GetCenter()


//------------------------------------------------------------------------------
template <typename Stream>
void geo::TPCGeo::PrintTPCInfo(
  Stream&& out,
  std::string indent /* = "" */,
  unsigned int verbosity /* = 1 */
) const {
  
  //----------------------------------------------------------------------------
  out << "TPC " << std::string(ID());
  
  if (verbosity-- <= 0) return; // 0
  
  //----------------------------------------------------------------------------
  out
    << " (" << Width() << " x " << Height() << " x " << Length() << ") cm^3 at "
      << GetCenter<geo::Point_t>();
  
  if (verbosity-- <= 0) return; // 1
  
  //----------------------------------------------------------------------------
  
  out << "\n" << indent
    << "drift direction " << DriftDir<geo::Vector_t>()
      << " from cathode around " << GetCathodeCenter<geo::Point_t>()
      << " through " << DriftDistance() << " cm toward "
      << Nplanes() << " wire planes"
    ;
  
  if (verbosity-- <= 0) return; // 2
  
  //----------------------------------------------------------------------------
  out << "\n" << indent
    << "maximum wires on any plane: " << MaxWires();
    
  if (verbosity-- <= 0) return; // 3
  
  //----------------------------------------------------------------------------
  out << "\n" << indent
    << "active volume ("
      << ActiveWidth() << " x " << ActiveHeight() << " x " << ActiveLength()
      << ") cm^3, front face at " << GetFrontFaceCenter<geo::Point_t>()
      << " cm;"
    << "\n" << indent << "main directions:"
      << " width " << WidthDir<geo::Vector_t>()
      << " height " << HeightDir<geo::Vector_t>()
      << " length " << LengthDir<geo::Vector_t>()
    ;
    
  if (verbosity-- <= 0) return; // 4
  
  //----------------------------------------------------------------------------
  // print also the containing box
  geo::BoxBoundedGeo const& box = BoundingBox();
  out << "\n" << indent
    << "bounding box: " << box.Min() << " -- " << box.Max();
  
//  if (verbosity-- <= 0) return; // 5
  
  //----------------------------------------------------------------------------
  // print also the active box
  geo::BoxBoundedGeo const& activeBox = ActiveBoundingBox();
  out << "\n" << indent
    << "active volume box: " << activeBox.Min() << " -- "  << activeBox.Max();
  
//  if (verbosity-- <= 0) return; // 6
  
  //----------------------------------------------------------------------------
} // geo::TPCGeo::PrintTPCInfo()

#endif // LARCOREALG_GEOMETRY_TPCGEO_H