GeometryBase.cxx

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \file    GeometryBase.cxx
/// \brief
/// \version $Id: GeometryBase.cxx,v 1.12 2012-12-19 05:34:46 gsdavies Exp $
/// \author  messier@indiana.edu
////////////////////////////////////////////////////////////////////////

#include "GeometryObjects/GeometryBase.h"

#include <sys/stat.h>
#include <sys/types.h>

#include <cassert>
#include <iostream>
#include <stdio.h>

//ROOT includes
#include "TGeoManager.h"
#include "TGeoNode.h"
#include "TStopwatch.h"
#include "TVirtualGeoPainter.h"

// Framework includes
#include "cetlib_except/exception.h"
#include "cetlib/search_path.h"
#include "messagefacility/MessageLogger/MessageLogger.h"
#include "CLHEP/Random/RandFlat.h"

// NOvA includes
#include "GeometryObjects/CellGeo.h"
#include "GeometryObjects/Geo.h"
#include "GeometryObjects/PlaneGeo.h"
#include "Utilities/func/MathUtil.h"

namespace geo
{
  //....................................................................
  static bool plane_sort(const PlaneGeo* p1, const PlaneGeo* p2)
  {
    double xyz1[3], xyz2[3];
    p1->Cell(0)->GetCenter(xyz1);
    p2->Cell(0)->GetCenter(xyz2);
    return xyz1[2]<xyz2[2];
  }

  //......................................................................
  GeometryBase::GeometryBase(const Params& params)
    : fParams(params)
    , fDetId(novadaq::cnv::kUNKNOWN_DET)
    , fRunNumber(0)
  {
    setDetectorID(novadaq::cnv::kUNKNOWN_DET);

    // If they want to force FCL usage only, they better have supplied a file
    assert(!(fParams.ForceUseFCLOnly() && fParams.GDMLFile().empty()));
    if(fParams.GDMLFile().empty()){
      // That's fine. Normal operating mode. We'll get GDML information in
      // beginRun().
      return;
    }

    //NOTE: Setting fROOTFile to the same thing as GDMLFile, in the future
    //this should be removed but for now EventGeneratorBase has a print
    //statement requiring the root file name
    fROOTFile = fParams.GDMLFile();

    fGDMLFromFCL = fParams.GDMLFile();

    // constructor decides if initialized value is a path or an environment variable
    cet::search_path sp("FW_SEARCH_PATH");
  
    if ( !sp.find_file(fParams.GDMLFile(), fGDMLFile) ) {
      throw cet::exception("GeometryFileLoad")
        << "cannot find GDML file " << fGDMLFile << " bail ungracefully\n"
        << __FILE__ << ":" << __LINE__ << "\n";
    }
    mf::LogWarning("LoadNewGeometry") << "loading new geometry files\n"
                                      << fGDMLFile << "\n";
  
    this->LoadGeometryFile(fGDMLFile);
  }

  //......................................................................
  GeometryBase::~GeometryBase()
  {
    for (unsigned int i=0; i<fPlanes.size(); ++i) {
      if (fPlanes[i]) { delete fPlanes[i]; fPlanes[i] = 0; }
    }

    //clean up an existing temporary file
    if (fGDMLFile.find("TEMP_GEO_GDML") != std::string::npos){
      this->RemoveTmpFile(fGDMLFile);
    }
    

  }
  
  //....................................................................
  ///
  /// Method to make a temporary file holding Detector geometry (gdml) file
  ///
  int GeometryBase::MakeTmpFile(std::string gdmlInfo)
  {
    int fStat = -1; //file status variable, -1 indicates error
    FILE *tmpFile; //file variable
    
    //hold temporary gdml file name
    std::string tmpFileName;
    
    //Use mkstemp to make a temporary file, mkstemp takes in a template
    //character array as an argument.  The last 6 characters must be XXXXXX
    //mkstemp will then attempt to open a file with that template and the 
    //last 6 characters will be changed to be unique.  If the file cannot
    //be successfully opened an error code will be returned.
    //If a file is properly opened it will be unique.  File is set to be
    //opened in /temp.
    char tmp[256];
    sprintf(tmp,"%s/TEMP_GEO_GDML_XXXXXX", fParams.StoreTempGeo().c_str());
    //check to see if file opens properly
    if ( ((fStat = mkstemp(tmp)) == -1) ||
           ((tmpFile = fopen(tmp, "w+")) == NULL) ) {
      
      if (fStat != -1){
        fclose(tmpFile);
        remove(tmp);
        std::cerr<<"Error making temporary file "
            <<tmp<<std::endl;
      }
      //if file could not be created properly exit and return -1
      std::cerr<<"Failed to produce temporary geometry file"<<std::endl;
      return -1;
    }
    
    //this means a unique file has been properly created and opened
    fclose(tmpFile); //properly close file
    
    //add temp file name to string with proper .gdml extension
    tmpFileName += tmp;
    tmpFileName += ".gdml";
    
    //now attempt to open new file with .gdml extension
    FILE *gdmlFile;
    gdmlFile = fopen(tmpFileName.c_str(), "w+");
    if (gdmlFile == NULL){
      //file was not opened properly, remove temp files
      //and exit with a status -1
      remove(tmp);
      remove(tmpFileName.c_str());
      std::cerr<<"Error making temporary file "
               <<tmpFileName<<std::endl;
      return -1;
    }
    
    //temporary gdml file was properly opened, now put our gdml info in file
    fputs (gdmlInfo.c_str(), gdmlFile);
    
    //now remove initial temp file without the .gdml extension
    remove(tmp);
    
    //close temporary gdml file
    fclose(gdmlFile);
    
    //now store name of temporary gdml file in GDMLFile
    //this way the temp file can also be accessed by GEANT4,
    //which loads the geometry directly
    fGDMLFile = tmpFileName;
    
    //return 0 indicating successful temporary file creation
    return 0;
  }
  
  //......................................................................
  ///
  /// Method to remove temporary file that holds detector geometry file
  void GeometryBase::RemoveTmpFile(std::string fileName){
    remove(fileName.c_str());
  }

  //......................................................................
  std::string GeometryBase::ExtractGDML() const
  {
    return ExtractGDML(fGDMLFile, true);
  }
  
  //......................................................................
  std::string GeometryBase::ExtractGDML(std::string fname, bool fullpath)
  {
    if(!fullpath) fname = FindGDMLFile(fname);

    std::string gdmlInfo;

    std::ifstream geoFile(fname.c_str());
    std::string line;
    
    if (geoFile.is_open()){
      while (geoFile.good()){
        getline(geoFile, line);
        gdmlInfo += line;
        gdmlInfo += '\n';
      }
      geoFile.close();
    }
    
   return gdmlInfo;
  }

  //......................................................................
  std::string GeometryBase::FindGDMLFile(std::string fname)
  {
    cet::search_path sp("FW_SEARCH_PATH");
    size_t lastSlash   = fname.find_last_of("/");
    std::string basename = fname.substr(lastSlash+1);

    fname = "Geometry/gdml/"+basename;

    std::string ret;
    if(!sp.find_file(fname, ret)){
      throw cet::exception("FindGDMLFile")
        << "cannot find GDML file " << fname << " bail ungracefully\n"
        << __FILE__ << ":" << __LINE__ << "\n";
    }
    return ret;
  }

