SensitiveDetector/src/SensitiveDetector.cc

0001 #include "SimG4Core/SensitiveDetector/interface/SensitiveDetector.h"
0002
0003 #include "FWCore/Utilities/interface/isFinite.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005
0006 #include "G4SDManager.hh"
0007 #include "G4Step.hh"
0008 #include "G4Track.hh"
0009 #include "G4StepPoint.hh"
0010 #include "G4Transform3D.hh"
0011 #include "G4LogicalVolumeStore.hh"
0012 #include "G4TouchableHistory.hh"
0013 #include "G4VUserTrackInformation.hh"
0014
0015 #include <sstream>
0016
0017 SensitiveDetector::SensitiveDetector(const std::string& iname,
0018                                      const SensitiveDetectorCatalog& clg,
0019                                      bool calo,
0020                                      const std::string& newcollname)
0021     : G4VSensitiveDetector(iname), m_isCalo(calo) {
0022   // for CMS hits
0023   m_namesOfSD.push_back(iname);
0024
0025   // Geant4 hit collection
0026   collectionName.insert(iname);
0027   if (!newcollname.empty())
0028     collectionName.insert(newcollname);
0029
0030   // register sensitive detector
0031   G4SDManager* SDman = G4SDManager::GetSDMpointer();
0032   SDman->AddNewDetector(this);
0033
0034   const std::vector<std::string_view>& lvNames = clg.logicalNames(iname);
0035   std::stringstream ss;
0036   for (auto& lvname : lvNames) {
0037     this->AssignSD({lvname.data(), lvname.size()});
0038     ss << " " << lvname << "\n";
0039   }
0040   if (newcollname.empty())
0041     ss << " Collection " << iname;
0042   else
0043     ss << " Collections " << iname << " and " << newcollname;
0044   edm::LogVerbatim("SensitiveDetector") << " <" << iname << "> : Assigns SD to " << lvNames.size() << " LVs "
0045                                         << ss.str();
0046 }
0047
0048 SensitiveDetector::~SensitiveDetector() {}
0049
0050 void SensitiveDetector::Initialize(G4HCofThisEvent* eventHC) {}
0051
0052 void SensitiveDetector::EndOfEvent(G4HCofThisEvent* eventHC) {}
0053
0054 void SensitiveDetector::AssignSD(const std::string& vname) {
0055   G4LogicalVolumeStore* theStore = G4LogicalVolumeStore::GetInstance();
0056   for (auto& lv : *theStore) {
0057     if (vname == lv->GetName()) {
0058       lv->SetSensitiveDetector(this);
0059     }
0060   }
0061 }
0062
0063 Local3DPoint SensitiveDetector::InitialStepPosition(const G4Step* step, coordinates cd) const {
0064   const G4StepPoint* preStepPoint = step->GetPreStepPoint();
0065   const G4ThreeVector& globalCoordinates = preStepPoint->GetPosition();
0066   if (cd == WorldCoordinates) {
0067     return ConvertToLocal3DPoint(globalCoordinates);
0068   }
0069   const G4TouchableHistory* theTouchable = static_cast<const G4TouchableHistory*>(preStepPoint->GetTouchable());
0070   const G4ThreeVector localCoordinates =
0071       theTouchable->GetHistory()->GetTopTransform().TransformPoint(globalCoordinates);
0072   return ConvertToLocal3DPoint(localCoordinates);
0073 }
0074
0075 Local3DPoint SensitiveDetector::FinalStepPosition(const G4Step* step, coordinates cd) const {
0076   const G4StepPoint* postStepPoint = step->GetPostStepPoint();
0077   const G4ThreeVector& globalCoordinates = postStepPoint->GetPosition();
0078   if (cd == WorldCoordinates) {
0079     return ConvertToLocal3DPoint(globalCoordinates);
0080   }
0081   const G4StepPoint* preStepPoint = step->GetPreStepPoint();
0082   const G4ThreeVector localCoordinates =
0083       preStepPoint->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(globalCoordinates);
0084   return ConvertToLocal3DPoint(localCoordinates);
0085 }
0086
0087 Local3DPoint SensitiveDetector::LocalPreStepPosition(const G4Step* step) const {
0088   const G4StepPoint* preStepPoint = step->GetPreStepPoint();
0089   G4ThreeVector localCoordinates =
0090       preStepPoint->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(preStepPoint->GetPosition());
0091   return ConvertToLocal3DPoint(localCoordinates);
0092 }
0093
0094 Local3DPoint SensitiveDetector::LocalPostStepPosition(const G4Step* step) const {
0095   const G4ThreeVector& globalCoordinates = step->GetPostStepPoint()->GetPosition();
0096   G4ThreeVector localCoordinates =
0097       step->GetPreStepPoint()->GetTouchable()->GetHistory()->GetTopTransform().TransformPoint(globalCoordinates);
0098   return ConvertToLocal3DPoint(localCoordinates);
0099 }
0100
0101 TrackInformation* SensitiveDetector::cmsTrackInformation(const G4Track* aTrack) {
0102   TrackInformation* info = (TrackInformation*)(aTrack->GetUserInformation());
0103   if (nullptr == info) {
0104     edm::LogWarning("SensitiveDetector") << " no TrackInformation available for trackID= " << aTrack->GetTrackID()
0105                                          << " inside SD " << GetName();
0106     G4Exception(
0107         "SensitiveDetector::cmsTrackInformation()", "sd01", FatalException, "cannot handle hits without trackinfo");
0108   }
0109   return info;
0110 }
0111
0112 void SensitiveDetector::setNames(const std::vector<std::string>& hnames) {
0113   m_namesOfSD.clear();
0114   m_namesOfSD = hnames;
0115 }
0116
0117 void SensitiveDetector::NaNTrap(const G4Step* aStep) const {
0118   G4Track* currentTrk = aStep->GetTrack();
0119   double ekin = currentTrk->GetKineticEnergy();
0120   if (ekin < 0.0) {
0121     const G4VPhysicalVolume* pCurrentVol = aStep->GetPreStepPoint()->GetPhysicalVolume();
0122     edm::LogWarning("SensitiveDetector") << "Negative kinetic energy Ekin(MeV)=" << ekin / CLHEP::MeV << " of "
0123                                          << currentTrk->GetDefinition()->GetParticleName()
0124                                          << " trackID= " << currentTrk->GetTrackID() << " inside "
0125                                          << pCurrentVol->GetName();
0126     currentTrk->SetKineticEnergy(0.0);
0127   }
0128   const G4ThreeVector& currentPos = currentTrk->GetPosition();
0129   double xyz = currentPos.x() + currentPos.y() + currentPos.z();
0130   const G4ThreeVector& currentMom = currentTrk->GetMomentum();
0131   xyz += currentMom.x() + currentMom.y() + currentMom.z();
0132
0133   if (edm::isNotFinite(xyz)) {
0134     const G4VPhysicalVolume* pCurrentVol = aStep->GetPreStepPoint()->GetPhysicalVolume();
0135     edm::LogWarning("SensitiveDetector") << "NaN detected for trackID= " << currentTrk->GetTrackID() << " inside "
0136                                          << pCurrentVol->GetName();
0137     G4Exception("SensitiveDetector::NaNTrap()", "sd01", FatalException, "corrupted event or step");
0138   }
0139 }