TouchableToDetectorElementFast.cc

Go to the documentation of this file.

#include "TouchableToDetectorElementFast.h"

#include "G4TouchableHistory.hh"

#include "DetDesc/ILVolume.h"
#include "DetDesc/IPVolume.h"
#include "DetDesc/IDetectorElement.h"
#include "DetDesc/IGeometryInfo.h"

#include "GaudiKernel/DataObject.h"
#include "GaudiKernel/IRegistry.h"
#include "GaudiKernel/RegistryEntry.h"
#include "GaudiKernel/DataObject.h"

#include <vector>

TouchableToDetectorElementFast::TouchableToDetectorElementFast(
                          const std::string& type,
                          const std::string& name,
                          const IInterface* parent)
       : GaudiTool(type,name,parent)
{
  declareInterface<ITouchableToDetectorElement>(this);
  ClearCache();  
}


StatusCode TouchableToDetectorElementFast::G4VolumeToDetDesc( 
        const G4VPhysicalVolume* inVol, 
        const IPVolume* &outVol )
{
  
  Relation* relation;
  StatusCode sc = GetRelation(inVol,relation);
  if(sc.isFailure()) {
    // This error generates too much verbiage.
    //err() << "Couldn't find IPVolume for " << inVol->GetName() << endreq;
    outVol = 0;
    return sc;
  }
  outVol = relation->mPhysVol;
  if(outVol==0) return StatusCode::FAILURE;

  return StatusCode::SUCCESS;
}


// ---------------------------------------------------------------------------------
template <class T> 
StatusCode TouchableToDetectorElementFast::FindObjectsInDirectory(
                  const std::string& dirname,
                  std::list<const T*>& list
                  )
{
  DataObject* d = 0;
  StatusCode sc = detSvc()->retrieveObject(dirname,d);

  if(sc.isFailure()) return sc;

  const T* t = dynamic_cast<const T*>(d);
  if(t) list.push_back(t);

  IRegistry* dr = d->registry();
  using namespace DataSvcHelpers;
  RegistryEntry* dre = dynamic_cast<RegistryEntry*>(dr);
  if (!dre) {
    err() << "Failed to get RegistryEntry on DataObject" << d << endreq;
    return StatusCode::FAILURE;
  }

  RegistryEntry::Iterator it = dre->begin(), done = dre->end();
  for (; it != done; ++it) {
    std::string id = (*it)->identifier();
    StatusCode sc = FindObjectsInDirectory<T>(id,list);
    if(sc.isFailure()) return sc;
  }

  return StatusCode::SUCCESS;  
}


