GeometryCore.h

Go to the documentation of this file.

#ifndef LARCOREALG_GEOMETRY_GEOMETRYCORE_H
#define LARCOREALG_GEOMETRY_GEOMETRYCORE_H

// LArSoft libraries
#include "larcorealg/CoreUtils/RealComparisons.h"
#include "larcorealg/Geometry/AuxDetGeo.h"
#include "larcorealg/Geometry/AuxDetSensitiveGeo.h"
#include "larcorealg/Geometry/BoxBoundedGeo.h"
#include "larcorealg/Geometry/ChannelMapAlg.h"
#include "larcorealg/Geometry/CryostatGeo.h"
#include "larcorealg/Geometry/GeoObjectSorter.h"
#include "larcorealg/Geometry/GeometryBuilder.h"
#include "larcorealg/Geometry/GeometryData.h"
#include "larcorealg/Geometry/GeometryDataContainers.h" // geo::TPCDataContainer
#include "larcorealg/Geometry/OpDetGeo.h"
#include "larcorealg/Geometry/PlaneGeo.h"
#include "larcorealg/Geometry/ReadoutDataContainers.h" // readout::ROPDataContainer
#include "larcorealg/Geometry/TPCGeo.h"
#include "larcorealg/Geometry/WireGeo.h"
#include "larcorealg/Geometry/details/geometry_iterator_maker.h"
#include "larcorealg/Geometry/details/geometry_iterators.h"
#include "larcorealg/Geometry/fwd.h"
#include "larcorealg/Geometry/geo_vectors_utils.h"       // geo::vect namespace
#include "larcoreobj/SimpleTypesAndConstants/RawTypes.h" // raw::ChannelID_t
#include "larcoreobj/SimpleTypesAndConstants/geo_types.h"
#include "larcoreobj/SimpleTypesAndConstants/geo_vectors.h"
#include "larcoreobj/SimpleTypesAndConstants/readout_types.h"

// Framework and infrastructure libraries
#include "fhiclcpp/ParameterSet.h"

// C/C++ standard libraries
#include <cstddef>  // size_t
#include <iterator> // std::forward_iterator_tag
#include <memory>   // std::shared_ptr<>
#include <set>
#include <string>
#include <type_traits> // std::is_base_of<>
#include <utility>
#include <vector>

// ROOT class prototypes
class TGeoManager;
class TGeoNode;
class TGeoVolume;
class TGeoMaterial;

namespace geo {

  // BEGIN Geometry group ------------------------------------------------------

  //
  // GeometryCore
  //

  class GeometryCore {
  public:
    template <typename Point>
    struct Segment : public std::pair<Point, Point> {

      // use the base class constructors
      using std::pair<Point, Point>::pair;

      Point const& start() const { return this->first; }
      Point& start() { return this->first; }

      Point const& end() const { return this->second; }
      Point& end() { return this->second; }

    }; // struct Segment_t

    using Segment_t = Segment<Point_t>;

    using CryostatList_t = GeometryData_t::CryostatList_t;
    using AuxDetList_t = GeometryData_t::AuxDetList_t;

    static constexpr std::size_t MaxWireDepthInGDML = 20U;

    explicit GeometryCore(fhicl::ParameterSet const& pset);

    // this object is not copiable nor moveable (see also issue #14384);
    // currently, auxiliary detectors are stored as bare pointers,
    // which prevents trivial copy or move.
    GeometryCore(GeometryCore const&) = delete;
    GeometryCore(GeometryCore&&) = delete;
    GeometryCore& operator=(GeometryCore const&) = delete;
    GeometryCore& operator=(GeometryCore&&) = delete;

    double DefaultWiggle() const { return fPositionWiggle; }

    std::string const& ROOTFile() const { return fROOTfile; }

    std::string const& GDMLFile() const { return fGDMLfile; }

    // BEGIN Detector information

    //
    // global features
    //
    std::string const& DetectorName() const { return fDetectorName; }

    //
    // position
    //

    TGeoVolume const* WorldVolume() const;

    void WorldBox(double* xlo, double* xhi, double* ylo, double* yhi, double* zlo, double* zhi)
      const;

    BoxBoundedGeo WorldBox() const;

    //
    Length_t SurfaceY() const { return fSurfaceY; }

    //
    // object description and information
    //

    TGeoManager* ROOTGeoManager() const;

    const std::string GetWorldVolumeName() const;

    BoxBoundedGeo DetectorEnclosureBox(std::string const& name = "volDetEnclosure") const;


    std::string VolumeName(Point_t const& point) const;

    std::vector<TGeoNode const*> FindAllVolumes(std::set<std::string> const& vol_names) const;

    std::vector<std::vector<TGeoNode const*>> FindAllVolumePaths(
      std::set<std::string> const& vol_names) const;

    TGeoMaterial const* Material(Point_t const& point) const;

    std::string MaterialName(Point_t const& point) const;

