ROOTGeomAnalyzer.cxx

Go to the documentation of this file.

//____________________________________________________________________________
/*
 Copyright (c) 2003-2020, The GENIE Collaboration
 For the full text of the license visit http://copyright.genie-mc.org

 Anselmo Meregaglia <anselmo.meregaglia \at cern.ch>
 ETH Zurich

 Costas Andreopoulos <constantinos.andreopoulos \at cern.ch>
 University of Liverpool & STFC Rutherford Appleton Laboratory

 Robert Hatcher <rhatcher \at fnal.gov>
 Fermilab
*/
//____________________________________________________________________________

#include <cassert>
#include <cstdlib>
#include <iomanip>
#include <set>

#include <TGeoVolume.h>
#include <TGeoManager.h>
#include <TGeoShape.h>
#include <TGeoMedium.h>
#include <TGeoMaterial.h>
#include <TGeoMatrix.h>
#include <TGeoNode.h>
#include <TObjArray.h>
#include <TLorentzVector.h>
#include <TList.h>
#include <TSystem.h>
#include <TMath.h>
#include <TPolyMarker3D.h>
#include <TGeoBBox.h>

#include "Framework/Conventions/GBuild.h"
#include "Framework/Conventions/Units.h"
#include "Framework/Conventions/Controls.h"
#include "Tools/Geometry/PathSegmentList.h"
#include "Framework/EventGen/PathLengthList.h"
#include "Framework/EventGen/GFluxI.h"
#include "Tools/Geometry/ROOTGeomAnalyzer.h"
#include "Tools/Geometry/GeomVolSelectorI.h"
#include "Framework/Messenger/Messenger.h"
#include "Framework/Numerical/RandomGen.h"
#include "Framework/ParticleData/PDGCodeList.h"
#include "Framework/ParticleData/PDGLibrary.h"
#include "Framework/Utils/PrintUtils.h"

using namespace genie;
using namespace genie::geometry;
using namespace genie::controls;

//#define RWH_DEBUG
//#define RWH_DEBUG_2
//#define RWH_COUNTVOLS

#ifdef RWH_COUNTVOLS
// keep some statistics about how many volumes traversed for each box face
long int mxsegments = 0; //rwh
long int curface = 0; //rwh
long int nswims[6]  = { 0, 0, 0, 0, 0, 0}; //rwh
long int nnever[6]  = { 0, 0, 0, 0, 0, 0}; //rwh
double   dnvols[6]  = { 0, 0, 0, 0, 0, 0}; //rwh
double   dnvols2[6] = { 0, 0, 0, 0, 0, 0}; //rwh
bool  accum_vol_stat = false;
#endif

//___________________________________________________________________________
ROOTGeomAnalyzer::ROOTGeomAnalyzer(string geometry_filename)
  : GeomAnalyzerI()
{
///
/// Constructor from a geometry file
///
  LOG("GROOTGeom", pDEBUG)
    << "ROOTGeomAnalyzer ctor \"" << geometry_filename << "\"";
  this->Initialize();
  this->Load(geometry_filename);
}

//___________________________________________________________________________
ROOTGeomAnalyzer::ROOTGeomAnalyzer(TGeoManager * gm)
  : GeomAnalyzerI()
{
///
/// Constructor from a TGeoManager
///
  LOG("GROOTGeom", pDEBUG)
    << "ROOTGeomAnalyzer ctor passed TGeoManager*";
  this->Initialize();
  this->Load(gm);
}

//___________________________________________________________________________
ROOTGeomAnalyzer::~ROOTGeomAnalyzer()
{
  this->CleanUp();

  if ( fmxddist > 0 || fmxdstep > 0 )
    LOG("GROOTGeom",pNOTICE)
      << "ROOTGeomAnalyzer "
      << " mxddist " << fmxddist
      << " mxdstep " << fmxdstep;
}

//===========================================================================
// Geometry driver interface implementation:

//___________________________________________________________________________
const PDGCodeList & ROOTGeomAnalyzer::ListOfTargetNuclei(void)
{
  return *fCurrPDGCodeList;
}

//___________________________________________________________________________
const PathLengthList & ROOTGeomAnalyzer::ComputeMaxPathLengths(void)
{
/// Computes the maximum path lengths for all materials in the input
/// geometry. The computed path lengths are in SI units (kgr/m^2, if
/// density weighting is enabled)

  LOG("GROOTGeom", pNOTICE)
     << "Computing the maximum path lengths for all materials";

  if (!fGeometry) {
      LOG("GROOTGeom", pFATAL) << "No ROOT geometry is loaded!!";
      exit(1);
  }

  //-- initialize max path lengths
  fCurrMaxPathLengthList->SetAllToZero();

  //-- select maximum path length calculation method
  if ( fFlux ) {
    this->MaxPathLengthsFluxMethod();
    // clear any accumulated exposure accounted generated
    // while exploring the geometry
    fFlux->Clear("CycleHistory");
  } else {
    this->MaxPathLengthsBoxMethod();
  }

  return *fCurrMaxPathLengthList;
}

//___________________________________________________________________________
const PathLengthList & ROOTGeomAnalyzer::ComputePathLengths(
                          const TLorentzVector & x, const TLorentzVector & p)
{
/// Computes the path-length within each detector material for a
/// neutrino starting from point x (master coord) and travelling along
/// the direction of p (master coord).
/// The computed path lengths are in SI units (kgr/m^2, if density
/// weighting is enabled)

  //LOG("GROOTGeom", pDEBUG)
  //     << "Computing path-lengths for the input neutrino";

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
  LOG("GROOTGeom", pDEBUG)
       << "\nInput nu: 4p (GeV) = " << utils::print::P4AsShortString(&p)
       << ", 4x (m,s) = " << utils::print::X4AsString(&x);
#endif

  // if trimming configure with neutrino ray's info
  if ( fGeomVolSelector ) {
    fGeomVolSelector->SetCurrentRay(x,p);
    fGeomVolSelector->SetSI2Local(1/this->LengthUnits());
  }

  TVector3 udir = p.Vect().Unit(); // unit vector along direction
  TVector3 pos = x.Vect();         // initial position
  this->SI2Local(pos);             // SI -> curr geom units

  if (!fMasterToTopIsIdentity) {
    this->Master2Top(pos);         // transform position (master -> top)
    this->Master2TopDir(udir);     // transform direction (master -> top)
  }

  // reset current list of path-lengths
  fCurrPathLengthList->SetAllToZero();

  //loop over materials & compute the path-length
  vector<int>::iterator itr;
  for (itr=fCurrPDGCodeList->begin();itr!=fCurrPDGCodeList->end();itr++) {

    int pdgc = *itr;

    Double_t pl = this->ComputePathLengthPDG(pos,udir,pdgc);
    fCurrPathLengthList->AddPathLength(pdgc,pl);

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
    LOG("GROOTGeom", pINFO)
      <<"Calculated path length for material: " << pdgc << " = " << pl;
#endif

  } // loop over materials

  this->Local2SI(*fCurrPathLengthList); // curr geom units -> SI

  return *fCurrPathLengthList;
}

