// $Id: DemoGeometryTools.cc,v 1.10 2015/07/23 20:47:54 kelsey Exp $
////////////////////////////////////////////////////////////////////////
//  File:        DemoGeometryTools.cc                                 //
//                                                                    //
//  Description: Utility class to search GEANT4 geometry structure    //
//		 Adapted from Pedro Arce, GAMOS/GmGeometryUtils       //
//                                                                    //
//  Author:      Michael Kelsey (SLAC)                                //
//  Date:        30 September 2011                                    //
//                                                                    //
//  20120217  M. Kelsey -- Add interface to take name and copy no for //
//		physical volume search.                               //
//  20120402  M. Kelsey -- Add static delimiter character for parsing //
//		Use '|' instead of '/' to allow latter within names.  //
//  20140627  M. Kelsey -- Change fn arg "name" to avoid shadowing.   //
//  20150723  M. Kelsey -- Add functions to get touchable by position //
//  20160204  M. Kelsey -- Can't use typedef name as variable name!   //
//  20160301  M. Kelsey -- Add functions to map LVs to PVs, and to    //
//		get touchables using LV name.                         //
//  20190125  M. Kelsey -- Use local G4Navigator for touchable search //
//  20190128  Don't fill Navigator world volume in BuildDictionaries. //
//  20190614  Replace ::exit() with DemoAbortJob.                     //
//  20200507  M. Kelsey -- Convert to thread-local singleton for MT.  //
////////////////////////////////////////////////////////////////////////

#include "globals.hh"
#include "Demog4base/CDMSGeometryTools.hh"
#include "Demog4base/CDMSAbortJob.hh"
#include "Demog4base/CDMSTouchable.hh"
#include "G4LogicalVolume.hh"
#include "G4LogicalVolumeStore.hh"
#include "G4Material.hh"
#include "G4Navigator.hh"
#include "G4PhysicalVolumeStore.hh"
#include "G4RegionStore.hh"
#include "G4TouchableHistory.hh"
#include "G4Track.hh"
#include "G4TransportationManager.hh"
#include "G4UserLimits.hh"
#include "G4VPhysicalVolume.hh"
#include "G4VSensitiveDetector.hh"
#include "G4VSolid.hh"
#include "G4VisAttributes.hh"
#include <sstream>

extern "C" {
#include <fnmatch.h>	/* Does wildcard "glob-style" string match */
}


G4ThreadLocal DemoGeometryTools* CDMSGeometryTools::theInstance = 0;

const char DemoGeometryTools::delimiter = '|';	// Separates hierarchy names

// Constructor and destructor

DemoGeometryTools* CDMSGeometryTools::GetInstance() {
  if (!theInstance) theInstance = new DemoGeometryTools;
  return theInstance;
}

DemoGeometryTools::CDMSGeometryTools() 
  : verboseLevel(0), theTopPV(0), theNavigator(0), theTouchable(0) {
  BuildDictionaries();
}
 
DemoGeometryTools::~CDMSGeometryTools() {
  delete theNavigator; theNavigator=0;
  delete theTouchable; theTouchable=0;
}
 

// Report current geometry analysis

void DemoGeometryTools::DumpSummary(std::ostream& out) {
  out << " DemoGeometryTools Summary of Current GEANT4 Geometry " << G4endl;

  theTopPV = GetTopPV();
  if (theTopPV == 0) {
    out << " No volume created " << G4endl;
    return;
  }

  out << " Geometry built inside world volume: " << theTopPV->GetName()
      << G4endl;

  // Access geometry tables for information
  const G4PhysicalVolumeStore* pvs = G4PhysicalVolumeStore::GetInstance();
  const G4LogicalVolumeStore* lvs  = G4LogicalVolumeStore::GetInstance();
  const G4MaterialTable* matTab    = G4Material::GetMaterialTable();

  // Get number of solids (< # LV if several LV share a solid)
  vlv::const_iterator lvcite;
  std::set<G4VSolid*> theSolids;
  for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) {
    theSolids.insert((*lvcite)->GetSolid());
  }

  // Report counts of each type of volume
  out << " Number of G4VSolids: " << theSolids.size()
      << "\n Number of G4LogicalVolumes: " << lvs->size()
      << "\n Number of G4VPhysicalVolumes: " << pvs->size()
      << "\n Number of Touchables: " << NumberOfTouchables()
      << "\n Number of G4Materials: " << matTab->size() << G4endl;
}
 
void DemoGeometryTools::DumpG4LVList(std::ostream& out) {
  out << " @@@@@@@@@@@@@@@@ G4LogicalVolume List  " << G4endl;
  const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
  vlv::const_iterator lvcite;
  for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) {
    out << "LV: " << (*lvcite)->GetName() << G4endl;
  }

}

void DemoGeometryTools::DumpG4LVTree(std::ostream& out) {
  out << " @@@@@@@@@@@@@@@@ DUMPING G4LogicalVolume Tree  " << G4endl;
  G4LogicalVolume* lv = GetTopLV();
  DumpG4LVLeaf(lv, 0, out);
}


// Dump list of registered materials

void DemoGeometryTools::DumpMaterialList(std::ostream& out) {
  out << " @@@@@@@@@@@@@@@@ DUMPING G4Material List ";
  const G4MaterialTable* matTab = G4Material::GetMaterialTable();
  out << " with " << matTab->size() << " materials " << G4endl;
  std::vector<G4Material*>::const_iterator matite;
  for (matite = matTab->begin(); matite != matTab->end(); matite++) {
    out << "Material: " << (*matite) << G4endl;
  }

}


// Report current Logical Volume recursively

void 
DemoGeometryTools::DumpG4LVLeaf(G4LogicalVolume* lv, unsigned int leafDepth,
				std::ostream& out) {
  out << IndentLeaf(leafDepth) << " LV:(" << leafDepth << ") "
      << lv->GetName() << G4endl;

  unsigned int siz = lv->GetNoDaughters();
  //--- Count duplicated leaves only once for reporting
  std::set< G4LogicalVolume* > lvDaughters;
  for (unsigned int ii = 0; ii < siz; ii++) {
    lvDaughters.insert(lv->GetDaughter(ii)->GetLogicalVolume());
  }

  std::set< G4LogicalVolume* >::const_iterator cite;
  for (cite = lvDaughters.begin(); cite != lvDaughters.end(); cite++) {
    DumpG4LVLeaf(*cite, leafDepth+1, out);
  } 
}


// Report Logical/Physical volume structure recursively