    double TotalMass() const { return TotalMass(GetWorldVolumeName()); }
    double TotalMass(std::string vol) const;


    double MassBetweenPoints(Point_t const& p1, Point_t const& p2) const;

    template <typename Stream>
    void Print(Stream&& out, std::string indent = "  ") const;

    std::string Info(std::string indent = "  ") const;

    // END Detector information

    template <typename GeoID>
    GeoID GetBeginID() const
    {
      GeoID id;
      GetBeginID(id);
      return id;
    }

    template <typename GeoID>
    GeoID GetNextID(GeoID const& id) const
    {
      auto nextID(id);
      IncrementID(nextID);
      return nextID;
    }

    template <typename GeoID>
    GeoID GetEndID() const
    {
      GeoID id;
      GetEndID(id);
      return id;
    }

    template <typename GeoID, typename ContextID>
    GeoID GetBeginID(ContextID const& id) const;

    template <typename GeoID, typename ContextID>
    GeoID GetEndID(ContextID const& id) const;


    //
    // group features
    //


    unsigned int Ncryostats() const { return Cryostats().size(); }
    unsigned int NElements() const { return Ncryostats(); }
    unsigned int NSiblingElements(CryostatID const&) const { return Ncryostats(); }

    //
    // access
    //


    bool HasCryostat(CryostatID const& cryoid) const { return cryoid.Cryostat < Ncryostats(); }
    bool HasElement(CryostatID const& cryoid) const { return HasCryostat(cryoid); }


    static constexpr CryostatID cryostat_zero{0};
    CryostatGeo const& Cryostat(CryostatID const& cryoid = cryostat_zero) const;
    CryostatGeo const& GetElement(CryostatID const& cryoid) const { return Cryostat(cryoid); }


    CryostatGeo const* CryostatPtr(CryostatID const& cryoid) const
    {
      return HasCryostat(cryoid) ? &(Cryostats()[cryoid.Cryostat]) : nullptr;
    }
    CryostatGeo const* GetElementPtr(CryostatID const& cryoid) const { return CryostatPtr(cryoid); }

    CryostatGeo const* PositionToCryostatPtr(Point_t const& point) const;

    CryostatID PositionToCryostatID(Point_t const& point) const;


    CryostatGeo const& PositionToCryostat(Point_t const& point) const;

    //
    // iterators
    //

    void GetBeginID(CryostatID& id) const { id = CryostatID(0, HasCryostat(CryostatID(0))); }

    void GetEndID(CryostatID& id) const { id = CryostatID(Ncryostats(), false); }

    bool IncrementID(CryostatID& id) const; // inline implementation

    // FIXME: Need Doxygen comments
    template <typename T>
    details::begin_type<T> begin() const
    {
      return details::IteratorMaker<T>::create_begin(this);
    }

    template <typename T>
    details::end_type<T> end() const
    {
      return details::IteratorMaker<T>::create_end(this);
    }

    template <typename T>
    details::range_type<T> Iterate() const
    {
      return details::IteratorMaker<T>::create_range(this);
    }

    template <typename T, typename BaseID>
    details::begin_type<T> begin(BaseID const& id) const
    {
      return details::IteratorMaker<T>::create_begin(this, id);
    }

    template <typename T, typename BaseID>
    details::end_type<T> end(BaseID const& id) const
    {
      return details::IteratorMaker<T>::create_end(this, id);
    }

    template <typename T, typename ID>
    details::range_type<T> Iterate(ID const& id) const
    {
      return details::IteratorMaker<T>::create_range(this, id);
    }

    //
    // single object features
    //

    Length_t CryostatHalfWidth(CryostatID const& cid = cryostat_zero) const;

    Length_t CryostatHalfHeight(CryostatID const& cid = cryostat_zero) const;

    Length_t CryostatLength(CryostatID const& cid = cryostat_zero) const;

    //
    // object description
    //


    std::string GetCryostatVolumeName(CryostatID const& cid) const;



    //
    // group features
    //

    unsigned int MaxTPCs() const;

    unsigned int TotalNTPC() const;

    template <typename T>
    TPCDataContainer<T> makeTPCData() const
    {
      return {Ncryostats(), MaxTPCs()};
    }

    template <typename T>
    TPCDataContainer<T> makeTPCData(T const& defValue) const
    {
      return {Ncryostats(), MaxTPCs(), defValue};
    }


    unsigned int NTPC(CryostatID const& cryoid = cryostat_zero) const
    {
      CryostatGeo const* pCryo = GetElementPtr(cryoid);
      return pCryo ? pCryo->NElements() : 0;
    }
    unsigned int NElements(CryostatID const& cryoid) const { return NTPC(cryoid); }
    unsigned int NSiblingElements(TPCID const& tpcid) const { return NTPC(tpcid); }