  //......................................................................
  ///
  /// Method to load detector geometry file
  ///
  bool GeometryBase::LoadGeometryFile(std::string gdmlfile,
                                      novadaq::cnv::DetId detId) {

    int old_verbosity = gGeoManager->GetVerboseLevel();
    
    // TGeoManager is too verbose when loading geometry.
    // Make it quiet.
    gGeoManager->SetVerboseLevel(0);

    // Open the GDML file
    struct stat sb;
    if (gdmlfile.empty() || stat(gdmlfile.c_str(), &sb)!=0) {
      // failed to resolve the file name
      throw cet::exception("GeometryFileLoad")
        << "empty GDML string " << gdmlfile << " bail ungracefully\n"
        << __FILE__ << ":" << __LINE__ << "\n";
    }
    // clear the plane array and id map
    for (unsigned int i=0; i<fPlanes.size(); ++i) {
      if (fPlanes[i]) { delete fPlanes[i]; fPlanes[i] = 0; }
    }
    fPlanes.clear();
    fPlanesInMuonCatcher.clear();
    fIdMap.clear();
    CellGeo::clear();

    TGeoManager::Import(gdmlfile.c_str());

    this->SetDrawOptions();

    std::vector<const TGeoNode*> n(16);
    n[0] = gGeoManager->GetTopNode();
    this->FindPlanes(n, 0);
    sort(fPlanes.begin(), fPlanes.end(), plane_sort);
    sort(fPlanesInMuonCatcher.begin(), fPlanesInMuonCatcher.end(), plane_sort);

    this->SetDetectorSize();
    this->BuildMaps();

    if(detId == novadaq::cnv::kUNKNOWN_DET){
      // have to set the fDetId based on the input file then
      if(     gdmlfile.find("ndos") != std::string::npos) detId = novadaq::cnv::kNDOS;
      else if(gdmlfile.find("near") != std::string::npos) detId = novadaq::cnv::kNEARDET;
      else if(gdmlfile.find("far")  != std::string::npos) detId = novadaq::cnv::kFARDET;
      else if(gdmlfile.find("testbeam-1x1-2block") != std::string::npos) detId = novadaq::cnv::kTESTBEAM;
      else if(gdmlfile.find("testbeam-1x1-3block") != std::string::npos) detId = novadaq::cnv::kTESTBEAM;
      else if(gdmlfile.find("testbeam-2x2-2block") != std::string::npos) detId = novadaq::cnv::kTESTBEAM;
      else if(gdmlfile.find("testbeam-2x2-3block") != std::string::npos) detId = novadaq::cnv::kTESTBEAM;
    }
    setDetectorID(detId);

    if(fIsDetectorBigBoxUsed)
      this->setDetectorBigBox(fBigBoxRange);
    else{
      fDetectorBigBoxXHi =  fDetHalfWidth;
      fDetectorBigBoxXLo = -fDetHalfWidth;
      fDetectorBigBoxYHi =  fDetHalfHeight;
      fDetectorBigBoxYLo = -fDetHalfHeight;
      fDetectorBigBoxZHi =  fDetLength;
      fDetectorBigBoxZLo =  0.;
    }

    mf::LogInfo("LoadGeometryFile")
      << "Detector geometry loaded from " << gdmlfile;

    
    // Restore verbosity level before we leave because
    // we have to be polite and listen
    gGeoManager->SetVerboseLevel(old_verbosity);
    
    return true;
  }// end of LoadGeometryFile


  //......................................................................
  /// 
  /// Method to extract gdml file basename and store as an std::string
  ///
  std::string GeometryBase::FileBaseName() const {
    assert(fGDMLFile.find(".gdml") != std::string::npos);

    // to find the base name of the file
    // use string operations to find the last
    // instance of "/" in the fROOTFile string
    // and the last instance of ".", the base
    // string should be between those
    size_t lastSlash   = fGDMLFile.find_last_of("/");
    std::string almost = fGDMLFile.substr(lastSlash+1);

    // resize almost to drop the file ".gdml"
    almost.resize(almost.length() - 5);

    return almost;
  }// end of FileBaseName

  //......................................................................
  ///
  /// Return the plane and cell number for the current tracking position
  /// held by the gGeoManager
  ///
  /// \param planeid  : pointer to plane number, >=0 on success, -1 on failure
  /// \param cellid   : pointer to cell number, >=0 on success, -1 on failure
  /// \returns >0 if cell is found, <0 if not found.
  ///
  int GeometryBase::CurrentCell(int* planeid, int* cellid) const {
    if (! this->IdToCell(this->CurrentCellId(), planeid, cellid)) {
      *planeid = -1;
      *cellid  = -1;
      throw cet::exception("UNKNOWN_ID")
        << "unknown id " << this->CurrentCellId() << "\n"
        << __FILE__ << ":" << __LINE__ << "\n";
    }
    return 1;
  }// end of

  //......................................................................
  ///
  /// Return the unique cell id for the current cell -- zero if not in any cell
  ///
  const CellUniqueId GeometryBase::CurrentCellId() const {
    const std::string path = gGeoManager->GetPath();
    //std::cout << "path to cell is " << path << std::endl;

    // check that the path actually contains a cell
    if(path.find("Cell") == std::string::npos){
      MF_LOG_DEBUG("Geometry") << "no Cell in the path to the current volume:" << path << std::endl;
      return 0;
    }
    if(path.find("HAirBubble") != std::string::npos){
      MF_LOG_DEBUG("Geometry") << "This an air bubble, not a cell:" << path << std::endl;
      return 0;
    }

    // We used to return PathToUniqueId(path.c_str()) at this point, but as of
    // the root6 upgrade, the paths are all like
    // vPlaneV_0/vModuleV_0/vExtruV_0/vCellV_0 rather than
    // vPlaneV_3/vModuleV_1/vExtruV_4/vCellV_1, so use the list of parent
    // pointers to determine the unique ID instead. The code here is derived
    // from what TGeoManager::GetPath() does.
    TGeoNodeCache* cache = gGeoManager->GetCurrentNavigator()->GetCache();
    std::vector<const TGeoNode*> ns;
    for(int i = 0; i <= cache->GetLevel(); ++i){
      ns.push_back(((TGeoNode**)cache->GetBranch())[i]);
    }
    return NodesToUniqueId(ns, ns.size()-1);
  }// end of CurrentCellId

  //......................................................................
  ///
  /// Return the unique cel id for a specified location and direction in the detector
  ///
  /// \param x - x position coordinate (cm, world coordinates)
  /// \param y - y position coordinate (cm, world coordinates)
  /// \param z - z position coordinate (cm, world coordinates)
  /// \param dxds - dxds direction in which to explore a step
  /// \param dyds - dyds direction in which to explore a step
  /// \param dzds - dzds direction in which to explore a step
  /// \param step - step size to take
  ///
  /// For the sake of efficiency, this method does no checking to ensure
  /// that the location is in fact inside a cell. The user is expected
  /// to perform that check before using the position (x,y,z).
  ///
  /// Note further, that this uses the gGeoManager->InitTrack method and
  /// will disrupt TGeo-based tracking and should be called once any
  /// TGeo-based tracking is completed.
  ///
  const CellUniqueId GeometryBase::CellId(const double& x, const double& y, const double& z,
                                          double dxds, double dyds, double dzds,
                                          double step) const {
    // The strategy for finding the cell ID is to try the exact point
    // we've been handed. If that does not return a valid cell ID, then
    // try shifting the point around slightly to account for tracking
    // and edge finding uncertainties. If after trying that we still
    // don't have a valid ID, return 0.

    gGeoManager->InitTrack(x, y, z, dxds, dyds, dzds);
    gGeoManager->FindNextBoundary(1.E-6);
    CellUniqueId ID  = this->CurrentCellId();
    if (ID != 0) return ID;

    double       eps = 1.0E-1;
    const double s   = std::max(step, 1.e-3);

    // walk the step points backwards and forwards to see if we can find the
    // cell
    while(eps < 1.){

      const double eps_s = eps*s;

      //     std::cout << "step back pos " << x-dxds*eps*s
      //               << "," << y-dyds*eps*s << "," << z-dzds*eps*s
      //               << " " << eps << " " << step << std::endl;
      gGeoManager->InitTrack(x-dxds*eps_s, y-dyds*eps_s, z-dzds*eps_s, dxds, dyds, dzds);
      gGeoManager->FindNextBoundary(1.E-6);
      ID = this->CurrentCellId();
      if (ID != 0) return ID;

      //     std::cout << "step forward pos " << x+dxds*eps*s
      //               << "," << y+dyds*eps*s << "," << z+dzds*eps*s << std::endl;
      gGeoManager->InitTrack(x+dxds*eps_s, y+dyds*eps_s, z+dzds*eps_s, dxds, dyds, dzds);
      ID = this->CurrentCellId();
      if (ID != 0) return ID;

      eps += 1.e-1;
    }

    // No valid ID found
    return 0;
  }// end of CellId

  //......................................................................
  ///
  /// Method to obtain current detector material at given (x,y,z) coordinates
  ///
  TGeoMaterial* GeometryBase::Material(double x, double y, double z) const {
    const TGeoNode* node = gGeoManager->FindNode(x,y,z);
    return node->GetMedium()->GetMaterial();
  }// end of Material

  //......................................................................
  ///
  /// Is the material m an active or passive material?
  ///
  bool GeometryBase::IsActive(const TGeoMaterial* m) const {
    const std::string scint = "Scintillator";
    if (scint == m->GetName()) return true;
    return false;
  }

  //......................................................................
  unsigned int GeometryBase::NPlanes() const
  {
    assert(fDetId != novadaq::cnv::kUNKNOWN_DET);
    return fPlanes.size();
  }

  //......................................................................
  double GeometryBase::DetHalfWidth() const
  {
    assert(fDetId != novadaq::cnv::kUNKNOWN_DET);
    return fDetHalfWidth;
  }

  //......................................................................
  double GeometryBase::DetHalfHeight() const
  {
    assert(fDetId != novadaq::cnv::kUNKNOWN_DET);
    return fDetHalfHeight;
  }

  //......................................................................
  double GeometryBase::DetLength() const
  {
    assert(fDetId != novadaq::cnv::kUNKNOWN_DET);
    return fDetLength;
  }

  //......................................................................