void DemoGeometryTools::DumpG4PVLVTree(std::ostream& out, G4int verbosity) 
{
  //dumps in the following order: 
  //    1) a LV with details
  //    2) list of PVs daughters of this LV with details
  //    3) list of LVs daughters of this LV and for each go to 1)

  //----- Get top PV
  G4LogicalVolume* topLV = GetTopLV();
  
  //----- Dump this leave (it will recursively dump all the tree)
  DumpPV(GetTopPV(), 0, out, verbosity);
  DumpG4PVLVLeaf(topLV, 0, out);
  
  //----- Dump the touchables (it will recursively dump all the tree)
  //t  DumpTouch(theTopPV, 0, out);
}

// Report Logical/Physical volume relationships

void DemoGeometryTools::DumpG4PVLVLeaf(G4LogicalVolume* lv, unsigned int leafDepth, std::ostream& out, G4int verbosity) {
  //----- Dump this LV 
  DumpLV(lv, leafDepth, out, verbosity);

  //----- Get LV daughters from list of PV daughters
  mmlvpv lvpvDaughters;
  std::set< G4LogicalVolume* > lvDaughters;
  int NoDaughters = lv->GetNoDaughters();
  while((NoDaughters--)>0) {
    G4VPhysicalVolume* pvD = lv->GetDaughter(NoDaughters);
    lvpvDaughters.insert(mmlvpv::value_type(pvD->GetLogicalVolume(), pvD));
    lvDaughters.insert(pvD->GetLogicalVolume());
  }
 
  std::set< G4LogicalVolume* >::const_iterator scite;
  mmlvpv::const_iterator mmcite;

  //----- Dump daughters PV and LV
  for (scite = lvDaughters.begin(); scite != lvDaughters.end(); scite++) {
    std::pair< mmlvpv::iterator, mmlvpv::iterator > mmER =
      lvpvDaughters.equal_range(*scite);
    
    //----- Dump daughters PV of this LV
    for (mmcite = mmER.first ; mmcite != mmER.second; mmcite++) {
      DumpPV((*mmcite).second, leafDepth+1, out, verbosity);
    }
    
    //----- Dump daughters LV
    DumpG4PVLVLeaf(*scite, leafDepth+1, out, verbosity);
  }
}
 
// Report G4LogicalVolume data

void DemoGeometryTools::DumpLV(G4LogicalVolume* lv, unsigned int leafDepth,
			       std::ostream& out, G4int verbosity) {
  int dumpVerbose = (verbosity)%10;
  G4String spaces = IndentLeaf(leafDepth);

  //----- dump name and details
  if (dumpVerbose >= 1) {
    out << leafDepth << spaces << "$$$ LOGICAL VOLUME = " << lv->GetName();
  }
  if (dumpVerbose >= 2) {
    out << "  Solid: " << lv->GetSolid()->GetName() << "  MATERIAL: "
	<< lv->GetMaterial()->GetName();
  }
  out << G4endl;

  //----- dump solid
  if (dumpVerbose >= 3) {
    DumpSolid(lv->GetSolid(), leafDepth, out);
  }

  //----- dump LV info 
  if (dumpVerbose >= 3) {
    //--- Visualisation attributes
    const G4VisAttributes* fVA = lv->GetVisAttributes();
    if (fVA) { 
      out << spaces << "  VISUALISATION ATTRIBUTES: " << "\n"
          << spaces << "    IsVisible " << fVA->IsVisible() << "\n"
          << spaces << "    IsDaughtersInvisible " << fVA->IsDaughtersInvisible() << "\n"
          << spaces << "    Colour " << fVA->GetColour() << "\n"
          << spaces << "    LineStyle " << fVA->GetLineStyle() << "\n"
          << spaces << "    LineWidth " << fVA->GetLineWidth() << "\n"
          << spaces << "    IsForceDrawingStyle " << fVA->IsForceDrawingStyle() << "\n"
          << spaces << "    ForcedDrawingStyle " << fVA->GetForcedDrawingStyle() << G4endl;
    }    
    
    //--- User Limits
    G4UserLimits *fUL = lv->GetUserLimits();
    G4Track dummy;
    if (fUL) {
      out << spaces << "    MaxAllowedStep " << fUL->GetMaxAllowedStep(dummy) << "\n"
          << spaces << "    UserMaxTrackLength " << fUL->GetUserMaxTrackLength(dummy) << "\n"
          << spaces << "    UserMaxTime " << fUL->GetUserMaxTime(dummy) << "\n"
          << spaces << "    UserMinEkine " << fUL->GetUserMinEkine(dummy) << "\n"
          << spaces << "    UserMinRange " << fUL->GetUserMinRange(dummy) << G4endl;
    }
    
    //--- other LV info
    if (lv->GetSensitiveDetector()) {
      out << spaces << "  IS SENSITIVE DETECTOR " << G4endl;
    }
    if (lv->GetFieldManager()) {
      out << spaces << "  FIELD ON " << G4endl;
    }

    out << spaces << "    Quality for optimisation, average number of voxels"
	<< " to be spent per content " << lv->GetSmartless() << G4endl;

    if (lv->GetFastSimulationManager()) {
      out << spaces << "   Logical Volume is an envelope for a"
	  << " FastSimulationManager " << G4endl;
    }
    out << spaces << "     Weight used in the event biasing technique = "
	<< lv->GetBiasWeight() << G4endl;
  }
}


// Report G4VPhysicalVolume data

void DemoGeometryTools::DumpPV(G4VPhysicalVolume* pv, unsigned int leafDepth, std::ostream& out, G4int verbosity) {
  if (pv == 0) return;
  int dumpVerbose = (verbosity/10)%10;
  G4String spaces = IndentLeaf(leafDepth);

  //----- PV info
  if (dumpVerbose >= 1) {
    G4String mother = "World";
    if (pv->GetMotherLogical()) mother = pv->GetMotherLogical()->GetName();
    out << leafDepth << spaces << "### PHYSICAL VOLUME = " << pv->GetName() 
	<< " Copy No " << pv->GetCopyNo() << " in " << mother;
    if (dumpVerbose == 1) out << G4endl;
  }
  
  if (dumpVerbose >= 2) {
    out << " at " << pv->GetTranslation();
  }

  if (!pv->IsReplicated()) {
    if (pv->GetRotation() == 0) {
      if (dumpVerbose >= 2) {
	out << " with no rotation" << G4endl;
      }
    } else {
      if (dumpVerbose >= 2) {
	//just rotation name 
	out  << " with rotation" << G4endl;
      } else if (dumpVerbose >= 3) {
	out  << " with rotation " << *(pv->GetRotation()) << G4endl;
      }
    }
  } else {
    if (dumpVerbose >= 3) {
      out << spaces << "    It is replica: " << G4endl;
      EAxis axis; G4int nReplicas; G4double width; G4double offset; G4bool consuming;
      pv->GetReplicationData(axis, nReplicas, width, offset, consuming);
      out << spaces << "     axis " << axis << G4endl
	  << spaces << "     nReplicas " << nReplicas << G4endl;
      if (pv->GetParameterisation() != 0) {
	out << spaces << "    It is parameterisation " << G4endl;
      } else { 
	out << spaces << "     width " << width << G4endl
	    << spaces << "     offset " << offset << G4endl
	    << spaces << "     consuming" <<  consuming << G4endl;
      }
      if (pv->GetParameterisation() != 0) {
	out << spaces << "    It is parameterisation " << G4endl;
      }
      
    }
  }
}