    //
    // access
    //
    bool HasTPC(TPCID const& tpcid) const
    {
      CryostatGeo const* pCryo = CryostatPtr(tpcid);
      return pCryo ? pCryo->HasTPC(tpcid) : false;
    }

    bool HasElement(TPCID const& tpcid) const { return HasTPC(tpcid); }


    static constexpr TPCID tpc_zero{cryostat_zero, 0};
    TPCGeo const& TPC(TPCID const& tpcid = tpc_zero) const { return Cryostat(tpcid).TPC(tpcid); }
    TPCGeo const& GetElement(TPCID const& tpcid) const { return TPC(tpcid); }


    TPCGeo const* TPCPtr(TPCID const& tpcid) const
    {
      CryostatGeo const* pCryo = CryostatPtr(tpcid);
      return pCryo ? pCryo->TPCPtr(tpcid) : nullptr;
    }
    TPCGeo const* GetElementPtr(TPCID const& tpcid) const { return TPCPtr(tpcid); }


    TPCID FindTPCAtPosition(Point_t const& point) const;

    TPCGeo const* PositionToTPCptr(Point_t const& point) const;


    TPCGeo const& PositionToTPC(Point_t const& point) const;

    TPCID PositionToTPCID(Point_t const& point) const;


    void GetBeginID(TPCID& id) const
    {
      GetBeginID(id.asCryostatID());
      id.TPC = 0;
    }

    void GetEndID(TPCID& id) const;

    bool IncrementID(TPCID& id) const; // inline implementation

    TPCID GetBeginTPCID(CryostatID const& id) const { return {id, 0}; }

    TPCID GetEndTPCID(CryostatID const& id) const;

    //
    // single object features
    //


    Length_t DetHalfWidth(TPCID const& tpcid = tpc_zero) const;


    Length_t DetHalfHeight(TPCID const& tpcid = tpc_zero) const;


    Length_t DetLength(TPCID const& tpcid = tpc_zero) const;


    Point_t GetTPCFrontFaceCenter(TPCID const& tpcid) const
    {
      return TPC(tpcid).GetFrontFaceCenter();
    }

    //
    // object description
    //


    std::string GetLArTPCVolumeName(TPCID const& tpcid = tpc_zero) const;



    //
    // group features
    //

    unsigned int MaxPlanes() const;

    template <typename T>
    PlaneDataContainer<T> makePlaneData() const
    {
      return {Ncryostats(), MaxTPCs(), MaxPlanes()};
    }

    template <typename T>
    PlaneDataContainer<T> makePlaneData(T const& defValue) const
    {
      return {Ncryostats(), MaxTPCs(), MaxPlanes(), defValue};
    }


    unsigned int Nplanes(TPCID const& tpcid = tpc_zero) const
    {
      TPCGeo const* pTPC = GetElementPtr(tpcid);
      return pTPC ? pTPC->NElements() : 0;
    }
    unsigned int NElements(TPCID const& tpcid) const { return Nplanes(tpcid); }
    unsigned int NSiblingElements(PlaneID const& planeid) const { return Nplanes(planeid); }

    unsigned int Nviews() const;

    //
    // access
    //


    bool HasPlane(PlaneID const& planeid) const
    {
      TPCGeo const* pTPC = TPCPtr(planeid);
      return pTPC ? pTPC->HasPlane(planeid) : false;
    }
    bool HasElement(PlaneID const& planeid) const { return HasPlane(planeid); }


    PlaneGeo const& Plane(PlaneID const& planeid) const { return TPC(planeid).Plane(planeid); }
    PlaneGeo const& GetElement(PlaneID const& planeid) const { return Plane(planeid); }


    PlaneGeo const* PlanePtr(PlaneID const& planeid) const
    {
      TPCGeo const* pTPC = TPCPtr(planeid);
      return pTPC ? pTPC->PlanePtr(planeid) : nullptr;
    }
    PlaneGeo const* GetElementPtr(PlaneID const& planeid) const { return PlanePtr(planeid); }

    //
    // iterators
    //

    void GetBeginID(PlaneID& id) const
    {
      GetBeginID(id.asTPCID());
      id.Plane = 0;
    }

    void GetEndID(PlaneID& id) const;

    bool IncrementID(PlaneID& id) const; // inline implementation

    PlaneID GetBeginPlaneID(CryostatID const& id) const { return {GetBeginTPCID(id), 0}; }

    PlaneID GetEndPlaneID(CryostatID const& id) const;

    PlaneID GetBeginPlaneID(TPCID const& id) const { return {id, 0}; }

    PlaneID GetEndPlaneID(TPCID const& id) const;

    //
    // single object features
    //


    Length_t PlanePitch(TPCID const& tpcid,
                        PlaneID::PlaneID_t p1 = 0,
                        PlaneID::PlaneID_t p2 = 1) const;
    Length_t PlanePitch(PlaneID const& pid1, PlaneID const& pid2) const;

    View_t View(PlaneID const& pid) const;

