Merge pull request cms-sw#40127 from civanch/fix_SD_crash

Protection for negative kinetic energy in Geant4 tracking
jfernan2 · Nov 23, 2022 · d2bfd2c · d2bfd2c
2 parents 0bcb76f + e56853c
commit d2bfd2c
Show file tree

Hide file tree

Showing 3 changed files with 21 additions and 2 deletions.
diff --git a/SimG4Core/Application/src/ExceptionHandler.cc b/SimG4Core/Application/src/ExceptionHandler.cc
@@ -33,7 +33,7 @@ bool ExceptionHandler::Notify(const char* exceptionOrigin,
   // part of exception happens outside tracking loop
   if (nullptr != track) {
     ekin = track->GetKineticEnergy();
-    message << "\n"
+    message << "\n CMS info: "
             << "TrackID=" << track->GetTrackID() << " ParentID=" << track->GetParentID() << "  "
             << track->GetParticleDefinition()->GetParticleName() << "; Ekin(MeV)=" << ekin / CLHEP::MeV
             << "; time(ns)=" << track->GetGlobalTime() / CLHEP::ns << "; status=" << track->GetTrackStatus()

diff --git a/SimG4Core/Application/src/SteppingAction.cc b/SimG4Core/Application/src/SteppingAction.cc
@@ -98,6 +98,16 @@ void SteppingAction::UserSteppingAction(const G4Step* aStep) {
   G4Track* theTrack = aStep->GetTrack();
   TrackStatus tstat = (theTrack->GetTrackStatus() == fAlive) ? sAlive : sKilledByProcess;
 
+  if (theTrack->GetKineticEnergy() < 0.0) {
+    if (nWarnings < 5) {
+      ++nWarnings;
+      edm::LogWarning("SimG4CoreApplication")
+          << "Track #" << theTrack->GetTrackID() << " " << theTrack->GetDefinition()->GetParticleName()
+          << " Ekin(MeV)= " << theTrack->GetKineticEnergy() / MeV;
+    }
+    theTrack->SetKineticEnergy(0.0);
+  }
+
   const G4StepPoint* preStep = aStep->GetPreStepPoint();
   const G4StepPoint* postStep = aStep->GetPostStepPoint();
 

diff --git a/SimG4Core/SensitiveDetector/src/SensitiveDetector.cc b/SimG4Core/SensitiveDetector/src/SensitiveDetector.cc
@@ -105,7 +105,16 @@ void SensitiveDetector::setNames(const std::vector<std::string>& hnames) {
 }
 
 void SensitiveDetector::NaNTrap(const G4Step* aStep) const {
-  const G4Track* currentTrk = aStep->GetTrack();
+  G4Track* currentTrk = aStep->GetTrack();
+  double ekin = currentTrk->GetKineticEnergy();
+  if (ekin < 0.0) {
+    const G4VPhysicalVolume* pCurrentVol = aStep->GetPreStepPoint()->GetPhysicalVolume();
+    edm::LogWarning("SensitiveDetector") << "Negative kinetic energy Ekin(MeV)=" << ekin / CLHEP::MeV << " of "
+                                         << currentTrk->GetDefinition()->GetParticleName()
+                                         << " trackID= " << currentTrk->GetTrackID() << " inside "
+                                         << pCurrentVol->GetName();
+    currentTrk->SetKineticEnergy(0.0);
+  }
   const G4ThreeVector& currentPos = currentTrk->GetPosition();
   double xyz = currentPos.x() + currentPos.y() + currentPos.z();
   const G4ThreeVector& currentMom = currentTrk->GetMomentum();