//-----------------------------------------------------------------------
void DemoGeometryTools::DumpSolid(G4VSolid* sol, unsigned int leafDepth, std::ostream& out) {
  out << IndentLeaf(leafDepth) << *(sol) << G4endl;
}


//-----------------------------------------------------------------------
G4VPhysicalVolume* DemoGeometryTools::GetTopPV() {
  return G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume();
}


//-----------------------------------------------------------------------
G4LogicalVolume* DemoGeometryTools::GetTopLV(G4bool bGeomInit) {
  G4LogicalVolume* topLV ;
  theTopPV = GetTopPV();

  if (verboseLevel) G4cout << "DemoGeometryTools::GetTopLV:"
			   << " topPV " << theTopPV << G4endl;

  if (bGeomInit && theTopPV) {
    topLV = theTopPV->GetLogicalVolume();
  } else {
    G4PhysicalVolumeStore* pvs = G4PhysicalVolumeStore::GetInstance();
    vpv::iterator cite;
    for (cite = pvs->begin(); cite != pvs->end(); cite++) {
      if ((*cite)->GetMotherLogical() == 0) break;
    }
    topLV = (*cite)->GetLogicalVolume();
  }

  return topLV;
}


// Search for (unique) material with specified name

G4Material* DemoGeometryTools::GetMaterial(const G4String& mname) {
  G4Material* mate = 0;
  const G4MaterialTable* matTab = G4Material::GetMaterialTable();
  std::vector<G4Material*>::const_iterator matite;
  for (matite = matTab->begin(); matite != matTab->end(); matite++) {
    if (WildcardMatch(mname, (*matite)->GetName())) {
      //    if ((*matite)->GetName() == mname) {
      if (mate != 0) {
	G4Exception("DemoGeometryTools::GetMaterial","WARNING",JustWarning,("More than one material with same name: "+mname + " The first will be taken").c_str());
      }
      mate = *matite;
    }
  }

  return mate;
}

std::vector<G4Material*> DemoGeometryTools::GetMaterials(const G4String& mname) {
  std::vector<G4Material*> mates;
  const G4MaterialTable* matTab = G4Material::GetMaterialTable();
  std::vector<G4Material*>::const_iterator matite;
  for (matite = matTab->begin(); matite != matTab->end(); matite++) {
    if (WildcardMatch(mname, (*matite)->GetName())) {
      //    if ((*matite)->GetName() == mname) {
      mates.push_back(*matite);
    }
  }

  return mates;
}


//-----------------------------------------------------------------------
int DemoGeometryTools::NumberOfTouchables() {
  return 0;		// FIXME:  Currently disabled

  int nTouch = 0;
  G4LogicalVolume* lv = GetTopLV();
  
  Add1Touchable(lv, nTouch, 1);
  return nTouch;
}


//---------------------------------------------------------------------------
void DemoGeometryTools::Add1Touchable(G4LogicalVolume* lv, int& nTouch,
				      int nReplicas) {
  int siz = lv->GetNoDaughters();
  for (int ii = 0; ii < siz; ii++) {
    G4int newnReplicas = 1;
    if (lv->GetDaughter(ii)->IsReplicated()) {
      EAxis axis;
      G4double width;
      G4double offset;
      G4bool consuming;
      lv->GetDaughter(ii)->GetReplicationData(axis, newnReplicas, width, offset, consuming);
    }

    nReplicas *= newnReplicas;
    nTouch += nReplicas;

    if (verboseLevel>1) {
      G4cout << " Add1Touchable " << nTouch << " nReplicas " << nReplicas
	     << " newnReplicas " << newnReplicas << " "
	     << lv->GetDaughter(ii)->GetName() << G4endl;
    }

    Add1Touchable(lv->GetDaughter(ii)->GetLogicalVolume(), nTouch, nReplicas);
  } 
}


// Build hierarchical structure of geometry for all physical volumes

