DDRec/src/CellIDPositionConverter.cpp

0001 //==========================================================================
0002 //  AIDA Detector description implementation
0003 //--------------------------------------------------------------------------
0004 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0005 // All rights reserved.
0006 //
0007 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0008 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0009 //
0010 // Author     : F.Gaede
0011 //
0012 //==========================================================================
0013
0014 #include <DDRec/CellIDPositionConverter.h>
0015
0016 #include <DD4hep/Detector.h>
0017 #include <DD4hep/detail/VolumeManagerInterna.h>
0018
0019 #include <TGeoManager.h>
0020
0021 namespace dd4hep {
0022   namespace rec {
0023
0024     using std::set;
0025
0026     const VolumeManagerContext*
0027     CellIDPositionConverter::findContext(const CellID& cellID) const {
0028       return _volumeManager.lookupContext( cellID ) ;
0029     }
0030
0031
0032     Position CellIDPositionConverter::position(const CellID& cell) const {
0033
0034       // untill we have the alignment map object, we return the nominal position
0035
0036       return positionNominal( cell ) ;
0037     }
0038
0039     Position CellIDPositionConverter::positionNominal(const CellID& cell) const {
0040
0041       double l[3], e[3], g[3];
0042
0043       const VolumeManagerContext* context = findContext( cell ) ;
0044
0045       if( context == NULL)
0046     return Position() ;
0047
0048       DetElement det = context->element ;
0049
0050
0051 #if 0  // this uses the deprecated VolumeManagerContext::placement
0052
0053       PlacedVolume pv = context->placement ;
0054
0055       if( ! pv.volume().isSensitive() )
0056     return Position() ;
0057
0058       SensitiveDetector sd = pv.volume().sensitiveDetector();
0059       Readout r = sd.readout() ;
0060
0061 #else
0062
0063       // use a recursive search for the Readout
0064
0065       Readout r = findReadout( det ) ;
0066
0067 #endif
0068
0069       Segmentation seg = r.segmentation() ;
0070       Position local = seg.position(cell);
0071
0072       local.GetCoordinates(l);
0073
0074       const TGeoMatrix& volToElement = context->toElement();
0075       volToElement.LocalToMaster(l, e);
0076
0077       const TGeoMatrix& elementToGlobal = det.nominal().worldTransformation();
0078       elementToGlobal.LocalToMaster(e, g);
0079
0080
0081       return Position(g[0], g[1], g[2]);
0082     }
0083
0084
0085
0086
0087     CellID CellIDPositionConverter::cellID(const Position& global) const {
0088
0089       CellID result(0) ;
0090
0091       TGeoManager *geoManager = _description->world().volume()->GetGeoManager() ;
0092
0093       PlacedVolume pv = geoManager->FindNode( global.x() , global.y() , global.z() ) ;
0094
0095       if(  pv.isValid() && pv.volume().isSensitive() ) {
0096
0097     TGeoHMatrix*  m = geoManager->GetCurrentMatrix() ;
0098
0099     double g[3], l[3] ;
0100     global.GetCoordinates( g ) ;
0101     m->MasterToLocal( g, l );
0102
0103     SensitiveDetector sd = pv.volume().sensitiveDetector();
0104     Readout r = sd.readout() ;
0105
0106     // collect all volIDs for the current path
0107     PlacedVolume::VolIDs volIDs ;
0108     volIDs.insert( std::end(volIDs), std::begin(pv.volIDs()), std::end(pv.volIDs())) ;
0109
0110     TGeoPhysicalNode pN( geoManager->GetPath() ) ;
0111
0112     unsigned motherCount = 0 ;
0113
0114     while( pN.GetMother( motherCount ) != NULL   ){
0115
0116         PlacedVolume mPv = pN.GetMother( motherCount++ ) ;
0117
0118         if( mPv.isValid() &&  pN.GetMother( motherCount ) != NULL )  // world has no volIDs
0119           volIDs.insert( std::end(volIDs), std::begin(mPv.volIDs()), std::end(mPv.volIDs())) ;
0120     }
0121
0122     VolumeID volIDPVs = r.idSpec().encode( volIDs ) ;
0123
0124     result = r.segmentation().cellID( Position( l[0], l[1], l[2] ) , global, volIDPVs  );
0125       }
0126
0127       return result ;
0128     }
0129
0130     // CellID CellIDPositionConverter::cellID(const Position& global) const {
0131
0132     //   CellID result(0) ;
0133
0134     //   DetElement motherDet = _description->world()  ; // could also start from an arbitrary DetElement here !?
0135
0136     //   DetElement det = findDetElement( global , motherDet ) ;
0137
0138     //   if( ! det.isValid() )
0139     //  return result ;
0140
0141     //   double g[3], e[3] , l[3] ;
0142     //   global.GetCoordinates( g ) ;
0143     //   det.nominal().worldTransformation().MasterToLocal( g, e );
0144
0145     //   PlacedVolume::VolIDs volIDs ;
0146
0147     //   PlacedVolume pv = findPlacement( Position( e[0], e[1] , e[2] ) , det.placement() , l , volIDs ) ;
0148
0149     //   TGeoManager *geoManager = det.volume()->GetGeoManager() ;
0150     //   // TGeoManager *geoManager = _description->world().volume()->GetGeoManager() ;
0151
0152     //   PlacedVolume pv1 = geoManager->FindNode( global.x() , global.y() , global.z() ) ;
0153
0154     //   // std::cout << " -> TGM : " << pv1.name() << "  valid: " << pv1.isValid() << " sensitive: "
0155     //   //         <<  (pv1.isValid() ? pv1.volume().isSensitive() : false ) << std::endl ;
0156     //   // std::cout << " -> FG :  " << pv.name() << "  valid: " << pv.isValid() << " sensitive: "
0157     //   //         <<  (pv.isValid() ? pv.volume().isSensitive() : false ) << std::endl ;
0158
0159
0160
0161     //  if(  pv.isValid() && pv.volume().isSensitive() ) {
0162
0163     //  SensitiveDetector sd = pv.volume().sensitiveDetector();
0164     //  Readout r = sd.readout() ;
0165
0166     //  VolumeID volIDElement = det.volumeID() ;
0167     //  // add the placed volumes volIDs:
0168     //  VolumeID volIDPVs = r.idSpec().encode( volIDs ) ;
0169
0170     //  result = r.segmentation().cellID( Position( l[0], l[1], l[2] ) , global, ( volIDElement | volIDPVs  ) );
0171
0172     //   // } else {
0173     //   //     std::cout << " *** ERROR : found no sensitive Placement " << pv.name()
0174     //   //           << "  for point " << global << " try with TGeoManager ... " << std::endl ;
0175
0176     //   //     TGeoManager *geoManager = det.volume()->GetGeoManager() ;
0177     //   //     // TGeoManager *geoManager = _description->world().volume()->GetGeoManager() ;
0178
0179     //   //     PlacedVolume p = geoManager->FindNode( global.x() , global.y() , global.z() ) ;
0180
0181     //   //     std::cout << " -> found: " << p.name() << "  valid: " << p.isValid() << " sensitive: "
0182     //   //           <<  (p.isValid() ? p.volume().isSensitive() : false ) << std::endl ;
0183
0184     //   }
0185
0186     //   return result ;
0187     // }
0188
0189
0190     namespace {
0191
0192       bool containsPoint( const DetElement& det, const Position& global ) {
0193
0194     if( det.volume().isValid() and det.volume().solid().isValid() ) {
0195
0196       double g[3], l[3] ;
0197       global.GetCoordinates( g ) ;
0198
0199       det.nominal().worldTransformation().MasterToLocal( g, l );
0200
0201       return det.volume().solid()->Contains( l ) ;
0202     }
0203
0204     return false ;
0205       }
0206
0207     }
0208
0209     DetElement CellIDPositionConverter::findDetElement(const Position& global,
0210                                const DetElement& d) const {
0211
0212       DetElement det = ( d.isValid() ? d : _description->world() ) ;
0213
0214       //      std::cout << " --- " << global << det.name() << std::endl ;
0215
0216       if( containsPoint( det, global ) ) {
0217
0218     if( det.children().empty() )  // no children -> we are done
0219       return det ;
0220
0221     // see if we have a child DetElement that contains the point ...
0222
0223     DetElement result ;
0224     for( const auto& it : det.children() ){
0225       // std::cout << " -   " << global << it.second.name()  << " " << containsPoint( it.second , global )
0226       //        << " nChild: " << it.second.children().size() << " isValid: " <<  it.second.isValid()
0227       //        << std::endl ;
0228
0229       if( containsPoint( it.second , global ) ){
0230         result = it.second ;
0231         break ;
0232       }
0233     }
0234     if( result.isValid() ){  // ... yes, we have
0235
0236       if( result.children().empty() )  // no more children -> done
0237         return result ;
0238       else
0239         return findDetElement( global, result ) ; // keep searching in children
0240     }
0241
0242       }
0243
0244       // point not contained
0245       return DetElement() ;
0246     }
0247
0248     PlacedVolume CellIDPositionConverter::findPlacement(const Position& pos, const  PlacedVolume& pv , double locPos[3], PlacedVolume::VolIDs& volIDs) const {
0249
0250
0251       double l[3] ;
0252       pos.GetCoordinates( l ) ;
0253
0254       //      std::cout << " --- " << pos << " " << pv.name() << " loc: (" << locPos[0] << "," << locPos[1] << "," << locPos[2] << ")" << std::endl ;
0255
0256       if( pv.volume().solid()->Contains( l ) ) {
0257
0258     // copy the volIDs
0259     volIDs.insert( std::end(volIDs), std::begin(pv.volIDs()), std::end(pv.volIDs()));
0260
0261     int ndau = pv->GetNdaughters() ;
0262
0263     if( ndau == 0 ) // no daughter volumes -> we are done
0264       return pv ;
0265
0266     // see if we have a daughter volume that contains the point ...
0267
0268     PlacedVolume result ;
0269     for (int i = 0 ; i < ndau; ++i) {
0270
0271       PlacedVolume pvDau = pv->GetDaughter( i );
0272       pvDau->MasterToLocal( l , locPos ) ;  // transform point to daughter's local frame
0273
0274       // std::cout << "   - " << pos << " " << pvDau.name()
0275       //        << " loc: (" << locPos[0] << "," << locPos[1] << "," << locPos[2] << ")"
0276       //        <<  pvDau.volume().solid()->Contains( locPos )
0277       //        << " ndau: " << pvDau->GetNdaughters()
0278       //        << std::endl ;
0279
0280
0281       if( pvDau.volume().solid()->Contains( locPos ) ) { // point is contained in daughter node
0282         result = pvDau ;
0283         volIDs.insert( std::end(volIDs), std::begin(pvDau.volIDs()), std::end(pvDau.volIDs()) );
0284         break ;
0285       }
0286     }
0287
0288     if( result.isValid() ){  // ... yes, we have
0289
0290       if( result->GetNdaughters() == 0 ){  // no more children -> done
0291         return result ;
0292       } else
0293         return findPlacement( Position( locPos[0], locPos[1] , locPos[2]  ), result , locPos , volIDs) ; // keep searching in daughter volumes
0294     }
0295
0296       }
0297
0298       return  PlacedVolume() ;
0299     }
0300
0301
0302     Readout CellIDPositionConverter::findReadout(const DetElement& det) const {
0303
0304       // first check if top level is a sensitive detector
0305       if (det.volume().isValid() and det.volume().isSensitive()) {
0306     SensitiveDetector sd = det.volume().sensitiveDetector();
0307     if (sd.isValid() and sd.readout().isValid()) {
0308       return sd.readout();
0309     }
0310       }
0311       // if not, return the first sensitive daughter volume's readout
0312
0313       Readout r = findReadout( det.placement() ) ;
0314       if( r.isValid() )
0315     return r ;
0316
0317       // nothing found !?
0318       return Readout();
0319     }
0320
0321     Readout CellIDPositionConverter::findReadout(const PlacedVolume& pv) const {
0322
0323       // first check if we are in a sensitive volume
0324       if( pv.volume().isSensitive() ){
0325     SensitiveDetector sd = pv.volume().sensitiveDetector();
0326         if (sd.isValid() and sd.readout().isValid()) {
0327       return sd.readout();
0328     }
0329       }
0330
0331       for (Int_t idau = 0, ndau = pv->GetNdaughters(); idau < ndau; ++idau) {
0332
0333     PlacedVolume dpv = pv->GetDaughter(idau);
0334     Readout r = findReadout( dpv ) ;
0335     if( r.isValid() )
0336       return r ;
0337       }
0338
0339       return Readout() ;
0340     }
0341
0342     std::vector<double> CellIDPositionConverter::cellDimensions(const CellID& cell) const {
0343       auto context = findContext( cell ) ;
0344       if( context == nullptr ) return { };
0345       dd4hep::Readout r  = findReadout( context->element ) ;
0346       dd4hep::Segmentation seg = r.segmentation() ;
0347       return seg.cellDimensions( cell );
0348     }
0349
0350
0351
0352
0353   } /* namespace rec */
0354 } /* namespace dd4hep */