  ///
  /// A typical y-position value at the surface (where earth meets air)
  /// for this detector site
  ///
  /// \returns typical y position at surface in units of cm
  ///
  /// \throws geo::Exception if detector ID is not set properly
  ///
  /// \todo: Only far det number is reasonable accurate. Near and IPND
  /// need some thought
  ///
  double GeometryBase::SurfaceY() const {

    switch(fDetId){
    case novadaq::cnv::kFARDET:
      return  9.94E2;  // far : surface ~10m  above det. center

    case novadaq::cnv::kNDOS:
      return -1.966E2; // NDOS: surface ~2m   below det. center

    case novadaq::cnv::kNEARDET:
      return 130.0E2;  // near: surface ~130m above det. center

    case novadaq::cnv::kTESTBEAM:
      return -1.966E2; // NDOS: surface ~2m   below det. center

    default:
      throw cet::exception("BadGeoConfig")
        << "Bad geometry configuration: " << (int)fDetId << "\n"
        << __FILE__ << ":" << __LINE__ << "\n";
    }
  }// end of SurfaceY

  //---------------------------------------------------------------------------
  TVector3 GeometryBase::NuMIBeamDirection() const
  {
    // These are obtained by plotting rec.mc.nu.p.fP.fZ/rec.mc.nu.p.fE etc from
    // CAF files and fitting the resulting distribution to a gaussian.
    // A true vertex inside of the detector is require as well as the NuE group FA reconstruction
    // and containment cuts.
    // CAF Files used:
    //   NearDet:
    //    prod_caf_S15-05-22a_nd_genie_fhc_nonswap_ndnewpos
    //   FarDet:
    //    prod_caf_S15-05-22a_fd_genie_fhc_nonswap_fdfluxv08 
    //   Unknown for NDOS

    switch(fDetId){
    case novadaq::cnv::kFARDET:   return TVector3(-6.833e-05,
                                                 +6.388e-02,
                                                 +9.980e-1).Unit();
    case novadaq::cnv::kNDOS:     return TVector3(-3.845e-4,
                                                 +6.517e-2,
                                                 +9.967e-1).Unit();
    case novadaq::cnv::kNEARDET:  return TVector3(-8.424e-04,
                                                 -6.174e-02,
                                                 +9.981e-1).Unit();
    case novadaq::cnv::kTESTBEAM: return TVector3(0,
                                                  0,
                                                  1).Unit();
    default:
      throw cet::exception("BadGeoConfig")
        << "Bad geometry configuration: " << fDetId << "\n"
        << __FILE__ << ":" << __LINE__ << "\n";
    }
  }

  //---------------------------------------------------------------------------
  TVector3 GeometryBase::BoosterBeamDirection() const
  {
    // No simulation available for NDOS Booster. These values copied from the Cana
    // package, where the comment read:
    //
    // NDOS BNB direction from Zelimir Djurcic <zdjurcic@hep.anl.gov> 2/26/2011,
    // private communication.
    //
    // ND BNB direction from Ryan Murphy <rwmurphy@indiana.edu> 10/14/2016
    // Used same method as described in NuMIBeamDirection()
    // CAF File used:
    //   NearDet:
    //   /nova/ana/users/rwmurphy/bnb_mc_all_hadd_caf.caf.root

    switch(fDetId){
    case novadaq::cnv::kNDOS: return TVector3(-0.367, 0.012, 0.930).Unit();

    case novadaq::cnv::kNEARDET: return TVector3(-2.961e-01, -1.206e-01,+9.475e-01).Unit();

    default:
      throw cet::exception("BadGeoConfig")
        << "Bad geometry configuration: " << fDetId << "\n"
        << __FILE__ << ":" << __LINE__ << "\n";
    }
  }

  ///
  /// Return the ranges of x,y and z for the "world volume" that the
  /// entire geometry lives in. If any pointers are 0, then those
  /// coordinates are ignored.
  ///
  /// \param xlo : On return, lower bound on x positions
  /// \param xhi : On return, upper bound on x positions
  /// \param ylo : On return, lower bound on y positions
  /// \param yhi : On return, upper bound on y positions
  /// \param zlo : On return, lower bound on z positions
  /// \param zhi : On return, upper bound on z positions
  ///
  void GeometryBase::WorldBox(double* xlo, double* xhi,
                              double* ylo, double* yhi,
                              double* zlo, double* zhi) const  {
    const TGeoShape* s = gGeoManager->GetVolume("vWorld")->GetShape();
    assert(s);

    if (xlo || xhi) {
      double x1, x2;
      s->GetAxisRange(1,x1,x2); if (xlo) *xlo = x1; if (xhi) *xhi = x2;
    }
    if (ylo || yhi) {
      double y1, y2;
      s->GetAxisRange(2,y1,y2); if (ylo) *ylo = y1; if (yhi) *yhi = y2;
    }
    if (zlo || zhi) {
      double z1, z2;
      s->GetAxisRange(3,z1,z2); if (zlo) *zlo = z1; if (zhi) *zhi = z2;
    }
  }// end of WorldBox

  ///
  /// Return the ranges of x,y and z for the "Detector volume" that the
  /// entire geometry lives in. If any pointers are 0, then those
  /// coordinates are ignored.
  ///
  /// \param xlo : On return, lower bound on x positions
  /// \param xhi : On return, upper bound on x positions
  /// \param ylo : On return, lower bound on y positions
  /// \param yhi : On return, upper bound on y positions
  /// \param zlo : On return, lower bound on z positions
  /// \param zhi : On return, upper bound on z positions
  ///
  void GeometryBase::DetectorEnclosureBox(double* xlo, double* xhi,
                                          double* ylo, double* yhi,
                                          double* zlo, double* zhi) const  {
    const TGeoShape* s = gGeoManager->GetVolume("vDetEnclosure")->GetShape();
    assert(s);

    if (xlo || xhi) {
      double x1, x2;
      s->GetAxisRange(1,x1,x2); if (xlo) *xlo = x1; if (xhi) *xhi = x2;
    }
    if (ylo || yhi) {
      double y1, y2;
      s->GetAxisRange(2,y1,y2); if (ylo) *ylo = y1; if (yhi) *yhi = y2;
    }
    if (zlo || zhi) {
      double z1, z2;
      s->GetAxisRange(3,z1,z2); if (zlo) *zlo = z1; if (zhi) *zhi = z2;
    }
  }// end ofDetectorEnclosureBox

  ///---------------------------------------------------------------------------
  /// Return the ranges of x,y and z for the "Detector volume" that the
  /// entire geometry lives in. If any pointers are 0, then those
  /// coordinates are ignored.
  /// \param xlo : On return, lower bound on x positions
  /// \param xhi : On return, upper bound on x positions
  /// \param ylo : On return, lower bound on y positions
  /// \param yhi : On return, upper bound on y positions
  /// \param zlo : On return, lower bound on z positions
  /// \param zhi : On return, upper bound on z positions
  ///
  void GeometryBase::DetectorBigBox(double* xlo, double* xhi,
                                    double* ylo, double* yhi,
                                    double* zlo, double* zhi) const {
    //we don't want to return without passing along the big box coordinates
    //note the big box is set by default to just be the detector size
    //if 'isDetectorBigBox' is false
    //if(!this->isDetectorBigBoxUsed())return;
    *xlo = fDetectorBigBoxXLo;
    *xhi = fDetectorBigBoxXHi;
    *ylo = fDetectorBigBoxYLo;
    *yhi = fDetectorBigBoxYHi;
    *zlo = fDetectorBigBoxZLo;
    *zhi = fDetectorBigBoxZHi;
  }// end of DetectorBigBox

  //......................................................................
  bool GeometryBase::IntersectsDetector(double *xyz, double* dxyz) const {
    return geo::IntersectsBox(xyz, dxyz,
                              -fDetHalfWidth, fDetHalfWidth,
                              -fDetHalfHeight, fDetHalfHeight,
                              0., fDetLength);
  }// end of IntersectsDetector

  //......................................................................
  bool GeometryBase::IntersectsBigBox(double *xyz, double* dxyz) const  {
    return geo::IntersectsBox(xyz, dxyz,
                              fDetectorBigBoxXLo, fDetectorBigBoxXHi,
                              fDetectorBigBoxYLo, fDetectorBigBoxYHi,
                              fDetectorBigBoxZLo, fDetectorBigBoxZHi);
  }// end of IntersectsBigBox

  //......................................................................

  void GeometryBase::SetDrawOptions() { }

  //......................................................................

