OpDetGeo.h

Go to the documentation of this file.

#ifndef LARCOREALG_GEOMETRY_OPDETGEO_H
#define LARCOREALG_GEOMETRY_OPDETGEO_H

// LArSoft libraries
#include "larcorealg/CoreUtils/RealComparisons.h"
#include "larcorealg/Geometry/LocalTransformationGeo.h"
#include "larcorealg/Geometry/TransformationMatrix.h"
#include "larcoreobj/SimpleTypesAndConstants/geo_optical_vectors.h"
#include "larcoreobj/SimpleTypesAndConstants/geo_types.h" // geo::OpDetID
#include "larcoreobj/SimpleTypesAndConstants/geo_vectors.h"

// ROOT libraries
#include "TClass.h"
#include "TGeoBBox.h"
#include "TGeoMatrix.h" // TGeoHMatrix
#include "TGeoSphere.h"
#include "TGeoTube.h"

// C/C++ standard libraries
#include <algorithm> // std::minmax()
#include <array>
#include <cassert>
#include <string>
#include <type_traits> // std::decay_t(), std::is_base_of_v
#include <typeinfo>    // typeid()
#include <vector>

// forward declarations
class TGeoNode;

namespace geo {

  class OpDetGeo {
  public:

    using LocalPoint_t = geo::OpticalPoint_t;

    using LocalVector_t = geo::OpticalVector_t;


    OpDetGeo(TGeoNode const& node, geo::TransformationMatrix&& trans);

    geo::OpDetID const& ID() const { return fID; }

    geo::Point_t const& GetCenter() const { return fCenter; }
    double RMin() const;
    double RMax() const;
    double HalfL() const;
    double HalfW() const;
    double HalfH() const;
    double Length() const { return 2.0 * HalfL(); }
    double Width() const { return 2.0 * HalfW(); }
    double Height() const { return 2.0 * HalfH(); }
    double ThetaZ() const;             
    double ThetaZ(bool degrees) const; 

    double CosThetaFromNormal(geo::Point_t const& point) const;

    double DistanceToPoint(geo::Point_t const& point) const;


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


    const TGeoNode* Node() const { return fOpDetNode; }

    // --- BEGIN -- detector shape ---------------------------------------------

    TGeoShape const* Shape() const { return Node()->GetVolume()->GetShape(); }

    template <typename ShapeObj>
    bool isShape() const;

    template <typename ShapeObj>
    bool isShapeLike() const;

    bool isTube() const { return isShapeLike<TGeoTube>(); }

    bool isBar() const { return isShape<TGeoBBox>(); }

    bool isSphere() const { return isShape<TGeoSphere>(); }

    // --- END -- detector shape -----------------------------------------------

    void UpdateAfterSorting(geo::OpDetID opdetid);

    template <typename Stream>
    void PrintOpDetInfo(Stream&& out, std::string indent = "", unsigned int verbosity = 0) const;

    std::string OpDetInfo(std::string indent = "", unsigned int verbosity = 0) const;

    static constexpr unsigned int MaxVerbosity = 2;

  private:
    using LocalTransformation_t =
      geo::LocalTransformationGeo<ROOT::Math::Transform3D, LocalPoint_t, LocalVector_t>;

    LocalTransformation_t fTrans; 
    const TGeoNode* fOpDetNode;   
    geo::Point_t fCenter;         

    geo::OpDetID fID; 

    TGeoTube const* asTube() const { return dynamic_cast<TGeoTube const*>(Shape()); }

    TGeoSphere const* asSphere() const { return dynamic_cast<TGeoSphere const*>(Shape()); }

    TGeoBBox const* asBox() const { return dynamic_cast<TGeoBBox const*>(Shape()); }

  }; // class OpDetGeo

} // namespace geo

//------------------------------------------------------------------------------
//--- template implementation
//---

template <typename ShapeObj>
bool geo::OpDetGeo::isShape() const
{
  static_assert(std::is_base_of_v<TGeoShape, std::decay_t<ShapeObj>>);

  // C++ understanding of the business instead of ROOT's (no strong reason)
  TGeoShape const* shape = Shape(); // needed to convince Clang 7 I really mean it
  return typeid(*shape) == typeid(std::decay_t<ShapeObj>);

} // geo::OpDetGeo::isShape()

//------------------------------------------------------------------------------
template <typename ShapeObj>
bool geo::OpDetGeo::isShapeLike() const
{
  static_assert(std::is_base_of_v<TGeoShape, std::decay_t<ShapeObj>>);

  // C++ understanding of the business instead of ROOT's (no strong reason)
  return dynamic_cast<std::decay_t<ShapeObj> const*>(Shape()) != nullptr;

} // geo::OpDetGeo::isShapeLike()

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

  lar::util::RealComparisons<double> cmp(1e-5);

  //----------------------------------------------------------------------------
  out << "optical detector " << ID() << " centered at " << GetCenter() << " cm";

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

  //----------------------------------------------------------------------------
  if (isTube()) {
    out << ", radius: " << RMax() << " cm";
    if (cmp.nonZero(RMin())) out << " (inner: " << RMin() << " cm)";
    out << ", length: " << Length() << " cm";
  }
  else if (isBar()) {
    out << ", bar size " << Width() << " x " << Height() << " x " << Length() << " cm";
  }
  else if (TGeoSphere const* sphere = asSphere(); sphere) {
    assert(isSphere());
    auto const [th1, th2] = std::minmax({sphere->GetTheta1(), sphere->GetTheta2()});
    out << ", ";
    // some information out of the interface
    if (cmp.zero(th1) && cmp.equal(th2, 180.0))
      out << "spherical";
    else if ((cmp.zero(th1) && cmp.equal(th2, 90.0)) ||
             (cmp.equal(th1, 90.0) && cmp.equal(th2, 180.0))) {
      out << "hemispherical";
    }
    else
      out << "spherical portion (" << th1 << " -> " << th2 << " degree)";
    out << " with external radius " << RMax() << " cm";
  }
  else
    out << ", shape: '" << Shape()->IsA()->GetName() << "'";

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

  //----------------------------------------------------------------------------
  out << ", theta(z): " << ThetaZ() << " rad";

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

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

} // geo::OpDetGeo::PrintOpDetInfo()

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

#endif // LARCOREALG_GEOMETRY_OPDETGEO_H