const vtch& DemoGeometryTools::GetTouchables(const G4String& tname) {
  if (verboseLevel) 
    G4cout << " DemoGeometryTools::GetTouchables " << tname << G4endl;

  // If name already requested, return result from map
  if (theTouchables.find(tname) != theTouchables.end())
    return theTouchables[tname];

  if (verboseLevel>1)
    G4cout << " Constructing touchable(s) for named volume(s) ..." << G4endl;

  vpsi ancestorsRequested = ExtractAncestorsRequested(tname);
  CollectAncestorLines(tname);

  // Create entry in buffer of previously requested touchables 
  // Construct DemoTouchable's for each line of PV's ancestors 
  vtch& vtouch = theTouchables[tname];
  vtouch.clear();

  G4int asiz = theAncestorLines.size();
  if (verboseLevel>1) 
    G4cout << " Got total of " << asiz << " ancestors" << G4endl;

  // Loop to each line of PV's ancestors
  for (G4int ii = 0; ii < asiz; ii++) {
    vpv pvlist = theAncestorLines[ii];
    G4int siz =  pvlist.size();
    std::vector<G4int> vi;

    // build list of copyNo's of each ancestor:

    // If it is replica of parameterisation 1 PV has many
    // touchables. Therefore if copyNo is -1, this means that you want
    // all the touchables of the replica, so you have to expand it
    for (G4int jj = 0; jj < siz; jj++) {
      G4int copyn;
      G4int copyNoAncestor = -1;

      // Simple placement: the copyNo you are requesting is the
      // PV's copyNo (already checked in AddParentPV method)
      if (!pvlist[jj]->IsReplicated()) copyNoAncestor = pvlist[jj]->GetCopyNo(); 
      else {
	// replica/parameterisation: PV::copyNo has no sense,
	// requested copyNo is actually meaning the number of the
	// G4Touchable, not the G4PhysicalVolume
	if (jj < G4int(ancestorsRequested.size()))
	  copyNoAncestor = ancestorsRequested[jj].second;
      }

      // here -1 means that it is a replica and you have not requested
      // any copy number of it. If it were a simple placement and you
      // have not requested any copy number, then you would have build
      // an ancestor line for each PV
      if (copyNoAncestor == -1) {
	EAxis axis;
	G4int nReplicas;
	G4double width;
	G4double offset;
	G4bool consuming;
	pvlist[jj]->GetReplicationData(axis, nReplicas, width, offset, consuming);

	copyn = -10-nReplicas; 	// Special for use in method ExpandCopyNolist
      } else {
	copyn = copyNoAncestor;
      }

      vi.push_back(copyn);
    }

    // Now expand list of copyNo's, if needed
    std::vector< std::vector<G4int> > vvi;
    ExpandCopyNoList(vvi, vi, 0);

    // Loop to the expanded list of copyNo's of this ancestor line
    G4int sizv = vvi.size();
    for (G4int kk = 0; kk < sizv; kk++) {
      if (verboseLevel>2) {      		// dump ancestor line
       	G4cout << ii << "theAncestorLine: ";
	for (G4int jj = 0; jj < G4int(pvlist.size()); jj++) {
	  G4cout << " " << pvlist[jj]->GetName()  << " " << pvlist[jj]->GetCopyNo()
		 << "=" << vvi[kk][jj] << " " << pvlist[jj] ;
	}
	G4cout << G4endl;
      }

      // Build DemoTouchable: from vpv you will get the PV info,
      // the list of copyNo's used will be vvi[kk], not the copyNo's of
      // the PV's, because the copy No has no meaning for replicas

      //build the long name to find if it exists already:
      G4String longName = "";
      G4String ancestorName;
      G4int nAnces = pvlist.size();
      for (G4int ll = 0; ll < nAnces-1; ll++) {
	G4VPhysicalVolume* pv = pvlist[ll];
	G4int copyNoAncestor =
	  pv->IsReplicated() ? vvi[kk][ll] : pv->GetCopyNo();

	ancestorName = delimiter + AppendCopyNo(pv->GetName(), copyNoAncestor);
	longName.prepend(ancestorName);
      }

      DemoTouchable* to = new CDMSTouchable(pvlist, vvi[kk], verboseLevel);

      if (verboseLevel>2) 
	G4cout << longName << " touchable name " << to->GetLongName() << G4endl;

      // Add it to the list of touchables to return
      vtouch.push_back(to);
    }	// Loop over PVs in ancestor list
  }	// Loop over ancestorLines

  if (verboseLevel>2) G4cout << " vtouch size " << vtouch.size() << G4endl;

  return vtouch;
}

//-----------------------------------------------------------------------
const vtch& DemoGeometryTools::GetTouchablesOfLV(const G4String& lvname) {
  if (verboseLevel)
    G4cout << " DemoGeometryTools::GetTouchablesOfLV " << lvname << G4endl;

  G4String lvident = lvname+":LV";	// For unique lookup of touchables

  // If name already requested, return result from map
  if (theTouchables.find(lvident) != theTouchables.end())
    return theTouchables[lvident];

  if (verboseLevel>1)
    G4cout << " Constructing touchables for logical volume ..." << G4endl;

  vtch& lvTouchables = theTouchables[lvident];

  vpv pvlist = GetPhysicalVolumesOfLV(lvname);	// All placements of LV
  if (pvlist.empty()) return lvTouchables;

  if (verboseLevel>1)
    G4cout << " ... found " << pvlist.size() << " placements ..." << G4endl;

  G4String pvname;
  for (size_t ii=0; ii<pvlist.size(); ii++) {
    G4VPhysicalVolume* ipv = pvlist[ii];
    pvname = ipv->GetName();
    if (!ipv->IsReplicated()) pvname = AppendCopyNo(pvname, ipv->GetCopyNo());

    const vtch& pvTouch = GetTouchables(pvname);
    lvTouchables.insert(lvTouchables.end(), pvTouch.begin(), pvTouch.end());
  }

  return lvTouchables;
}

//-----------------------------------------------------------------------
DemoTouchable* CDMSGeometryTools::GetTouchable(const G4ThreeVector& pos) {
  if (!theNavigator) RebuildDictionaries();	// This should never happen
  if (!theNavigator->GetWorldVolume()) theNavigator->SetWorldVolume(GetTopPV());

  if (!theTouchable) theTouchable = new G4TouchableHistory;
  theNavigator->LocateGlobalPointAndUpdateTouchable(pos, theTouchable, false);

  return new DemoTouchable(theTouchable);
}

//-----------------------------------------------------------------------
std::vector<G4VPhysicalVolume*> 
DemoGeometryTools::GetPhysicalVolumeCopy(const G4String& pvname, G4int copyNo) {
  std::stringstream svname;
  svname << pvname;
  if (copyNo>=0) svname << ":" << copyNo;
  return GetPhysicalVolumes(svname.str());
}

//-----------------------------------------------------------------------
std::vector<G4VPhysicalVolume*> 
DemoGeometryTools::GetPhysicalVolumes(const G4String& pvname, G4int nVols) {
  vpv vvolu;
  std::string::size_type ial = pvname.rfind(":");

  if (verboseLevel) 
    G4cout << " DemoGeometryTools::GetPhysicalVolumes volname " << pvname
	   << " " << nVols << G4endl;

  G4String volname = "";
  G4int volcopy = 0;
  if (ial != std::string::npos) {
    std::string::size_type ial2 = pvname.rfind(":",ial-1);
    if (ial2 != std::string::npos) {
      G4cerr << "DemoGeometryTools::GetPhysicalVolumes: name corresponds to"
	     << " a touchable " << pvname;
    } else { 
      volname = pvname.substr(0, ial);
      volcopy = atoi(pvname.substr(ial+1, pvname.length()).c_str());
    }
  } else {
    volname = pvname;
    volcopy = -1;
  }

  G4PhysicalVolumeStore* pvs = G4PhysicalVolumeStore::GetInstance();
  vpv::iterator citepv;
  for (citepv = pvs->begin(); citepv != pvs->end(); citepv++) {
    if (verboseLevel>1) {
      G4cout << " DemoGeometryTools::GetPhysicalVolumes  volname " << volname
	     << " copy " << volcopy << " TRY " << (*citepv)->GetName()
	     << " " << (*citepv)->GetCopyNo() << " "
	     << WildcardMatch(volname, (*citepv)->GetName()) << G4endl;
    }

    if (WildcardMatch(volname, (*citepv)->GetName()) 
	&& ((*citepv)->GetCopyNo() == volcopy || -1 == volcopy)) {

      if (verboseLevel>2)
	G4cout << " DemoGeometryTools::GetPhysicalVolumes found " << G4endl;

      vvolu.push_back(*citepv);
    }
  }

  //----- Check that at least one volume was found
  if (nVols != -1 && G4int(vvolu.size()) != nVols) {
    std::stringstream message;
    message << "Wrong number of physical volumes: found " << vvolu.size()
	    << ", requested " << nVols;

    DemoAbortJob("CDMSGeometryTools::GetPhysicalVolumes", message.str());
  }

  return vvolu;
}