StatusCode TouchableToDetectorElementFast::GetBestDetectorElement(
                          const G4TouchableHistory* inHistory,  // The current particle location
                          const std::vector<std::string>& inPaths,// Name(s) of specific detElements, or paths to be searched
                          const IDetectorElement* &outElement,  // output: The best element (may be zero!) 
                          int& outCompatibility )              // output: the goodness of the match (higher is better)
{
  // First, look at the touchable volume and get the top-level directory, which should be our current DetectorElement.
  G4VPhysicalVolume* g4top = inHistory->GetVolume(inHistory->GetHistoryDepth()-1);
  if(!g4top) {
    err() << "Got null G4VPhysicalVolume pointer in the bottom of G4TouchableHistory. Huh???" << endreq;
    return StatusCode::FAILURE;
  }


  if( mWorldElementName != g4top->GetName() ) {
    // Uh oh!  Our cache is bad. Clear it.
    ClearCache();
    // Set the new one.
    if(!exist<IDetectorElement>(detSvc(),g4top->GetName()),false) {
      err() << "Can't find detector element at the bottom of the G4TouchableHistory, name = " << g4top->GetName() << endreq;
      return StatusCode::FAILURE;
    }

    mWorldElement = getDet<IDetectorElement>(g4top->GetName());    
    mWorldElementName = mWorldElement->name();
    if(!mWorldElement) {
      err() << "Unable to find world element. G4 top is " << g4top->GetName() << endreq;
    }
  }

  // Ok, now build up the total relative path from the world element to the top of the TouchableHistory.
  ILVolume::ReplicaPath touchablePath;
  for(int ind=inHistory->GetHistoryDepth()-2; ind>=0;  --ind){
    G4VPhysicalVolume* g4pv = inHistory->GetVolume(ind);
    // Look it up.
    Relation* rel;
    StatusCode sc = GetRelation(g4pv,rel);
    if(sc.isFailure()){
      err() << "Failure from GetRelation(" << g4pv->GetName() << ")" << endreq;
      return sc;
    }

    // Tack this onto the end of the list.
    touchablePath.insert(touchablePath.end(),rel->mRpath.begin(),rel->mRpath.end());
    
  } // End loop over history
  

  // So, now we have our TouchablePath. Let's now look at each detelem.


  // Is this a new search, or does it match an old one?
  if(mLastSearchPaths == inPaths) {
    // Do nothing! Our cache is still good!
  } else {
     mElementList.clear();
    // Search each path, looking for things.
    for(unsigned int i = 0; i< inPaths.size(); ++i) {
      info() << "Looking for structures in path " << inPaths[i] << endreq;
      StatusCode sc = FindObjectsInDirectory<IDetectorElement>(inPaths[i],mElementList);
      if(sc.isFailure()) {
        warning() << "Couldn't resolve DetectorElement search path " <<inPaths[i] << endreq;
      }
    }
    if (inPaths.size()) {
        // w/out this message a quiet, long running G4 looks like we are still taking up time
        info() << "Finished search structures in " << inPaths.size() << " paths" << endreq;
    }
    mLastSearchPaths = inPaths;    
  }

  // Default: none is best:
  // At the moment, the best possible match is to the top Element, which has only basic compatibility.
   outCompatibility = 0x7fffffff; // Maximum possible int.
   outElement = mWorldElement;
  
  // Now compare each one to our current.
  ElementList_t::iterator it = mElementList.begin();
  Relation* rel;
  while( it != mElementList.end() ) {
    if(*it==0) {
      err() << "Found a zero DetElem in our list. How did that happen?"<< endreq;
      ++it; continue;
    }
    StatusCode sc = GetRelation(*it,rel);
    if(sc.isFailure()) { ++it; continue; }
    
    // Is it supported by our volume?
    if(rel->isNull()){
      // delete this entry.. it's doing us no good.
      it = mElementList.erase(it);
    } else {
      ILVolume::ReplicaPath &elementPath = rel->mRpath;
    
      int compat = Compatability(touchablePath,elementPath);
      if((compat>=0) && (compat<outCompatibility)) {
        outElement = *it;
        outCompatibility = compat;
      }
      ++it;
    }
  }

  return StatusCode::SUCCESS; 
}


int  TouchableToDetectorElementFast::Compatability(const ILVolume::ReplicaPath& inPlace, 
                                                   const ILVolume::ReplicaPath& inContainer)
{
  // Check to see that the container is at least as specific as the place. 
  // If not, the container specifies a deeper detector element than the place
  // e.g. This condition catches container=PMT, place=rock
  unsigned int sizeContainer = inContainer.size();
  unsigned int sizePlace     = inPlace.size();

  if(sizeContainer > sizePlace) return -1;

  // Check through to see if it matches at every level.
  ILVolume::ReplicaPath::const_iterator container_it = inContainer.begin();
  ILVolume::ReplicaPath::const_iterator place_it     = inPlace.begin();
  ILVolume::ReplicaPath::const_iterator container_end = inContainer.end();
  while(container_it != container_end) {
    if(*container_it != *place_it) return -1;
    container_it++;
    place_it++;
  }
  
  // Ok, it checks at every level of the container. The goodness of the match is the specificity of the container.
  return sizePlace - sizeContainer;
}



