In This Package:

TouchableToDetectorElement.cc

Go to the documentation of this file.

#include "TouchableToDetectorElement.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>

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


StatusCode TouchableToDetectorElement::G4VolumeToDetDesc( 
        const G4VPhysicalVolume* inVol, 
        const IPVolume* &outVol )        
{
  ILVolume::ReplicaPath rp;
  return G4VolumeToDetDesc(inVol,outVol,rp);
}

StatusCode TouchableToDetectorElement::G4VolumeToDetDesc( 
        const G4VPhysicalVolume* inVol, 
        const IPVolume* &outVol,
        ILVolume::ReplicaPath& outPath )
{
  
  if(!inVol) {
    return StatusCode::FAILURE;
  }
  
  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)) {
    err() << "Can't find lvolume at path " << lname << " when parsing physvol " << inVol->GetName() << endreq;
    return StatusCode::FAILURE;
  }    
  lv= getDet<ILVolume>(lname);
  
  // 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:
    outPath.push_back(replica);      
    
    // Advance lvolume
    lv = pv->lvolume();
  }

  outVol = pv;

  return StatusCode::SUCCESS;
}


// ---------------------------------------------------------------------------------
template <class T> 
StatusCode TouchableToDetectorElement::FindObjectsInDirectory(
                  const std::string& dirname,
                  std::vector<const T*>& list
                  )
{
  info() << "Looking for structures in " << dirname << endreq;
  if (!existDet<DataObject>(dirname)) {
    err() << "No " << dirname << " in TDS" << endreq;
    return StatusCode::FAILURE;
  }
  const DataObject *d = getDet<DataObject>(dirname);

  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 TouchableToDetectorElement::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)
{
  TouchableToDetectorElement::SpecifiedVolume g4place; 
  StatusCode sc = GetSpecifiedVolume(inHistory, g4place);
  if(sc.isFailure()) return sc;

  // At the moment, the best possible match is to the top Element, which has only basic compatibility.
  outCompatibility = 0x7fffffff; // Maximum possible int.
  outElement = g4place.mTop;

  // Find a list we want to compare to.
  //debug() << "Getting list of elements" << endreq;
  std::vector<const IDetectorElement*> elements;
  
  // Search each path, looking for things.
  for(unsigned int i = 0; i< inPaths.size(); ++i) {
    sc = FindObjectsInDirectory<IDetectorElement>(inPaths[i],elements);
    if(sc.isFailure()) return sc;
  }

  // Now compare each one to our current.
  for(unsigned int i=0; i<elements.size(); ++i) {
    TouchableToDetectorElement::SpecifiedVolume trialVolume;
    sc = GetSpecifiedVolume(elements[i], g4place.mTop, trialVolume);
    if(sc.isFailure()) continue;

    int compat = Compatability(g4place,trialVolume);
    if((compat >= 0) && (compat<outCompatibility)) {
      outElement = elements[i];
      outCompatibility = compat;
    }
  }

  return StatusCode::SUCCESS; 
}

StatusCode TouchableToDetectorElement::GetSpecifiedVolume(
            const G4TouchableHistory* inHistory, 
            TouchableToDetectorElement::SpecifiedVolume& outVol)
{
  //debug() << "Building SpecVol from G4TouchableHistory" << endreq;
  if(!inHistory) {
    err() << "No good TouchableHistory provided." << endreq;
    return StatusCode::FAILURE;
  }

  // Try to find the topmost level, which should be a DetectorElement.
  G4VPhysicalVolume* top = inHistory->GetVolume(inHistory->GetHistoryDepth()-1);
  if(!top) {
    err() << "Got null G4VPhysicalVolume pointer in the bottom of G4TouchableHistory. Huh???" << endreq;
    return StatusCode::FAILURE;
  }
  if(!exist<IDetectorElement>(detSvc(),top->GetName()),false) {
    err() << "Can't find detector element at the bottom of the G4TouchableHistory, name = " << top->GetName() << endreq;
    return StatusCode::FAILURE;
  }

  outVol.mTop = getDet<IDetectorElement>(top->GetName());
  outVol.mRpath.clear();

   // Ok, now parse the rest of the history.
   // This could in principle be sped up by caching the resulting sub-SpecfiedVolumes and storing keyed by
   // G4VPhysicalVolume pointer. But that would result in a fair number of things laying around, and it's
   // not clear if the speed is important yet.

  for(int ind=inHistory->GetHistoryDepth()-2; ind>=0;  --ind){
    G4VPhysicalVolume* g4pv = inHistory->GetVolume(ind);
     // Now, parse the volume name.
    ILVolume::ReplicaPath subpath;
    const IPVolume* pv;
    StatusCode sc = G4VolumeToDetDesc(g4pv,pv,subpath);
    if(sc.isFailure()) return sc;

    // Tack this onto the end of the list.
    outVol.mRpath.insert(outVol.mRpath.end(),subpath.begin(),subpath.end());
    
  } // End loop over history
  
  return StatusCode::SUCCESS;
}

StatusCode TouchableToDetectorElement::GetSpecifiedVolume(const IDetectorElement* inElem,   const IDetectorElement* inTopElem, TouchableToDetectorElement::SpecifiedVolume& outVol)
{
  const IGeometryInfo* topGeo = inTopElem->geometry();
  const IGeometryInfo* trialGeo = inElem->geometry();

  // Try to build up a full replica path.
  outVol.mRpath.clear();
  outVol.mTop = 0;
  
  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;
    outVol.mRpath.insert(outVol.mRpath.begin(),rpath.begin(),rpath.end());
    if(0==curGeo) break;
  }
  if(curGeo==topGeo) outVol.mTop = inTopElem;
  else outVol.mTop = 0; // 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.
  
  // debug() << "Built SpecVol from " << inElem->name() << ":" << endreq;
  // if(outVol.mTop)
  //    debug() << outVol.mTop->name();
  // for(unsigned int i=0;i<outVol.mRpath.size();i++) debug() << " [" << outVol.mRpath[i] << "]";
  // debug() << endreq;
  
  return StatusCode::SUCCESS; 
}

int   TouchableToDetectorElement::Compatability(const TouchableToDetectorElement::SpecifiedVolume& inPlace, const TouchableToDetectorElement::SpecifiedVolume& inContainer)
{
  // Check to ensure that they share a top volume. This will fail for DetElements at, for instance, different
  // sites than the simulated site.
  if(inPlace.mTop != inContainer.mTop) return -1; // They don't match at the top.
  
  // 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
  if(inContainer.mRpath.size() > inPlace.mRpath.size()) return -1;

  // Check through to see if it matches at every level.
  for(unsigned int i=0; i< inContainer.mRpath.size(); ++i ){
    if(inContainer.mRpath[i] != inPlace.mRpath[i]) return -1; // There's a mismatch, throw this one out.
  }
  
  // Ok, it checks at every level of the container. 
  // The goodness of the match is the specificity of the place minus the specificity of the container.
  // This is zero for a perfect match.
  return inPlace.mRpath.size() - inContainer.mRpath.size();
}

---