//___________________________________________________________________________
const TVector3 & ROOTGeomAnalyzer::GenerateVertex(
              const TLorentzVector & x, const TLorentzVector & p, int tgtpdg)
{
/// Generates a random vertex, within the detector material with the input
/// PDG code, for a neutrino starting from point x (master coord) and
/// travelling along the direction of p (master coord).

  LOG("GROOTGeom", pNOTICE)
       << "Generating vtx in material: " << tgtpdg
       << " along the input neutrino direction";

  int nretry = 0;
  retry:  // goto label in case of abject failure
  nretry++;

  // reset current interaction vertex
  fCurrVertex->SetXYZ(0.,0.,0.);

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
  LOG("GROOTGeom", pDEBUG)
       << "\nInput nu: 4p (GeV) = " << utils::print::P4AsShortString(&p)
       << ", 4x (m,s) = " << utils::print::X4AsString(&x);
#endif

  if (!fGeometry) {
      LOG("GROOTGeom", pFATAL) << "No ROOT geometry is loaded!!";
      exit(1);
  }

  // calculate the max path length for the selected material starting from
  // x and looking along the direction of p
  TVector3 udir = p.Vect().Unit();
  TVector3 pos = x.Vect();
  this->SI2Local(pos);           // SI -> curr geom units

  if (!fMasterToTopIsIdentity) {
    this->Master2Top(pos);       // transform position (master -> top)
    this->Master2TopDir(udir);   // transform direction (master -> top)
  }

  double maxwgt_dist = this->ComputePathLengthPDG(pos,udir,tgtpdg);
  if ( maxwgt_dist <= 0 ) {
    LOG("GROOTGeom", pERROR)
     << "The current trajectory does not cross the selected material!!";
    return *fCurrVertex;
  }

  // generate random number between 0 and max_dist
  RandomGen * rnd = RandomGen::Instance();
  double genwgt_dist(maxwgt_dist * rnd->RndGeom().Rndm());

  LOG("GROOTGeom", pINFO)
    << "Swim mass: Top Vol dir = " << utils::print::P3AsString(&udir)
    << ", pos = " << utils::print::Vec3AsString(&pos);
  LOG("GROOTGeom", pINFO)
     << "Max {L x Density x Weight} given (init,dir) = " << maxwgt_dist;
  LOG("GROOTGeom", pINFO)
       << "Generated 'distance' in selected material = " << genwgt_dist;
#ifdef RWH_DEBUG
  if ( ( fDebugFlags & 0x01 ) ) {
    fCurrPathSegmentList->SetDoCrossCheck(true);       //RWH
    LOG("GROOTGeom", pINFO) << *fCurrPathSegmentList;  //RWH
    double mxddist = 0, mxdstep = 0;
    fCurrPathSegmentList->CrossCheck(mxddist,mxdstep);
    fmxddist = TMath::Max(fmxddist,mxddist);
    fmxdstep = TMath::Max(fmxdstep,mxdstep);
  }
#endif

  // compute the pdg weight for each material just once, then use a stl map
  PathSegmentList::MaterialMap_t wgtmap;
  PathSegmentList::MaterialMapCItr_t mitr     =
    fCurrPathSegmentList->GetMatStepSumMap().begin();
  PathSegmentList::MaterialMapCItr_t mitr_end =
    fCurrPathSegmentList->GetMatStepSumMap().end();
  // loop over map to get tgt weight for each material (once)
  // steps outside the geometry may have no assigned material
  for ( ; mitr != mitr_end; ++mitr ) {
    const TGeoMaterial* mat = mitr->first;
    double wgt = ( mat ) ? this->GetWeight(mat,tgtpdg) : 0;
    wgtmap[mat] = wgt;
#ifdef RWH_DEBUG
    if ( ( fDebugFlags & 0x02 ) ) {
      LOG("GROOTGeom", pINFO)
        << " wgtmap[" << mat->GetName() << "] pdg " << tgtpdg << " wgt " << Form("%.6f",wgt);
    }
#endif
  }

  // walk down the path to pick the vertex
  const genie::geometry::PathSegmentList::PathSegmentV_t& segments =
    fCurrPathSegmentList->GetPathSegmentV();
  genie::geometry::PathSegmentList::PathSegVCItr_t sitr;
  double walked = 0;
  for ( sitr = segments.begin(); sitr != segments.end(); ++sitr) {
    const genie::geometry::PathSegment& seg = *sitr;
    const TGeoMaterial* mat = seg.fMaterial;
    double trimmed_step = seg.GetSummedStepRange();
    double wgtstep = trimmed_step * wgtmap[mat];
    double beyond = walked + wgtstep;
#ifdef RWH_DEBUG
    if ( ( fDebugFlags & 0x04 ) ) {
      LOG("GROOTGeom", pINFO)
        << " beyond " << beyond << " genwgt_dist " << genwgt_dist
        << " trimmed_step " << trimmed_step << " wgtstep " << wgtstep;
    }
#endif
    if ( beyond > genwgt_dist ) {
      // the end of this segment is beyond our generation point
#ifdef RWH_DEBUG
      if ( ( fDebugFlags & 0x08 ) ) {
        LOG("GROOTGeom", pINFO)
          << "Choose vertex pos walked=" << walked
          << " beyond=" << beyond
          << " wgtstep " << wgtstep
          << " ( " << trimmed_step << "*" << wgtmap[mat] << ")"
          << " look for " << genwgt_dist
          << " in " << seg.fVolume->GetName() << " "
          << mat->GetName();
      }
#endif
      // choose a vertex in this segment (possibly multiple steps)
      double frac = ( genwgt_dist - walked ) / wgtstep;
      if ( frac > 1.0 ) {
        LOG("GROOTGeom", pWARN)
          << "Hey, frac = " << frac << " ( > 1.0 ) "
          << genwgt_dist << " " << walked << " " << wgtstep;
      }
      pos = seg.GetPosition(frac);
      fGeometry -> SetCurrentPoint (pos[0],pos[1],pos[2]);
      fGeometry -> FindNode();
      LOG("GROOTGeom", pINFO)
        << "Choose vertex position in " << seg.fVolume->GetName() << " "
         << utils::print::Vec3AsString(&pos);
      break;
    }
    walked = beyond;
  }

  LOG("GROOTGeom", pNOTICE)
     << "The vertex was placed in volume: "
     << fGeometry->GetCurrentVolume()->GetName()
     << ", path: " << fGeometry->GetPath();

  // warn for any volume overshoots
  bool ok = this->FindMaterialInCurrentVol(tgtpdg);
  if (!ok) {
    LOG("GROOTGeom", pWARN)
       << "Geometry volume was probably overshot";
    LOG("GROOTGeom", pWARN)
       << "No material with code = " << tgtpdg << " could be found at genwgt_dist="
       << genwgt_dist << " (maxwgt_dist=" << maxwgt_dist << ")";
    if ( nretry < 10 ) {
      LOG("GROOTGeom", pWARN)
        << "retry placing vertex";
      goto retry;  // yeah, I know! MyBad.
    }
  }

  if (!fMasterToTopIsIdentity) {
     this->Top2Master(pos); // transform position (top -> master)
  }

  this->Local2SI(pos);   // curr geom units -> SI

  fCurrVertex->SetXYZ(pos[0],pos[1],pos[2]);

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
  LOG("GROOTGeom", pDEBUG)
      << "Vtx (m) = " << utils::print::Vec3AsString(&pos);
#endif

  return *fCurrVertex;
}

//===========================================================================
// Driver configuration methods:

//___________________________________________________________________________
void ROOTGeomAnalyzer::SetLengthUnits(double u)
{
/// Use the units of the input geometry,
///    e.g. SetLengthUnits(genie::units::centimeter)
/// GENIE uses the physical system of units (hbar=c=1) almost throughtout
/// so everything is expressed in GeV but when analyzing detector geometries
/// use meters. Setting input geometry units will allow the code to compute
/// the conversion factor.
/// As input, use one of the constants in $GENIE/src/Conventions/Units.h

  fLengthScale = u/units::meter;
  LOG("GROOTGeom", pNOTICE)
     << "Geometry length units scale factor (geom units -> m): "
     << fLengthScale;
}

