In This Package:

TH2DE.cc

Go to the documentation of this file.

#include "TH2DE.h"

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

#include "G4TouchableHistory.hh"
#include "G4VPhysicalVolume.hh"

#include <sstream>

using namespace std;

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

TH2DE::~TH2DE()
{
    this->ClearCache();
}

static void split(vector<string>& out, const string& src, char del)
{
    out.clear();
    string::size_type pos = 0, siz = src.size();
    while (pos < siz) {
        string::size_type found = src.find(del,pos);
        if (found == string::npos) {
            out.push_back(src.substr(pos,siz-pos));
            break;
        }
        out.push_back(src.substr(pos,found-pos));
        pos = found+1;
    }
}

static string dumpstr(const G4TouchableHistory* g4hist)
{
    stringstream ss;
    int siz = g4hist->GetHistoryDepth();
    for (int ind=0; ind < siz; ++ind) {
        ss << g4hist->GetVolume(ind)->GetName() << "\n";
    }
    ss << ends;
    return ss.str();
}

StatusCode TH2DE::GetBestDetectorElement(const G4TouchableHistory* inHistory,
                                        const std::vector<std::string>& /*notused*/,
                                        const IDetectorElement* &outElement,
                                        int& outCompatibility)
{
    if (!inHistory->GetHistoryDepth()) {
        warning() << "TH2DE::GetBestDetectorElement given an empty history" << endreq;
        return StatusCode::FAILURE;
    }

    const IDetectorElement* de = this->CheckCache(inHistory);
    if (de) {
        outElement = de;
        return StatusCode::SUCCESS;
    }
        
    TouchableHistory_t th;      // for cache
    NameHistory_t name_history; // always read backwards, [0] = daughter, [n] = ancestors
    const int depth = inHistory->GetHistoryDepth();
    for (int ind=0; ind < depth; ++ind) {
        G4VPhysicalVolume* g4pv = inHistory->GetVolume(ind);
        th.push_back(g4pv);

        string full_name = g4pv->GetName();

        verbose() << ind << ": " << full_name << endreq;

        vector<string> names;
        split(names,full_name,'#');

#if 0
        for (size_t iname=0; iname<names.size(); ++iname) {
            debug() << iname << ": [" << names[iname] << "]";
        }
        debug() << endreq;
#endif

        if (depth-ind == 1) { // have /dd/Structure top level DE
            if (names.size() != 1) {
                warning() << "got unknown type at top of history: " << full_name << ", " <<names.size() << endreq;
                return StatusCode::FAILURE;
            }

            DataObject* obj = 0;
            StatusCode sc = detSvc()->retrieveObject(names[0],obj);
            if (sc.isFailure()) {
                warning() << "failed to get DetectorElement: " << names[0] << endreq;
                return StatusCode::FAILURE;
            }
            de = dynamic_cast<const IDetectorElement*>(obj);
            if (!de) {
                warning() << "failed to dynamic_cast<DetectorElement*>: " << names[0] << endreq;
                return StatusCode::FAILURE;
            }
            break;
        }
            

        if (names.size() > 2) {   // have pv1#pv2#pv3 style path
            DataObject* obj = 0;
            StatusCode sc = detSvc()->retrieveObject(names[0],obj);
            if (sc.isFailure()) {
                warning() << "failed to get LVolume: " << names[0] << endreq;
                return StatusCode::FAILURE;
            }
            const ILVolume* lv = dynamic_cast<const ILVolume*>(obj);
            if (!lv) {
                warning() << "failed to dynamic_cast<LVolume*>: " << names[0] << endreq;
                return StatusCode::FAILURE;
            }
            const IPVolume* pv = (*lv)[names[1]];
            lv = pv->lvolume();
            NameHistory_t reverse;
            for (string::size_type iname = 2; iname < names.size(); ++iname) {
                pv = (*lv)[names[iname]];
                reverse.push_back(LvPvPair_t(lv->name(),pv->name()));
                lv = pv->lvolume();
            }
            name_history.insert(name_history.end(),reverse.rbegin(),reverse.rend());
        }

        name_history.push_back(LvPvPair_t(names[0],names[1]));

    } // loop over history

    if (!de) {
        const int depth = inHistory->GetHistoryDepth();
        warning() << "failed to find top level supporting DetectorElement, history has depth of " << depth << "\n";
        for (int ind=0; ind < depth; ++ind) {
            G4VPhysicalVolume* g4pv = inHistory->GetVolume(ind);
            warning() << "\n\t(" << ind << ") " << g4pv->GetName();
        }
        warning() << endreq;
        return StatusCode::FAILURE;
    }

    de = this->FindDE(de,name_history);
    if (!de) {
        warning() << "failed to find DetectorElement for TouchableHistory:\n" << dumpstr(inHistory) << endreq;
        return StatusCode::FAILURE;
    }

    m_THcache[th] = de;
    outElement = de;
    outCompatibility = name_history.size();
    return StatusCode::SUCCESS;
}

int TH2DE::InHistory(const IDetectorElement* de, const NameHistory_t& name_history)
{
    const IGeometryInfo* gi = de->geometry();
    if (!gi->hasSupport()) return -1;

    const ILVolume::ReplicaPath& rpath = gi->supportPath();
    const IGeometryInfo* support_gi = gi->supportIGeometryInfo();
    const ILVolume* lv = support_gi->lvolume();

    verbose() << "InHistory de=" << de->name() << endreq;

    // Walk the DE's support and match against name_history;
    size_t index = name_history.size();
    for (size_t ind = 0; index && ind < rpath.size(); ++ind) {
        IPVolume* pv = lv->pvolumes()[rpath[ind]];

        --index;
        const LvPvPair_t& check = name_history[index];

        verbose() << "("<<index<<") lvpv=" << lv->name() << "#" << pv->name() 
                  << " =?= " 
                  << check.first << "#" << check.second << endreq;

        if (lv->name() != check.first) return -1;
        if (pv->name() != check.second) return -1;

        lv = pv->lvolume();
    }
    return index;
}

