TPCGeo.h

Go to the documentation of this file.

#ifndef LARCOREALG_GEOMETRY_TPCGEO_H
#define LARCOREALG_GEOMETRY_TPCGEO_H

// LArSoft libraries
#include "larcorealg/Geometry/BoxBoundedGeo.h"
#include "larcorealg/Geometry/GeoObjectSorter.h"
#include "larcorealg/Geometry/LocalTransformationGeo.h"
#include "larcorealg/Geometry/PlaneGeo.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 "TGeoMatrix.h"
#include "TGeoVolume.h"

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

class TGeoNode;

namespace geo {

  //......................................................................
  class TPCGeo : public BoxBoundedGeo {
  public:
    using ID_t = TPCID;
    using PlaneCollection_t = std::vector<geo::PlaneGeo>;
    using GeoNodePath_t = geo::WireGeo::GeoNodePath_t;

    using ElementIteratorBox = PlaneCollection_t const&;


    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(TGeoNode const& node, geo::TransformationMatrix&& trans, PlaneCollection_t&& planes);


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

    decltype(auto) WidthDir() const { return fWidthDir; }

    decltype(auto) HeightDir() const { return fHeightDir; }

    decltype(auto) LengthDir() const { return fLengthDir; }

    DriftDirection_t DriftDirection() const { return fDriftDirection; }

    Vector_t DriftDir() const { return fDriftDir; }

    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;

    // @{
    ElementIteratorBox IterateElements() const;
    ElementIteratorBox IteratePlanes() const { return IterateElements(); }
    // @}

    std::set<geo::View_t> Views() const;



    Point_t GetCenter() const { return toWorldCoords(LocalPoint_t{}); }

    Point_t GetActiveVolumeCenter() const { return fActiveCenter; }

    Point_t GetCathodeCenter() const { return GetCathodeCenterImpl(); }

    geo::Point_t GetFrontFaceCenter() const { return GetFrontFaceCenterImpl(); }

    geo::BoxBoundedGeo const& BoundingBox() const { return *this; }

    geo::BoxBoundedGeo const& ActiveBoundingBox() const { return fActiveBox; }

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


    Vector_t RefWidthDir() const { return ReferencePlane().WidthDir(); }


    Vector_t RefDepthDir() const { return ReferencePlane().DepthDir(); }


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


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


    DecomposedVector_t DecomposePoint(geo::Point_t const& point) const
    {
      return ReferencePlane().DecomposePointWidthDepth(point);
    }


    Point_t ProjectionReferencePoint() const { return ReferencePlane().GetCenter(); }


    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 Projection(geo::Vector_t const& v) const
    {
      return ReferencePlane().VectorWidthDepthProjection(v);
    }
    Projection_t VectorProjection(geo::Vector_t const& v) const { return Projection(v); }


    Point_t ComposePoint(DecomposedVector_t const& decomp) const
    {
      return ReferencePlane().ComposePoint(decomp);
    }


    Point_t ComposePoint(double distance, Projection_t const& proj) const
    {
      return ReferencePlane().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);
    }


    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;

    std::string TPCInfo(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<ROOT::Math::Transform3D, LocalPoint_t, LocalVector_t>;

    LocalTransformation_t fTrans; 

    PlaneCollection_t 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 UpdatePlaneCache();

    void InitTPCBoundaries();

    void SortPlanes(std::vector<geo::PlaneGeo>&) const;

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

//------------------------------------------------------------------------------
//--- template implementation
//---
//------------------------------------------------------------------------------
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();

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

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

  out << "\n"
      << indent << "drift direction " << DriftDir() << " from cathode around " << GetCathodeCenter()
      << " 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() << " cm;"
      << "\n"
      << indent << "main directions:"
      << " width " << WidthDir() << " height " << HeightDir() << " length " << LengthDir();

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

  //----------------------------------------------------------------------------
  // print also the containing box
  geo::BoxBoundedGeo const& box = BoundingBox();
  out << "\n" << indent << "bounding box: " << box.Min() << " -- " << box.Max();

  //----------------------------------------------------------------------------
  // print also the active box
  geo::BoxBoundedGeo const& activeBox = ActiveBoundingBox();
  out << "\n" << indent << "active volume box: " << activeBox.Min() << " -- " << activeBox.Max();

  //----------------------------------------------------------------------------
} // geo::TPCGeo::PrintTPCInfo()

#endif // LARCOREALG_GEOMETRY_TPCGEO_H