StatusCode TouchableToDetectorElementFast::GetRelation(const G4VPhysicalVolume* inVol, Relation* &outRelation )
{
 
  if(!inVol) {
    return StatusCode::FAILURE;
  }
  
  G4ToRelationMap_t::iterator it = mG4ToRelationMap.find(inVol);
  if(it != mG4ToRelationMap.end()) {
    outRelation = &(it->second);
    return StatusCode::SUCCESS;
  }
  
  // Create the new lookup.
  Relation rel;
  
  // Ok, we failed to find it. Time to do some actual work.
  std::string pvpath = inVol->GetName();
  
  const ILVolume* lv = 0;
  const IPVolume* pv = 0;
   // This is a top-level call
   // Look for an explicit logical volume name.
  std::string lname = pvpath;
  std::string::size_type p = pvpath.find_first_of('#');
  if(p!= std::string::npos) {
    lname = pvpath.substr(0,p);
    pvpath.erase(0,p+1); // strip off this part.
  } 

  // Find the logical volume.
  if(!exist<ILVolume>(detSvc(),lname,false)) {
    // It's possible that we've hit the top (geant) volume, which is actually a DetectorElement.
    if(exist<IDetectorElement>(detSvc(),lname,false)) {
      // Ok, this is a special case.
      IDetectorElement* de = getDet<IDetectorElement>(lname);
      IGeometryInfo* geo = de->geometry();
      assert(geo);
      if(geo->hasSupport()) {
        // This may be the top geant volume, but it's supported in DetDesc space. 
        // So, find the support, and work down.
        ILVolume::ReplicaPath rpath;
        IGeometryInfo* nextgeo;
        geo->location(nextgeo, rpath);
        // take off one level.
        assert(rpath.size()>0);
        ILVolume::ReplicaType last = rpath.back();
        rpath.pop_back();
        lv = nextgeo->lvolume(nextgeo,rpath);
        pv = (*lv)[last];
        // output.
        rel.mRpath = ILVolume::ReplicaPath(); // this case isn't valid in the overall tree
        rel.mLogVol = lv;
        rel.mPhysVol = pv;
        it =  mG4ToRelationMap.insert(G4ToRelationMap_t::value_type(inVol,rel)).first;
        outRelation = &(it->second);
        return StatusCode::SUCCESS;
      }
      // Ok, this looks like something strange.
      if(geo->hasLVolume()) {
         lv = de->geometry()->lvolume();
         // just try the first physvol. this is correct for the genuine top volume.
         pv = (*lv)[0];
         rel.mRpath = ILVolume::ReplicaPath(); // this case isn't valid in the overall tree
         rel.mLogVol = lv;
         rel.mPhysVol = pv;
         it =  mG4ToRelationMap.insert(G4ToRelationMap_t::value_type(inVol,rel)).first;
         outRelation = &(it->second);
         return StatusCode::SUCCESS;
       }

    }
    // Wierd.. it's not an lvolume or a detelement, or the detelement couldn't be evaluated.
    err() << "Can't find lvolume at path " << lname << endreq;    
    rel.mLogVol = 0;
    rel.mPhysVol = 0;
    rel.mRpath = ILVolume::ReplicaPath();
    it =  mG4ToRelationMap.insert(G4ToRelationMap_t::value_type(inVol,rel)).first;
    return StatusCode::FAILURE;
  } 
  lv= getDet<ILVolume>(lname);
  // This is the correct supporting volume.
  rel.mLogVol = lv;
  
  // Keep eating the remaining string until we have fully parsed the remaining logical/physical path sequence.
  while(pvpath.size() > 0) {
    // Make sure the current logical volume is good.
    if(!lv) {
      err() << "Can't find lvolume." << endreq;
      return StatusCode::FAILURE;
    }

    // parse out the next bit before a hash mark:
    std::string pname = pvpath;
    std::string::size_type p = pvpath.find_first_of('#');
    if(p != std::string::npos ) {
      pname.erase(p,std::string::npos);
      pvpath.erase(0,p+1); // Ok, take this stuff off the front.
    } else {
      pvpath.erase();
    }
  
    ILVolume::ReplicaType replica = 0;
    pv = 0;
    for ( ILVolume::ReplicaType i = 0; i< lv->noPVolumes(); ++i ){
      const IPVolume* apv = lv->pvolume(i);
      if(!apv) {
        err() << "Bad link to IPVolume from " << lv->name() << endreq;
        return StatusCode::FAILURE;
      }
      // Look for a match.
      if(apv->name() == pname){ 
        pv = apv; 
        replica = i; 
        break;
      }
    } // End of loop over physical volumes in lv
    
    if(!pv) {
      err() << "Couldn't find a match for physical volume with name " << pname << " in lvolume name " << lv->name() << endreq;
      return StatusCode::FAILURE;
    }
    // Finally, this is what we were after:
    rel.mRpath.push_back(replica);      
    
    // Advance lvolume
    lv = pv->lvolume();
  }

  // The physical volume
  rel.mPhysVol = pv;
  
  it =  mG4ToRelationMap.insert(G4ToRelationMap_t::value_type(inVol,rel)).first;
  outRelation = &(it->second);
  
  return StatusCode::SUCCESS;
  
  
}

