GeometryBase.h

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \file  GeometryBase.h
/// \brief Encapsulate the geometry of one entire detector (near, far,
///        ndos)
///
/// \version $Id: GeometryBase.h,v 1.6 2012-12-05 02:24:17 bckhouse Exp $
/// \author  messier@indiana.edu
////////////////////////////////////////////////////////////////////////
#ifndef GEO_GEOMETRYBASE_H
#define GEO_GEOMETRYBASE_H

#include "GeometryObjects/CellUniqueId.h"
#include "GeometryObjects/CoordinateTransformation.h"
#include "GeometryObjects/OfflineChan.h"
#include "GeometryObjects/PlaneGeo.h"
#include "NovaDAQConventions/DAQConventions.h"

#include "fhiclcpp/types/Atom.h"
#include "fhiclcpp/types/Table.h"

#include <map>
#include <set>
#include <string>
#include <vector>

#include "TVector3.h"

class TGeoNode;
class TGeoMaterial;
class TGeoManager;

namespace CLHEP{class HepRandomEngine;}

namespace geo {
  class PlaneGeo;
  class CellGeo;
  /// Return codes which flag errors and the like
  enum _return_codes {
    kPLANE_NOT_FOUND = 999000
  };

  typedef struct{
    OfflineChan chan;  ///< Channel passed through by line
    TVector3 entry;         ///< Entry point into cell
    TVector3 exit;          ///< Exit point from cell
  } CellOnLine;
  
  /// The geometry of one entire detector (near, far, ipnd)
  class GeometryBase {
  public:
    template<class T> using Atom = fhicl::Atom<T>;
    template<class T> using Table = fhicl::Table<T>;
    using Comment = fhicl::Comment;
    using Name = fhicl::Name;

    struct DetectorParams
    {
      Atom<bool> BigBoxUsed{Name("BigBoxUsed")};
      Atom<double> BigBoxRange{Name("BigBoxRange"), Comment("in cm")};
      Atom<double> FiducialVolumeXHi{Name("FiducialVolumeXHi")};
      Atom<double> FiducialVolumeXLo{Name("FiducialVolumeXLo")};
      Atom<double> FiducialVolumeYHi{Name("FiducialVolumeYHi")};
      Atom<double> FiducialVolumeYLo{Name("FiducialVolumeYLo")};
      Atom<double> FiducialVolumeZHi{Name("FiducialVolumeZHi")};
      Atom<double> FiducialVolumeZLo{Name("FiducialVolumeZLo")};
    };

    struct Params
    {
      Atom<std::string> GDMLFile{Name("GDML")};
      Atom<bool> ForceUseFCLOnly{Name("ForceUseFCLOnly"), Comment("Force geometry to only use the value set by the fcl file"), false};
      Atom<std::string> StoreTempGeo{Name("StoreTempGeo"), Comment("Folder to store temporary geometry files in")};

      Table<DetectorParams> ndos{Name("ndos")};
      Table<DetectorParams> nd{Name("nd")};
      Table<DetectorParams> fd{Name("fd")};
      Table<DetectorParams> tb{Name("tb")};
    };

    GeometryBase(const Params& params);
    virtual ~GeometryBase();

    // Information about planes and cells
    const PlaneGeo*               Plane(unsigned int i)           const;
    const std::set<unsigned int>& GetPlanesByView(View_t v=kXorY) const;
    const unsigned int            NextPlaneInView(unsigned int p,
                                                  int d=+1)       const;
    const unsigned int            NextPlaneOtherView(unsigned int p,
                                                     int d=+1)    const;
    const unsigned int            FirstPlaneInMuonCatcher()       const;
    double CellTpos(unsigned int ip, unsigned int ic, double w=0.0) const;
    void   CellInfo(unsigned int ip, unsigned int ic,
                    View_t* view=0, double* pos=0, double* dpos=0) const;

    const CellGeo*  IdToCell(const CellUniqueId& id,
                            int* iplane, int* icell) const;

    const PlaneGeo* IdToPlane(const CellUniqueId& id,
                              int* iplane) const;

    /// Same as IdToCell and IdToPlane, but doesn't calculate CellGeo (PlaneGeo),
    /// just the global plane number or cell number
    int  getPlaneID(const CellUniqueId& id) const;
    int  getCellID (const CellUniqueId& id) const;
    bool getPlaneAndCellID(const CellUniqueId& id, int& plane, int& cell) const;
    int  getPlaneID(const double& x, const double& y, const double& z) const                              {return getPlaneID(CellId(x,y,z));}
    int  getCellID (const double& x, const double& y, const double& z) const                              {return getCellID(CellId(x,y,z));}
    bool getPlaneAndCellID(const double& x, const double& y, const double& z, int& plane, int& cell) const{return getPlaneAndCellID(CellId(x,y,z), plane, cell);}

    
    void CountCellsOnLine(double X1, double Y1, double Z1,
                          double X2, double Y2, double Z2,
                          std::vector<OfflineChan>& Xhitsonline,
                          std::vector<OfflineChan>& Yhitsonline);