    SigType_t SignalType(PlaneID const& pid) const;



    //
    // group features
    //


    unsigned int Nwires(PlaneID const& planeid) const
    {
      PlaneGeo const* pPlane = GetElementPtr(planeid);
      return pPlane ? pPlane->NElements() : 0;
    }
    unsigned int NElements(PlaneID const& planeid) const { return Nwires(planeid); }
    unsigned int NSiblingElements(WireID const& wireid) const { return Nwires(wireid); }

    unsigned int MaxWires() const;


    //
    // access
    //


    bool HasWire(WireID const& wireid) const
    {
      PlaneGeo const* pPlane = PlanePtr(wireid);
      return pPlane ? pPlane->HasWire(wireid) : false;
    }
    bool HasElement(WireID const& wireid) const { return HasWire(wireid); }


    WireGeo const* WirePtr(WireID const& wireid) const
    {
      PlaneGeo const* pPlane = PlanePtr(wireid);
      return pPlane ? pPlane->WirePtr(wireid) : nullptr;
    } // WirePtr()
    WireGeo const* GetElementPtr(WireID const& wireid) const { return WirePtr(wireid); }


    WireGeo const& Wire(WireID const& wireid) const { return Plane(wireid).Wire(wireid); }
    WireGeo const& WireIDToWireGeo(WireID const& wireid) const { return Wire(wireid); }
    WireGeo const& GetElement(WireID const& wireid) const { return Wire(wireid); }

    //
    // iterators
    //

    void GetBeginID(WireID& id) const
    {
      GetBeginID(id.asPlaneID());
      id.Wire = 0;
    }

    void GetEndID(WireID& id) const;

    bool IncrementID(WireID& id) const; // inline implementation

    WireID GetBeginWireID(CryostatID const& id) const { return {GetBeginPlaneID(id), 0}; }

    WireID GetEndWireID(CryostatID const& id) const;

    WireID GetBeginWireID(TPCID const& id) const { return {PlaneID(id, 0), 0}; }

    WireID GetEndWireID(TPCID const& id) const;

    WireID GetBeginWireID(PlaneID const& id) const { return {id, 0}; }

    WireID GetEndWireID(PlaneID const& id) const;

    //
    // single object features
    //


    static constexpr PlaneID plane_zero{tpc_zero, 0};
    Length_t WirePitch(PlaneID const& planeid = plane_zero) const;

    Length_t WirePitch(View_t view) const;


    double WireAngleToVertical(View_t view, TPCID const& tpcid) const;



    //
    // simple geometry queries
    //

    void WireEndPoints(WireID const& wireid, double* xyzStart, double* xyzEnd) const;


    Segment<Point_t> WireEndPoints(WireID const& wireID) const;


    //
    // closest wire
    //

    WireID NearestWireID(Point_t const& point, PlaneID const& planeid) const;


    Length_t WireCoordinate(Point_t const& pos, PlaneID const& planeid) const;

    //
    // wire intersections
    //


    bool WireIDsIntersect(WireID const& wid1, WireID const& wid2, Point_t& intersection) const;


    bool WireIDsIntersect(WireID const& wid1,
                          WireID const& wid2,
                          WireIDIntersection& widIntersect) const;

    bool IntersectionPoint(WireID const& wid1, WireID const& wid2, double& y, double& z) const;

    PlaneID ThirdPlane(PlaneID const& pid1, PlaneID const& pid2) const;

    double ThirdPlaneSlope(PlaneID const& pid1,
                           double slope1,
                           PlaneID const& pid2,
                           double slope2,
                           PlaneID const& output_plane) const;

    double ThirdPlaneSlope(PlaneID const& pid1,
                           double slope1,
                           PlaneID const& pid2,
                           double slope2) const;


    double ThirdPlaneSlope(PlaneID::PlaneID_t plane1,
                           double slope1,
                           PlaneID::PlaneID_t plane2,
                           double slope2,
                           TPCID const& tpcid) const
    {
      return ThirdPlaneSlope(PlaneID(tpcid, plane1), slope1, PlaneID(tpcid, plane2), slope2);
    }

    double ThirdPlane_dTdW(PlaneID const& pid1,
                           double slope1,
                           PlaneID const& pid2,
                           double slope2,
                           PlaneID const& output_plane) const;

    double ThirdPlane_dTdW(PlaneID const& pid1,
                           double slope1,
                           PlaneID const& pid2,
                           double slope2) const;

    static double ComputeThirdPlaneSlope(double angle1,
                                         double slope1,
                                         double angle2,
                                         double slope2,
                                         double angle_target);

    static double ComputeThirdPlane_dTdW(double angle1,
                                         double pitch1,
                                         double dTdW1,
                                         double angle2,
                                         double pitch2,
                                         double dTdW2,
                                         double angle_target,
                                         double pitch_target);



    //
    // group features
    //

    unsigned int NOpDets() const;