//___________________________________________________________________________
void ROOTGeomAnalyzer::SetDensityUnits(double u)
{
/// Like SetLengthUnits, but for density (default units = kgr/m3)

  fDensityScale = u / (units::kilogram / units::meter3);
  LOG("GROOTGeom", pNOTICE)
    << "Geometry density units scale factor (geom units -> kgr/m3): "
    << fDensityScale;
}

//___________________________________________________________________________
void ROOTGeomAnalyzer::SetMaxPlSafetyFactor(double sf)
{
/// Set a factor that can multiply the computed max path lengths.
/// The maximum path lengths are computed by performing an MC scanning of
/// the input geometry. If you configure the scanner with a low number of
/// points or rays you might understimate the path lengths, so you might
/// want to 'inflate' them a little bit using this method.
/// Do not set this number too high, because the max interaction probability
/// will be grossly overestimated and you would need lots of attempts before
/// getting a flux neutrino to interact...

  fMaxPlSafetyFactor = sf;

  LOG("GROOTGeom", pNOTICE)
    << "Max path length safety factor: " << fMaxPlSafetyFactor;
}

//___________________________________________________________________________
void ROOTGeomAnalyzer::SetMixtureWeightsSum(double sum)
{
/// Set it to x, if the relative weight proportions of elements in a mixture
/// add up to x (eg x=1, 100, etc). Set it to a negative value to explicitly
/// compute the correct weight normalization.

  fMixtWghtSum = sum;
}

//___________________________________________________________________________
void ROOTGeomAnalyzer::SetTopVolName(string name)
{
/// Set the name of the top volume.
/// This driver would ask the TGeoManager::GetTopVolume() for the top volume.
/// Use this method for changing this if for example you want to set a smaller
/// volume as the top one so as to generate events only in a specific part of
/// your detector.

  if (name.size() == 0) return;

  fTopVolumeName = name;
  LOG("GROOTGeom",pNOTICE) << "Geometry Top Volume name: " << fTopVolumeName;

  TGeoVolume * gvol = fGeometry->GetVolume(fTopVolumeName.c_str());
  if (!gvol) {
     LOG("GROOTGeom",pWARN) << "Could not find volume: " << name.c_str();
     LOG("GROOTGeom",pWARN) << "Will not change the current top volume";
     fTopVolumeName = "";
     return;
  }

  // Get a matrix connecting coordinates of master and top volumes.
  // The matrix will be used for transforming the coordinates of incoming
  // flux neutrinos & generated interaction vertices.
  // This is needed (in case that the input top volume != master volume)
  // because ROOT always sets the coordinate system origin at the centre of
  // the specified top volume (whereas GENIE assumes that the global reference
  // frame is that of the master volume)

  TGeoIterator next(fGeometry->GetMasterVolume());
  TGeoNode *node;
  TString nodeName, volNameStr;
  const char* volName = fTopVolumeName.c_str();
  while ((node = next())) {
    nodeName = node->GetVolume()->GetName();
    if (nodeName == volName) {
      if (fMasterToTop) delete fMasterToTop;
      fMasterToTop = new TGeoHMatrix(*next.GetCurrentMatrix());
      fMasterToTopIsIdentity = fMasterToTop->IsIdentity();
      break;
    }
  }

  // set volume name
  fTopVolume = gvol;
  fGeometry->SetTopVolume(fTopVolume);
}

//===========================================================================
// Geometry/Unit transforms:

//___________________________________________________________________________
void ROOTGeomAnalyzer::Local2SI(PathLengthList & pl) const
{
/// convert path lengths from current geometry units to SI units
///

  double scaling_factor = this->LengthUnits();
  if (this->WeightWithDensity()) { scaling_factor *= this->DensityUnits(); }

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
  LOG("GROOTGeom", pDEBUG)
    << "Scaling path-lengths from local units -> meters "
    << ((this->WeightWithDensity()) ? "* kgr/m^3" : "")
    << " - scale = " << scaling_factor;
#endif

  PathLengthList::iterator pliter;
  for(pliter = pl.begin(); pliter != pl.end(); ++pliter)
  {
    int pdgc = pliter->first;
    pl.ScalePathLength(pdgc, scaling_factor);
  }
}

//___________________________________________________________________________
void ROOTGeomAnalyzer::Local2SI(TVector3 & vec) const
{
/// convert position vector from current geometry units to SI units
///

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
    LOG("GROOTGeom", pDEBUG)
      << "Position (loc): " << utils::print::Vec3AsString(&vec);
#endif

    vec *= (this->LengthUnits());

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
    LOG("GROOTGeom", pDEBUG)
      << "Position (SI): " << utils::print::Vec3AsString(&vec);
#endif
}

//___________________________________________________________________________
void ROOTGeomAnalyzer::SI2Local(TVector3 & vec) const
{
/// convert position vector from SI units to current geometry units
///
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
    LOG("GROOTGeom", pDEBUG)
      << "Position (SI): " << utils::print::Vec3AsString(&vec);
#endif

    vec *= (1./this->LengthUnits());

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
    LOG("GROOTGeom", pDEBUG)
      << "Position (loc): " << utils::print::Vec3AsString(&vec);
#endif
}

//___________________________________________________________________________
void ROOTGeomAnalyzer::Master2Top(TVector3 & vec) const
{
/// transform the input position vector from the master volume coordinate
/// system to the specified top volume coordinate system, but not
/// change the units.

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
    LOG("GROOTGeom", pDEBUG)
      << "Position (coord:master): " << utils::print::Vec3AsString(&vec);
#endif

    Double_t mast[3], top[3];
    vec.GetXYZ(mast);
    fMasterToTop->MasterToLocal(mast, top);
    vec.SetXYZ(top[0], top[1], top[2]);

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
    LOG("GROOTGeom", pDEBUG)
      << "Position (coord:top): " << utils::print::Vec3AsString(&vec);
#endif
}

//___________________________________________________________________________
void ROOTGeomAnalyzer::Master2TopDir(TVector3 & vec) const
{
/// transform the input direction vector from the master volume coordinate
/// system to the specified top volume coordinate system, but not
/// change the units.

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
    LOG("GROOTGeom", pDEBUG)
      << "Direction (coord:master): " << utils::print::Vec3AsString(&vec);
#endif

    Double_t mast[3], top[3];
    vec.GetXYZ(mast);
    fMasterToTop->MasterToLocalVect(mast, top);
    vec.SetXYZ(top[0], top[1], top[2]);

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
    LOG("GROOTGeom", pDEBUG)
      << "Direction (coord:top): " << utils::print::Vec3AsString(&vec);
#endif
}

//___________________________________________________________________________
void ROOTGeomAnalyzer::Top2Master(TVector3 & vec) const
{
/// transform the input position vector from the specified top volume
/// coordinate system to the master volume coordinate system, but not
/// change the units.

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
    LOG("GROOTGeom", pDEBUG)
      << "Position (coord:top): " << utils::print::Vec3AsString(&vec);
#endif

    Double_t mast[3], top[3];
    vec.GetXYZ(top);
    fMasterToTop->LocalToMaster(top, mast);
    vec.SetXYZ(mast[0], mast[1], mast[2]);

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
    LOG("GROOTGeom", pDEBUG)
      << "Position (coord:master): " << utils::print::Vec3AsString(&vec);
#endif
}