    void CountCellsOnLine(TVector3 r1, TVector3 r2,
                          std::vector<OfflineChan>& Xhitsonline,
                          std::vector<OfflineChan>& Yhitsonline);

    /// \brief Make list of cells in each view traversed by a line segment and the entry and exit points into each cell
    void CountCellsOnLine(double X1, double Y1, double Z1,
                          double X2, double Y2, double Z2,
                          std::vector<CellOnLine>& Xhitsonline,
                          std::vector<CellOnLine>& Yhitsonline);

    void CountCellsOnLine(TVector3 r1, TVector3 r2,
                          std::vector<CellOnLine>& Xhitsonline,
                          std::vector<CellOnLine>& Yhitsonline);


    /// \brief Make list of cells in each view traversed by a line segment
    ///
    /// Unlike \ref CountCellsOnLine, which uses the full geometry, here we
    /// make some assumptions about cuboidal cells and the effects of
    /// alignment. For this loss in accuracy we gain a lot of speed.
    void CountCellsOnLineFast(double x1, double y1, double z1,
                              double x2, double y2, double z2,
                              std::vector<OfflineChan>& Xhitsonline,
                              std::vector<OfflineChan>& Yhitsonline);

    void CountCellsOnLineFast(TVector3 r1, TVector3 r2,
                              std::vector<OfflineChan>& Xhitsonline,
                              std::vector<OfflineChan>& Yhitsonline);

    /// 2D version
    void CountCellsOnLineFast(geo::View_t view,
                              double v1, double z1,
                              double v2, double z2,
                              std::vector<OfflineChan>& hitsonline);

    // Size of detector, information about "world" frames
    unsigned int NPlanes()       const;
    double       DetHalfWidth()  const;
    double       DetHalfHeight() const;
    double       DetLength()     const;
    double       SurfaceY()      const;

    /// \brief Direction of neutrinos from the NuMI beam (unit vector)
    ///
    /// Because of the varying location of pion decays in the decay pipe, this
    /// can only be an average value in the ND and NDOS.
    TVector3     NuMIBeamDirection() const;

    /// Direction of neutrinos from the Booster beam (unit vector)
    TVector3     BoosterBeamDirection() const;

    ///get distance from local z position in cell to apd in cm, including pigtail
    double DistToElectronics(double localz, const CellGeo& cell) const;

    /// Return length of fiber in cm from end of cell to apd
    double GetPigtail(const CellUniqueId& id) const;

    //  Boxes
    void         WorldBox(double* xlo_cm, double* xhi_cm,
                          double* ylo_cm, double* yhi_cm,
                          double* zlo_cm, double* zhi_cm) const;

    void         DetectorEnclosureBox(double* xlo_cm, double* xhi_cm,
                                      double* ylo_cm, double* yhi_cm,
                                      double* zlo_cm, double* zhi_cm) const;

    void         DetectorBigBox(double* xlo_cm, double* xhi_cm,
                                double* ylo_cm, double* yhi_cm,
                                double* zlo_cm, double* zhi_cm) const;

    bool         IntersectsDetector(double* xyz_cm, double* dxyz) const;
    bool         IntersectsBigBox(double* xyz_cm, double* dxyz) const;
    // DetectorBigBox stuff
    /// Do we use the Detector Big Box
    bool         isDetectorBigBoxUsed() const {return fIsDetectorBigBoxUsed;}
    /// Is the particle inside the detector Big Box?
    bool         isInsideDetectorBigBox(const double x_cm, const double y_cm, const double z_cm) const;

    /// Is the particle inside the detector Fiducial Volume?
    bool         isInsideFiducialVolume(const double x_cm, const double y_cm, const double z_cm) const;
    bool         isInsideFiducialVolume(const double* xyz_dimensionless) const;

    void FiducialBox(TVector3& r0, TVector3& r1) const;

    // Useful for tracking and swimming through a TGeometry
    int                 CurrentCell(int* ip, int* ic)           const;
    const CellUniqueId  CurrentCellId()                         const;
    const CellUniqueId  CellId(const double& x, const double& y, const double& z,
                               double dxds=0., double dyds=0.,
                               double dzds=1.,double step=0.01) const;
    TGeoMaterial* Material(double x, double y, double z)        const;
    bool          IsActive(const TGeoMaterial* m)               const;
    double        TotalMass(const char *volume="vDetEnclosure") const;
    double        MassBetweenPoints(double *p1, double *p2)     const;
    
    void          MaterialsBetweenPoints(const double *p1,
                                         const double *p2,
                                         std::vector<double>& ds,
                                         std::vector<const TGeoMaterial*>& mat) const;

    void          MaterialsBetweenPoints(TVector3 v1,
                                         TVector3 v2,
                                         std::vector<double>& ds,
                                         std::vector<const TGeoMaterial*>& mat) const;
    