    //
    // access
    //
    OpDetGeo const& OpDetGeoFromOpChannel(unsigned int OpChannel) const;

    OpDetGeo const& OpDetGeoFromOpDet(unsigned int OpDet) const;


    unsigned int GetClosestOpDet(Point_t const& point) const;

    //
    // object description
    //

    std::string OpDetGeoName(CryostatID const& cid = cryostat_zero) const;




    //
    // group features
    //

    unsigned int NAuxDets() const { return AuxDets().size(); }

    unsigned int NAuxDetSensitive(size_t const& aid) const;

    //
    // access
    //

    AuxDetGeo const& AuxDet(unsigned int const ad = 0) const;

    unsigned int FindAuxDetAtPosition(Point_t const& point, double tolerance = 0) const;

    void FindAuxDetSensitiveAtPosition(Point_t const& point,
                                       std::size_t& adg,
                                       std::size_t& sv,
                                       double tolerance = 0) const;

    AuxDetGeo const& PositionToAuxDet(Point_t const& point,
                                      unsigned int& ad,
                                      double tolerance = 0) const;

    const AuxDetSensitiveGeo& PositionToAuxDetSensitive(Point_t const& point,
                                                        size_t& ad,
                                                        size_t& sv,
                                                        double tolerance = 0) const;

    const AuxDetGeo& ChannelToAuxDet(
      std::string const& auxDetName,
      uint32_t const& channel) const; // return the AuxDetGeo for the given detector
                                      // name and channel

    const AuxDetSensitiveGeo& ChannelToAuxDetSensitive(
      std::string const& auxDetName,
      uint32_t const& channel) const; // return the AuxDetSensitiveGeo for the given



    //
    // group features
    //

    unsigned int Nchannels() const;

    unsigned int Nchannels(readout::ROPID const& ropid) const;

    std::vector<raw::ChannelID_t> ChannelsInTPCs() const;
    //
    std::set<View_t> const& Views() const { return allViews; }

    //
    // access
    //

    bool HasChannel(raw::ChannelID_t channel) const;


    raw::ChannelID_t PlaneWireToChannel(WireID const& wireid) const;

    //
    // single object features
    //

    SigType_t SignalType(raw::ChannelID_t const channel) const;

    View_t View(raw::ChannelID_t const channel) const;

    std::vector<WireID> ChannelToWire(raw::ChannelID_t const channel) const;

    readout::ROPID ChannelToROP(raw::ChannelID_t channel) const;

    //
    // geometry queries
    //

    raw::ChannelID_t NearestChannel(Point_t const& worldLoc, PlaneID const& planeid) const;

    bool ChannelsIntersect(raw::ChannelID_t c1, raw::ChannelID_t c2, double& y, double& z) const;



    //
    // group features
    //


    unsigned int NTPCsets(geo::CryostatID const& cryoid) const;
    unsigned int NSiblingElements(readout::TPCsetID const& tpcsetid) const
    {
      return NTPCsets(tpcsetid);
    }

    unsigned int MaxTPCsets() const;

    template <typename T>
    readout::TPCsetDataContainer<T> makeTPCsetData() const
    {
      return {Ncryostats(), MaxTPCsets()};
    }

    template <typename T>
    readout::TPCsetDataContainer<T> makeTPCsetData(T const& defValue) const
    {
      return {Ncryostats(), MaxTPCsets(), defValue};
    }

    //
    // access
    //
    bool HasTPCset(readout::TPCsetID const& tpcsetid) const;

    bool HasElement(readout::TPCsetID const& tpcsetid) const { return HasTPCset(tpcsetid); }

    readout::TPCsetID FindTPCsetAtPosition(Point_t const& worldLoc) const;

    //
    // mapping
    //
    readout::TPCsetID TPCtoTPCset(TPCID const& tpcid) const;

    std::vector<TPCID> TPCsetToTPCs(readout::TPCsetID const& tpcsetid) const;


    void GetBeginID(readout::TPCsetID& id) const
    {
      GetBeginID(id.asCryostatID());
      id.TPCset = 0;
    }

    void GetEndID(readout::TPCsetID& id) const
    {
      GetEndID(id.asCryostatID());
      id.TPCset = 0;
    }

    bool IncrementID(readout::TPCsetID& id) const; // inline implementation

    readout::TPCsetID GetBeginTPCsetID(CryostatID const& id) const { return {id, 0}; }

    readout::TPCsetID GetEndTPCsetID(CryostatID const& id) const { return {id.Cryostat + 1, 0}; }



    //
    // group features
    //


    unsigned int NROPs(readout::TPCsetID const& tpcsetid) const;
    unsigned int NSiblingElements(readout::ROPID const& ropid) const { return NROPs(ropid); }

    unsigned int MaxROPs() const;

    template <typename T>
    readout::ROPDataContainer<T> makeROPdata() const
    {
      return {Ncryostats(), MaxTPCsets(), MaxROPs()};
    }

