DDG4/src/Geant4VolumeManager.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     : M.Frank
0011 //
0012 //==========================================================================
0013
0014 // Framework include files
0015 #include <DD4hep/Printout.h>
0016 #include <DD4hep/Volumes.h>
0017 #include <DD4hep/DetElement.h>
0018 #include <DD4hep/DetectorTools.h>
0019 #include <DD4hep/VolumeManager.h>
0020 #include <DD4hep/detail/VolumeManagerInterna.h>
0021 #include <DDG4/Geant4VolumeManager.h>
0022 #include <DDG4/Geant4TouchableHandler.h>
0023 #include <DDG4/Geant4Mapping.h>
0024
0025 // Geant4 include files
0026 #include <G4VTouchable.hh>
0027 #include <G4LogicalVolume.hh>
0028 #include <G4VPhysicalVolume.hh>
0029
0030 // C/C++ include files
0031 #include <set>
0032 #include <stdexcept>
0033 #include <string>
0034 #include <utility>
0035 #include <vector>
0036
0037 //#define VOLMGR_HAVE_DEBUG_INFO  1
0038
0039 #ifdef VOLMGR_HAVE_DEBUG_INFO
0040 /// Namespace for the AIDA detector description toolkit
0041 namespace dd4hep {
0042   /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
0043   namespace sim {
0044
0045     class Geant4GeometryInfo::DebugInfo  {
0046     public:
0047       typedef std::vector<const G4VPhysicalVolume*>  Geant4PlacementPath;
0048       std::map<Geant4PlacementPath, Placement>       g4Paths;
0049     };
0050
0051   }    // End namespace sim
0052 }      // End namespace dd4hep
0053 #endif
0054
0055
0056 /// Namespace for the AIDA detector description toolkit
0057 namespace dd4hep {
0058   /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
0059   namespace sim {
0060     /// Print volume ID
0061     void _print_volumeid(const std::string& tag, const IDDescriptor& iddesc, DDSegmentation::VolumeID volID);
0062
0063   }    // End namespace sim
0064 }      // End namespace dd4hep
0065
0066 using namespace dd4hep::sim;
0067 using namespace dd4hep;
0068
0069 #include <DDG4/Geant4AssemblyVolume.h>
0070 using VolIDDescriptor = std::pair<VolumeID,std::vector<std::pair<const BitFieldElement*, VolumeID> > >;
0071
0072 namespace  {
0073
0074   /// Helper class to populate the Geant4 volume manager
0075   /**
0076    *  \author  M.Frank
0077    *  \version 1.0
0078    *  \ingroup DD4HEP_SIMULATION
0079    */
0080   struct Populator {
0081
0082     typedef std::vector<const TGeoNode*> Chain;
0083     // typedef std::map<VolumeID,Geant4TouchableHandler::Geant4PlacementPath> Registries;
0084     typedef std::set<VolumeID> Registries;
0085
0086     /// Reference to the Detector instance
0087     const Detector&     m_detDesc;
0088     /// Set of already added entries
0089     Registries          m_entries;
0090     /// Reference to Geant4 translation information
0091     Geant4GeometryInfo& m_geo;
0092     /// Debug flag for population
0093     long                m_debug { 0 };
0094
0095     /// Default constructor
0096     Populator(const Detector& description, Geant4GeometryInfo& g, long dbg)
0097       : m_detDesc(description), m_geo(g), m_debug(dbg)
0098     {
0099 #ifdef VOLMGR_HAVE_DEBUG_INFO
0100       if ( nullptr == g.g4DebugInfo )  {
0101         g.g4DebugInfo = new Geant4GeometryInfo::DebugInfo();
0102       }
0103 #endif
0104     }
0105     /// Default destructor
0106     ~Populator()  {
0107 #ifdef VOLMGR_HAVE_DEBUG_INFO
0108       if ( g.g4DebugInfo )  {
0109         delete g.g4DebugInfo;
0110         g.g4DebugInfo = nullptr;
0111       }
0112 #endif
0113     }
0114
0115
0116     /// Populate the Volume manager
0117     void populate(DetElement e)  {
0118       const DetElement::Children& children = e.children();
0119       m_entries.clear();
0120       for( const auto& i : children )  {
0121         DetElement   de = i.second;
0122         PlacedVolume pv = de.placement();
0123         if( pv.isValid() )  {
0124           Chain chain;
0125           SensitiveDetector sd;
0126           PlacedVolume::VolIDs ids;
0127           m_entries.clear();
0128           chain.emplace_back(m_detDesc.world().placement().ptr());
0129           scanPhysicalVolume(pv.ptr(), std::move(ids), sd, chain);
0130           continue;
0131         }
0132         printout(WARNING, "Geant4VolumeManager",
0133                  "++ Detector element %s of type %s has no placement.",
0134                  de.name(), de.type().c_str());
0135       }
0136       /// Needed to compute the cellID of parameterized volumes
0137       for( const auto& pv : m_geo.g4Placements )  {
0138         if( pv.second->IsParameterised() )
0139           m_geo.g4Parameterised[pv.second] = pv.first;
0140         if( pv.second->IsReplicated() )
0141           m_geo.g4Replicated[pv.second] = pv.first;
0142       }
0143       m_entries.clear();
0144     }
0145
0146     /// Scan a single physical volume and look for sensitive elements below
0147     void scanPhysicalVolume(const TGeoNode* node, PlacedVolume::VolIDs ids, SensitiveDetector& sd, Chain& chain) {
0148       PlacedVolume pv = node;
0149       Volume       vol = pv.volume();
0150       PlacedVolume::VolIDs pv_ids = pv.volIDs();
0151
0152       chain.emplace_back(node);
0153       ids.PlacedVolume::VolIDs::Base::insert(ids.end(), pv_ids.begin(), pv_ids.end());
0154       if( vol.isSensitive() )  {
0155         sd = vol.sensitiveDetector();
0156         if( sd.readout().isValid() )  {
0157           add_entry(sd, node, ids, chain);
0158         }
0159         else  {
0160           printout(WARNING, "Geant4VolumeManager",
0161                    "populate: Strange constellation volume %s is sensitive, but has no readout! sd:%p", pv.volume().name(),
0162                    sd.ptr());
0163         }
0164       }
0165       for( Int_t idau = 0, ndau = node->GetNdaughters(); idau < ndau; ++idau )  {
0166         TGeoNode* daughter = node->GetDaughter(idau);
0167         PlacedVolume placement(daughter);
0168         if( placement.data() ) {
0169           scanPhysicalVolume(daughter, ids, sd, chain);
0170         }
0171       }
0172       chain.pop_back();
0173     }
0174
0175     void add_entry(SensitiveDetector sd, const TGeoNode* n, const PlacedVolume::VolIDs& ids, const Chain& nodes) {
0176       Chain        control;
0177       Volume       vol;
0178       Readout      rdout           = sd.readout();
0179       IDDescriptor iddesc          = rdout.idSpec();
0180       VolumeID     code            = iddesc.encode(ids);
0181       PrintLevel   print_action    = (m_debug&Geant4VolumeManager::PRINT_ACTION) ? ALWAYS : m_geo.printLevel;
0182       PrintLevel   print_chain     = (m_debug&Geant4VolumeManager::PRINT_CHAIN)  ? ALWAYS : m_geo.printLevel;
0183       PrintLevel   print_res       = (m_debug&Geant4VolumeManager::PRINT_RESULT) ? ALWAYS : m_geo.printLevel;
0184       bool         print_nodes     = (m_debug&Geant4VolumeManager::PRINT_NODES)  ? true : false;
0185       Geant4TouchableHandler::Geant4PlacementPath path;
0186       Registries::const_iterator i = m_entries.find(code);
0187
0188       printout(print_action,"Geant4VolumeManager","+++ Add path:%s vid:%016X",
0189                detail::tools::placementPath(nodes, false).c_str(), code);
0190
0191       if( i == m_entries.end() ) {
0192         path.reserve(nodes.size());
0193         for( Chain::const_reverse_iterator k = nodes.rbegin(), kend=nodes.rend(); k != kend; ++k ) {
0194           const TGeoNode* node = *(k);
0195           auto g4pit = m_geo.g4Placements.find(node);
0196           if( g4pit != m_geo.g4Placements.end() )  {
0197             G4VPhysicalVolume* phys = g4pit->second;
0198             if( phys->IsParameterised() )  {
0199               PlacedVolume pv(n);
0200               PlacedVolumeExtension* ext = pv.data();
0201               if( nullptr == ext->params->field )  {
0202                 ext->params->field = iddesc.field(ext->volIDs.at(0).first);
0203               }
0204             }
0205             path.emplace_back(phys);
0206             printout(print_chain, "Geant4VolumeManager",
0207                      "+++     Chain: Node OK: %s [%s]", node->GetName(), phys->GetName().c_str());
0208             continue;
0209           }
0210           control.insert(control.begin(),node);
0211           vol = Volume(node->GetVolume());
0212           auto iVolImp = m_geo.g4VolumeImprints.find(vol);
0213           if ( iVolImp != m_geo.g4VolumeImprints.end() )  {
0214             for(const auto& imp : iVolImp->second )  {
0215               const auto& c = imp.first;
0216               if ( c.size() <= control.size() && control == c )  {
0217                 path.emplace_back(imp.second);
0218                 printout(print_chain, "Geant4VolumeManager", "+++     Chain: Node OK: %s %s -> %s",
0219                          node->GetName(), detail::tools::placementPath(c,false).c_str(),
0220                          imp.second->GetName().c_str());
0221                 control.clear();
0222                 break;
0223               }
0224             }
0225           }
0226         }
0227         if ( control.empty() )  {
0228           printout(print_res, "Geant4VolumeManager", "+++     Encoded Volume  IDs:%s",
0229                    detail::tools::toString(iddesc,ids,code).c_str());
0230           path.erase(path.begin()+path.size()-1);
0231           printout(print_res, "Geant4VolumeManager", "+++     Map %016X to Geant4 Path:%s",
0232                    (void*)code, Geant4TouchableHandler::placementPath(path).c_str());
0233           auto hash = detail::hash64(&path[0], path.size()*sizeof(path[0]));
0234           bool missing_hash_path = m_geo.g4Paths.find(hash) == m_geo.g4Paths.end();
0235 #ifdef VOLMGR_HAVE_DEBUG_INFO
0236           {
0237             bool missing_real_path = m_geo.g4DebugInfo->g4Paths.find(path) == m_geo.g4DebugInfo->g4Paths.end();
0238             if ( missing_real_path != missing_hash_path )   {
0239               if ( !path.empty() )
0240                 printout(ERROR,"Geant4VolumeManager"," New   G4 path: %s", Geant4TouchableHandler::placementPath(path).c_str());
0241               if ( !nodes.empty() )
0242                 printout(ERROR,"Geant4VolumeManager","     TGeo path: %s", detail::tools::placementPath(nodes,false).c_str());
0243               printout(ERROR,"Geant4VolumeManager",  " Offend.VolIDs: %s", detail::tools::toString(iddesc,ids,code).c_str());
0244             }
0245             if ( missing_real_path ) {
0246               Geant4GeometryInfo::PlacementFlags opt;
0247               opt.flags.parametrised = path.front()->IsParameterised() ? 1 : 0;
0248               opt.flags.replicated   = path.front()->IsReplicated()    ? 1 : 0;
0249               m_geo.g4DebugInfo->g4Paths[path] = { code, opt.value };
0250             }
0251           }
0252 #endif
0253           if ( missing_hash_path ) {
0254             Geant4GeometryInfo::PlacementFlags opt;
0255             opt.flags.parametrised = path.front()->IsParameterised() ? 1 : 0;
0256             opt.flags.replicated   = path.front()->IsReplicated()    ? 1 : 0;
0257             m_geo.g4Paths[hash]    = { code, opt.value };
0258             if( m_debug&Geant4VolumeManager::PRINT_VOLIDS )  {
0259               std::string idstr = iddesc.str(code);
0260               printout(ALWAYS, "Geant4VolumeManager",
0261                        "+++     Decoded Volume  IDs: %016llX -> %s", code, idstr.c_str());
0262             }
0263             m_entries.emplace(code);
0264             return;
0265           }
0266           /// This is a normal case for parametrized volumes and no error
0267           if ( !path.empty() && (path.front()->IsParameterised() || path.front()->IsReplicated()) )  {
0268             return;
0269           }
0270           printout(ERROR, "Geant4VolumeManager", "populate: Severe error: Duplicated Geant4 path!!!! %s %s",
0271                    " [THIS SHOULD NEVER HAPPEN]", Geant4TouchableHandler::placementPath(path).c_str());
0272           goto Err;
0273         }
0274         printout(INFO, "Geant4VolumeManager", "Control block has still %d entries:%s",
0275                  int(control.size()), detail::tools::placementPath(control,true).c_str());
0276         print_nodes = true;
0277         goto Err;
0278       }
0279       else  {
0280         /// This is a normal case for parametrized volumes and no error
0281         if ( !path.empty() && (path.front()->IsParameterised() || path.front()->IsReplicated()) )  {
0282           return;
0283         }
0284       }
0285       printout(ERROR, "Geant4VolumeManager", "populate: Severe error: Duplicated Volume entry: 0x%X"
0286                " [THIS SHOULD NEVER HAPPEN]", code);
0287
0288     Err:
0289       if ( i != m_entries.end() )
0290         printout( ERROR,"Geant4VolumeManager"," Known G4 path: %s", Geant4TouchableHandler::placementPath(path).c_str() );
0291       if ( !path.empty() )
0292         printout( ERROR,"Geant4VolumeManager"," New   G4 path: %s", Geant4TouchableHandler::placementPath(path).c_str() );
0293       if ( !nodes.empty() )  {
0294         printout( ERROR,"Geant4VolumeManager","     TGeo path: %s", detail::tools::placementPath(nodes,false).c_str() );
0295         if( print_nodes )  {
0296           std::string node_path;
0297           for( std::size_t in=0; in<nodes.size(); ++in )  {
0298             PlacedVolume pv(nodes[in]);
0299             node_path += "/";
0300             node_path += pv.name();
0301             printout( ALWAYS,"Geant4VolumeManager", "     TGeo Node[%ld]:  %s [%p]  Volids: '%s'",
0302                       in, node_path.c_str(), (void*)nodes[in], pv.volIDs().str().c_str() );
0303           }
0304         }
0305       }
0306       printout( ERROR,"Geant4VolumeManager",  " Offend.VolIDs: %s", detail::tools::toString(iddesc,ids,code).c_str() );
0307       throw std::runtime_error("Failed to populate Geant4 volume manager!");
0308     }
0309   };
0310 }
0311
0312 /// Initializing constructor. The tree will automatically be built if possible
0313 Geant4VolumeManager::Geant4VolumeManager(const Detector& description, Geant4GeometryInfo* info, long debug)
0314   : Handle<Geant4GeometryInfo>(info)  {
0315   if( info && info->valid )  {
0316     if( !info->has_volmgr )  {
0317       Populator p(description, *info, debug);
0318       printout( ALWAYS, "Geant4VolumeManager", "+++ Populating Geant4 volume manager.");
0319       p.populate(description.world());
0320       printout( ALWAYS, "Geant4VolumeManager",
0321                 "+++ Geant4 volume manager populated with %ld sensitive path entries.",
0322                 info->g4Paths.size() );
0323       if( debug&PRINT_ENTRIES )  {
0324         int count = 0;
0325         VolumeManager volmgr = description.volumeManager();
0326         for( auto it=info->g4Paths.begin(); it != info->g4Paths.end(); ++it, ++count )  {
0327           VolumeID volid = it->second.volumeID;
0328           VolumeManagerContext* context = volmgr.lookupContext(volid);
0329           if( context )  {
0330             std::string  path = context->element.path();
0331             PlacedVolume plac = context->volumePlacement();
0332             SensitiveDetector sens = plac.volume().sensitiveDetector();
0333             std::string       idstr = sens.idSpec().str(volid);
0334             printout( ALWAYS, "Geant4VolumeManager", "%8d:  %016X %s -> %s",
0335                       count, volid, path.c_str(), plac.name());
0336             printout(ALWAYS, "Geant4VolumeManager", "%8s   %16s %s", "", "", idstr.c_str());
0337           }
0338           else  {
0339             printout( ERROR, "Geant4VolumeManager",
0340                       "Missing volume manager entry: volume ID %016X", volid);
0341           }
0342         }
0343       }
0344       info->has_volmgr = true;
0345     }
0346     return;
0347   }
0348   except("Geant4VolumeManager", "Attempt populate from invalid Geant4 geometry info [Invalid-Info]");
0349 }
0350
0351 /// Helper: Generate placement path from touchable object
0352 std::vector<const G4VPhysicalVolume*>
0353 Geant4VolumeManager::placementPath(const G4VTouchable* touchable, bool exception) const  {
0354   Geant4TouchableHandler handler(touchable);
0355   return handler.placementPath(exception);
0356 }
0357
0358 /// Check the validity of the information before accessing it.
0359 bool Geant4VolumeManager::checkValidity() const  {
0360   if( !isValid() )  {
0361     except("Geant4VolumeManager", "Attempt to use invalid Geant4 volume manager [Invalid-Handle]");
0362   }
0363   else if( !ptr()->valid )  {
0364     except("Geant4VolumeManager", "Attempt to use invalid Geant4 geometry info [Invalid-Info]");
0365   }
0366   return true;
0367 }
0368
0369 namespace  {
0370   std::string debug_status(const Geant4VolumeManager* mgr)  {
0371     char text[256];
0372     auto* p = mgr->ptr();
0373     if ( p )  {
0374       ::snprintf(text, sizeof(text), "==> #path entries: %ld valid: %s has_volmgr: %s",
0375                  p->g4Paths.size(), yes_no(p->valid), yes_no(p->has_volmgr));
0376       return { text };
0377     }
0378     return { "Invalid handle to Geant4GeometryInfo" };
0379   }
0380 }
0381
0382 /// Access CELLID by Geant4 touchable object
0383 VolumeID Geant4VolumeManager::volumeID(const G4VTouchable* touchable) const  {
0384   Geant4TouchableHandler handler(touchable);
0385   std::vector<const G4VPhysicalVolume*> path = handler.placementPath();
0386   if( !isValid() )  {
0387     printout(INFO, "Geant4VolumeManager", "+++   INVALID Geant4VolumeManager handle.");
0388     return NonExisting;
0389   }
0390   else if( !ptr()->valid )  {
0391     printout(INFO, "Geant4VolumeManager", "+++   INVALID Geant4VolumeManager [Not initialized]");
0392     return NonExisting;
0393   }
0394   else if( path.empty() )  {
0395     printout(INFO, "Geant4VolumeManager", "+++   EMPTY volume Geant4 Path: %s",
0396              Geant4TouchableHandler::placementPath(path).c_str());
0397     return NonExisting;
0398   }
0399   else  {
0400     uint64_t hash = detail::hash64(&path[0], sizeof(path[0])*path.size());
0401     auto i = ptr()->g4Paths.find(hash);
0402     if( i != ptr()->g4Paths.end() )  {
0403       const auto& e = (*i).second;
0404       VolumeID volid = e.volumeID;
0405       /// No parametrization or replication.
0406       if( e.flags == 0 )  {
0407         return volid;
0408       }
0409       const auto& paramterised = ptr()->g4Parameterised;
0410       const auto& replicated   = ptr()->g4Replicated;
0411       /// This is incredibly slow .... but what can I do ? Need a better idea.
0412       for( std::size_t j=0; j < path.size(); ++j )  {
0413         const auto* phys = path[j];
0414         if( phys->IsParameterised() )  {
0415           int copy_no = touchable->GetCopyNumber(j);
0416           const auto it = paramterised.find(phys);
0417           if( it != paramterised.end() )  {
0418             //printout(INFO,"Geant4VolumeManager",
0419             //         "Copy number:   %ld  <--> %ld", copy_no, long(phys->GetCopyNo()));
0420             const auto* field = (*it).second.data()->params->field;
0421             volid |= IDDescriptor::encode(field, copy_no);
0422             continue;
0423           }
0424           except("Geant4VolumeManager",
0425                  "Error  Geant4VolumeManager::volumeID(const G4VTouchable* touchable)");
0426         }
0427         else if( phys->IsReplicated() )   {
0428           int copy_no = touchable->GetCopyNumber(j);
0429           const auto it = replicated.find(phys);
0430           if( it != replicated.end() )  {
0431             const auto* field = (*it).second.data()->params->field;
0432             volid |= IDDescriptor::encode(field, copy_no);
0433             continue;
0434           }
0435           except("Geant4VolumeManager",
0436                  "Error  Geant4VolumeManager::volumeID(const G4VTouchable* touchable)");
0437         }
0438       }
0439       return volid;
0440     }
0441     if( !path[0] )  {
0442       printout(INFO, "Geant4VolumeManager", "+++   Bad Geant4 volume path: \'%s\' [invalid path] %s",
0443                Geant4TouchableHandler::placementPath(path).c_str(), debug_status(this).c_str());
0444       return InvalidPath;
0445     }
0446     else if( !path[0]->GetLogicalVolume()->GetSensitiveDetector() )  {
0447       printout(DEBUG, "Geant4VolumeManager", "+++   Bad Geant4 volume path: \'%s\' [insensitive] %s",
0448                Geant4TouchableHandler::placementPath(path).c_str(), debug_status(this).c_str());
0449       return Insensitive;
0450     }
0451     printout(INFO, "Geant4VolumeManager",
0452              "+++   Bad Geant4 volume path: \'%s\' [missing entry] %s",
0453              Geant4TouchableHandler::placementPath(path).c_str(), debug_status(this).c_str());
0454     return NonExisting;
0455   }
0456   printout(INFO, "Geant4VolumeManager", "+++   Bad Geant4 volume path: \'%s\' %s",
0457            Geant4TouchableHandler::placementPath(path).c_str(), yes_no(path.empty()));
0458   return NonExisting;
0459 }
0460
0461 /// Access fully decoded volume fields  by placement path
0462 void Geant4VolumeManager::volumeDescriptor(const std::vector<const G4VPhysicalVolume*>& path,
0463                                            VolIDDescriptor& vol_desc) const
0464 {
0465   vol_desc.second.clear();
0466   vol_desc.first = NonExisting;
0467   if( !path.empty() && checkValidity() )  {
0468     auto hash = detail::hash64(&path[0], sizeof(path[0])*path.size());
0469     auto i = ptr()->g4Paths.find(hash);
0470     if( i != ptr()->g4Paths.end() )  {
0471       VolumeID vid = (*i).second.volumeID;
0472       G4LogicalVolume* lvol = path[0]->GetLogicalVolume();
0473       if( lvol->GetSensitiveDetector() ) {
0474         const auto* node = path[0];
0475         const auto& pm = ptr()->g4Placements;
0476         for( const auto& ipm : pm )  {
0477           if ( ipm.second == node )  {
0478             PlacedVolume      pv  = ipm.first;
0479             SensitiveDetector sd  = pv.volume().sensitiveDetector();
0480             IDDescriptor      dsc = sd.readout().idSpec();
0481             vol_desc.first = vid;
0482             dsc.decodeFields(vid, vol_desc.second);
0483             return;
0484           }
0485         }
0486       }
0487       vol_desc.first = Insensitive;
0488       return;
0489     }
0490     if( !path[0] )
0491       vol_desc.first = InvalidPath;
0492     else if( !path[0]->GetLogicalVolume()->GetSensitiveDetector() )
0493       vol_desc.first = Insensitive;
0494     else
0495       vol_desc.first = NonExisting;
0496   }
0497 }
0498
0499 /// Access fully decoded volume fields by Geant4 touchable object
0500 void Geant4VolumeManager::volumeDescriptor(const G4VTouchable* touchable,
0501                                            VolIDDescriptor&    vol_desc)  const  {
0502   volumeDescriptor(placementPath(touchable), vol_desc);
0503 }