diff --git a/DDRec/include/DDRec/CellIDPositionConverter.h b/DDRec/include/DDRec/CellIDPositionConverter.h
index 92365505e672aa6a68e69b0173a5104b2e931cbd..2f27e27382778970e3dc6a249e78d307e5d5deab 100644
--- a/DDRec/include/DDRec/CellIDPositionConverter.h
+++ b/DDRec/include/DDRec/CellIDPositionConverter.h
@@ -33,7 +33,7 @@ namespace DD4hep {
public:
/// The constructor - takes the main lcdd object.
- CellIDPositionConverter( Geometry::LCDD& lcdd ) {
+ CellIDPositionConverter( Geometry::LCDD& lcdd ) : _lcdd( &lcdd ) {
_volumeManager = Geometry::VolumeManager::getVolumeManager(lcdd);
}
@@ -61,6 +61,21 @@ namespace DD4hep {
const Geometry::VolumeManagerContext* findContext(const CellID& cellID) const;
+
+ /** Find the DetElement that contains the given point - if no DetElement is found, an
+ * invalid DetElement is returned. Uses the optionally given DetElement as start for the search.
+ */
+ Geometry::DetElement findDetElement(const Geometry::Position& global,
+ const Geometry::DetElement& det=Geometry::DetElement() ) const;
+
+
+ /** Find the lowest daughter Placement in the Placement that
+ * contains the point (in the coordinate system of the mother placement) the local coordintates
+ * in this placement are also returned.
+ */
+ Geometry::PlacedVolume findPlacement(const Geometry::Position& point, const Geometry::PlacedVolume& mother, double locPos[3]) const ;
+
+
/** Find the readout object for the given DetElement. If the DetElement is sensitive the corresondig
* Readout is returned, else a recursive search in the daughter volumes (nodes) of this DetElement's
* volume is performed and the first Readout object is returned.
@@ -77,6 +92,7 @@ namespace DD4hep {
protected:
Geometry::VolumeManager _volumeManager{} ;
+ const Geometry::LCDD* _lcdd ;
};
diff --git a/DDRec/src/CellIDPositionConverter.cpp b/DDRec/src/CellIDPositionConverter.cpp
index 21d3a324e3eee58dd7512e4e6e2d5114c0e77b34..4951548d7352bbc4a31490887e43905af4b45642 100644
--- a/DDRec/src/CellIDPositionConverter.cpp
+++ b/DDRec/src/CellIDPositionConverter.cpp
@@ -46,7 +46,7 @@ namespace DD4hep {
DetElement det = context->element ;
-#if 0 // this uses the deprected VolumeManagerContext::placement
+#if 0 // this uses the deprecated VolumeManagerContext::placement
PlacedVolume pv = context->placement ;
@@ -85,11 +85,148 @@ namespace DD4hep {
CellID CellIDPositionConverter::cellID(const Position& global) const {
- //FIXME: how to best implement this ?
+ 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 pv = findPlacement( Position( e[0], e[1] , e[2] ) , det.placement() , l ) ;
+
+ if( pv.isValid() && pv.volume().isSensitive() ) {
+
+ Geometry::SensitiveDetector sd = pv.volume().sensitiveDetector();
+ Readout r = sd.readout() ;
+ VolumeID volID = det.volumeID() ;
+
+ result = r.segmentation().cellID( Position( l[0], l[1], l[2] ) , global, volID );
+
+ } else {
+
+ std::cout << " *** ERROR : found non-sensitive Placement " << pv.name()
+ << " for point " << global << std::endl ;
+ }
+
+
+ return result ;
+ }
- return 0 ;
+
+ namespace {
+
+ bool containsPoint( const DetElement& det, const Geometry::Position& global ) {
+
+ if( det.volume().isValid() and det.volume().solid().isValid() ) {
+
+ double g[3], l[3] ;
+ global.GetCoordinates( g ) ;
+
+ det.nominal().worldTransformation().MasterToLocal( g, l );
+
+ return det.volume().solid()->Contains( l ) ;
+ }
+
+ return false ;
+ }
+
}
+ DetElement CellIDPositionConverter::findDetElement(const Geometry::Position& global,
+ const DetElement& d) const {
+
+ DetElement det = ( d.isValid() ? d : _lcdd->world() ) ;
+
+ // std::cout << " --- " << global << det.name() << std::endl ;
+
+ if( containsPoint( det, global ) ) {
+
+ if( det.children().empty() ) // no children -> we are done
+ return det ;
+
+ // see if we have a child DetElement that contains the point ...
+
+ DetElement result ;
+ for( auto it : det.children() ){
+ // std::cout << " - " << global << it.second.name() << " " << containsPoint( it.second , global )
+ // << " nChild: " << it.second.children().size() << " isValid: " << it.second.isValid()
+ // << std::endl ;
+
+ if( containsPoint( it.second , global ) ){
+ result = it.second ;
+ break ;
+ }
+ }
+ if( result.isValid() ){ // ... yes, we have
+
+ if( result.children().empty() ) // no more children -> done
+ return result ;
+ else
+ return findDetElement( global, result ) ; // keep searching in children
+ }
+
+ }
+
+ // point not contained
+ return DetElement() ;
+ }
+
+ Geometry::PlacedVolume CellIDPositionConverter::findPlacement(const Geometry::Position& pos, const Geometry::PlacedVolume& pv , double locPos[3]) const {
+
+
+ double l[3] ;
+ pos.GetCoordinates( l ) ;
+
+ std::cout << " --- " << pos << " " << pv.name() << " loc: (" << locPos[0] << "," << locPos[1] << "," << locPos[2] << ")" << std::endl ;
+
+ if( pv.volume().solid()->Contains( l ) ) {
+
+ int ndau = pv->GetNdaughters() ;
+
+ if( ndau == 0 ) // no daughter volumes -> we are done
+ return pv ;
+
+ // see if we have a daughter volume that contains the point ...
+
+ PlacedVolume result ;
+ for (int i = 0 ; i < ndau; ++i) {
+
+ PlacedVolume pvDau = pv->GetDaughter( i );
+ pvDau->MasterToLocal( l , locPos ) ; // transform point to daughter's local frame
+
+ std::cout << " - " << pos << " " << pvDau.name()
+ << " loc: (" << locPos[0] << "," << locPos[1] << "," << locPos[2] << ")"
+ << pvDau.volume().solid()->Contains( locPos )
+ << " ndau: " << pvDau->GetNdaughters()
+ << std::endl ;
+
+
+ if( pvDau.volume().solid()->Contains( locPos ) ) { // point is contained in daughter node
+ result = pvDau ;
+ break ;
+ }
+ }
+
+ if( result.isValid() ){ // ... yes, we have
+
+ if( result->GetNdaughters() == 0 ){ // no more children -> done
+ return result ;
+ } else
+ return findPlacement( Position( locPos[0], locPos[1] , locPos[2] ), result , locPos ) ; // keep searching in daughter volumes
+ }
+
+ }
+
+ return PlacedVolume() ;
+ }
+
+
Readout CellIDPositionConverter::findReadout(const Geometry::DetElement& det) const {
@@ -120,13 +257,9 @@ namespace DD4hep {
}
}
- // traverse the daughter nodes:
- const TGeoNode* node = pv.ptr();
-
- for (Int_t idau = 0, ndau = node->GetNdaughters(); idau < ndau; ++idau) {
+ for (Int_t idau = 0, ndau = pv->GetNdaughters(); idau < ndau; ++idau) {
- TGeoNode* daughter = node->GetDaughter(idau);
- PlacedVolume dpv( daughter );
+ PlacedVolume dpv = pv->GetDaughter(idau);
Readout r = findReadout( dpv ) ;
if( r.isValid() )
return r ;
diff --git a/UtilityApps/src/test_cellid_position_converter.cpp b/UtilityApps/src/test_cellid_position_converter.cpp
index 2ec12ded2fb8123975f798123f5e95d48a1a278d..0f54788559f26cef5c83a5a42ec66311f6d47fbc 100644
--- a/UtilityApps/src/test_cellid_position_converter.cpp
+++ b/UtilityApps/src/test_cellid_position_converter.cpp
@@ -40,6 +40,8 @@ static DDTest test( "cellid_position_converter" ) ;
//=============================================================================
const double epsilon = dd4hep::micrometer ;
+const int maxHit = 10 ;
+
double dist( const Position& p0, const Position& p1 ){
Position p2 = p1 - p0 ;
@@ -77,7 +79,7 @@ int main_wrapper(int argc, char** argv ){
// use only hits from these collections
std::set< std::string > subset = {} ;
- //{"BeamCalCollection" } ; //EcalBarrelCollection" } ; //"HcalBarrelRegCollection"} ;
+ //{"BeamCalCollection" } ; //{"EcalBarrelCollection" } ; //{"HcalBarrelRegCollection"} ;
// ignore all hits from these collections
@@ -107,9 +109,11 @@ int main_wrapper(int argc, char** argv ){
std::string cellIDEcoding = col->getParameters().getStringVal("CellIDEncoding") ;
- DD4hep::BitField64 idDecoder( cellIDEcoding ) ;
+ DD4hep::BitField64 idDecoder0( cellIDEcoding ) ;
+ DD4hep::BitField64 idDecoder1( cellIDEcoding ) ;
- int nHit = col->getNumberOfElements() ;
+ int nHit = std::min( col->getNumberOfElements(), maxHit ) ;
+
for(int i=0 ; i< nHit ; ++i){
@@ -119,25 +123,34 @@ int main_wrapper(int argc, char** argv ){
DD4hep::long64 id1 = sHit->getCellID1() ;
DD4hep::long64 id = ( id1 << 32 | id0 ) ;
- idDecoder.setValue( id ) ;
+ idDecoder0.setValue( id ) ;
Position point( sHit->getPosition()[0]* dd4hep::mm , sHit->getPosition()[1]* dd4hep::mm , sHit->getPosition()[2]* dd4hep::mm ) ;
// ====== test cellID to position and positio to cellID conversion ================================
+ DetElement det = idposConv.findDetElement( point ) ;
- // CellID idFromDecoder = idposConv.cellID( point ) ;
+ CellID idFromDecoder = idposConv.cellID( point ) ;
- Position pointFromDecoder = idposConv.position( id ) ;
+ idDecoder1.setValue( idFromDecoder ) ;
+
+
+ std::stringstream sst ;
+ sst << " compare ids: " << det.name() << " " << idDecoder0.valueString() << " - " << idDecoder1.valueString() ;
- // test( idFromDecoder , id , " compare ids: " ) ;
+ test( id, idFromDecoder, sst.str() ) ;
+
+
+ Position pointFromDecoder = idposConv.position( id ) ;
double d = dist(pointFromDecoder, point) ;
- std::stringstream sst ;
- sst << " dist " << d << " ( " << pointFromDecoder << " ) - ( " << point << " )" ;
+ std::stringstream sst1 ;
+ sst1 << " dist " << d << " ( " << point << " ) - ( " << pointFromDecoder << " ) - detElement: "
+ << det.name() ;
- test( d < epsilon , true , sst.str() ) ;
+ test( d < epsilon , true , sst1.str() ) ;
}
}