//-----------------------------------------------------------------------
std::vector<G4String> 
DemoGeometryTools::GetPhysicalVolumeNames(const G4String& pvname) {
  std::vector<G4String> vvolu;
  std::string::size_type ial = pvname.rfind(":");

  if (verboseLevel)
    G4cout << " DemoGeometryTools::GetPhysicalVolumes volname " << ial
	   << " " << pvname << G4endl;

  G4String volname = "";
  G4int volcopy = 0;
  if (ial != std::string::npos) {
    std::string::size_type ial2 = pvname.rfind(":",ial-1);
    if (ial2 != std::string::npos) {
      G4cerr << "DemoGeometryTools::GetPhysicalVolumes: name corresponds to"
	     << " a touchable " << pvname;
    }else { 
      volname = pvname.substr(0, ial);
      volcopy = atoi(pvname.substr(ial+1, pvname.length()).c_str());
    }
  } else {
    volname = pvname;
    volcopy = -1;
  }

  G4PhysicalVolumeStore* pvs = G4PhysicalVolumeStore::GetInstance();
  vpv::iterator citepv;
  for (citepv = pvs->begin(); citepv != pvs->end(); citepv++) {
    if (verboseLevel>1) {
      G4cout << " DemoGeometryTools::GetPhysicalVolumeNames volname " << volname
	     << " copy " << volcopy << " TRY " << (*citepv)->GetName()
	     << " " << (*citepv)->GetCopyNo() << " "
	     << WildcardMatch(volname, (*citepv)->GetName()) << G4endl;
    }

    if (WildcardMatch(volname, (*citepv)->GetName()) 
	&& ((*citepv)->GetCopyNo() == volcopy || -1 == volcopy)) {

      if (verboseLevel>2) 
	G4cout << " DemoGeometryTools::GetPhysicalVolumeNames found " << G4endl;

      G4String pvName = (*citepv)->GetName();
      if ((*citepv)->IsReplicated()) {
	EAxis axis;
	G4int nReplicas;
	G4double width;
	G4double offset;
	G4bool consuming;
	(*citepv)->GetReplicationData(axis, nReplicas, width, offset, consuming);
	for (G4int ii = 0; ii < nReplicas; ii++) {
	  vvolu.push_back(AppendCopyNo(pvName,ii));
	}
      } else {
	vvolu.push_back(AppendCopyNo(pvName,(*citepv)->GetCopyNo()));
      }
    }
  }

  return vvolu;
}


//-----------------------------------------------------------------------
std::vector<G4LogicalVolume*> 
DemoGeometryTools::GetLogicalVolumes(const G4String& lvname, G4int nVols) {
  if (verboseLevel) 
    G4cout << " DemoGeometryTools::GetLogicalVolumes " << lvname << G4endl;

  vlv vvolu;
  G4int ial = lvname.rfind(":");
  if (ial != -1) {
    G4cerr << "DemoGeometryTools::GetLogicalVolumes: '" << lvname
	   << "' corresponds to a touchable or physcal volume" << G4endl;
    return vvolu;
  }

  G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
  vlv::iterator citelv;
  for (citelv = lvs->begin(); citelv != lvs->end(); citelv++) {
    if (WildcardMatch(lvname,(*citelv)->GetName())) {
      vvolu.push_back(*citelv);
    }
  }

  if (nVols != -1 && G4int(vvolu.size()) != nVols) {
    std::stringstream message;
    message << "Wrong number of logical volumes: found " << vvolu.size()
	    << ", requested " << nVols;

    DemoAbortJob("CDMSGeometryTools::GetLogicalVolumes", message.str());
  } 

  return vvolu;
}


//-----------------------------------------------------------------------
std::vector<G4VPhysicalVolume*> 
DemoGeometryTools::GetPhysicalVolumesOfLV(const G4String& lvname) {
  if (verboseLevel) 
    G4cout << " DemoGeometryTools::GetPhysicalVolumesOfLV " << lvname << G4endl;

  vlv lvlist = GetLogicalVolumes(lvname);
  return GetPhysicalVolumes(lvlist);
}

//-----------------------------------------------------------------------
// NOTE: Adapted from http://stackoverflow.com/questions/771453/
namespace {
  template <typename pair_t> struct second_t {
    typename pair_t::second_type operator()(const pair_t& p) { return p.second; }
  };

  template <typename map_t> second_t<typename map_t::value_type>
  get_second(const map_t&) { return second_t<typename map_t::value_type>(); }
}

//-----------------------------------------------------------------------
std::vector<G4VPhysicalVolume*> 
DemoGeometryTools::GetPhysicalVolumes(G4LogicalVolume* lv) {
  if (verboseLevel) {
    G4cout << " DemoGeometryTools::GetPhysicalVolumes "
	   << (lv?lv->GetName():"NULL") << G4endl;
  }

  vpv pvlist;			// Output buffer
  if (!lv) return pvlist;

  std::pair<mmlvpv::iterator, mmlvpv::iterator> lprange = theLVPVTree.equal_range(lv);

  std::transform(lprange.first, lprange.second, std::back_inserter(pvlist),
		 get_second(theLVPVTree));

  return pvlist;
}

//-----------------------------------------------------------------------
std::vector<G4VPhysicalVolume*> 
DemoGeometryTools::GetPhysicalVolumes(const vlv& lvlist) {
  if (verboseLevel) 
    G4cout << " DemoGeometryTools::GetPhysicalVolumes (LVlist)" << G4endl;

  vpv pvlist;		// Output buffer
  if (lvlist.empty()) return pvlist;

  vpv pvoflv;
  for (size_t ilv=0; ilv<lvlist.size(); ilv++) {
    pvoflv = GetPhysicalVolumes(lvlist[ilv]);
    pvlist.insert(pvlist.end(), pvoflv.begin(), pvoflv.end());
  }

  return pvlist;
}