    bool          calculateMassesLong(const unsigned int number_of_points,
                                      CLHEP::HepRandomEngine& engine) const;
    double        DEdge(const double* x0, const double* dx,
                        double *exitp=0)                        const;
    double        ClosestApproach(unsigned int    ip,
                                  unsigned int    ic,
                                  const TVector3& vtx,
                                  const TVector3& dir,
                                  double          dx=0.0,
                                  double          dy=0.0,
                                  double*         w=0,
                                  double*         s=0,
                                  TVector3*       point_in_cell=0,
                                  TVector3*       point_on_track=0) const;
    std::string   GDMLFile()       const { return fGDMLFile;  }
    std::string   ROOTFile()       const { return fROOTFile; }
    std::string   GDMLFromFCL()    const { return fGDMLFromFCL; }
    std::string   FileBaseName()   const;
    int           Run()            const { return fRunNumber; }
    TGeoManager*  ROOTGeoManager() const;

    /// Prefer \ref ds::DetectorService::DetId() instead
    novadaq::cnv::DetId DetId() const { return fDetId; }

    CoordinateTransformation* getCoordinateTransformation() const { return (CoordinateTransformation*)&fCoordinateTransformation;}
    /// Extract contents from \ref fGDMLFile and return as a string
    std::string ExtractGDML() const;
    /// \brief Extract contents from \a fname and return as a string
    ///
    /// \param fullpath If false, search for the file in FW_SEARCH_PATH
    static std::string ExtractGDML(std::string fname, bool fullpath);

  protected:

    typedef std::unordered_map<CellUniqueId, OfflineChan> UIDMap;
    geo::OfflineChan getPlaneCellMap(const CellUniqueId& id) const;

    /// Method to set DetectorID
    virtual void setDetectorID(novadaq::cnv::DetId);

    void       setDetectorBigBox(double detector_bigbox_range);

    void FindPlanes(std::vector<const TGeoNode*>& n, unsigned int depth, bool inMuonCatcher = false);
    void MakePlane (std::vector<const TGeoNode*>& n, unsigned int depth, bool inMuonCatcher);
    void SetDrawOptions();
    void BuildMaps();
    void SetDetectorSize();

    /// In case of unknown ID, guesses from the filename
    bool LoadGeometryFile(std::string gdmlfile, novadaq::cnv::DetId det_id = novadaq::cnv::kUNKNOWN_DET);
    
    int MakeTmpFile(std::string gdmlInfo);
    void RemoveTmpFile(std::string fileName);

    /// Search for Geometry/gdml/fname in FW_SEARCH_PATH, return full path
    static std::string FindGDMLFile(std::string fname);

  protected:
    Params fParams;

    std::string                  fGDMLFile;    ///< gdml file holding the geometry
    std::string                  fROOTFile;    ///< root file holding the geometry
    std::string                  fGDMLFromFCL; ///< keep track of original fcl parameter

    novadaq::cnv::DetId          fDetId;           ///< id for the detector being used


    /// Coordinate Transformation class
    CoordinateTransformation fCoordinateTransformation;

    std::vector<PlaneGeo*> fPlanes;      ///< The detector planes
    std::vector<PlaneGeo*> fPlanesInMuonCatcher; ///< Same pointers as fPlanes
    UIDMap                 fIdMap;       ///< Unique ID -> Plane,Cell

    std::set<unsigned int> fAllPlanes;   ///< List of all planes
    std::set<unsigned int> fXplanes;     ///< List of X measuring planes
    std::set<unsigned int> fYplanes;     ///< List of Y measuring planes

    double fDetLength;                   ///< Detector length (cm)
    double fDetHalfHeight;               ///< Detector 1/2 height (cm)
    double fDetHalfWidth;                ///< Detector 1/2 width  (cm)

    // Detector Big Box is the box surrounding the detector enclosure
    // It is usually defined by adding 15m to each if the enclosure dimension,
    // but could be redefined in the detectorbigbox.xml
    // This is used in Geant - if the particle leaves it, stop tracking it.
    // As well as in the Cosmic generator - if the cosmic ray particle doesn't intersect it,
    // don't put it in the stack.
    // DetectorBigBox dimensions in cm
    bool                   fIsDetectorBigBoxUsed;  ///< Do we need to use the BigBox cut?
    double                 fBigBoxRange;           ///< Range of big box

    /// Dimensions of the DetectorBigBox in cm
    double                 fDetectorBigBoxXHi;
    double                 fDetectorBigBoxXLo;
    double                 fDetectorBigBoxYHi;
    double                 fDetectorBigBoxYLo;
    double                 fDetectorBigBoxZHi;
    double                 fDetectorBigBoxZLo;

    fhicl::ParameterSet fConfigNDOS, fConfigND, fConfigFD, fConfigTB;

    double                 fFiducialVolumeXHi;
    double                 fFiducialVolumeXLo;
    double                 fFiducialVolumeYHi;
    double                 fFiducialVolumeYLo;
    double                 fFiducialVolumeZHi;
    double                 fFiducialVolumeZLo;
    
    std::string            fStoreTempGeo; ///< folder to store temp geo file
    
    unsigned int           fRunNumber;             ///< Run number of configuration
  };
}
#endif // GEO_GEOMETRYBASE_H
////////////////////////////////////////////////////////////////////////