//___________________________________________________________________________
void ROOTGeomAnalyzer::Top2MasterDir(TVector3 & vec) const
{
/// transform the input direction vector from the specified top volume
/// coordinate system to the master volume coordinate system.

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
    LOG("GROOTGeom", pDEBUG)
      << "Direction (coord:top): " << utils::print::Vec3AsString(&vec);
#endif

    Double_t mast[3], top[3];
    vec.GetXYZ(top);
    fMasterToTop->LocalToMasterVect(top, mast);
    vec.SetXYZ(mast[0], mast[1], mast[2]);

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
    LOG("GROOTGeom", pDEBUG)
      << "Direction (coord:master): " << utils::print::Vec3AsString(&vec);
#endif
}

//===========================================================================
// Private methods:

//___________________________________________________________________________
void ROOTGeomAnalyzer::Initialize(void)
{
  LOG("GROOTGeom", pNOTICE)
                << "Initializing ROOT geometry driver & setting defaults";

  fCurrMaxPathLengthList = 0;
  fCurrPathLengthList    = 0;
  fCurrPathSegmentList   = 0;
  fGeomVolSelector       = 0;
  fCurrPDGCodeList       = 0;
  fTopVolume             = 0;
  fTopVolumeName         = "";
  fKeepSegPath           = false;

  // some defaults:
  this -> SetScannerNPoints    (200);
  this -> SetScannerNRays      (200);
  this -> SetScannerNParticles (10000);
  this -> SetScannerFlux       (0);
  this -> SetMaxPlSafetyFactor (1.1);
  this -> SetLengthUnits       (genie::units::meter);
  this -> SetDensityUnits      (genie::units::kilogram/genie::units::meter3);
  this -> SetWeightWithDensity (true);
  this -> SetMixtureWeightsSum (-1.);

  fMasterToTopIsIdentity = true;

  fmxddist = 0;
  fmxdstep = 0;
  fDebugFlags = 0;
}

//___________________________________________________________________________
void ROOTGeomAnalyzer::CleanUp(void)
{
  LOG("GROOTGeom", pNOTICE) << "Cleaning up...";

  if ( fCurrPathSegmentList   ) delete fCurrPathSegmentList;
  if ( fCurrPathLengthList    ) delete fCurrPathLengthList;
  if ( fCurrMaxPathLengthList ) delete fCurrMaxPathLengthList;
  if ( fCurrPDGCodeList       ) delete fCurrPDGCodeList;
  if ( fMasterToTop           ) delete fMasterToTop;
}

//___________________________________________________________________________
void ROOTGeomAnalyzer::Load(string filename)
{
/// Load the detector geometry from the input ROOT file
///
  LOG("GROOTGeom", pNOTICE) << "Loading geometry from: " << filename;

  bool is_accessible = ! (gSystem->AccessPathName( filename.c_str() ));
  if (!is_accessible) {
     LOG("GROOTGeom", pFATAL)
       << "The ROOT geometry doesn't exist! Initialization failed!";
     exit(1);
  }

// ROOT versions 6.16 - 6.24 defaulted to kG4Units [ugh]
#if ROOT_VERSION_CODE >= ROOT_VERSION(6,16,0)
  if (TGeoManager::GetDefaultUnits() != TGeoManager::kRootUnits) {
    TGeoManager::LockDefaultUnits(false);
    TGeoManager::SetDefaultUnits(TGeoManager::kRootUnits);
    TGeoManager::LockDefaultUnits(true);
  }
#endif

  TGeoManager * gm = TGeoManager::Import(filename.c_str());

  this->Load(gm);
}

//___________________________________________________________________________
void ROOTGeomAnalyzer::Load(TGeoManager * gm)
{
/// Load the detector geometry from the input TGeoManager

  LOG("GROOTGeom", pNOTICE)
         << "A TGeoManager is being loaded to the geometry driver";
  fGeometry = gm;

  if (!fGeometry) {
    LOG("GROOTGeom", pFATAL) << "Null TGeoManager! Aborting";
  }
  assert(fGeometry);

  this->BuildListOfTargetNuclei();

  const PDGCodeList & pdglist = this->ListOfTargetNuclei();

  fTopVolume             = 0;
  fCurrPathSegmentList   = new PathSegmentList();
  fCurrPathLengthList    = new PathLengthList(pdglist);
  fCurrMaxPathLengthList = new PathLengthList(pdglist);
  fCurrVertex            = new TVector3(0.,0.,0.);

  // ask geometry manager for its top volume
  fTopVolume = fGeometry->GetTopVolume();
  if (!fTopVolume) {
      LOG("GROOTGeom", pFATAL) << "Could not get top volume!!!";
  }
  assert(fTopVolume);

  // load matrix (identity) of top volume
  fMasterToTop = new TGeoHMatrix(*fGeometry->GetCurrentMatrix());
  fMasterToTopIsIdentity = true;

//#define PRINT_MATERIALS
#ifdef PRINT_MATERIALS
  fGeometry->GetListOfMaterials()->Print();
  fGeometry->GetListOfMedia()->Print();
#endif

}

//___________________________________________________________________________
void ROOTGeomAnalyzer::BuildListOfTargetNuclei(void)
{
/// Determine possible target PDG codes.
/// Note:  If one is using a top volume other than the master level
/// then the final list might include PDG codes that will never
/// be seen during swimming through the volumes if those code are only
/// found in materials outside the top volume.

  fCurrPDGCodeList = new PDGCodeList;

  if (!fGeometry) {
    LOG("GROOTGeom", pFATAL) << "No ROOT geometry is loaded!!";
    exit(1);
  }

  TObjArray * volume_list = fGeometry->GetListOfVolumes();
  if (!volume_list) {
     LOG("GROOTGeom", pERROR)
        << "Null list of geometry volumes. Can not find build target list!";
     return;
  }

  std::set<Int_t> seen_mat; // list of materials we've already handled
  std::vector<TGeoVolume*> volvec; //RWH

  int numVol = volume_list->GetEntries();
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
  LOG("GROOTGeom", pDEBUG) << "Number of volumes found: " << numVol;
#endif

  for (int ivol = 0; ivol < numVol; ivol++) {
      TGeoVolume * volume = dynamic_cast <TGeoVolume *>(volume_list->At(ivol));
      if (!volume) {
         LOG("GROOTGeom", pWARN)
           << "Got a null geometry volume!! Skiping current list element";
         continue;
      }
      TGeoMaterial * mat = volume->GetMedium()->GetMaterial();

      // shortcut if we've already seen this material
      Int_t mat_indx = mat->GetIndex();
      if ( seen_mat.find(mat_indx) != seen_mat.end() ) continue;
      seen_mat.insert(mat_indx);
      volvec.push_back(volume); //RWH

      if (mat->IsMixture()) {
         TGeoMixture * mixt = dynamic_cast <TGeoMixture*> (mat);
         int Nelements = mixt->GetNelements();
         for (int i=0; i<Nelements; i++) {
            int ion_pdgc = this->GetTargetPdgCode(mixt,i);
            fCurrPDGCodeList->push_back(ion_pdgc);
         }
      } else {
          int ion_pdgc = this->GetTargetPdgCode(mat);
          fCurrPDGCodeList->push_back(ion_pdgc);
      }
  }
  // sort the list
  // we don't calculate this list but once per geometry and a sorted
  // list is easier to read so this doesn't cost much
  std::sort(fCurrPDGCodeList->begin(),fCurrPDGCodeList->end());

}