    template <typename T>
    readout::ROPDataContainer<T> makeROPdata(T const& defValue) const
    {
      return {Ncryostats(), MaxTPCsets(), MaxROPs(), defValue};
    }

    //
    // access
    //
    bool HasROP(readout::ROPID const& ropid) const;

    bool HasElement(readout::ROPID const& ropid) const { return HasROP(ropid); }

    //
    // mapping
    //
    readout::ROPID WirePlaneToROP(PlaneID const& planeid) const;

    std::vector<PlaneID> ROPtoWirePlanes(readout::ROPID const& ropid) const;

    std::vector<TPCID> ROPtoTPCs(readout::ROPID const& ropid) const;

    raw::ChannelID_t FirstChannelInROP(readout::ROPID const& ropid) const;


    void GetBeginID(readout::ROPID& id) const
    {
      GetBeginID(id.asTPCsetID());
      id.ROP = 0;
    }

    void GetEndID(readout::ROPID& id) const
    {
      GetEndID(id.asTPCsetID());
      id.ROP = 0;
    }

    bool IncrementID(readout::ROPID& id) const; // inline implementation

    readout::ROPID GetBeginROPID(CryostatID const& id) const { return {GetBeginTPCsetID(id), 0}; }

    readout::ROPID GetEndROPID(CryostatID const& id) const { return {GetEndTPCsetID(id), 0}; }

    readout::ROPID GetBeginROPID(readout::TPCsetID const& id) const { return {id, 0}; }

    readout::ROPID GetEndROPID(readout::TPCsetID const& id) const { return {GetNextID(id), 0}; }

    View_t View(readout::ROPID const& ropid) const;

    SigType_t SignalType(readout::ROPID const& ropid) const;



    //
    // group features
    //

    unsigned int NOpChannels() const;

    unsigned int MaxOpChannel() const;

    // Number of hardware channels for a given optical detector
    unsigned int NOpHardwareChannels(int opDet) const;

    //
    // access
    //

    bool IsValidOpChannel(int opChannel) const;

    unsigned int OpChannel(int detNum, int hardwareChannel) const;

    unsigned int OpDetFromOpChannel(int opChannel) const;

    unsigned int HardwareChannelFromOpChannel(int opChannel) const;

    unsigned int OpDetFromCryo(unsigned int o, unsigned int c) const;


    //
    // unsorted methods
    //


    void LoadGeometryFile(std::string gdmlfile,
                          std::string rootfile,
                          GeometryBuilder& builder,
                          bool bForceReload = false);

    void LoadGeometryFile(std::string gdmlfile, std::string rootfile, bool bForceReload = false);

    void ApplyChannelMap(std::unique_ptr<ChannelMapAlg> pChannelMap);

  private:
    CryostatList_t& Cryostats() { return fGeoData.cryostats; }
    CryostatList_t const& Cryostats() const { return fGeoData.cryostats; }

    AuxDetList_t& AuxDets() { return fGeoData.auxDets; }
    AuxDetList_t const& AuxDets() const { return fGeoData.auxDets; }

    GeometryData_t fGeoData; 

    double fSurfaceY;          
    std::string fDetectorName; 
    std::string fGDMLfile;     
    std::string fROOTfile;     
    double fMinWireZDist;      
    double fPositionWiggle;    

    fhicl::ParameterSet fBuilderParameters;
    std::unique_ptr<const ChannelMapAlg> fChannelMapAlg;

    // cached values
    std::set<View_t> allViews; 

    std::vector<TGeoNode const*> FindDetectorEnclosure(
      std::string const& name = "volDetEnclosure") const;

    bool FindFirstVolume(std::string const& name, std::vector<const TGeoNode*>& path) const;

    void BuildGeometry(GeometryBuilder& builder);

    bool WireIDIntersectionCheck(const WireID& wid1, const WireID& wid2) const;

    void SortGeometry(GeoObjectSorter const& sorter);

    void UpdateAfterSorting();

    void ClearGeometry();

  }; // class GeometryCore

  // END Geometry group --------------------------------------------------------

} // namespace geo

//******************************************************************************
//*** inline implementation
//***
inline bool geo::GeometryCore::IncrementID(CryostatID& id) const
{
  ++id.Cryostat;
  if (id) id.isValid = HasCryostat(id); // if invalid already, it stays so
  return bool(id);
}

inline bool geo::GeometryCore::IncrementID(TPCID& id) const
{
  unsigned int const nTPCsInCryo = NTPC(id);
  if (++id.TPC < nTPCsInCryo) return bool(id); // if was invalid, it stays so
  // no more TPCs in this cryostat
  id.TPC = 0;
  return IncrementID(id.asCryostatID()); // also sets validity
}

