DDG4/src/Geant4SensDetAction.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/Primitives.h>
0017 #include <DD4hep/InstanceCount.h>
0018
0019 #include <DDG4/Geant4Kernel.h>
0020 #include <DDG4/Geant4Mapping.h>
0021 #include <DDG4/Geant4StepHandler.h>
0022 #include <DDG4/Geant4SensDetAction.h>
0023 #include <DDG4/Geant4VolumeManager.h>
0024 #include <DDG4/Geant4MonteCarloTruth.h>
0025
0026 // Geant4 include files
0027 #include <G4Step.hh>
0028 #include <G4SDManager.hh>
0029 #include <G4VSensitiveDetector.hh>
0030
0031 // C/C++ include files
0032 #include <stdexcept>
0033
0034 #ifdef DD4HEP_USE_GEANT4_UNITS
0035 #define MM_2_CM 1.0
0036 #else
0037 #define MM_2_CM 0.1
0038 #endif
0039
0040 using namespace dd4hep::sim;
0041
0042 #if 0
0043 namespace {
0044   Geant4ActionSD* _getSensitiveDetector(const std::string& name) {
0045     G4SDManager* mgr = G4SDManager::GetSDMpointer();
0046     G4VSensitiveDetector* sd = mgr->FindSensitiveDetector(name);
0047     if (0 == sd) {
0048       dd4hep::except("Geant4Sensitive", "DDG4: You requested to configure actions "
0049                      "for the sensitive detector %s,\nDDG4: which is not known to Geant4. "
0050                      "Are you sure you already converted the geometry?", name.c_str());
0051     }
0052     Geant4ActionSD* action_sd = dynamic_cast<Geant4ActionSD*>(sd);
0053     if (0 == action_sd) {
0054       throw dd4hep::except("Geant4Sensitive", "DDG4: You may only configure actions "
0055                            "for sensitive detectors of type Geant4ActionSD.\n"
0056                            "DDG4: The sensitive detector of %s is of type %s, which is incompatible.", name.c_str(),
0057                            typeName(typeid(*sd)).c_str());
0058     }
0059     return action_sd;
0060   }
0061 }
0062 #endif
0063
0064 /// Standard action constructor
0065 Geant4ActionSD::Geant4ActionSD(const std::string& nam)
0066   : Geant4Action(0, nam) {
0067   InstanceCount::increment(this);
0068 }
0069
0070 /// Default destructor
0071 Geant4ActionSD::~Geant4ActionSD() {
0072   InstanceCount::decrement(this);
0073 }
0074
0075 /// Standard constructor
0076 Geant4Filter::Geant4Filter(Geant4Context* ctxt, const std::string& nam)
0077   : Geant4Action(ctxt, nam) {
0078   InstanceCount::increment(this);
0079 }
0080
0081 /// Standard destructor
0082 Geant4Filter::~Geant4Filter() {
0083   InstanceCount::decrement(this);
0084 }
0085
0086 /// Filter action. Return true if hits should be processed
0087 bool Geant4Filter::operator()(const G4Step*) const {
0088   return true;
0089 }
0090
0091 /// Filter action. Return true if hits should be processed
0092 bool Geant4Filter::operator()(const Geant4FastSimSpot*) const {
0093   except("The filter action %s does not support the GFLASH/FastSim interface for Geant4.", c_name());
0094   return false;
0095 }
0096
0097 /// Constructor. The detector element is identified by the name
0098 Geant4Sensitive::Geant4Sensitive(Geant4Context* ctxt, const std::string& nam, DetElement det, Detector& det_ref)
0099   : Geant4Action(ctxt, nam), m_detDesc(det_ref), m_detector(det)
0100 {
0101   InstanceCount::increment(this);
0102   if (!det.isValid()) {
0103     except("DDG4: Detector elemnt for %s is invalid.", nam.c_str());
0104   }
0105   declareProperty("HitCreationMode", m_hitCreationMode = SIMPLE_MODE);
0106   m_sequence     = context()->kernel().sensitiveAction(m_detector.name());
0107   m_sensitive    = m_detDesc.sensitiveDetector(det.name());
0108   m_readout      = m_sensitive.readout();
0109   m_segmentation = m_readout.segmentation();
0110 }
0111
0112 /// Standard destructor
0113 Geant4Sensitive::~Geant4Sensitive() {
0114   m_filters(&Geant4Filter::release);
0115   m_filters.clear();
0116   InstanceCount::decrement(this);
0117 }
0118
0119 /// Add an actor responding to all callbacks. Sequence takes ownership.
0120 void Geant4Sensitive::adoptFilter(Geant4Action* action)   {
0121   Geant4Filter* filter = dynamic_cast<Geant4Filter*>(action);
0122   adopt(filter);
0123 }
0124
0125 /// Add an actor responding to all callbacks. Sequence takes ownership.
0126 void Geant4Sensitive::adopt(Geant4Filter* filter) {
0127   if (filter) {
0128     filter->addRef();
0129     m_filters.add(filter);
0130     return;
0131   }
0132   except("Attempt to add invalid sensitive filter!");
0133 }
0134
0135 /// Add an actor responding to all callbacks. Sequence takes ownership.
0136 void Geant4Sensitive::adoptFilter_front(Geant4Action* action)   {
0137   Geant4Filter* filter = dynamic_cast<Geant4Filter*>(action);
0138   adopt_front(filter);
0139 }
0140
0141 /// Add an actor responding to all callbacks. Sequence takes ownership.
0142 void Geant4Sensitive::adopt_front(Geant4Filter* filter) {
0143   if (filter) {
0144     filter->addRef();
0145     m_filters.add_front(filter);
0146     return;
0147   }
0148   except("Attempt to add invalid sensitive filter!");
0149 }
0150
0151 /// Callback before hit processing starts. Invoke all filters.
0152 bool Geant4Sensitive::accept(const G4Step* step) const {
0153   bool (Geant4Filter::*filter)(const G4Step*) const = &Geant4Filter::operator();
0154   bool result = m_filters.filter(filter, step);
0155   return result;
0156 }
0157
0158 /// GFLASH/FastSim interface: Callback before hit processing starts. Invoke all filters.
0159 bool Geant4Sensitive::accept(const Geant4FastSimSpot* spot) const {
0160   bool (Geant4Filter::*filter)(const Geant4FastSimSpot*) const = &Geant4Filter::operator();
0161   bool result = m_filters.filter(filter, spot);
0162   return result;
0163 }
0164
0165 /// Access to the sensitive detector object
0166 void Geant4Sensitive::setDetector(Geant4ActionSD* sens_det) {
0167   m_sensitiveDetector = sens_det;
0168 }
0169
0170 /// Access to the sensitive detector object
0171 Geant4ActionSD& Geant4Sensitive::detector() const {
0172   if (m_sensitiveDetector)
0173     return *m_sensitiveDetector;
0174   //m_sensitiveDetector = _getSensitiveDetector(m_detector.name());
0175   //if (  m_sensitiveDetector ) return *m_sensitiveDetector;
0176   except("DDG4: The sensitive detector for action %s was not properly configured.", name().c_str());
0177   throw std::runtime_error("Geant4Sensitive::detector");
0178 }
0179
0180 /// Access to the hosting sequence
0181 Geant4SensDetActionSequence& Geant4Sensitive::sequence() const {
0182   return *m_sequence;
0183 }
0184
0185 /// Access the detector desciption object
0186 dd4hep::Detector& Geant4Sensitive::detectorDescription() const {
0187   return m_detDesc;
0188 }
0189
0190 /// Access HitCollection container names
0191 const std::string& Geant4Sensitive::hitCollectionName(std::size_t which) const {
0192   return sequence().hitCollectionName(which);
0193 }
0194
0195 /// Retrieve the hits collection associated with this detector by its serial number
0196 Geant4HitCollection* Geant4Sensitive::collection(std::size_t which) {
0197   return sequence().collection(which);
0198 }
0199
0200 /// Retrieve the hits collection associated with this detector by its collection identifier
0201 Geant4HitCollection* Geant4Sensitive::collectionByID(std::size_t id) {
0202   return sequence().collectionByID(id);
0203 }
0204
0205 /// Define collections created by this sensitivie action object
0206 void Geant4Sensitive::defineCollections() {
0207 }
0208
0209 /// Method invoked at the begining of each event.
0210 void Geant4Sensitive::begin(G4HCofThisEvent* /* HCE */) {
0211 }
0212
0213 /// Method invoked at the end of each event.
0214 void Geant4Sensitive::end(G4HCofThisEvent* /* HCE */) {
0215 }
0216
0217 /// G4VSensitiveDetector interface: Method for generating hit(s) using the information of G4Step object.
0218 bool Geant4Sensitive::process(const G4Step* /* step */, G4TouchableHistory* /* history */) {
0219   return false;
0220 }
0221
0222 /// GFLASH/FastSim interface: Method for generating hit(s) using the information of the Geant4FastSimSpot object.
0223 bool Geant4Sensitive::processFastSim(const Geant4FastSimSpot* /* spot */, G4TouchableHistory* /* history */) {
0224   except("The sensitive action %s does not support the GFLASH/FastSim interface for Geant4.", c_name());
0225   return false;
0226 }
0227
0228 /// Method is invoked if the event abortion is occured.
0229 void Geant4Sensitive::clear(G4HCofThisEvent* /* HCE */) {
0230 }
0231
0232 /// Mark the track to be kept for MC truth propagation during hit processing
0233 void Geant4Sensitive::mark(const G4Track* track) const  {
0234   Geant4MonteCarloTruth* truth = context()->event().extension<Geant4MonteCarloTruth>(false);
0235   if ( truth ) truth->mark(track);
0236 }
0237
0238 /// Mark the track of this step to be kept for MC truth propagation during hit processing
0239 void Geant4Sensitive::mark(const G4Step* step) const  {
0240   Geant4MonteCarloTruth* truth = context()->event().extension<Geant4MonteCarloTruth>(false);
0241   if ( truth ) truth->mark(step);
0242 }
0243
0244 /// Returns the volumeID of the sensitive volume corresponding to the step
0245 long long int Geant4Sensitive::volumeID(const G4Step* step) {
0246   Geant4StepHandler stepH(step);
0247   Geant4VolumeManager volMgr = Geant4Mapping::instance().volumeManager();
0248   VolumeID id = volMgr.volumeID(stepH.preTouchable());
0249   return id;
0250 }
0251
0252 /// Returns the volumeID of the sensitive volume corresponding to the touchable history
0253 long long int Geant4Sensitive::volumeID(const G4VTouchable* touchable) {
0254   Geant4VolumeManager volMgr = Geant4Mapping::instance().volumeManager();
0255   VolumeID id = volMgr.volumeID(touchable);
0256   return id;
0257 }
0258
0259 /// Returns the cellID(volumeID+local coordinate encoding) of the sensitive volume corresponding to the step
0260 long long int Geant4Sensitive::cellID(const G4Step* step) {
0261   Geant4StepHandler h(step);
0262   Geant4VolumeManager volMgr = Geant4Mapping::instance().volumeManager();
0263   VolumeID volID = volMgr.volumeID(h.preTouchable());
0264   if ( m_segmentation.isValid() )  {
0265     std::exception_ptr eptr;
0266     G4ThreeVector global = 0.5 * (h.prePosG4()+h.postPosG4());
0267     G4ThreeVector local  = h.preTouchable()->GetHistory()->GetTopTransform().TransformPoint(global);
0268     Position loc(local.x()*MM_2_CM, local.y()*MM_2_CM, local.z()*MM_2_CM);
0269     Position glob(global.x()*MM_2_CM, global.y()*MM_2_CM, global.z()*MM_2_CM);
0270     try  {
0271       VolumeID cID = m_segmentation.cellID(loc, glob, volID);
0272       return cID;
0273     }
0274     catch(const std::exception& e)   {
0275       eptr = std::current_exception();
0276       error("cellID: failed to access segmentation for VolumeID: %016lX [%ld]  [%s]", volID, volID, e.what());
0277       error("....... G4-local:   (%f, %f, %f) G4-global:   (%f, %f, %f)",
0278         local.x(), local.y(), local.z(), global.x(), global.y(), global.z());
0279       error("....... TGeo-local: (%f, %f, %f) TGeo-global: (%f, %f, %f)",
0280         loc.x(), loc.y(), loc.z(), glob.x(), glob.y(), glob.z());
0281       error("....... Pre-step:  %s  SD: %s", h.volName(h.pre), h.sdName(h.pre).c_str());
0282       if ( h.post )
0283         error("....... Post-step: %s  SD: %s", h.volName(h.post), h.sdName(h.post).c_str());
0284       std::rethrow_exception(std::move(eptr));
0285     }
0286   }
0287   return volID;
0288 }
0289
0290 /// Returns the cellID(volumeID+local coordinate encoding) of the sensitive volume corresponding to the touchable history
0291 long long int Geant4Sensitive::cellID(const G4VTouchable* touchable, const G4ThreeVector& global) {
0292   Geant4VolumeManager volMgr = Geant4Mapping::instance().volumeManager();
0293   VolumeID volID = volMgr.volumeID(touchable);
0294   if ( m_segmentation.isValid() )  {
0295     std::exception_ptr eptr;
0296     G4ThreeVector local  = touchable->GetHistory()->GetTopTransform().TransformPoint(global);
0297     Position loc (local.x()*MM_2_CM, local.y()*MM_2_CM, local.z()*MM_2_CM);
0298     Position glob(global.x()*MM_2_CM, global.y()*MM_2_CM, global.z()*MM_2_CM);
0299     try  {
0300       VolumeID cID = m_segmentation.cellID(loc, glob, volID);
0301       return cID;
0302     }
0303     catch(const std::exception& e)   {
0304       auto* pvol = touchable->GetVolume();
0305       auto* vol = pvol->GetLogicalVolume();
0306       auto* sd = vol->GetSensitiveDetector();
0307       eptr = std::current_exception();
0308       error("cellID: failed to access segmentation for VolumeID: %016lX [%ld]  [%s]", volID, volID, e.what());
0309       error("....... G4-local:   (%f, %f, %f) G4-global:   (%f, %f, %f)",
0310         local.x(), local.y(), local.z(), global.x(), global.y(), global.z());
0311       error("....... TGeo-local: (%f, %f, %f) TGeo-global: (%f, %f, %f)",
0312         loc.x(), loc.y(), loc.z(), glob.x(), glob.y(), glob.z());
0313       error("....... Touchable:  %s  SD: %s", vol->GetName().c_str(), sd ? sd->GetName().c_str() : "???");
0314       std::rethrow_exception(std::move(eptr));
0315     }
0316   }
0317   return volID;
0318 }
0319
0320 /// Standard constructor
0321 Geant4SensDetActionSequence::Geant4SensDetActionSequence(Geant4Context* ctxt, const std::string& nam)
0322   : Geant4Action(ctxt, nam), m_hce(0), m_detector(0)
0323 {
0324   m_needsControl = true;
0325   context()->sensitiveActions().insert(name(), this);
0326   /// Update the sensitive detector type, so that the proper instance is created
0327   m_sensitive = context()->detectorDescription().sensitiveDetector(nam);
0328   m_sensitiveType = m_sensitive.type();
0329   InstanceCount::increment(this);
0330 }
0331
0332 /// Default destructor
0333 Geant4SensDetActionSequence::~Geant4SensDetActionSequence() {
0334   m_filters(&Geant4Filter::release);
0335   m_actors(&Geant4Sensitive::release);
0336   m_filters.clear();
0337   m_actors.clear();
0338   InstanceCount::decrement(this);
0339 }
0340
0341 /// Set or update client context
0342 void Geant4SensDetActionSequence::updateContext(Geant4Context* ctxt)    {
0343   m_context = ctxt;
0344   m_actors.updateContext(ctxt);
0345   m_filters.updateContext(ctxt);
0346 }
0347
0348 /// Add an actor responding to all callbacks. Sequence takes ownership.
0349 void Geant4SensDetActionSequence::adoptFilter(Geant4Action* action)   {
0350   Geant4Filter* filter = dynamic_cast<Geant4Filter*>(action);
0351   adopt(filter);
0352 }
0353
0354 /// Add an actor responding to all callbacks
0355 void Geant4SensDetActionSequence::adopt(Geant4Sensitive* sensitive) {
0356   if (sensitive) {
0357     sensitive->addRef();
0358     m_actors.add(sensitive);
0359     return;
0360   }
0361   except("Attempt to add invalid sensitive actor!");
0362 }
0363
0364 /// Add an actor responding to all callbacks. Sequence takes ownership.
0365 void Geant4SensDetActionSequence::adopt(Geant4Filter* filter) {
0366   if (filter) {
0367     filter->addRef();
0368     m_filters.add(filter);
0369     return;
0370   }
0371   except("Attempt to add invalid sensitive filter!");
0372 }
0373
0374 /// Initialize the usage of a hit collection. Returns the collection identifier
0375 std::size_t Geant4SensDetActionSequence::defineCollection(Geant4Sensitive* owner, const std::string& collection_name, create_t func) {
0376   m_collections.emplace_back(collection_name, make_pair(owner,func));
0377   return m_collections.size() - 1;
0378 }
0379
0380 /// Called at construction time of the sensitive detector to declare all hit collections
0381 std::size_t Geant4SensDetActionSequence::Geant4SensDetActionSequence::defineCollections(Geant4ActionSD* sens_det) {
0382   std::size_t count = 0;
0383   m_detector = sens_det;
0384   m_actors(&Geant4Sensitive::setDetector, sens_det);
0385   m_actors(&Geant4Sensitive::defineCollections);
0386   for (HitCollections::const_iterator i = m_collections.begin(); i != m_collections.end(); ++i) {
0387     sens_det->defineCollection((*i).first);
0388     ++count;
0389   }
0390   return count;
0391 }
0392
0393 /// Access HitCollection container names
0394 const std::string& Geant4SensDetActionSequence::hitCollectionName(std::size_t which) const {
0395   if (which < m_collections.size()) {
0396     return m_collections[which].first;
0397   }
0398   static std::string blank = "";
0399   except("The collection name index for subdetector %s is out of range!", c_name());
0400   return blank;
0401 }
0402
0403 /// Retrieve the hits collection associated with this detector by its serial number
0404 Geant4HitCollection* Geant4SensDetActionSequence::collection(std::size_t which) const {
0405   if (which < m_collections.size()) {
0406     int hc_id = m_detector->GetCollectionID(which);
0407     Geant4HitCollection* c = (Geant4HitCollection*) m_hce->GetHC(hc_id);
0408     if (c)
0409       return c;
0410     except("The collection index for subdetector %s is wrong!", c_name());
0411   }
0412   except("The collection name index for subdetector %s is out of range!", c_name());
0413   return 0;
0414 }
0415
0416 /// Retrieve the hits collection associated with this detector by its collection identifier
0417 Geant4HitCollection* Geant4SensDetActionSequence::collectionByID(std::size_t id) const {
0418   Geant4HitCollection* c = (Geant4HitCollection*) m_hce->GetHC(id);
0419   if (c)
0420     return c;
0421   except("The collection index for subdetector %s is wrong!", c_name());
0422   return 0;
0423 }
0424
0425 /// Callback before hit processing starts. Invoke all filters.
0426 bool Geant4SensDetActionSequence::accept(const G4Step* step) const {
0427   bool (Geant4Filter::*filter)(const G4Step*) const = &Geant4Filter::operator();
0428   bool result = m_filters.filter(filter, step);
0429   return result;
0430 }
0431
0432 /// Callback before hit processing starts. Invoke all filters.
0433 bool Geant4SensDetActionSequence::accept(const Geant4FastSimSpot* spot) const {
0434   bool (Geant4Filter::*filter)(const Geant4FastSimSpot*) const = &Geant4Filter::operator();
0435   bool result = m_filters.filter(filter, spot);
0436   return result;
0437 }
0438
0439 /// G4VSensitiveDetector interface: Method for generating hit(s) using the information of G4Step object.
0440 bool Geant4SensDetActionSequence::process(const G4Step* step, G4TouchableHistory* history) {
0441   bool result = false;
0442   for (Geant4Sensitive* sensitive : m_actors)  {
0443     if ( sensitive->accept(step) )
0444       result |= sensitive->process(step, history);
0445   }
0446   m_process(step, history);
0447   return result;
0448 }
0449
0450 /// GFLASH/FastSim interface: Method for generating hit(s) using the information of the Geant4FastSimSpot object.
0451 bool Geant4SensDetActionSequence::processFastSim(const Geant4FastSimSpot* spot, G4TouchableHistory* history)  {
0452   bool result = false;
0453   for (Geant4Sensitive* sensitive : m_actors)  {
0454     if ( sensitive->accept(spot) )
0455       result |= sensitive->processFastSim(spot, history);
0456   }
0457   m_process(spot, history);
0458   return result;
0459 }
0460
0461 /** G4VSensitiveDetector interface: Method invoked at the begining of each event.
0462  *  The hits collection(s) created by this sensitive detector must
0463  *  be set to the G4HCofThisEvent object at one of these two methods.
0464  */
0465 void Geant4SensDetActionSequence::begin(G4HCofThisEvent* hce) {
0466   m_hce = hce;
0467   for (std::size_t count = 0; count < m_collections.size(); ++count) {
0468     const HitCollection& cr = m_collections[count];
0469     Geant4HitCollection* col = (*cr.second.second)(name(), cr.first, cr.second.first);
0470     int id = m_detector->GetCollectionID(count);
0471     m_hce->AddHitsCollection(id, col);
0472   }
0473   m_actors(&Geant4Sensitive::begin, m_hce);
0474   m_begin (m_hce);
0475 }
0476
0477 /// G4VSensitiveDetector interface: Method invoked at the end of each event.
0478 void Geant4SensDetActionSequence::end(G4HCofThisEvent* hce) {
0479   m_end(hce);
0480   m_actors(&Geant4Sensitive::end, hce);
0481   // G4HCofThisEvent must be availible until end-event. m_hce = 0;
0482 }
0483
0484 /// G4VSensitiveDetector interface: Method invoked if the event was aborted.
0485 /** Hits collections created but not being set to G4HCofThisEvent
0486  *  at the event should be deleted.
0487  *  Collection(s) which have already set to G4HCofThisEvent
0488  *  will be deleted automatically.
0489  */
0490 void Geant4SensDetActionSequence::clear() {
0491   m_clear (m_hce);
0492   m_actors(&Geant4Sensitive::clear, m_hce);
0493 }
0494
0495 /// Default destructor
0496 Geant4SensDetSequences::~Geant4SensDetSequences() {
0497   detail::releaseObjects(m_sequences);
0498   m_sequences.clear();
0499 }
0500
0501 /// Access sequence member by name
0502 Geant4SensDetActionSequence* Geant4SensDetSequences::operator[](const std::string& nam) const {
0503   std::string n = "SD_Seq_" + nam;
0504   Members::const_iterator i = m_sequences.find(n);
0505   if (i != m_sequences.end())
0506     return (*i).second;
0507   except("Attempt to access undefined SensDetActionSequence: %s ", nam.c_str());
0508   return nullptr;
0509 }
0510
0511 /// Access sequence member by name
0512 Geant4SensDetActionSequence* Geant4SensDetSequences::find(const std::string& name) const {
0513   std::string nam = "SD_Seq_" + name;
0514   Members::const_iterator i = m_sequences.find(nam);
0515   if (i != m_sequences.end())
0516     return (*i).second;
0517   return 0;
0518 }
0519
0520 /// Insert sequence member
0521 void Geant4SensDetSequences::insert(const std::string& name, Geant4SensDetActionSequence* seq) {
0522   if (seq) {
0523     std::string nam = "SD_Seq_" + name;
0524     seq->addRef();
0525     m_sequences[nam] = seq;
0526     return;
0527   }
0528   except("Attempt to add invalid sensitive sequence with name:%s", name.c_str());
0529 }
0530
0531 /// Clear the sequence list
0532 void Geant4SensDetSequences::clear()   {
0533   m_sequences.clear();
0534 }