  ///
  /// Recursively search for planes in the geometry
  ///
  /// \param n     - Array holding the path of nodes traversed
  /// \param depth - Depth of the search to this point
  /// \param inMuonCatcher - internal flag, recursed into the muon catcher?
  /// \throws geo::Exception is maximum depth is exceeded.
  ///
  void GeometryBase::FindPlanes(std::vector<const TGeoNode*>& n,
                                unsigned int depth,
                                bool inMuonCatcher)
  {
    const char* nm = n[depth]->GetName();
    if (strncmp(nm,"vPlane", 6)==0) {
      this->MakePlane(n, depth, inMuonCatcher);
      return;
    }

    if(strncmp(nm, "vBlockMuon", 10) == 0) inMuonCatcher = true;

    // Explore the next layer down
    unsigned int deeper = depth+1;
    if (deeper>=n.size()) {
      throw cet::exception("BAD_NODE")
        << "Exceeded maximum depth (" << n.size() << ") " << deeper << "\n"
        << __FILE__ << ":" << __LINE__ << "\n";
    }
    const TGeoVolume* v = n[depth]->GetVolume();
    int nd = v->GetNdaughters();
    for (int i=0; i<nd; ++i) {
      n[deeper] = v->GetNode(i);
      this->FindPlanes(n, deeper, inMuonCatcher);
    }
  }// end of FindPlanes

  //........................................................................
  void GeometryBase::MakePlane(std::vector<const TGeoNode*>& n,
                               unsigned int depth,
                               bool inMuonCatcher)
  {
    fPlanes.push_back(new PlaneGeo(n, depth));

    if(inMuonCatcher) fPlanesInMuonCatcher.push_back(fPlanes.back());
  }

  //......................................................................
  ///
  /// Return the geometry description of the ith plane in the detector.
  ///
  /// \param i : input plane number, starting from 0
  /// \returns plane geometry for ith plane
  ///
  /// \throws geo::Exception if iplane is outside allowed range
  ///
  const PlaneGeo* GeometryBase::Plane(unsigned int i) const   {
    if (i>=fPlanes.size()) {
      throw cet::exception("BAD_PLANE_NUMBER")
        << "Plane index " << i << " out of range " << fPlanes.size() << "\n"
        << __FILE__ << ":" << __LINE__ << "\n";
      return 0;
    }
    return fPlanes[i];
  }//end of Plane

  //........................................................................
  ///
  /// Return the transverse position of the cell in which ever view its
  /// in
  ///
  /// \param ip - plane number
  /// \param ic - cell number
  /// \param w  - optional position along the local cell z-axis
  ///
  /// \returns - transverse position in units of cm
  ///
  /// \throws - geo::Exception is plane, cell, or view are out of range
  ///
  double GeometryBase::CellTpos(unsigned int ip, unsigned int ic, double w) const {
    if (ip>=this->NPlanes()) {
      throw cet::exception("BAD_PLANE_NUMBER")
        << " iplane " << ip << " [0," << this->NPlanes() << ")\n"
        << __FILE__ << ":" << __LINE__ << "\n";
    }
    const PlaneGeo* p = this->Plane(ip);
    if (p) {
      if (ic>=p->Ncells()) {
        throw cet::exception("BAD_CELL_NUMBER")
          << " iplane " << ic << " [0," << p->Ncells() << ")\n"
          << __FILE__ << ":" << __LINE__ << "\n";
      }
      const CellGeo* c = p->Cell(ic);
      if (c) {
        double pos[3];
        c->GetCenter(pos,w);
        if (p->View() == kX) { return pos[0]; }
        if (p->View() == kY) { return pos[1]; }
        throw cet::exception("BAD_PLANE_VIEW")
          << "view " << p->View() << "\n"
          << __FILE__ << ":" << __LINE__ << "\n";
        return 0.0;
      }
    }
    return 0.0;
  }// end of CellTpos

  //......................................................................
  ///
  /// Given a plane and cell number, return which view the cell
  /// measures, its nominal center position, and nominal size about that
  /// center position.
  ///
  /// \param ip - plane number
  /// \param ic - cell number
  /// \param view - on return, the view this cell measures in
  /// \param pos  - the nominal (xyz) center position of this cell (cm)
  /// \param dpos - half-dimensions of the cell in xyz (cm)
  /// \throw geo::Exception
  ///
  void GeometryBase::CellInfo(unsigned int ip, unsigned int ic,
                              View_t* view, double* pos, double* dpos) const  {
    if (ip>=this->NPlanes()) {
      throw cet::exception("BAD_PLANE_NUMBER")
        << " iplane " << ip << " [0," << this->NPlanes() << ")\n"
        << __FILE__ << ":" << __LINE__ << "\n";
    }
    const PlaneGeo* p = this->Plane(ip);
    if (p) {
      if (view) *view = p->View();

      if (ic>=p->Ncells()) {
        throw cet::exception("BAD_CELL_NUMBER")
          << " iplane " << ic << " [0," << p->Ncells() << ")\n"
          << __FILE__ << ":" << __LINE__ << "\n";
      }
      const CellGeo* c = p->Cell(ic);
      if (pos) c->GetCenter(pos);

      if (dpos) {
        if (p->View() == kX) {
          dpos[0] = c->HalfW();
          dpos[1] = c->HalfL();
          dpos[2] = c->HalfD();
        }
        else if (p->View() == kY) {
          dpos[0] = c->HalfL();
          dpos[1] = c->HalfW();
          dpos[2] = c->HalfD();
        }
      }
    }
  }// end of CellInfo

  //......................................................................
  ///
  /// Given a unique cell identifier, look up the plane and cell
  /// numbers.
  ///
  /// \param id : Cell identifier
  /// \param iplane : pointer to return plane number
  /// \param icell  : pointer to return cell number
  ///
  /// \returns cell geometry description
  /// \throws  geo::Exception if id is not found
  ///
  const CellGeo* GeometryBase::IdToCell(const CellUniqueId& id,
                                        int* iplane,
                                        int* icell) const {
    const OfflineChan offline_chan = getPlaneCellMap(id);
    *iplane = (int)offline_chan.Plane();
    *icell  = (int)offline_chan.Cell();
    return this->Plane(*iplane)->Cell(*icell);
  }// end of IdToCell

  //......................................................................
  ///
  /// Given a unique cell identifier, look up the plane number.
  ///
  /// \param id : Cell identifier
  /// \param iplane : optional, pointer to return plane number
  ///
  /// \returns Plane geometry description
  /// \throws  geo::Exception if id is not found
  ///

  const PlaneGeo* GeometryBase::IdToPlane(const CellUniqueId& id, int* iplane) const {
    *iplane = getPlaneCellMap(id).Plane();
    return this->Plane(*iplane);
  }// end of IdToCell
  //......................................................................
  int  GeometryBase::getPlaneID(const CellUniqueId& id) const{
    return getPlaneCellMap(id).Plane();
  }
  //......................................................................
  int  GeometryBase::getCellID (const CellUniqueId& id) const{
    return getPlaneCellMap(id).Cell();
  }
  //......................................................................
  bool GeometryBase::getPlaneAndCellID(const CellUniqueId& id, int& plane, int& cell) const{
    const OfflineChan offline_chan = getPlaneCellMap(id);
    plane = (int)offline_chan.Plane();
    cell  = (int)offline_chan.Cell();
    return true;
  }
  //......................................................................
  OfflineChan GeometryBase::getPlaneCellMap(const CellUniqueId& id) const{

    UIDMap::const_iterator itr = fIdMap.find(id);
    if (itr == fIdMap.end()) {
      throw cet::exception("BAD_CELL_UNIQUE_ID")
        << "CellUniqueId " << id << "\n"
        << __FILE__ << ":" << __LINE__ << "\n";
    }

    return itr->second;
  }
  
  //......................................................................
  /// Return list of planes which measure the requested projection
  ///
  /// \param v : X/Y or All (see View_t typedef in PlaneGeo.h)
  /// \returns : a set of plane indices satisfying the request
  /// \throws : geo::Exception if initialization fails
  ///
  const std::set<unsigned int>& GeometryBase::GetPlanesByView(View_t v) const {
    // Return the user's choice, fall back on all planes
    if      (v==kX) return fXplanes;
    else if (v==kY) return fYplanes;
    return fAllPlanes;
  }// end of GetPlanesByView

  //......................................................................
  /// Return the index of the next plane in a particular view
  ///
  /// \param p1 : Index of current plane
  /// \param d  : Direction to step in (>0 next down stream, <0 upstream)
  ///
  /// \returns Index of next plane in same view. If no plane can be
  /// found (eg. reach end of detector) return kPLANE_NOT_FOUND.
  ///
  /// \throws : Asserts that p1 is in valid range
  ///
  const unsigned int GeometryBase::NextPlaneInView(unsigned int p1, int d) const {
    assert(p1<fPlanes.size());

    // March along from p1 in the direction d until we find a match
    int s = (d>0 ? 1 : -1);
    if (s<0 && p1==0)                return kPLANE_NOT_FOUND;
    if (s>0 && p1==fPlanes.size()-1) return kPLANE_NOT_FOUND;
    for (unsigned int p2=p1+s; 1; p2+=s) {
      if (fPlanes[p1]->View() == fPlanes[p2]->View()) return p2;
      if (p2==fPlanes.size()-1) return kPLANE_NOT_FOUND;
      if (p2==0)                return kPLANE_NOT_FOUND;
    }
  }// end of NextPlaneInView