//___________________________________________________________________________
int ROOTGeomAnalyzer::GetTargetPdgCode(const TGeoMaterial * const m) const
{
  int A = TMath::Nint(m->GetA());
  int Z = TMath::Nint(m->GetZ());

  int pdgc = pdg::IonPdgCode(A,Z);

  return pdgc;
}

//___________________________________________________________________________
int ROOTGeomAnalyzer::GetTargetPdgCode(
                     const TGeoMixture * const m, int ielement) const
{
  int A = TMath::Nint(m->GetAmixt()[ielement]);
  int Z = TMath::Nint(m->GetZmixt()[ielement]);

  int pdgc = pdg::IonPdgCode(A,Z);

  return pdgc;
}

//___________________________________________________________________________
double ROOTGeomAnalyzer::GetWeight(const TGeoMaterial * mat, int pdgc)
{
/// Get the weight of the input material.
/// Return the weight only if the material's pdg code matches the input code.
/// If the material is found to be a mixture, call the corresponding method
/// for mixtures.
/// Weight is in the curr geom density units.

  if (!mat) {
    LOG("GROOTGeom", pERROR) << "Null input material. Return weight = 0.";
    return 0;
  }

  bool exists = fCurrPDGCodeList->ExistsInPDGCodeList(pdgc);
  if (!exists) {
    LOG("GROOTGeom", pERROR) << "Target doesn't exist. Return weight = 0.";
    return 0;
  }

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
  LOG("GROOTGeom",pDEBUG)
       << "Curr. material: A/Z = " << mat->GetA() << " / " << mat->GetZ();
#endif

  // if the input material is a mixture, get a the sum of weights for
  // all matching elements
  double weight = 0.;
  if (mat->IsMixture()) {
    const TGeoMixture * mixt = dynamic_cast <const TGeoMixture*> (mat);
    if (!mixt) {
     LOG("GROOTGeom", pERROR) << "Null input mixture. Return weight = 0.";
     return 0;
    }
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
    LOG("GROOTGeom", pDEBUG)
      << "Material : " << mat->GetName()
          << " is a mixture with " << mixt->GetNelements() << " elements";
#endif
    // loop over elements & sum weights of matching elements
    weight = this->GetWeight(mixt,pdgc);
    return weight;
  } // is mixture?

  // pure material
  int ion_pdgc = this->GetTargetPdgCode(mat);
  if (ion_pdgc != pdgc) return 0.;

  if (this->WeightWithDensity())
    weight = mat->GetDensity(); // material density (curr geom units)
  else
    weight = 1.0;

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
  LOG("GROOTGeom", pDEBUG)
       << "Weight[mat:" << mat->GetName() << "] = " << weight;
#endif

  return weight;
}

//___________________________________________________________________________
double ROOTGeomAnalyzer::GetWeight(const TGeoMixture * mixt, int pdgc)
{
/// Loop over the mixture elements, find the one matching the input pdgc
/// and  return its weight.
/// Weight is in the curr geom density units.

  double weight = 0;

  int nm = 0;
  for (int i = 0; i < mixt->GetNelements(); i++) {
     double dw = (this->GetWeight(mixt,i,pdgc));
     if (dw>0) nm++;
     weight += dw;
  }

  if (nm>1) {
     for (int j = 0; j < mixt->GetNelements(); j++) {
           LOG("GROOTGeom", pWARN)
              << "[" << j << "] Z = " << mixt->GetZmixt()[j]
              << ", A = " << mixt->GetAmixt()[j]
              << " (pdgc = " << this->GetTargetPdgCode(mixt,j)
              << "), w = " << mixt->GetWmixt()[j];
     }
     LOG("GROOTGeom", pERROR)
        << "Material pdgc = " << pdgc << " appears " << nm
        << " times (>1) in mixture = " << mixt->GetName();
     LOG("GROOTGeom", pFATAL)
        << "Your geometry must be incorrect - Aborting";
     exit(1);
  }

  // if we are not weighting with the density then the weight=1 if the pdg
  // code was matched for any element of this mixture
  if ( !this->WeightWithDensity() && weight>0. ) weight=1.0;

  return weight;
}

//___________________________________________________________________________
double ROOTGeomAnalyzer::GetWeight(const TGeoMixture* mixt, int ielement, int pdgc)
{
/// Get the weight of the input ith element of the input material.
/// Return the weight only if the element's pdg code matches the input code.
/// Weight is in the curr geom density units.

//  int ion_pdgc = this->GetTargetPdgCode(mixt->GetElement(ielement));
  int ion_pdgc = this->GetTargetPdgCode(mixt, ielement);
  if (ion_pdgc != pdgc) return 0.;

  double d = mixt->GetDensity();         // mixture density (curr geom units)
  double w = mixt->GetWmixt()[ielement]; // relative proportion by mass

  double wtot = this->MixtureWeightsSum();

  // <0 forces explicit calculation of relative proportion normalization
  if (wtot < 0) {
    wtot = 0;
    for (int i = 0; i < mixt->GetNelements(); i++) {
      wtot += (mixt->GetWmixt()[i]);
    }
  }
  assert(wtot>0);

  w /= wtot;
  double weight = d*w;

  return weight;
}

//___________________________________________________________________________
void ROOTGeomAnalyzer::MaxPathLengthsFluxMethod(void)
{
/// Use the input flux driver to generate "rays", and then follow them through
/// the detector and figure out the maximum path length for each material

  LOG("GROOTGeom", pNOTICE)
               << "Computing the maximum path lengths using the FLUX method";

  int iparticle = 0;
  PathLengthList::const_iterator pl_iter;

  const int nparticles = abs(this->ScannerNParticles());

  // if # scanner particles is negative, this signals that the user
  // desires to force rays to have the maximum energy (useful if the
  // volume considered changes size with neutrino energy)
  bool rescale_e = (this->ScannerNParticles() < 0 );
  double emax = fFlux->MaxEnergy();
  if ( rescale_e ) {
    LOG("GROOTGeom", pNOTICE)
      << "max path lengths with FLUX method forcing Enu=" << emax;
  }

  while (iparticle < nparticles ) {

    bool ok = fFlux->GenerateNext();
    if (!ok) {
       LOG("GROOTGeom", pWARN) << "Couldn't generate a flux neutrino";
       continue;
    }

    TLorentzVector   nup4  = fFlux->Momentum();
    if ( rescale_e ) {
      double ecurr = nup4.E();
      if ( ecurr > 0 ) nup4 *= (emax/ecurr);
    }
    const TLorentzVector & nux4  = fFlux->Position();

    //LOG("GMCJDriver", pNOTICE)
    //   << "\n [-] Generated flux neutrino: "
    //   << "\n  |----o 4-momentum : " << utils::print::P4AsString(&nup4)
    //   << "\n  |----o 4-position : " << utils::print::X4AsString(&nux4);

    const PathLengthList & pl = this->ComputePathLengths(nux4, nup4);

    bool enters = false;

    for (pl_iter = pl.begin(); pl_iter != pl.end(); ++pl_iter) {
       int    pdgc       = pl_iter->first;
       double pathlength = pl_iter->second;

       if ( pathlength > 0 ) {
          pathlength *= (this->MaxPlSafetyFactor());

          pathlength = TMath::Max(pathlength, fCurrMaxPathLengthList->PathLength(pdgc));
          fCurrMaxPathLengthList->SetPathLength(pdgc,pathlength);
          enters = true;
       }
    }
    if (enters) iparticle++;
  }
}