inline bool geo::GeometryCore::IncrementID(PlaneID& id) const
{
  // this implementation is non-optimal, in that the cryostat lookup is
  // performed both here and, potentially, in IncrementID(TPCID)
  unsigned int const nPlanesInTPC = Nplanes(id);
  if (++id.Plane < nPlanesInTPC) return bool(id); // if was invalid, stays so
  // no more planes in this TPCs
  id.Plane = 0;
  return IncrementID(id.asTPCID()); // also sets validity
}

inline bool geo::GeometryCore::IncrementID(WireID& id) const
{
  // this implementation is non-optimal, in that the TPC lookup is
  // performed both here and, potentially, in IncrementID(PlaneID)
  unsigned int const nWiresInPlane = Nwires(id);
  if (++id.Wire < nWiresInPlane) return bool(id); // if was invalid, stays so
  // no more wires in this plane
  id.Wire = 0;
  return IncrementID(id.asPlaneID()); // also sets validity
}

inline bool geo::GeometryCore::IncrementID(readout::TPCsetID& id) const
{
  unsigned int const nTPCsetsInCryo = NTPCsets(id);
  if (++id.TPCset < nTPCsetsInCryo) return bool(id); // if was invalid, it stays so
  // no more TPC sets in this cryostat
  id.TPCset = 0;
  return IncrementID(id.asCryostatID()); // also sets validity
}

inline bool geo::GeometryCore::IncrementID(readout::ROPID& id) const
{
  // this implementation is non-optimal, in that the cryostat lookup is
  // performed both here and, potentially, in IncrementID(TPCsetID)
  unsigned int const nROPinTPC = NROPs(id);
  if (++id.ROP < nROPinTPC) return bool(id); // if was invalid, stays so
  // no more readout planes in this TPC set
  id.ROP = 0;
  return IncrementID(id.asTPCsetID()); // also sets validity
}

//******************************************************************************
//***  template implementation
//***
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
// template member function specializations
namespace geo {

  // TPCID
  template <>
  inline TPCID GeometryCore::GetBeginID<TPCID, CryostatID>(CryostatID const& id) const
  {
    return GetBeginTPCID(id);
  }

  template <>
  inline TPCID GeometryCore::GetEndID<TPCID, CryostatID>(CryostatID const& id) const
  {
    return GetEndTPCID(id);
  }

  // PlaneID
  template <>
  inline PlaneID GeometryCore::GetBeginID<PlaneID, CryostatID>(CryostatID const& id) const
  {
    return GetBeginPlaneID(id);
  }

  template <>
  inline PlaneID GeometryCore::GetBeginID<PlaneID, TPCID>(TPCID const& id) const
  {
    return GetBeginPlaneID(id);
  }

  template <>
  inline PlaneID GeometryCore::GetEndID<PlaneID, CryostatID>(CryostatID const& id) const
  {
    return GetEndPlaneID(id);
  }

  template <>
  inline PlaneID GeometryCore::GetEndID<PlaneID, TPCID>(TPCID const& id) const
  {
    return GetEndPlaneID(id);
  }

  // WireID
  template <>
  inline WireID GeometryCore::GetBeginID<WireID, CryostatID>(CryostatID const& id) const
  {
    return GetBeginWireID(id);
  }

  template <>
  inline WireID GeometryCore::GetBeginID<WireID, TPCID>(TPCID const& id) const
  {
    return GetBeginWireID(id);
  }

  template <>
  inline WireID GeometryCore::GetBeginID<WireID, PlaneID>(PlaneID const& id) const
  {
    return GetBeginWireID(id);
  }

  template <>
  inline WireID GeometryCore::GetEndID<WireID, CryostatID>(CryostatID const& id) const
  {
    return GetEndWireID(id);
  }

  template <>
  inline WireID GeometryCore::GetEndID<WireID, TPCID>(TPCID const& id) const
  {
    return GetEndWireID(id);
  }

  template <>
  inline WireID GeometryCore::GetEndID<WireID, PlaneID>(PlaneID const& id) const
  {
    return GetEndWireID(id);
  }

  // TPCsetID
  template <>
  inline readout::TPCsetID GeometryCore::GetBeginID<readout::TPCsetID, CryostatID>(
    CryostatID const& id) const
  {
    return GetBeginTPCsetID(id);
  }

  template <>
  inline readout::TPCsetID GeometryCore::GetEndID<readout::TPCsetID, CryostatID>(
    CryostatID const& id) const
  {
    return GetEndTPCsetID(id);
  }

  // ROPID
  template <>
  inline readout::ROPID GeometryCore::GetBeginID<readout::ROPID, CryostatID>(
    CryostatID const& id) const
  {
    return GetBeginROPID(id);
  }

  template <>
  inline readout::ROPID GeometryCore::GetEndID<readout::ROPID, CryostatID>(
    CryostatID const& id) const
  {
    return GetEndROPID(id);
  }

  template <>
  inline readout::ROPID GeometryCore::GetBeginID<readout::ROPID, readout::TPCsetID>(
    readout::TPCsetID const& id) const
  {
    return GetBeginROPID(id);
  }