  //......................................................................
  /// Return the index of the next plane in a particular view
  ///
  /// \param p1 : Index of current plane
  /// \param d  : Direction to step in (>0 next down stream, <0 upstream)
  ///
  /// \returns index of next plane in same view.
  ///
  /// \throws : Asserts that p1 is in valid range
  ///
  const unsigned int GeometryBase::NextPlaneOtherView(unsigned int p1, int d) const
  {
    assert(p1<fPlanes.size());

    // March along from p1 in the direction d until we find a match
    int s = (d>0 ? 1 : -1);
    if (s<0 && p1==0)                return kPLANE_NOT_FOUND;
    if (s>0 && p1==fPlanes.size()-1) return kPLANE_NOT_FOUND;
    for (unsigned int p2=p1+s; 1; p2+=s) {
      if (fPlanes[p1]->View() != fPlanes[p2]->View()) return p2;
      if (p2==fPlanes.size()) return kPLANE_NOT_FOUND;
      if (p2==0)              return kPLANE_NOT_FOUND;
    }
  }// end of NextPlaneOtherView

  //......................................................................
  /// \brief Returns the index of the first plane contained in the muon catcher
  ///
  /// Or kPLANE_NOT_FOUND if there is no muon catcher
  const unsigned int GeometryBase::FirstPlaneInMuonCatcher() const
  {
    if(fPlanesInMuonCatcher.empty()){
      mf::LogWarning("GeometryBase") << "Requesting muon catcher info for a detector without one";
      return kPLANE_NOT_FOUND;
    }

    auto it = std::find(fPlanes.begin(), fPlanes.end(), fPlanesInMuonCatcher[0]);

    assert(it != fPlanes.end());
    return it-fPlanes.begin();
  }

  //......................................................................
  ///
  /// Build maps used for quick access to the geometry
  ///
  void GeometryBase::BuildMaps ()
  {
    fAllPlanes.clear();
    fXplanes.clear();
    fYplanes.clear();

    // Sets which list planes by view
    for (unsigned int i=0; i<fPlanes.size(); ++i) {
      fAllPlanes.insert(i);
      if      (fPlanes[i]->View()==kX) fXplanes.insert(i);
      else if (fPlanes[i]->View()==kY) fYplanes.insert(i);
      else {
        throw cet::exception("BAD_GEO_CONFIG")
          << "Bad plane view" << fPlanes[i]->View() << ", plane " << i << "\n"
          << __FILE__ << ":" << __LINE__ << "\n";
      }
    }

    // Unique Id -> cell map
    for (unsigned int i=0; i<fPlanes.size(); ++i) {
      const PlaneGeo* p = fPlanes[i];
      for (unsigned int j=0; j<p->Ncells(); ++j) {
        // Assert that the cell ids are in fact unique
        if (fIdMap.find(this->Plane(i)->Cell(j)->Id())!=fIdMap.end()) {
          throw cet::exception("BAD_CELL_UNIQUE_ID")
            << "Duplicate Cell UIDs" << this->Plane(i)->Cell(j)->Id()
            << ", plane " << i << ", cell " << j << "\n"
            << __FILE__ << ":" << __LINE__ << "\n";
        }
        OfflineChan c(i,j);
        fIdMap[this->Plane(i)->Cell(j)->Id()] = c;
      }
    }
  }// end of BuildMaps

  //......................................................................
  ///
  /// Return the total mass of the detector
  ///
  ///
  double GeometryBase::TotalMass(const char *volume) const
  {
    double mass = 0.;

    // the detector resides in the vDetEnclosure volume so start there
    const TGeoVolume *enclosure = gGeoManager->FindVolumeFast(volume);

    int old_verbosity = gGeoManager->GetVerboseLevel();
    gGeoManager->SetVerboseLevel(0);
    // loop over the daughters in the volume and let ROOT calculate the
    // mass in kg for you for the daughter nodes that have vBlock in
    // their names this only goes one level down, so no worries of
    // double counting non-superblock blocks
    //
    for(int d = 0; d < enclosure->GetNdaughters(); ++d){
      const TGeoNode* node = enclosure->GetNode(d);
      const char* nm = node->GetName();
      if(strncmp(nm,"vSuperBlock",11) == 0){
        std::string name = nm;
        name.erase(name.find('_'));
        mass += gGeoManager->FindVolumeFast(name.c_str())->Weight();
      }
    }
    gGeoManager->SetVerboseLevel(old_verbosity);
    return mass;
  }// end of TotalMass


  //---------------------------------------------------------------------------
  ///
  /// Method to calculate masses of the detector
  ///
  bool GeometryBase::calculateMassesLong(const unsigned int number_of_points,
                                         CLHEP::HepRandomEngine& engine) const {
    
    std::map <std::string, double> density;         // Map of densities of all materials in the detector. Dimensionless. Mapped by material names.
    std::map <std::string, double> volume;          // Map of volumes of all materials in the detector. Dimensionless. Mapped by material names.
    std::map <std::string, double> mass;            // Map of total masses of all materials the detector. Dimensionless. Mapped by material names.
    std::map <std::string, double> volume_fiducial; // Map of volumes of all materials in the fiducial region of the detector. Dimensionless. Mapped by material names.
    std::map <std::string, double> mass_fiducial;   // Map of total masses of all materials in the fiducial region of the detector. Dimensionless. Mapped by material names.
    
    CLHEP::RandFlat   flat(engine);
    
    TStopwatch timer;
    timer.Start();
    
    // Effective cell volume for a point
    const double cell_volume = (2.0*fDetHalfWidth*CLHEP::cm)*(2.0*fDetHalfHeight*CLHEP::cm)*(fDetLength*CLHEP::cm) / double(number_of_points);

    // Loop over points
    for(unsigned int ipoint = 0; ipoint < number_of_points; ++ipoint){

      // Set the progress bar
      ROOTGeoManager()->GetGeomPainter()->OpProgress("Calculating masses", ipoint, number_of_points-1, &timer, kFALSE);

      // X, Y, Z of the point in cm
      const double x = flat.fire(-fDetHalfWidth , fDetHalfWidth ) * CLHEP::cm;
      const double y = flat.fire(-fDetHalfHeight, fDetHalfHeight) * CLHEP::cm;
      const double z = flat.fire(0.0            , fDetLength    ) * CLHEP::cm;
      const double current_position[3] = {x,y,z};
      
      // Get the material of the point
      const TGeoMaterial* current_material = Material(x/CLHEP::cm,
                                                      y/CLHEP::cm,
                                                      z/CLHEP::cm
                                                      );
      // Get the material name
      const std::string current_material_name = current_material->GetName();

      // See if we have not defined the material
      if(density.find(current_material_name) == density.end()){
        density        [current_material_name] = current_material->GetDensity() * (CLHEP::g/CLHEP::cm3); //http://indra.in2p3.fr/KaliVedaDoc/1.8.1/src/KVedaLossMaterial.h.html
        // Initialize these to zero
        volume         [current_material_name] = 0.0;
        volume_fiducial[current_material_name] = 0.0;
        mass           [current_material_name] = 0.0;
        mass_fiducial  [current_material_name] = 0.0;
      }// end of searching current_material_name

      volume[current_material_name] += cell_volume;
      mass  [current_material_name] += cell_volume * density[current_material_name];
      
      // See if the point is inside fiducial volume
      if(isInsideFiducialVolume(current_position)){
        volume_fiducial[current_material_name] += cell_volume;
        mass_fiducial  [current_material_name] += cell_volume * density[current_material_name];
      }// end of checking if the point is inside fiducial volume
      
    }// end of looping over points
    
    std::cout<<"Calculated masses, volumes etc.:\n";
    
    // Printing the mass table for each detector material
    for(std::map<std::string, double>::iterator iter = density.begin(); iter != density.end(); ++iter)  {
      
      const std::string key_str = iter->first;
      
      std::cout<<std::setw(9)<<"\t\tMaterial="<<std::setw(20)<<key_str
               <<std::setw(10)<<" density="<<std::setw(10)<<density[key_str]/CLHEP::kg*CLHEP::m3<<std::setw(6)<<"kg/m3"
               <<std::setw(9)<<" volume="<<std::setw(10)<<volume[key_str]/CLHEP::m3<<std::setw(2)<<"m3"
               <<std::setw(7)<<" mass="<<std::setw(10)<<mass[key_str] / CLHEP::kg<<std::setw(2)<<"kg"
               <<std::setw(17)<<" volume_fiducial="<<std::setw(10)<<volume_fiducial[key_str] / CLHEP::m3<<std::setw(1)<<"m3"
               <<std::setw(15)<<" mass_fiducial="<<std::setw(10)<<mass_fiducial[key_str] / CLHEP::kg<<std::setw(3)<<"kg\n";
      
    }// end of printing the mass table for each detector material

    return true;
  }
  