//___________________________________________________________________________
void ROOTGeomAnalyzer::MaxPathLengthsBoxMethod(void)
{
/// Generate points in the geometry's bounding box and for each point
/// generate random rays, follow them through the detector and figure out
/// the maximum path length for each material

  LOG("GROOTGeom", pNOTICE)
    << "Computing the maximum path lengths using the BOX method";
#ifdef RWH_COUNTVOLS
  accum_vol_stat = true;
#endif

  int  iparticle = 0;
  bool ok = true;
  TLorentzVector nux4;
  TLorentzVector nup4;

  PathLengthList::const_iterator pl_iter;

  while ( (ok = this->GenBoxRay(iparticle++,nux4,nup4)) ) {

    //LOG("GMCJDriver", pNOTICE)
    //  << "\n [-] Generated flux neutrino: "
    //  << "\n  |----o 4-momentum : " << utils::print::P4AsString(&nup4)
    //  << "\n  |----o 4-position : " << utils::print::X4AsString(&nux4);

    const PathLengthList & pllst = this->ComputePathLengths(nux4, nup4);

    for (pl_iter = pllst.begin(); pl_iter != pllst.end(); ++pl_iter) {
       int    pdgc = pl_iter->first;
       double pl   = pl_iter->second;

       if (pl>0) {
          pl *= (this->MaxPlSafetyFactor());

          pl = TMath::Max(pl, fCurrMaxPathLengthList->PathLength(pdgc));
          fCurrMaxPathLengthList->SetPathLength(pdgc,pl);
       }
    }
  }

  // print out the results
  LOG("GROOTGeom", pDEBUG)
    << "DensWeight \"" << (fDensWeight?"true":"false")
    << "\" MixtWghtSum " << fMixtWghtSum;
  LOG("GROOTGeom", pDEBUG) << "CurrMaxPathLengthList: "
    << *fCurrMaxPathLengthList;

#ifdef RWH_COUNTVOLS
  // rwh
  // print statistics for average,rms of number of volumes seen for
  // various rays for each face
  for (int j = 0; j < 6; ++j ) {
    long int ns = nswims[j];
    double   x  = dnvols[j];
    double   x2 = dnvols2[j];
    if ( ns == 0 ) ns = 1;
    double avg = x / (double)ns;
    double rms = TMath::Sqrt((x2/(double)ns) - avg*avg);
    LOG("GROOTGeom", pNOTICE)
      << "RWH: nswim after BOX face " << j << " is " << ns
      << " avg " << avg << " rms " << rms
      << " never " << nnever[j];
  }
  LOG("GROOTGeom", pNOTICE)
    << "RWH: Max PathSegmentList size " << mxsegments;
  accum_vol_stat = false;
#endif

}

//___________________________________________________________________________
bool ROOTGeomAnalyzer::GenBoxRay(int indx, TLorentzVector& x4, TLorentzVector& p4)
{
/// Generate points in the geometry's bounding box and for each point generate
/// random rays -- a pseudo-flux -- in master coordinates and SI units

  firay++;
  fnewpnt = false;

  // first time through ... special case
  if ( indx == 0 )           { fiface = 0; fipoint = 0; firay = 0; fnewpnt = true; }

  if ( firay   >= fNRays   ) {             fipoint++;   firay = 0; fnewpnt = true; }
  if ( fipoint >= fNPoints ) { fiface++;   fipoint = 0; firay = 0; fnewpnt = true; }
  if ( fiface  >= 3 ) {
    LOG("GROOTGeom",pINFO) << "Box surface scanned: " << indx << " points/rays";
    return false; // signal end
  }

  if ( indx == 0 ) {
    // get geometry's bounding box
    //LOG("GROOTGeom", pNOTICE) << "Getting a TGeoBBox enclosing the detector";
    TGeoShape * TS  = fTopVolume->GetShape();
    TGeoBBox *  box = (TGeoBBox *)TS;
    //get box origin and dimensions (in the same units as the geometry)
    fdx = box->GetDX(); // half-length
    fdy = box->GetDY(); // half-length
    fdz = box->GetDZ(); // half-length
    fox = (box->GetOrigin())[0];
    foy = (box->GetOrigin())[1];
    foz = (box->GetOrigin())[2];

    LOG("GROOTGeom",pNOTICE)
      << "Box size (GU)   :"
      << " x = " << 2*fdx << ", y = " << 2*fdy << ", z = " << 2*fdz;
    LOG("GROOTGeom",pNOTICE)
      << "Box origin (GU) :"
      << " x = " << fox   << ", y = " << foy   << ", z = " <<   foz;
    LOG("GROOTGeom",pNOTICE)
      << "Will generate [" << fNPoints << "] random points / box surface";
    LOG("GROOTGeom",pNOTICE)
      << "Will generate [" << fNRays   << "] rays / point";

#ifdef VALIDATE_CORNERS
    // RWH: test that we know the coordinate transforms for the box corners
    for (int sz = -1; sz <= +1; ++sz) {
      for (int sy = -1; sy <= +1; ++sy) {
        for (int sx = -1; sx <= +1; ++sx) {
          if (sx == 0 || sy == 0 || sz == 0 ) continue;
          TVector3& pos = fGenBoxRayPos;
          pos.SetXYZ(fox+(sx*fdx),foy+(sy*fdy),foz+(sz*fdz));
          TVector3 master(fGenBoxRayPos);
          this->Top2Master(master);   // transform position (top -> master)
          this->Local2SI(master);
          TVector3 pos2(master);
          this->Master2Top(pos2);
          this->SI2Local(pos2);
          LOG("GROOTGeom", pNOTICE)
            << "TopVol corner "
            << " [" << pos[0] << "," << pos[1] << "," << pos[2] << "] "
            << "Master corner "
            << " [" << master[0] << "," << master[1] << "," << master[2] << "] "
            << " top again"
            << " [" << pos2[0] << "," << pos2[1] << "," << pos2[2] << "] ";
        }
      }
    }
#endif

  }

  RandomGen * rnd = RandomGen::Instance();

  double eps = 0.0; //1.0e-12; // tweak to make sure we're inside the box

  switch ( fiface ) {
  case 0:
    {
      //top:
      if ( firay == fNRays-1 && fipoint == fNPoints-1 )
        LOG("GROOTGeom",pINFO) << "Box surface scanned: [TOP]";
      fGenBoxRayDir.SetXYZ(-0.5+rnd->RndGeom().Rndm(),
                           -rnd->RndGeom().Rndm(),
                           -0.5+rnd->RndGeom().Rndm());
      if ( fnewpnt )
        fGenBoxRayPos.SetXYZ(fox-fdx+eps+2*(fdx-eps)*rnd->RndGeom().Rndm(),
                             foy+fdy-eps,
                             foz-fdz+eps+2*(fdz-eps)*rnd->RndGeom().Rndm());
      break;
    }
  case 1:
    {
      //left:
      if ( firay == fNRays-1 && fipoint == fNPoints-1 )
        LOG("GROOTGeom",pINFO) << "Box surface scanned: [LEFT]";
      fGenBoxRayDir.SetXYZ(rnd->RndGeom().Rndm(),
                           -0.5+rnd->RndGeom().Rndm(),
                           -0.5+rnd->RndGeom().Rndm());
      if ( fnewpnt )
        fGenBoxRayPos.SetXYZ(fox-fdx+eps,
                             foy-fdy+eps+2*(fdy-eps)*rnd->RndGeom().Rndm(),
                             foz-fdz+eps+2*(fdz-eps)*rnd->RndGeom().Rndm());
      break;
    }
  case 2:
    {
      //front: (really what I, RWH, would call the back)
      if ( firay == fNRays-1 && fipoint == fNPoints-1 )
        LOG("GROOTGeom",pINFO) << "Box surface scanned: [FRONT]";
      fGenBoxRayDir.SetXYZ(-0.5+rnd->RndGeom().Rndm(),
                           -0.5+rnd->RndGeom().Rndm(),
                           -rnd->RndGeom().Rndm());
      if ( fnewpnt )
        fGenBoxRayPos.SetXYZ(fox-fdx+eps+2*(fdx-eps)*rnd->RndGeom().Rndm(),
                             foy-fdy+eps+2*(fdy-eps)*rnd->RndGeom().Rndm(),
                             foz+fdz+eps);
      break;
    }
  }
/*
    //bottom:
      pos.SetXYZ(ox-dx+2*dx*rnd->RndGeom().Rndm(),
                 oy-dy,
                 oz-dz+2*dz*rnd->RndGeom().Rndm());
      dir.SetXYZ(-0.5+rnd->RndGeom().Rndm(),
                 rnd->RndGeom().Rndm(),
                 -0.5+rnd->RndGeom().Rndm());
    //right:
      pos.SetXYZ(ox+dx,
                 oy-dy+2*dy*rnd->RndGeom().Rndm(),
                 oz-dz+2*dz*rnd->RndGeom().Rndm());
      dir.SetXYZ(-rnd->RndGeom().Rndm(),
                 -0.5+rnd->RndGeom().Rndm(),
                 -0.5+rnd->RndGeom().Rndm());
    //back:
      pos.SetXYZ(ox-dx+2*dx*rnd->RndGeom().Rndm(),
                 oy-dy+2*dy*rnd->RndGeom().Rndm(),
                 oz-dz);
      dir.SetXYZ(-0.5+rnd->RndGeom().Rndm(),
                 -0.5+rnd->RndGeom().Rndm(),
                  rnd->RndGeom().Rndm());
*/

#ifdef RWH_COUNTVOLS
  curface = fiface;
#endif

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
  if ( fnewpnt )
    LOG("GROOTGeom", pNOTICE)
      << "GenBoxRay(topvol) "
      << " iface " << fiface << " ipoint " << fipoint << " iray " << firay
      << " newpnt " << (fnewpnt?"true":"false")
      << " x " << utils::print::Vec3AsString(&fGenBoxRayPos)
      << " p " << utils::print::P3AsString(&fGenBoxRayDir);
#endif

  if ( fnewpnt ) {
    if ( ! fMasterToTopIsIdentity) {
      this->Top2Master(fGenBoxRayPos); // transform position (top -> master)
    }
    this->Local2SI(fGenBoxRayPos);
  }
  this->Top2MasterDir(fGenBoxRayDir);  // transform direction (top -> master)

  x4.SetVect(fGenBoxRayPos);
  p4.SetVect(fGenBoxRayDir.Unit());

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
  LOG("GROOTGeom", pNOTICE)
    << "GenBoxRay(master) "
    << " iface " << fiface << " ipoint " << fipoint << " iray " << firay
    << " newpnt " << (fnewpnt?"true":"false")
    << " x " << utils::print::Vec3AsString(&fGenBoxRayPos)
    << " p " << utils::print::P3AsString(&fGenBoxRayDir);
#endif

  return true;
}