  template <>
  inline readout::ROPID GeometryCore::GetEndID<readout::ROPID, readout::TPCsetID>(
    readout::TPCsetID const& id) const
  {
    return GetEndROPID(id);
  }

} // namespace geo

inline geo::GeometryCore::Segment<geo::Point_t> geo::GeometryCore::WireEndPoints(
  WireID const& wireid) const
{
  WireGeo const& wire = Wire(wireid);
  return {wire.GetStart(), wire.GetEnd()};
}

//------------------------------------------------------------------------------
template <typename Stream>
void geo::GeometryCore::Print(Stream&& out, std::string indent /* = "  " */) const
{

  out << "Detector " << DetectorName() << " has " << Ncryostats() << " cryostats and " << NAuxDets()
      << " auxiliary detectors:";

  auto const& detEnclosureBox = DetectorEnclosureBox();
  out << "\n"
      << indent << "Detector enclosure: " << detEnclosureBox.Min() << " -- "
      << detEnclosureBox.Max() << " cm => ( " << detEnclosureBox.SizeX() << " x "
      << detEnclosureBox.SizeY() << " x " << detEnclosureBox.SizeZ() << " ) cm^3";

  for (auto const& cryostat : Iterate<CryostatGeo>()) {
    out << "\n" << indent;
    cryostat.PrintCryostatInfo(std::forward<Stream>(out), indent + "  ", cryostat.MaxVerbosity);

    const unsigned int nTPCs = cryostat.NTPC();
    for (unsigned int t = 0; t < nTPCs; ++t) {
      const TPCGeo& tpc = cryostat.TPC(t);

      out << "\n" << indent << "  ";
      tpc.PrintTPCInfo(std::forward<Stream>(out), indent + "    ", tpc.MaxVerbosity);

      const unsigned int nPlanes = tpc.Nplanes();
      for (unsigned int p = 0; p < nPlanes; ++p) {
        const PlaneGeo& plane = tpc.Plane(p);
        const unsigned int nWires = plane.Nwires();

        out << "\n" << indent << "    ";
        plane.PrintPlaneInfo(std::forward<Stream>(out), indent + "      ", plane.MaxVerbosity);
        SigType_t const sigType = SignalType(plane.ID());
        out << "\n"
            << indent << "      "
            << "signal type: " << SignalTypeName(sigType) << " (" << static_cast<int>(sigType)
            << ")";

        for (unsigned int w = 0; w < nWires; ++w) {
          const WireGeo& wire = plane.Wire(w);
          WireID wireID(plane.ID(), w);

          // the wire should be aligned on z axis, half on each side of 0,
          // in its local frame
          out << "\n" << indent << "      " << wireID << " ";
          wire.PrintWireInfo(std::forward<Stream>(out), indent + "      ", wire.MaxVerbosity);
        } // for wire
      }   // for plane
    }     // for TPC

    unsigned int nOpDets = cryostat.NOpDet();
    for (unsigned int iOpDet = 0; iOpDet < nOpDets; ++iOpDet) {
      OpDetGeo const& opDet = cryostat.OpDet(iOpDet);
      out << "\n" << indent << "  [OpDet #" << iOpDet << "] ";
      opDet.PrintOpDetInfo(std::forward<Stream>(out), indent + "  ", opDet.MaxVerbosity);
    } // for
  }   // for cryostat

  unsigned int const nAuxDets = NAuxDets();
  for (unsigned int iDet = 0; iDet < nAuxDets; ++iDet) {
    AuxDetGeo const& auxDet = AuxDet(iDet);

    out << "\n" << indent << "[#" << iDet << "] ";
    auxDet.PrintAuxDetInfo(std::forward<Stream>(out), indent + "  ", auxDet.MaxVerbosity);

    unsigned int const nSensitive = auxDet.NSensitiveVolume();
    switch (nSensitive) {
    case 0: break;
    case 1: {
      AuxDetSensitiveGeo const& auxDetS = auxDet.SensitiveVolume(0U);
      out << "\n" << indent << "  ";
      auxDetS.PrintAuxDetInfo(std::forward<Stream>(out), indent + "    ", auxDetS.MaxVerbosity);
      break;
    }
    default:
      for (unsigned int iSens = 0; iSens < nSensitive; ++iSens) {
        out << "\n" << indent << "[#" << iSens << "] ";
        AuxDetSensitiveGeo const& auxDetS = auxDet.SensitiveVolume(iSens);
        auxDetS.PrintAuxDetInfo(std::forward<Stream>(out), indent + "    ", auxDetS.MaxVerbosity);
      } // for
      break;
    } // if sensitive detectors

  } // for auxiliary detector

  out << '\n';

} // geo::GeometryCore::Print()

#endif // LARCOREALG_GEOMETRY_GEOMETRYCORE_H