diff --git a/DDRec/include/DDRec/CellIDPositionConverter.h b/DDRec/include/DDRec/CellIDPositionConverter.h index c0e378cd2ecc4b5e18b3e3079532293329c7d7f6..17965823d579c22fb3446dcac8a7cf39c1dba6e4 100644 --- a/DDRec/include/DDRec/CellIDPositionConverter.h +++ b/DDRec/include/DDRec/CellIDPositionConverter.h @@ -9,7 +9,6 @@ #include <set> #include <string> -class TGeoManager; namespace DD4hep { namespace DDRec { @@ -40,9 +39,6 @@ namespace DD4hep { /// Destructor virtual ~CellIDPositionConverter(){} ; - /// returns the global cell ID from a given global position - CellID cellID(const Geometry::Position& global) const; - /** Return the nominal global position for a given cellID of a sensitive volume. * No Alignment corrections are applied. * If no sensitive volume is found, (0,0,0) is returned. @@ -55,6 +51,15 @@ namespace DD4hep { */ Geometry::Position position(const CellID& cellID) const; + + /** Return the global cellID for the given global position. + * Note: this call is rather slow - only use it when really needed ! + * + */ + CellID cellID(const Geometry::Position& global) const; + + + /** Find the context with DetElement, placements etc for a given cellID of a sensitive volume. * Returns NULL if not found (e.g. if the cellID does not correspond to a sensitive volume). */ diff --git a/DDRec/src/CellIDPositionConverter.cpp b/DDRec/src/CellIDPositionConverter.cpp index eef501f60906b6413dec229609d51edc32388b11..7af0665599d077d3a57eda187c1ff43cda5e4f0c 100644 --- a/DDRec/src/CellIDPositionConverter.cpp +++ b/DDRec/src/CellIDPositionConverter.cpp @@ -4,6 +4,8 @@ #include "DD4hep/LCDD.h" #include "DD4hep/objects/VolumeManagerInterna.h" +#include "TGeoManager.h" + namespace DD4hep { namespace DDRec { @@ -83,43 +85,107 @@ namespace DD4hep { CellID CellIDPositionConverter::cellID(const Position& global) const { - + CellID result(0) ; - DetElement motherDet = _lcdd->world() ; // could also start from an arbitrary DetElement here !? - - DetElement det = findDetElement( global , motherDet ) ; - - if( ! det.isValid() ) - return result ; - - double g[3], e[3] , l[3] ; - global.GetCoordinates( g ) ; - det.nominal().worldTransformation().MasterToLocal( g, e ); + TGeoManager *geoManager = _lcdd->world().volume()->GetGeoManager() ; - PlacedVolume::VolIDs volIDs ; - - PlacedVolume pv = findPlacement( Position( e[0], e[1] , e[2] ) , det.placement() , l , volIDs ) ; + PlacedVolume pv = geoManager->FindNode( global.x() , global.y() , global.z() ) ; if( pv.isValid() && pv.volume().isSensitive() ) { - + + TGeoHMatrix* m = geoManager->GetCurrentMatrix() ; + + double g[3], l[3] ; + global.GetCoordinates( g ) ; + m->MasterToLocal( g, l ); + Geometry::SensitiveDetector sd = pv.volume().sensitiveDetector(); Readout r = sd.readout() ; - VolumeID volIDElement = det.volumeID() ; - // add the placed volumes volIDs: - VolumeID volIDPVs = r.idSpec().encode( volIDs ) ; + // collect all volIDs for the current path + PlacedVolume::VolIDs volIDs ; + volIDs.insert( std::end(volIDs), std::begin(pv.volIDs()), std::end(pv.volIDs())) ; + + TGeoPhysicalNode pN( geoManager->GetPath() ) ; - result = r.segmentation().cellID( Position( l[0], l[1], l[2] ) , global, ( volIDElement | volIDPVs ) ); + unsigned motherCount = 0 ; + + while( pN.GetMother( motherCount ) != NULL ){ + + PlacedVolume mPv = pN.GetMother( motherCount++ ) ; + + if( mPv.isValid() && pN.GetMother( motherCount ) != NULL ) // world has no volIDs + volIDs.insert( std::end(volIDs), std::begin(mPv.volIDs()), std::end(mPv.volIDs())) ; + } - // } else { - // std::cout << " *** ERROR : found non-sensitive Placement " << pv.name() - // << " for point " << global << std::endl ; + VolumeID volIDPVs = r.idSpec().encode( volIDs ) ; + + result = r.segmentation().cellID( Position( l[0], l[1], l[2] ) , global, volIDPVs ); } return result ; } + // CellID CellIDPositionConverter::cellID(const Position& global) const { + + // CellID result(0) ; + + // DetElement motherDet = _lcdd->world() ; // could also start from an arbitrary DetElement here !? + + // DetElement det = findDetElement( global , motherDet ) ; + + // if( ! det.isValid() ) + // return result ; + + // double g[3], e[3] , l[3] ; + // global.GetCoordinates( g ) ; + // det.nominal().worldTransformation().MasterToLocal( g, e ); + + // PlacedVolume::VolIDs volIDs ; + + // PlacedVolume pv = findPlacement( Position( e[0], e[1] , e[2] ) , det.placement() , l , volIDs ) ; + + // TGeoManager *geoManager = det.volume()->GetGeoManager() ; + // // TGeoManager *geoManager = _lcdd->world().volume()->GetGeoManager() ; + + // PlacedVolume pv1 = geoManager->FindNode( global.x() , global.y() , global.z() ) ; + + // // std::cout << " -> TGM : " << pv1.name() << " valid: " << pv1.isValid() << " sensitive: " + // // << (pv1.isValid() ? pv1.volume().isSensitive() : false ) << std::endl ; + // // std::cout << " -> FG : " << pv.name() << " valid: " << pv.isValid() << " sensitive: " + // // << (pv.isValid() ? pv.volume().isSensitive() : false ) << std::endl ; + + + + // if( pv.isValid() && pv.volume().isSensitive() ) { + + // Geometry::SensitiveDetector sd = pv.volume().sensitiveDetector(); + // Readout r = sd.readout() ; + + // VolumeID volIDElement = det.volumeID() ; + // // add the placed volumes volIDs: + // VolumeID volIDPVs = r.idSpec().encode( volIDs ) ; + + // result = r.segmentation().cellID( Position( l[0], l[1], l[2] ) , global, ( volIDElement | volIDPVs ) ); + + // // } else { + // // std::cout << " *** ERROR : found no sensitive Placement " << pv.name() + // // << " for point " << global << " try with TGeoManager ... " << std::endl ; + + // // TGeoManager *geoManager = det.volume()->GetGeoManager() ; + // // // TGeoManager *geoManager = _lcdd->world().volume()->GetGeoManager() ; + + // // PlacedVolume p = geoManager->FindNode( global.x() , global.y() , global.z() ) ; + + // // std::cout << " -> found: " << p.name() << " valid: " << p.isValid() << " sensitive: " + // // << (p.isValid() ? p.volume().isSensitive() : false ) << std::endl ; + + // } + + // return result ; + // } + namespace { @@ -232,8 +298,7 @@ namespace DD4hep { return PlacedVolume() ; } - - + Readout CellIDPositionConverter::findReadout(const Geometry::DetElement& det) const { // first check if top level is a sensitive detector