//________________________________________________________________________
double ROOTGeomAnalyzer::ComputePathLengthPDG(
                  const TVector3 & r0, const TVector3 & udir, int pdgc)
{
/// Compute the path length for the material with pdg-code = pdc, staring
/// from the input position r (top vol coord & units) and moving along the
/// direction of the unit vector udir (top vol coord).

  double pl = 0; // path-length (x density, if density-weighting is ON)

  this->SwimOnce(r0,udir);

  double step   = 0;
  double weight = 0;

  //  const TGeoVolume   * vol = 0;
  //  const TGeoMedium   * med = 0;
  const TGeoMaterial * mat = 0;

  // loop over independent materials, which is shorter or equal to # of volumes
  PathSegmentList::MaterialMapCItr_t itr     =
    fCurrPathSegmentList->GetMatStepSumMap().begin();
  PathSegmentList::MaterialMapCItr_t itr_end =
    fCurrPathSegmentList->GetMatStepSumMap().end();
  for ( ; itr != itr_end; ++itr ) {
    mat  = itr->first;
    if ( ! mat ) continue;  // segment outside geometry has no material
    step = itr->second;
    weight = this->GetWeight(mat,pdgc);
    pl += (step*weight);
  }

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
  LOG("GROOTGeom", pDEBUG)
       << "PathLength[" << pdgc << "] = " << pl << " in curr geom units";
#endif

  return pl;
}