  //......................................................................
  ///
  /// Return the column density between 2 points in gm/cm^2
  ///
  /// \param p1 : pointer to array holding xyz of first point in world
  /// coordinates
  ///
  /// \param p2 : pointer to array holding xyz of second point in world
  /// coordinates
  ///
  double GeometryBase::MassBetweenPoints(double *p1, double *p2) const
  {
    // The purpose of this method is to determine the column density
    // between the two points given.  Do that by starting at p1 and
    // stepping until you get to the node of p2.  calculate the distance
    // between the point just inside that node and p2 to get the last
    // bit of column density
    double columnD = 0.;

    // first initialize a track - get the direction cosines
    double length = util::pythag(p2[0]-p1[0],
                                 p2[1]-p1[1],
                                 p2[2]-p1[2]);
    double dxyz[3] = {(p2[0]-p1[0])/length,
                      (p2[1]-p1[1])/length,
                      (p2[2]-p1[2])/length};

    gGeoManager->InitTrack(p1,dxyz);
    gGeoManager->FindNextBoundary(1.E-6);
    
    // might be helpful to have a point to a TGeoNode
    TGeoNode *node = gGeoManager->GetCurrentNode();

    // check that the points are not in the same volume already.
    // if they are in different volumes, keep stepping until you
    // are in the same volume as the second point
    while(!gGeoManager->IsSameLocation(p2[0], p2[1], p2[2])){
      gGeoManager->FindNextBoundary();

      columnD +=
        gGeoManager->GetStep() *
        node->GetMedium()->GetMaterial()->GetDensity();

      // the act of stepping puts you in the next node and returns that
      // node
      node = gGeoManager->Step();
    } //end loop to get to volume of second point

    // now you are in the same volume as the last point, but not at that
    // point.  get the distance between the current point and the last
    // one
    const double *current = gGeoManager->GetCurrentPoint();
    length = util::pythag(p2[0]-current[0],
                          p2[1]-current[1],
                          p2[2]-current[2]);
    columnD += length*node->GetMedium()->GetMaterial()->GetDensity();

    return columnD;
  }// end of MassBetweenPoints

  //......................................................................

  void GeometryBase::MaterialsBetweenPoints
  (const double *p1,
   const double *p2,
   std::vector<double>& ds,
   std::vector<const TGeoMaterial*>& m) const
  {
    double d[3], dmag;
    dmag = util::pythag(p2[0]-p1[0],p2[1]-p1[1],p2[2]-p1[2]);
    d[0] = (p2[0]-p1[0])/dmag;
    d[1] = (p2[1]-p1[1])/dmag;
    d[2] = (p2[2]-p1[2])/dmag;
    
    gGeoManager->InitTrack(p1,d);
    gGeoManager->FindNextBoundary(1.E-6);
    TGeoNode *node = gGeoManager->GetCurrentNode();

    while (!gGeoManager->IsSameLocation(p2[0], p2[1], p2[2])) {
      gGeoManager->FindNextBoundary();
      
      // If we take a step and haven't gone anywhere, then we will be
      // stuck in a very long loop. So we will quit if the next step
      // takes us nowhere.
      double step = gGeoManager->GetStep();
      if(step < 1.0e-12) break;

      ds.push_back(step);
      m.push_back(node->GetMedium()->GetMaterial());
      
      node = gGeoManager->Step();
      if (node==0) {
        std::cerr << __FILE__ << ":" << __LINE__ 
                  << " Cannot continue tracking MaterialsBetweenPoints(). "
                  << "Null node pointer." << std::endl;
        return;
      }
    }
    //
    // Handle final step by hand
    //
    const double* p3 = gGeoManager->GetCurrentPoint();
    ds.push_back(util::pythag(p3[0]-p2[0],p3[1]-p2[1],p3[2]-p2[2]));
    m.push_back(node->GetMedium()->GetMaterial());
  }

  ///
  ///Overloaded MaterialsBetweenPoints to for using TVector3 instead of 
  ///arrays because rb::Track trajectory points are TVector3 and there 
  ///is a desire to track a 'track' through the geometry.  This will make
  ///that easier.  Function is really just a wrapper around the array version
  ///
  void GeometryBase::MaterialsBetweenPoints
  (TVector3 v1,
   TVector3 v2,
   std::vector<double>& ds,
   std::vector<const TGeoMaterial*>& m) const
  {
    double p1[3], p2[3];
    v1.GetXYZ(p1);
    v2.GetXYZ(p2);
    
    MaterialsBetweenPoints(p1, p2, ds, m);
  }
  
  ///
  /// Calculate the detector size based on the geometry. Requires a
  /// little bit of work, so do this once when the geometry is updated
  ///
  void GeometryBase::SetDetectorSize()
  {
    double dummy;

    // Compute the height and width based on the sizes of the vertical
    // and horizontal planes. Remember, in the local plane frame z goes
    // along the length of the modules. Hence the use of axis=3 to get
    // the detector height (from the vertical planes) and width (from
    // the horizontal planes).
    const TGeoVolume* vvp = gGeoManager->GetVolume("vPlaneV");
    const TGeoVolume* hvp = gGeoManager->GetVolume("vPlaneH");

    // Handle the case of the IPND where the volumes are names slightly
    // different accounting for the different widths of the modules.
    if (vvp==0) vvp = gGeoManager->GetVolume("vPlaneVNDOS");
    if (hvp==0) hvp = gGeoManager->GetVolume("vPlaneHNDOS");
    if (vvp==0 && hvp==0) {
      vvp = gGeoManager->GetVolume("vPlaneVND");
      hvp = gGeoManager->GetVolume("vPlaneHND");
    }

    /// If there no Vertical or horizontal planes, then we go get the size from the Detector Enclosure
    /// It used to be that once this condition is true, then it we would just throw an exception.
    /// However, now Giuseppe Ferone and Denis Perevalov are studying interactions in simple detector,
    /// which don't have planes in it. Thus, they don't want the code breaking.
    if (vvp==0 || hvp==0) {
      mf::LogWarning("SetDetectorSize") 
        << "Warning! BAD_GEO_CONFIG. Unable to find shapes to set detector size\n"
        << "Now getting the detector size from vDetEnclosure\n";
      
      /// Getting Detector enclosure volume
      const TGeoVolume* vdet_enclosure = gGeoManager->GetVolume("vDetEnclosure");
      
      /// Lo and Hi ranges of each of the axes. To be filled later
      double x_lo, x_hi;
      double y_lo, y_hi;
      double z_lo, z_hi;
      
      /// Fill them
      vdet_enclosure->GetShape()->GetAxisRange(1, x_lo, x_hi);
      vdet_enclosure->GetShape()->GetAxisRange(2, y_lo, y_hi);
      vdet_enclosure->GetShape()->GetAxisRange(3, z_lo, z_hi);
      
      /// Set detector size
      fDetHalfWidth  = x_hi;
      fDetHalfHeight = y_hi;
      fDetLength     = z_hi - z_lo;
      
      return;
    }
    
    const TGeoShape* vp = vvp->GetShape();
    const TGeoShape* hp = hvp->GetShape();
    vp->GetAxisRange(3,dummy,fDetHalfHeight);
    hp->GetAxisRange(3,dummy,fDetHalfWidth);
    
    // Remember -- in the local plane coordinate frame, the "depth" is
    // along y, hence axis=2. In the world coordinate this is along z.
    double vplanez1, vplanez2;
    vp->GetAxisRange(2,vplanez1,vplanez2);
    
    // Compute the detector length based on the center positions of
    // cells in the first and last plane. Adjust to account for the
    // plane widths. This ajustment works if the vertical and horizontal
    // planes have the same depth, or if the detector begins and ends
    // with vertical planes. As of the TDR both of these conditions was
    // true for all the NOvA detectors.
    double xyz1[3] = {0,0,0};
    double xyz2[3] = {0,0,0};
    this->Plane(0)->               Cell(0)->GetCenter(xyz1);
    this->Plane(fPlanes.size()-1)->Cell(0)->GetCenter(xyz2);
    fDetLength = (xyz2[2]-xyz1[2])+(vplanez2-vplanez1);
  }// end of SetDetectorSize
  
  //......................................................................
  
  void GeometryBase::setDetectorBigBox(double detector_bigbox_range) {
    fIsDetectorBigBoxUsed = true;

    double x1, x2;
    double y1, y2;
    double z1, z2;

    /// Get the Detector size
    this->DetectorEnclosureBox(&x1, &x2, &y1, &y2, &z1, &z2);

    /// Dimensions of the DetectorBigBox in cm
    fDetectorBigBoxXHi = x2 + detector_bigbox_range;
    fDetectorBigBoxXLo = x1 - detector_bigbox_range;
    fDetectorBigBoxYHi = y2 + detector_bigbox_range;
    fDetectorBigBoxYLo = y1 - detector_bigbox_range;
    fDetectorBigBoxZHi = z2 + detector_bigbox_range;
    fDetectorBigBoxZLo = z1 - detector_bigbox_range;

  }// end of setDetectorBigBox