const IDetectorElement* TH2DE::FindChildDE(const IDetectorElement* de, NameHistory_t& name_history)
{
    IDetectorElement::IDEContainer& children = de->childIDetectorElements();
    size_t nchildren = children.size();
    verbose() << "Finding children from " << nchildren << endreq;
    for (size_t ichild = 0; ichild < nchildren; ++ichild) {
        IDetectorElement* child = children[ichild];
        int index = this->InHistory(child,name_history);

        if (index<0) continue;

        verbose() << "Found child: " << child->name() << " at index " <<index<<" lv#pv=" << name_history[index].first << "#"<< name_history[index].second << endreq;

        // strip off used history
        while ((int)name_history.size() > index) {
            LvPvPair_t lvpv = name_history.back();
            verbose () << "\tpoping: index="<<index<<", size=" << name_history.size() << " [" << lvpv.first << "#" << lvpv.second << "]" << endreq;
            name_history.pop_back();
        }
        return child;
    }
    return 0;
}

const IDetectorElement* TH2DE::FindDE(const IDetectorElement* de, NameHistory_t& name_history)
{
    // If exhausted the NH then the current DE must be the one
    if (!name_history.size()) return de;


#if 0
    debug() << "FindDE: " << de->name() << endreq;
    for (size_t inh=0; inh<name_history.size(); ++inh) {
        debug() << inh <<": " << name_history[inh].first << "#" << name_history[inh].second << endreq;
    }
#endif

    LvPvPair_t lvpv(name_history.back().first,name_history.back().second);

    string de_lvname = de->geometry()->lvolumeName();
    if (de_lvname != lvpv.first) {
        warning() << "The given DE's LV does not match LV from top of history: "
                  << de_lvname <<" != "<< lvpv.first << endreq;
        return 0;
    }

    // Find immediate child that points into the history
    const IDetectorElement* child_de = this->FindChildDE(de,name_history);

    // If one found, recurse.
    if (child_de) return this->FindDE(child_de,name_history);

    // If we get here, we have reached the end of possible DEs.
    // Pop off the touchable history coresponding to current DE and return;
    name_history.pop_back();
    return de;
}

const IDetectorElement* TH2DE::CheckCache(const G4TouchableHistory* g4hist)
{
    TouchableHistory_t th;
    const int depth = g4hist->GetHistoryDepth();
    for (int ind=0; ind < depth; ++ind) {
        G4VPhysicalVolume* g4pv = g4hist->GetVolume(ind);
        th.push_back(g4pv);
    }
    return m_THcache[th];
}

StatusCode TH2DE::G4VolumeToDetDesc(const G4VPhysicalVolume* inVol,
                                    const IPVolume* &outVol)
{
    
    if (!inVol) {
      warning() << "TH2DE:  No PhysicalVolume detected!   " <<endreq;
      return StatusCode::FAILURE;
    }
    outVol = m_PVcache[inVol];
    if (outVol) return StatusCode::FAILURE;

    vector<string> names;
    split(names,inVol->GetName(),'#');

    DataObject* obj = 0;
    StatusCode sc = detSvc()->retrieveObject(names[0],obj);
    if (sc.isFailure()) return sc;

    if (names.size() == 1) {    // this PV is named after a detector element
        const IDetectorElement* de = dynamic_cast<const IDetectorElement*>(obj);
        const IGeometryInfo* gi = de->geometry();
        const ILVolume* lv_target = gi->lvolume();
        const ILVolume* lv = gi->supportIGeometryInfo()->lvolume();
        const ILVolume::ReplicaPath& rpath = gi->supportPath();
        const IPVolume* pv = 0;
        for (size_t ind=0; ind < rpath.size(); ++ind) {
            pv = lv->pvolumes()[rpath[ind]];
            if (!pv) return StatusCode::FAILURE;
            lv = pv->lvolume();
            if (!lv) return StatusCode::FAILURE;
        }
        if (lv != lv_target) {
            warning() << "Don't come full circle: " << lv_target->name() << " != " << lv->name() << endreq;
            return StatusCode::FAILURE;
        }
        outVol = pv;
    }
    else {
        const ILVolume* lv = dynamic_cast<const ILVolume*>(obj);
        const IPVolume* pv = (*lv)[names[1]];
        if (!pv) return StatusCode::FAILURE;

        if (names.size() > 2) {   // have pv1#pv2#pv3 style path
            lv = pv->lvolume();
            if (!lv) return StatusCode::FAILURE;

            for (string::size_type iname = 2; iname < names.size(); ++iname) {
                pv = (*lv)[names[iname]];
                if (!pv) return StatusCode::FAILURE;
                lv = pv->lvolume();
                if (!lv) return StatusCode::FAILURE;
            }
        }
        outVol = pv;
    }
    if (outVol) {
        m_PVcache[inVol] = outVol;
        return StatusCode::SUCCESS;
    }
    return StatusCode::FAILURE;
}
 
StatusCode TH2DE::ClearCache()
{
    m_THcache.clear();
    m_PVcache.clear();
    return StatusCode::SUCCESS;
}

---