//-----------------------------------------------------------------------
std::vector<G4Region*> DemoGeometryTools::GetRegions(const G4String& rname) {
  if (verboseLevel) 
    G4cout << " DemoGeometryTools::GetRegions " << rname << G4endl;

  std::vector<G4Region*> regions;
  G4RegionStore* regionStore = G4RegionStore::GetInstance();  
  for (unsigned int jj = 0; jj < regionStore->size(); jj++) {
    G4Region* reg = (*regionStore)[jj];
    if (WildcardMatch(rname,reg->GetName())) {
      regions.push_back(reg);
    }
  }
  
  return regions;
}

//-----------------------------------------------------------------------
void DemoGeometryTools::CollectAncestorLines(const G4String& tname) {
  if (verboseLevel)
    G4cout << " DemoGeometryTools::CollectAncestorLines " << tname << G4endl;

  // 'ancestorsRequested' will host the list of ancestors in pairs
  //   ancestor name - ancestor copy number 
  theAncestorLines.clear();
  vpsi ancestorsRequested = ExtractAncestorsRequested(tname);

  // True name of volume requested is the top name in the list of ancestors
  G4String volName = ancestorsRequested[0].first;

  // Get how many physical volumes are there of each ancestor:
  //   if there are 4 volumes A, each volume A is placed 2 times in
  //   volume B and volume B 3 times in volume C, there has to be
  //   constructed 4 X 2 X 3 DemoTouchable's

  // loop to all the PV's with the same name than volName
  G4PhysicalVolumeStore* pvs = G4PhysicalVolumeStore::GetInstance();
  vpv::iterator cite;
  for (cite = pvs->begin(); cite != pvs->end(); cite++) {
    G4VPhysicalVolume* pv0 = *cite;

    if (WildcardMatch(volName, pv0->GetName())) {
      if (verboseLevel>1) 
	G4cout << " Found matching volume pv0 " << pv0->GetName()
	       << " :" << pv0->GetCopyNo() << G4endl;

      // For each PV (assume PV is inside LV)
      // 1. get the LV corresponding to this PV
      // 2. get the LV parent
      // 3. get all PV corresponding to this LV

      // create vector of PV ancestors of this PV: each PV will
      // have associated a vector of its PV's ancestors. All the
      // vectors will be then stored in 'theAncestorLines'
      vpv pvlist;
      
      // check that it is the copy number we are looking for
      bool copyOK = CheckPVCopyNo(pv0, ancestorsRequested[0].second);
      if (copyOK) { 
	//-- add it to the vector of ancestors as the first one
	pvlist.push_back(pv0);
	//-- add parent (recursively all ancestors will be added)
	AddParentPV(pv0, pvlist, ancestorsRequested, 0);
      }
    }
  }
}

//-----------------------------------------------------------------------
vpsi DemoGeometryTools::ExtractAncestorsRequested(const G4String& aname) const {
  if (verboseLevel) {
    G4cout << " DemoGeometryTools::ExtractAncestorsRequested " << aname
	   << G4endl;
  }

  //---------- decompose name in the volume's and copyNo's of the ancestors:
  // |volA:2|volB:4|volC:3  =>  <volA,2>,  <volB,4>, <volC,3>
  G4String namet = aname;
  vpsi ancestorsRequested;
  //---------- loop until a delimiter is found in position 0, or not found
  for (;;) {  
    G4String shortName;
    G4int fnumb = namet.rfind(':');
    G4int fslash = namet.rfind(delimiter);

    if (verboseLevel>2)
      G4cout << namet << " test fnumb " << fnumb << " fslash " << fslash << G4endl;

    G4int numb;
    if (fnumb < fslash || fnumb == -1) {
      //----- name without number, then set number is -1
      shortName = namet.substr(fslash+1, namet.size()-fslash);
      numb = -1;
    } else {
      shortName = namet.substr(fslash+1, fnumb-fslash-1);
      numb = atoi((namet.substr(fnumb+1, namet.size()-1)).c_str());
      //---- check that copy number requested is not negative (or 0)
      if (numb < -1) {
	DemoAbortJob("CDMSGeometryTools::ExtractAncestorsRequested",
		     "Cannot ask for a copy number smaller than 0 "+aname);
      }
    }
    std::pair<G4String, G4int> psi(shortName, numb);

    if (verboseLevel>1)
      G4cout << " DemoGeometryTools::ExtractAncestorsRequested new node "
	     << psi.first << " : " << psi.second << G4endl;

    ancestorsRequested.push_back(psi);
    if (fslash == 0 || fslash == -1) break;
    namet = namet.substr(0, fslash);
  }

  return ancestorsRequested;
}