  //......................................................................
  bool GeometryBase::isInsideDetectorBigBox(const double x, const double y, const double z) const {
    if( x < fDetectorBigBoxXHi && x > fDetectorBigBoxXLo &&
        y < fDetectorBigBoxYHi && y > fDetectorBigBoxYLo &&
        z < fDetectorBigBoxZHi && z > fDetectorBigBoxZLo)
      return true;
    return false;
  }// end of GeometryBase::isInsideDetectorBigBox

  //......................................................................
  ///
  /// Is the particle inside the fiducial volume?
  ///
  /// \param x_cm : x coordinate in cm
  ///
  /// \param y_cm : y coordinate in cm
  ///
  /// \param z_cm : z coordinate in cm
  ///
  bool GeometryBase::isInsideFiducialVolume(const double x_cm, const double y_cm, const double z_cm) const {

    return (x_cm < fFiducialVolumeXHi && x_cm > fFiducialVolumeXLo &&
            y_cm < fFiducialVolumeYHi && y_cm > fFiducialVolumeYLo &&
            z_cm < fFiducialVolumeZHi && z_cm > fFiducialVolumeZLo);
  }// end of GeometryBase::isInsideFiducialVolume

  //......................................................................
  ///
  /// Is the particle inside the fiducial volume?
  ///
  /// \param xyz : pointer to array holding xyz in world coordinates
  ///
  bool GeometryBase::isInsideFiducialVolume(const double* xyz) const {

    return isInsideFiducialVolume(xyz[0]/CLHEP::cm,
                                  xyz[1]/CLHEP::cm,
                                  xyz[2]/CLHEP::cm);
  }// end of GeometryBase::isInsideFiducialVolume

  //......................................................................
  void GeometryBase::FiducialBox(TVector3& r0, TVector3& r1) const
  {
    r0 = TVector3(fFiducialVolumeXLo, fFiducialVolumeYLo, fFiducialVolumeZLo);
    r1 = TVector3(fFiducialVolumeXHi, fFiducialVolumeYHi, fFiducialVolumeZHi);
  }

  //----------------------------------------------------------------------
  ///
  /// Compute the distance from a point to the detector wall along a
  /// particular direction.
  ///
  /// \param x0    - Input, [x,y,z] start point (cm)
  /// \param dx    - Input, direction vector
  /// \param exitp - Ouput, Calculated [x,y,z] of exit point (cm)
  ///
  /// \returns distance to exit point in cm
  ///
  double GeometryBase::DEdge(const double* x0,
                             const double* dx,
                             double*       exitp) const
  {
    double xyzout[3];
    geo::ProjectToBoxEdge(x0, dx,
                          -fDetHalfWidth,  fDetHalfWidth,
                          -fDetHalfHeight, fDetHalfHeight,
                          0,               fDetLength,
                          xyzout);
    if (exitp!=0) {
      exitp[0] = xyzout[0];
      exitp[1] = xyzout[1];
      exitp[2] = xyzout[2];
    }
    return util::pythag(x0[0]-xyzout[0],
                        x0[1]-xyzout[1],
                        x0[2]-xyzout[2]);
  }// end of GeometryBase::DEdge

  //--------------------------------------------------------------------------
  /// \brief Find the distance of closest approach between a cell and a track
  ///
  /// @param ip  - plane number
  /// @param ic  - cell number
  /// @param vtx - track start point
  /// @param dir - track direction
  /// @param dx  - offset within the cell (in cm, dx=0 is center of cell)
  /// @param dy  - offset within the cell (in cm, dy=0 is center of cell)
  /// @param w   - Distance from read-out end where the closest approach occurs (cm)
  /// @param s   - Distance down the track where the closest approach occurs (cm)
  /// @param point_in_cell - Point inside the cell where the closest approach occurs
  /// @param point_on_track - Point on track where the closest approach occurs
  /// @returns The distance^2 of closest approach
  ///
  double GeometryBase::ClosestApproach(unsigned int    ip,
                                       unsigned int    ic,
                                       const TVector3& vtx,
                                       const TVector3& dir,
                                       double          dx,
                                       double          dy,
                                       double*         w,
                                       double*         s,
                                       TVector3*       point_in_cell,
                                       TVector3*       point_on_track) const
  {
    // ==
    // Find the line that runs down the cell at the request position
    //
    const PlaneGeo* plane = this->Plane(ip);
    View_t view = plane->View();;

    const CellGeo* cell = plane->Cell(ic);

    double p0[3], p1[3];
    double dz = cell->HalfL();
    cell->GetCenter(p0, +dz);
    cell->GetCenter(p1, -dz);
    
    // Shift to the ray requested
    p0[2] += dy;
    p1[2] += dy;
    if      (view==kX) { p0[0] += dx; p1[0] += dx; }
    else if (view==kY) { p0[1] += dx; p1[1] += dx; }

    TVector3 P0(p0[0], p0[1], p0[2]);
    TVector3 P1(p1[0], p1[1], p1[2]);
    // ==
    
    // Get two points on the track, the vertex and a point at vertex+dir
    TVector3 Q1(vtx); Q1 += dir;
    
    // Use the closest approach function to compute the rest
    return geo::ClosestApproach(P0,  P1, 
                                vtx, Q1, 
                                w,   s, 
                                point_in_cell,
                                point_on_track);
  }
  
  //--------------------------------------------------------------------------
  TGeoManager*  GeometryBase::ROOTGeoManager() const { return gGeoManager;}


  //--------------------------------------------------------------------------
  ///
  /// Method to set DetectorID
  /// 
  void GeometryBase::setDetectorID(novadaq::cnv::DetId det_id)
  {
    MF_LOG_DEBUG("Geometry")
      <<"GeometryBase::setDetectorID\n";

    fDetId = det_id;

    // Change coordinate transformation class as well
    fCoordinateTransformation.setDetectorAndBeam(fDetId, geo::BEAMTYPE_NUMI);
    //fCoordinateTransformation.print();

    // Work out which config to use
    DetectorParams config;
    switch(fDetId){
    case novadaq::cnv::kNDOS:     config = fParams.ndos(); break;
    case novadaq::cnv::kNEARDET:  config = fParams.nd();   break;
    case novadaq::cnv::kFARDET:   config = fParams.fd();   break;
    case novadaq::cnv::kTESTBEAM: config = fParams.tb();   break;
    default: return; // Probably unknown or testbeam. Nothing sane to set the
                     // below variables to.
    }

    // Set fiducial volume
    fFiducialVolumeXHi =  fDetHalfWidth  - config.FiducialVolumeXHi();
    fFiducialVolumeXLo = -fDetHalfWidth  + config.FiducialVolumeXLo();
    fFiducialVolumeYHi =  fDetHalfHeight - config.FiducialVolumeYHi();
    fFiducialVolumeYLo = -fDetHalfHeight + config.FiducialVolumeYLo();
    fFiducialVolumeZHi =  fDetLength     - config.FiducialVolumeZHi();
    fFiducialVolumeZLo =                   config.FiducialVolumeZLo();

    fIsDetectorBigBoxUsed = config.BigBoxUsed();
    fBigBoxRange = config.BigBoxRange();

  }// end of GeometryBase::setDetectorID


  //--------------------------------------------------------------------------
  ///
  /// Method to count the number of cells passed through by a line segment 
  /// between two points in the detector.
  ///
  void GeometryBase::CountCellsOnLine(double x1, double y1, double z1, double x2, double y2, double z2,
                                      std::vector<OfflineChan>& xhitsonline,
                                      std::vector<OfflineChan>& yhitsonline)
  {
    std::vector<CellOnLine> xhs, yhs;
    CountCellsOnLine(x1, y1, z1, x2, y2, z2, xhs, yhs);

    for(const CellOnLine& c: xhs) xhitsonline.push_back(c.chan);
    for(const CellOnLine& c: yhs) yhitsonline.push_back(c.chan);
  }  // end CountCellsOnLine

  //---------------------------------------------------------------------------
  ///
  /// Method to count the number of cells passed through by a line segment 
  /// between two points in the detector.
  ///
  void GeometryBase::CountCellsOnLine(TVector3 r1, TVector3 r2,
                                  std::vector<OfflineChan>& xhitsonline,
                                  std::vector<OfflineChan>& yhitsonline)
  {
    CountCellsOnLine(r1.X(), r1.Y(), r1.Z(), r2.X(), r2.Y(), r2.Z(),
                     xhitsonline, yhitsonline);
  }

