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