//-----------------------------------------------------------------------
void DemoGeometryTools::AddParentPV(const G4VPhysicalVolume* pv,
				    std::vector<G4VPhysicalVolume*>& pvlist,
				    vpsi& ancestorsRequested,
				    G4int ancestorLevel) {
  if (verboseLevel) {
    G4cout << " DemoGeometryTools::AddParentPV of " << pv->GetName()
	   << ":" << pv->GetCopyNo() << " ancestorsRequested.size() "
	   << ancestorsRequested.size() << " ancestorLevel "
	   << ancestorLevel << G4endl;
  }
  
  //----- 1. get the LV parent:
  // As a PV does not have any pointer to its parent (only if the
  // placement is done inside PV, instead of inside LV, this last
  // being the only option currently supported), we have to look for
  // its parent LV and then get the list of all the PV's of this LV.
  const G4LogicalVolume* lvParent = thePVLVTree[ const_cast<G4VPhysicalVolume*>(pv) ];

  if (verboseLevel>1)
    G4cout << " lvParent " << lvParent << " PVLV size " << thePVLVTree.size()
	   << G4endl;
  
  if (lvParent == 0) {	   //if top volume, store this chain of ancestors
    theAncestorLines.push_back(pvlist);
    return;
  }

  //----- get which is the parent name and copy number you are looking for 
  G4String parentName;
  G4int parentCopyNo;
  if (ancestorLevel+1 < int(ancestorsRequested.size())) {
    parentName = ancestorsRequested[ancestorLevel+1].first;
    parentCopyNo = ancestorsRequested[ancestorLevel+1].second;
  } else {
    //no parent defined, all will serve
    parentName = " ";
    parentCopyNo = -1;
  }
  
  //----- check that parent is the one asked for
  if (!WildcardMatch(parentName, lvParent->GetName()) && parentName != " ") return;
  //    if (lvParent->GetName() != parentName && parentName != " ") return;
  
  //----- 2.  get all PV corresponding to this LV 
  mmlvpv::const_iterator mmcite;
  std::pair<mmlvpv::iterator, mmlvpv::iterator> mmdi;
  mmdi = theLVPVTree.equal_range(const_cast<G4LogicalVolume*>(lvParent));
  //    G4cout << " theLVPVTree.equal_range " << ((*(mmdi.first)).first)->GetName() 
  //	   << " " << ((*(mmdi.second)).first)->GetName() << G4endl;
  if (mmdi.first == mmdi.second) { 
    //-second == theLVPVTree.end()) {
    //this should never happen if there is no bug in DetRep
    std::stringstream message;
    message << "For PV " << pv->GetName() << " LV parent not found ";
    if (lvParent) message << lvParent->GetName();

    DemoAbortJob("CDMSGeometryTools::AddParentPV", message.str());
  }
  
  if (verboseLevel>2) {
    G4int iim = 0;
    for (mmcite = mmdi.first; mmcite != mmdi.second; mmcite++) {
      iim++;
      G4cout << " LV parent name: " << lvParent->GetName()
	     << " : PV corresponding to LV parent = " << iim << G4endl;
    }
  }
  
  //t    G4int iim = 0;
  //---- Loop to all the PV's corresponding to LV parent
  for (mmcite = mmdi.first; mmcite != mmdi.second; mmcite++) {
    G4VPhysicalVolume* pv0 = (*mmcite).second;
    
    if (verboseLevel>2) {
      G4cout << " another PV with = LV as parent " << (*mmcite).first
	     << " " << pv0->GetName() << " " << pv0->GetCopyNo() << G4endl;
    }
    
    //--- check that it is the copy number that you want
    bool copyOK = CheckPVCopyNo(pv0, parentCopyNo);
    if (copyOK) {
      //--- add it to the list of ancestors
      pvlist.push_back(pv0);
      //--- look one level up
      AddParentPV(pv0, pvlist, ancestorsRequested, ancestorLevel+1);
      
      //--- when looking for next ancestors, vpv will be added a new
      //    element, but when looking to the next PV in this loop, the
      //    PV added should be deleted
      pvlist.pop_back();
    }
  }
}
 


//-----------------------------------------------------------------------
bool DemoGeometryTools::
CheckPVCopyNo(const G4VPhysicalVolume* pv, G4int copyNo) const {
  bool copyOK = false;

  if (copyNo == -1) copyOK = true;
  else {  //--- Simple placement: check that PV's copy number = copyNo
    if (!pv->IsReplicated()) { 
      copyOK = (copyNo == pv->GetCopyNo());

      if (verboseLevel>1) {
	G4cout << " DemoGeometryTools::CheckPVCopyNo " << copyNo << " =? "
	       << pv->GetCopyNo() << " OK= " << copyOK << G4endl; 
      }

      //--- Replica or parameterisation: PV::theCopyNo has no sense, a PV has many touchables, so we will take the PV always (later, when building the Cg4DrTouchable, we will get the touchable number copyNo), except if the copyNo we are asking for is bigger than the number of copies replicated
    } else {
      //    G4PVReplica* pvrep = static_cast<G4PVReplica*>(pv);
      EAxis axis;
      G4int nReplicas;
      G4double width;
      G4double offset;
      G4bool consuming;
      pv->GetReplicationData(axis, nReplicas, width, offset, consuming);
      copyOK = (copyNo < nReplicas);

      if (verboseLevel>1) {
	G4cout << " DemoGeometryTools::CheckPVCopyNo(replica) " << copyNo
	       << " =? " << nReplicas << " OK= " << copyOK << G4endl; 
      }
    }
  }

  return copyOK;
}



//-----------------------------------------------------------------------
void DemoGeometryTools::ExpandCopyNoList(std::vector<std::vector<G4int> >& vvi,
					 std::vector<G4int> vi, G4int level) {

  if (verboseLevel) 
    G4cout << " DemoGeometryTools::ExpandCopyNoList " << level << G4endl;

  G4int siz = vi.size();
  G4int ii, jj;
  std::vector< G4int > viNew = vi; 
  for (ii = level; ii < siz; ii++) {
    G4int iiNew;
    //---- if it is < -10, it means it has to be expanded: from 1 -13 0,
    //     you create three lists: 1 0 0, 1 1 0, 1 2 0 
    if (vi[ii] < -10) {
      iiNew = -vi[ii] - 10 - 1;
      for (jj = 0; jj < iiNew; jj++) {
	viNew[ii] = jj;
	ExpandCopyNoList(vvi, viNew, level+1);
      }
    } else {
      iiNew = vi[ii];
    }
    viNew[ii] = iiNew;
  }
  vvi.push_back(viNew);

  if (verboseLevel>1) {
    G4cout << vvi.size() << " expandList ";
    for (ii=0; ii < siz; ii ++) {
      G4cout << " " << viNew[ii];
    }
    G4cout << G4endl;
  }

}


//-----------------------------------------------------------------------
void DemoGeometryTools::RebuildDictionaries() {
  thePVs.clear();
  theLVPVTree.clear();
  thePVLVTree.clear();
  theTouchables.clear();
  theTopPV = 0;
  delete theNavigator; theNavigator=0;

  BuildDictionaries();
}
  

//-----------------------------------------------------------------------
void DemoGeometryTools::BuildDictionaries() {
  G4PhysicalVolumeStore* pvs = G4PhysicalVolumeStore::GetInstance();
  vpv::iterator citepv;
  //  G4cout << " pvs size " << pvs->size() << G4endl;
  for (citepv = pvs->begin(); citepv != pvs->end(); citepv++) {
    //build multimap of PV name - list of PV*'s 
    thePVs.insert(mmspv::value_type((*citepv)->GetName(), *citepv));
    //build multimap of LV* - list of its PV*'s
    theLVPVTree.insert(mmlvpv::value_type((*citepv)->GetLogicalVolume(), *citepv));
    //-    cout << " inserting theLVPVTree " << (*citepv)->GetLogicalVolume()->GetName() << " " << (*citepv)->GetName() << endl;
  }

  //build map of PV* - parent LV* : the only info in GEANT4 is from a LV, the list of its daughters PV's
  G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
  vlv::iterator citelv;
  //  G4cout << " lvs size " << lvs->size() << G4endl;

  G4int ii, siz;
  for (citelv = lvs->begin(); citelv != lvs->end(); citelv++) {
    siz = (*citelv)->GetNoDaughters();
    for (ii=0; ii< siz; ii++) {
      thePVLVTree[(*citelv)->GetDaughter(ii)] = *citelv;
    }
    //    cout << " G4LOGICALVOLUME " << (*citelv)->GetName() << " " << (*citelv)->GetMaterial()->GetName() << endl;
  }

  // Initialize local navigator with new geometry structure
  if (!theNavigator) theNavigator = new G4Navigator;
}