  //---------------------------------------------------------------------------
  ///
  /// Method to count the number of cells passed through by a line segment
  /// between two points in the detector and keep track of the entry and exit
  /// points of the line on the cells
  ///
  void GeometryBase::CountCellsOnLine(double x1, double y1, double z1, double x2, double y2, double z2,
                                      std::vector<CellOnLine>& xhitsonline, std::vector<CellOnLine>& yhitsonline)
  {

    // This function returns a vector of cells crossed by the line and the
    // entry and exit points of the line through the cells for each view. If
    // the end point is within a cell, the entry or exit point is set to be
    // the end of the line.
    //
    // WARNING: This function makes use of gGeoManager and is therefore
    //            accessing global variables!  Take caution to properly
    //            initialize these global variables when accessing them after
    //            using this function!

    // define initial position and normalized direction vector
    double xyz[3]  = {x1, y1, z1};
    double d       = util::pythag(x2-x1, y2-y1, z2-z1);
    double dxyz[3] = {(x2-x1)/d, (y2-y1)/d, (z2-z1)/d};

    gGeoManager->InitTrack(xyz, dxyz);
    const double* pos = gGeoManager->GetCurrentPoint();

    double dtemp = util::pythag(pos[0]-x1, pos[1]-y1, pos[2]-z1);

    bool enter = false;

    // We update the properties of this cell intersection as we go, and store
    // it one every exit.
    CellOnLine col;

    // create a list of positions in scintillator
    while(dtemp <= d) {
      // if the current material is scintillator, count that as a cell passed through
      if(strncmp("Scintillator", gGeoManager->GetCurrentNode()->GetVolume()->GetMaterial()->GetName(), 12) == 0) {
        enter = true;

        col.entry = TVector3(pos[0],pos[1],pos[2]);

        int plane, cell;
        IdToCell(CurrentCellId(), &plane, &cell);

        col.chan = OfflineChan(plane, cell);
      }

      // step forward to the next material boundary
      gGeoManager->FindNextBoundary();
      gGeoManager->Step(kTRUE,kTRUE);

      dtemp = util::pythag(pos[0]-x1, pos[1]-y1, pos[2]-z1);

      if(enter){
        enter = false;

        if(dtemp <= d){
          col.exit = TVector3(pos[0], pos[1], pos[2]);
        }
        else{
          col.exit = TVector3(x2, y2, z2);
        }

        const View_t view = Plane(col.chan.Plane())->View();
        if(view == geo::kX){ xhitsonline.push_back(col); }
        if(view == geo::kY){ yhitsonline.push_back(col); }

      }
    }

  }  // end CountCellsOnLine

  //---------------------------------------------------------------------------
  ///
  /// Method to count the number of cells passed through by a line segment 
  /// between two points in the detector and keep track of the entry and exit points
  /// of the line on the cells
  ///
  void GeometryBase::CountCellsOnLine(TVector3 r1, TVector3 r2,
                                      std::vector<CellOnLine>& Xhitsonline,
                                      std::vector<CellOnLine>& Yhitsonline)
  {
    CountCellsOnLine(r1.X(), r1.Y(), r1.Z(), r2.X(), r2.Y(), r2.Z(),
                     Xhitsonline, Yhitsonline);
  }

  //---------------------------------------------------------------------------
  ///
  /// Fast method to count the number of cells passed through by a line 
  /// segment between two points in the detector.
  ///
  void GeometryBase::CountCellsOnLineFast(geo::View_t view,
                                      double v1, double z1,
                                      double v2, double z2,
                                      std::vector<OfflineChan>& hitsonline)
  {
    assert(view == geo::kX || view == geo::kY);

    std::vector<OfflineChan> junk;
    if(view == geo::kX)
      CountCellsOnLineFast(v1, 0, z1, v2, 0, z2, hitsonline, junk);
    else
      CountCellsOnLineFast(0, v1, z1, 0, v2, z2, junk, hitsonline);
  }

  //---------------------------------------------------------------------------
  ///
  /// Fast method to count the number of cells passed through by a line
  /// segment between two points in the detector.
  ///
  void GeometryBase::CountCellsOnLineFast(double x1, double y1, double z1,
                                          double x2, double y2, double z2,
                                          std::vector<OfflineChan>& xhitsonline,
                                          std::vector<OfflineChan>& yhitsonline)
  {
    CountCellsOnLineFast(TVector3(x1, y1, z1), TVector3(x2, y2, z2),
                         xhitsonline, yhitsonline);
  }

  //---------------------------------------------------------------------------
  ///
  /// Fast method to count the number of cells passed through by a line
  /// segment between two points in the detector.
  ///
  void GeometryBase::CountCellsOnLineFast(TVector3 r1, TVector3 r2,
                                      std::vector<OfflineChan>& xhitsonline,
                                      std::vector<OfflineChan>& yhitsonline)
  {
    // Otherwise we would count when the ray was outside the detector in the
    // other view.
    ClampRayToDetectorVolume(&r1, &r2, this);

    // We assume downstream-going below
    if(r1.Z() > r2.Z()) std::swap(r1, r2);

    const double z1 = r1.Z();
    const double z2 = r2.Z();

    const unsigned int planeMax = fPlanes.size();
    for(unsigned int plane = 0; plane < planeMax; ++plane){
      const geo::PlaneGeo* geoplane = fPlanes[plane];
      const geo::CellGeo* geocell0 = geoplane->Cell(0);
      // Assuming here perfect alignment and identical cells across plane
      double xyzp[3];
      geocell0->GetCenter(xyzp);
      const double zc = xyzp[2];
      const double dz = geocell0->HalfD();

      if(z1 > zc+dz) continue;
      // Plane positions increase monotonically
      if(z2 < zc-dz) break;

      const geo::View_t view = geoplane->View();

      const double v1 = r1[view];
      const double v2 = r2[view];

      const double denom = z2-z1;
      const double dvdz = denom ? (v2-v1)/denom : 0;

      // Position entering the cell or starting
      const double zin = std::max(zc-dz, z1);
      double vlo = v1+(zin-z1)*dvdz;
      // Position exiting or ending
      const double zout = std::min(zc+dz, z2);
      double vhi = v1+(zout-z1)*dvdz;
      // Ensure they're from low to high instead
      if(vlo > vhi) std::swap(vlo, vhi);

      for(unsigned int cell = 0; cell < geoplane->Ncells(); ++cell){
        double xyz[3];
        geoplane->Cell(cell)->GetCenter(xyz);
        const double vc = xyz[view];
        const double dv = geoplane->Cell(cell)->HalfW();

        if(vlo > vc+dv) continue;
        // Cell positions increase monotonically
        if(vhi < vc-dv) break;

        if(view == geo::kX) xhitsonline.push_back(geo::OfflineChan(plane, cell));
        if(view == geo::kY) yhitsonline.push_back(geo::OfflineChan(plane, cell));
      } // end for cell
    } // end for plane
  }

  //......................................................................
  /// Return distance from a point in the cell to readout electronics, including pigtail
  /// \param localz : Z position in cm in local cell frame
  double GeometryBase::DistToElectronics(double localz, const CellGeo& cell) const
  {
    return (cell.HalfL() - localz + this->GetPigtail(cell.Id()));
  }

  //......................................................................
  double GeometryBase::GetPigtail(const CellUniqueId& id) const
  {
    int planeid,cellid;
    this->IdToCell(id, &planeid, &cellid);
    View_t view = (this->IdToPlane(id,&planeid))->View();
    //NOTE: Do we think 32 cells per module can ever change
    //This number can be found in DAQChannelMap but did not want to add that
    //dependency to geometry unless needed.
    const int kCellsPerModule = 32;
    cellid = cellid % kCellsPerModule;
    // In vertical planes, high cell numbers have longer fibres.
    // In horizontal planes it's the opposite, so flip it round.
    // For Test Beam, vertical planes actually follow the same pattern as
    // the horizontals, so flip both.
    if( view == geo::kY || this->DetId() == novadaq::cnv::kTESTBEAM ) cellid = kCellsPerModule-cellid-1;
    // This should really never happen, but just to be safe...
    if(cellid < 0 || cellid >= kCellsPerModule) return 100;
    // Email from Tom Chase 2011-04-29
    // NB: this isn't just a ~3.8cm change per cell. Every 8 cells something
    // different happens.
    const double kPigtails[kCellsPerModule] = {
      34.5738, 38.4379, 42.3020, 46.1660, 50.0301, 53.8942, 57.7582, 61.6223,
      64.7504, 68.6144, 72.4785, 76.3426, 80.2067, 84.0707, 87.9348, 91.0790,
      95.3301, 99.1941, 103.058, 106.922, 110.786, 114.650, 118.514, 122.379,
      125.507, 129.371, 133.235, 137.099, 140.963, 144.827, 148.691, 150.751
    };
    return kPigtails[cellid];
  }
  
} // end namespace geo
////////////////////////////////////////////////////////////////////////