//________________________________________________________________________
void ROOTGeomAnalyzer::SwimOnce(const TVector3 & r0, const TVector3 & udir)
{
/// Swim through the geometry from the from the input position
/// r0 (top vol coord & units) and moving along the direction of the
/// unit vector udir (topvol coord) to create a filled PathSegmentList

  int nvolswim = 0; //rwh

  if ( ! fCurrPathSegmentList ) fCurrPathSegmentList = new PathSegmentList();

  // don't swim if the current PathSegmentList is up-to-date
  if ( fCurrPathSegmentList->IsSameStart(r0,udir) ) return;

  // start fresh
  fCurrPathSegmentList->SetAllToZero();

  // set start info so next time we don't swim for the same ray
  fCurrPathSegmentList->SetStartInfo(r0,udir);

  PathSegment ps_curr;

  bool found_vol (false);
  bool keep_on   (true);

  double step    = 0;
  double raydist = 0;

  const TGeoVolume *   vol = 0;
  const TGeoMedium *   med = 0;
  const TGeoMaterial * mat = 0;

  // determining the geometry path is expensive, do it only if necessary
  bool selneedspath = ( fGeomVolSelector && fGeomVolSelector->GetNeedPath() );
  const bool fill_path = fKeepSegPath || selneedspath;

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
  LOG("GROOTGeom", pNOTICE)
    << "SwimOnce x [" << r0[0] << "," << r0[1] << "," << r0[2]
    << "] udir [" << udir[0] << "," << udir[1] << "," << udir[2];
#endif

  fGeometry -> SetCurrentDirection (udir[0],udir[1],udir[2]);
  fGeometry -> SetCurrentPoint     (r0[0],  r0[1],  r0[2]  );

  while (!found_vol || keep_on) {
     keep_on = true;

     fGeometry->FindNode();

     ps_curr.SetEnter( fGeometry->GetCurrentPoint() , raydist );
     vol = fGeometry->GetCurrentVolume();
     med = vol->GetMedium();
     mat = med->GetMaterial();
     ps_curr.SetGeo(vol,med,mat);
#ifdef PATHSEG_KEEP_PATH
     if (fill_path) ps_curr.SetPath(fGeometry->GetPath());
#endif

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
#ifdef DUMP_SWIM
       LOG("GROOTGeom", pDEBUG) << "Current volume: " << vol->GetName()
                             << " pos " << fGeometry->GetCurrentPoint()[0]
                             << " "     << fGeometry->GetCurrentPoint()[1]
                             << " "     << fGeometry->GetCurrentPoint()[2]
                             << " dir " << fGeometry->GetCurrentDirection()[0]
                             << " "     << fGeometry->GetCurrentDirection()[1]
                             << " "     << fGeometry->GetCurrentDirection()[2]
                             << "[path: " << fGeometry->GetPath() << "]";
#endif
#endif

     // find the start of top
     if (fGeometry->IsOutside() || !vol) {
        keep_on = false;
        if (found_vol) break;
        step = 0;
          this->StepToNextBoundary();
        //rwh//raydist += step;  // STNB doesn't actually "step"

#ifdef RWH_DEBUG
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
        LOG("GROOTGeom", pDEBUG) << "Outside ToNextBoundary step: " << step
                                 << " raydist: " << raydist;
#endif
#endif

        while (!fGeometry->IsEntering()) {
          step = this->Step();
          raydist += step;
#ifdef RWH_DEBUG
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
          LOG("GROOTGeom", pDEBUG)
              << "Stepping... [step size = " << step << "]";
          LOG("GROOTGeom", pDEBUG) << "Outside step: " << step
                                   << " raydist: " << raydist;
#endif
#endif
          if (this->WillNeverEnter(step)) {
#ifdef RWH_COUNTVOLS
            if ( accum_vol_stat ) {
              // this really shouldn't happen for the box exploration...
              // if coord transforms done right
              // could happen for neutrinos on a flux window
              nnever[curface]++;  //rwh
              if ( nnever[curface]%21 == 0 )
                LOG("GROOTGeom", pNOTICE)
                  << "curface " << curface << " " << nswims[curface]
                  << " never " << nnever[curface]
                  << " x [" << r0[0] << "," << r0[1] << "," << r0[2] << "] "
                  << " p [" << udir[0] << "," << udir[1] << "," << udir[2] << "]";
            }
#endif
            fCurrPathSegmentList->SetAllToZero();
            return;
          }
        } // finished while

        ps_curr.SetExit(fGeometry->GetCurrentPoint());
        ps_curr.SetStep(step);
        if ( ( fDebugFlags & 0x10 ) ) {
          // In general don't add the path segments from the start point to
          // the top volume (here for debug purposes)
          // Clear out the step range even if we keep it
          ps_curr.fStepRangeSet.clear();
          LOG("GROOTGeom", pNOTICE)
            << "debug: step towards top volume: " << ps_curr;
          fCurrPathSegmentList->AddSegment(ps_curr);
        }

     }  // outside or !vol

     if (keep_on) {
       if (!found_vol) found_vol = true;

       step   = this->StepUntilEntering();
       raydist += step;

       ps_curr.SetExit(fGeometry->GetCurrentPoint());
       ps_curr.SetStep(step);
       fCurrPathSegmentList->AddSegment(ps_curr);

       nvolswim++; //rwh

#ifdef DUMP_SWIM
       LOG("GROOTGeom", pDEBUG) << "Current volume: " << vol->GetName()
                                << " step " << step << " in " << mat->GetName();
#endif

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
       LOG("GROOTGeom", pDEBUG)
         << "Cur med.: " << med->GetName() << ", mat.: " << mat->GetName();
       LOG("GROOTGeom", pDEBUG)
         << "Step = " << step; // << ", weight = " << weight;
#endif
     } //keep on
  }

#ifdef RWH_COUNTVOLS
  if ( accum_vol_stat ) {
    nswims[curface]++;   //rwh
    dnvols[curface]  += (double)nvolswim;
    dnvols2[curface] += (double)nvolswim * (double)nvolswim;
    long int ns = fCurrPathSegmentList->size();
    if ( ns > mxsegments ) mxsegments = ns;
  }
#endif

//rwh:debug
#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
  LOG("GROOTGeom", pDEBUG)
    << "PathSegmentList size " << fCurrPathSegmentList->size();
#endif

#ifdef RWH_DEBUG_2
  if ( ( fDebugFlags & 0x20 ) ) {
    fCurrPathSegmentList->SetDoCrossCheck(true);       //RWH
    LOG("GROOTGeom", pNOTICE) << "Before trimming" << *fCurrPathSegmentList;
    double mxddist = 0, mxdstep = 0;
    fCurrPathSegmentList->CrossCheck(mxddist,mxdstep);
    fmxddist = TMath::Max(fmxddist,mxddist);
    fmxdstep = TMath::Max(fmxdstep,mxdstep);
  }
#endif

  // PathSegmentList trimming occurs here!
  if ( fGeomVolSelector ) {
    PathSegmentList* altlist =
      fGeomVolSelector->GenerateTrimmedList(fCurrPathSegmentList);
    std::swap(altlist,fCurrPathSegmentList);
    delete altlist;  // after swap delete original
  }

  fCurrPathSegmentList->FillMatStepSum();

#ifdef RWH_DEBUG_2
  if ( fGeomVolSelector) {
    // after FillMatStepSum() so one can see the summed mass
    if ( ( fDebugFlags & 0x40 ) ) {
      fCurrPathSegmentList->SetPrintVerbose(true);
      LOG("GROOTGeom", pNOTICE) << "After  trimming" << *fCurrPathSegmentList;
      fCurrPathSegmentList->SetPrintVerbose(false);
    }
  }
#endif


  return;
}

//___________________________________________________________________________
bool ROOTGeomAnalyzer::FindMaterialInCurrentVol(int tgtpdg)
{
  TGeoVolume * vol = fGeometry -> GetCurrentVolume();
  if(vol) {
    TGeoMaterial * mat = vol->GetMedium()->GetMaterial();
    if(mat->IsMixture()) {
      TGeoMixture * mixt = dynamic_cast <TGeoMixture*> (mat);
      for(int i = 0; i < mixt->GetNelements(); i++) {
         int pdg = this->GetTargetPdgCode(mixt, i);
         if(tgtpdg == pdg) return true;
      }
    } else {
       int pdg = this->GetTargetPdgCode(mat);
       if(tgtpdg == pdg) return true;
    }
  } else {
     LOG("GROOTGeom", pWARN) << "Current volume is null!";
     return false;
  }
  return false;
}
//___________________________________________________________________________
double ROOTGeomAnalyzer::StepToNextBoundary(void)
{
  fGeometry->FindNextBoundary();
  double step=fGeometry->GetStep();
  return step;
}
//___________________________________________________________________________
double ROOTGeomAnalyzer::Step(void)
{
  fGeometry->Step();
  double step=fGeometry->GetStep();
  return step;
}
//___________________________________________________________________________
double ROOTGeomAnalyzer::StepUntilEntering(void)
{
  this->StepToNextBoundary();  // doesn't actually step, so don't include in sum
  double step = 0; //

  while(!fGeometry->IsEntering()) {
    step += this->Step();
  }

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__

  bool isen = fGeometry->IsEntering();
  bool isob = fGeometry->IsOnBoundary();

  LOG("GROOTGeom",pDEBUG)
      << "IsEntering = "     << utils::print::BoolAsYNString(isen)
      << ", IsOnBoundary = " << utils::print::BoolAsYNString(isob)
      << ", Step = " << step;
#endif

  return step;
}
//___________________________________________________________________________
bool ROOTGeomAnalyzer::WillNeverEnter(double step)
{
// If the neutrino trajectory would never enter the detector, then the
// TGeoManager::GetStep returns the maximum step (1E30).
// Compare surrent step with max step and figure out whether the particle
// would never enter the detector

  if (step > 9.99E29) {

#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
     LOG("GROOTGeom", pINFO) << "Wow! Current step is dr = " << step;
     LOG("GROOTGeom", pINFO) << "This trajectory isn't entering the detector";
#endif
     return true;

  } else
    return false;
}

//___________________________________________________________________________