StatusCode TouchableToDetectorElementFast::GetRelation(const IDetectorElement* inElem, Relation* &outRelation )
{
  if(!inElem) return StatusCode::FAILURE;
  
  // Check the cache... is it there?
  DetElemToRelationMap_t::iterator it = mDetElemToRelationMap.find(inElem);
  if(it != mDetElemToRelationMap.end()) {
    outRelation = &(it->second);
    return StatusCode::SUCCESS;
  }
  
  // It's not in the cache, so we actually have to do some work.
  Relation rel;
  
  const IGeometryInfo* topGeo = mWorldElement->geometry();
  const IGeometryInfo* trialGeo = inElem->geometry();

  // Try to build up a full replica path.
  rel.mRpath.clear();
  
  const IGeometryInfo* curGeo = trialGeo;
  while(curGeo!=topGeo) {
    //std::string sstart;
    //curGeo->location(sstart, rpath);
    
    ILVolume::ReplicaPath rpath;
    IGeometryInfo* nextGeo;
    curGeo->location(nextGeo, rpath);
    curGeo = nextGeo;
    rel.mRpath.insert(rel.mRpath.begin(),rpath.begin(),rpath.end());
    if(0==curGeo) break;
  }
  if(curGeo==topGeo) {
    // We reached the top and it's a match to our world volume. Success!
    rel.mLogVol = curGeo->lvolume();
    // Navigate back to our physical volume.
    const ILVolume* lv = rel.mLogVol;
    const IPVolume* pv = 0;
    for(unsigned int i=0;i<rel.mRpath.size();i++) {
      if(lv)
        pv = (*lv)[rel.mRpath[i]];
      if(pv)
        lv = pv->lvolume();
    }
    rel.mPhysVol = pv;
  }
  else {
    //
    // This is an error, but I don't know how to fix it.. how do you go from GeometryInfo to a detector element, 
    // especially if it's an 'orphan' or a 'ghost'? Maybe it doesn't matter, since it isnt' used.
  }
  
  it = mDetElemToRelationMap.insert(DetElemToRelationMap_t::value_type(inElem,rel)).first;
  outRelation = &(it->second);

  return StatusCode::SUCCESS; 
  
}


StatusCode TouchableToDetectorElementFast::ClearCache()
{
  mWorldElement =0;
  mWorldElementName = "";
  mG4ToRelationMap.clear();
  mDetElemToRelationMap.clear();
  mLastSearchPaths.clear();
  mElementList.clear();
  return StatusCode::SUCCESS;
}