//-----------------------------------------------------------------------
std::set<G4String> DemoGeometryTools::GetAllSDTypes() {
 std::set<G4String> sdtypes;
 
 const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
 vlv::const_iterator lvcite;
 for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) {
   if ((*lvcite)->GetSensitiveDetector() != 0) {

     if (verboseLevel>1) {
       G4cout << " SD type " << (*lvcite)->GetSensitiveDetector()->GetPathName()
	      << G4endl;
     }

     G4String sdtype = (*lvcite)->GetSensitiveDetector()->GetPathName();
     sdtypes.insert(sdtype.substr(1,sdtype.length()-2));
   }
 }

 return sdtypes;

}

//-----------------------------------------------------------------------
G4double DemoGeometryTools::GetDistanceToOutNoDirection(const G4Track* aTrack) {
  // G4double safety = aTrack->GetStep()->GetPreStepPoint()->GetSafety();
  G4VPhysicalVolume* pv = aTrack->GetVolume();
  if (pv == 0) {
    G4cerr << " DemoGeometryTools::GetDistanceToOutNoDirection:  no PV associated to Track " << aTrack->GetPosition() << " particle" << aTrack->GetDefinition()->GetParticleName() << " energy " << aTrack->GetKineticEnergy() << G4endl;
    return  0.;
  }
  G4VSolid* solid = pv->GetLogicalVolume()->GetSolid();
  const G4VTouchable* touch = aTrack->GetTouchable();
  G4ThreeVector localPos = touch->GetHistory()->GetTopTransform().TransformPoint(aTrack->GetPosition());
  G4double d2out = solid->DistanceToOut(localPos);
  
  return d2out;
}


//-----------------------------------------------------------------------
G4double DemoGeometryTools::GetDistanceToOutDirection(const G4Track* aTrack) {  
  // G4double safety = aTrack->GetStep()->GetPreStepPoint()->GetSafety();
  G4VPhysicalVolume* pv = aTrack->GetVolume();
  if (pv == 0) {
    G4cerr << " DemoGeometryTools::GetDistanceToOutDirection:  no PV associated to Track " << aTrack->GetPosition() << " particle" << aTrack->GetDefinition()->GetParticleName() << " energy " << aTrack->GetKineticEnergy() << G4endl;
    return  0.;
  }
  G4VSolid* solid = pv->GetLogicalVolume()->GetSolid();
  const G4VTouchable* touch = aTrack->GetTouchable();
  G4ThreeVector localPos = touch->GetHistory()->GetTopTransform().TransformPoint(aTrack->GetPosition());
  G4ThreeVector localDir = touch->GetHistory()->GetTopTransform().TransformAxis(aTrack->GetMomentumDirection());

  G4double d2out = solid->DistanceToOut(localPos, localDir);
  
  return d2out;
}


//-----------------------------------------------------------------------
G4String DemoGeometryTools::GetG4TouchableLongName(G4TouchableHistory* touch) {
  G4String tname = "";

  for (G4int ii = touch->GetHistoryDepth(); ii >= 0; ii--) {
    G4VPhysicalVolume* pv = touch->GetVolume(ii);
    tname += delimiter + AppendCopyNo(pv->GetName(), pv->GetCopyNo());
  }

  return tname;
}


//-----------------------------------------------------------------------
std::vector<G4Material*> DemoGeometryTools::GetMaterialsUsed() {
  std::vector<G4Material*> mateVector;
  std::set<G4Material*> mateSet;

  //---- Assume all G4LogicalVolume's are part of the geometry tree
  const G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
  vlv::const_iterator lvcite;
  for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) {
    mateSet.insert((*lvcite)->GetMaterial());
  }

  std::set<G4Material*>::const_iterator ites;
  for (ites = mateSet.begin(); ites != mateSet.end(); ites++) {
    mateVector.push_back(*ites); 
  }

  return mateVector;
}


//--------------------------------------------------------------------
G4double DemoGeometryTools::GetCubicVolume(G4VPhysicalVolume* pv) {
  return pv->GetLogicalVolume()->GetSolid()->GetCubicVolume();
}


//------------------------------------------------------------------------
G4String DemoGeometryTools::BuildTouchableName(const G4ThreeVector& pos) {
  GetTouchable(pos);		// Fills local buffer 'theTouchable'
  return GetG4TouchableLongName(theTouchable);
}

// Extract short (PV) name from ancestor list at given level (0 = end leaf)
G4String DemoGeometryTools::ExtractShortName(const G4String& longname,
					     G4int level) const {
  if (verboseLevel) {
    G4cout << " DemoGeometryTools::ExtractShortName " << longname << " "
	   << level << G4endl;
  }

  vpsi ancestors = ExtractAncestorsRequested(longname);
  if (level < 0 || (size_t)level >= ancestors.size()) {
    if (verboseLevel) G4cerr << " no ancestor at level " << level << G4endl;
    return G4String("");
  }

  return AppendCopyNo(ancestors[level].first, ancestors[level].second);
}


// Attach copy number to end of input string

G4String DemoGeometryTools::AppendCopyNo(const G4String& pvname, G4int cn) const {
  std::ostringstream nameS(pvname, std::ios::ate);
  nameS << ":" << cn;
  return nameS.str();
}

G4int DemoGeometryTools::ExtractCopyNo(const G4String& pvname) const {
  G4int fnumb = pvname.rfind(':');
  if (fnumb == -1) return 0;		// No copy number included in name

  return atoi((pvname.substr(fnumb+1, pvname.size()-1)).c_str());
}


// Generate indention level for printing out tree structures

G4String DemoGeometryTools::IndentLeaf(unsigned int leafDepth) const {
  G4String spaces;
  return spaces.insert(0, leafDepth*2, ' ');
}


// Compare words, allowing for possible wildcards in first argument

bool DemoGeometryTools::
WildcardMatch(const G4String& glob, const G4String& word) const {
  // fnmatch() returns 0 for success, non-zero error code for failure
  return (0 == fnmatch(glob.c_str(), word.c_str(), 0));
}