fix memory corruption in scoring (#326)

This PR fixes a major memory corruption problem in the scoring:

Previously, the + 1 here:
```
      ::new (toucheableHandleStorage) G4TouchableHandle{::new (toucheableHistoryStorage + 1) G4TouchableHistory};
```

Was a too small offset. Thus, the `fPostG4TouchableHistoryHandle`
overlapped with the `fPreG4TouchableHistoryHandle` and overwrote its
values. This is now fixed. The error broke the physics validation (note:
the physics validation was done using the latest version of G4HepEm
which will be made available via #324):

Before this PR:

![Screenshot from 2024-12-11
16-55-26](https://github.com/user-attachments/assets/cce5bcc8-2fcb-4361-b683-6361b7ddd1e2)
(the absolute scale here is off due to a change in the number of
particles that was not incorporated into the plot, the relative error is
correct)
After this PR:
![Screenshot from 2024-12-13
15-47-10](https://github.com/user-attachments/assets/99afb009-85f3-4a7c-becf-e3586d313620)



Additionally, some quality of life improvements are added:

1. The way of accessing the correct ID was reduced from:
```
hitID = fScoringMap[aStep->GetPreStepPoint()->GetTouchable()->GetHistory()->GetTopVolume()->GetInstanceID()];
```
to
```
std::size_t hitID = fScoringMap[aStep->GetPreStepPoint()->GetPhysicalVolume()->GetInstanceID()];
```
2. some unused functions were removed
3. some doubled includes were removed
4. a new option of `--accumulated_events` was added to the
integrationBenchmark example. Using that option, the edep of all events
is written to a single line in the output.csv.

examples/AsyncExample/include/SensitiveDetector.hh

@@ -56,9 +56,7 @@ public:
   /// Process energy deposit from the full simulation.
   virtual G4bool ProcessHits(G4Step *aStep, G4TouchableHistory *aROhist) final;
 
-  virtual G4bool ProcessHits(int hitID, double energy) final;
-
-  SimpleHit *RetrieveAndSetupHit(G4TouchableHistory *aTouchable);
+  SimpleHit *RetrieveAndSetupHit(std::size_t hitID);
 
   std::vector<G4LogicalVolume *> fSensitiveLogicalVolumes;
   /// Physical Volumes where we want to score

examples/AsyncExample/src/SensitiveDetector.cc

@@ -26,10 +26,8 @@
 // ********************************************************************
 //
 #include "SensitiveDetector.hh"
-#include "SimpleHit.hh"
 
 #include "G4HCofThisEvent.hh"
-#include "G4TouchableHistory.hh"
 #include "G4Step.hh"
 #include "G4SDManager.hh"
 
@@ -68,9 +66,9 @@ G4bool SensitiveDetector::ProcessHits(G4Step *aStep, G4TouchableHistory *)
   G4double edep = aStep->GetTotalEnergyDeposit();
   if (edep == 0.) return true;
 
-  G4TouchableHistory *aTouchable = (G4TouchableHistory *)(aStep->GetPreStepPoint()->GetTouchable());
+  std::size_t hitID = fScoringMap[aStep->GetPreStepPoint()->GetPhysicalVolume()->GetInstanceID()];
 
-  auto hit = RetrieveAndSetupHit(aTouchable);
+  auto hit = RetrieveAndSetupHit(hitID);
 
   // Add energy deposit from G4Step
   hit->AddEdep(edep);
@@ -85,22 +83,6 @@ G4bool SensitiveDetector::ProcessHits(G4Step *aStep, G4TouchableHistory *)
 
   return true;
 }
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-G4bool SensitiveDetector::ProcessHits(int hitID, double energy)
-{
-
-  if (energy == 0.) return true;
-
-  auto hit = (*fHitsCollection)[hitID];
-
-  // Add energy deposit from G4FastHit
-  hit->AddEdep(energy);
-  // set type to fast sim
-  hit->SetType(1);
-
-  return true;
-}
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
@@ -117,4 +99,4 @@ SimpleHit *SensitiveDetector::RetrieveAndSetupHit(G4TouchableHistory *aTouchable
   }
 
   return (*fHitsCollection)[hitID];
-}
+}

examples/Example1/include/SensitiveDetector.hh

@@ -57,7 +57,7 @@ public:
   /// Process energy deposit from the full simulation.
   virtual G4bool ProcessHits(G4Step *aStep, G4TouchableHistory *aROhist) final;
 
-  SimpleHit *RetrieveAndSetupHit(G4TouchableHistory *aTouchable);
+  SimpleHit *RetrieveAndSetupHit(std::size_t hitID);
 
   std::vector<G4LogicalVolume *> fSensitiveLogicalVolumes;
   /// Physical Volumes where we want to score

examples/Example1/src/SensitiveDetector.cc

@@ -26,10 +26,8 @@
 // ********************************************************************
 //
 #include "SensitiveDetector.hh"
-#include "SimpleHit.hh"
 
 #include "G4HCofThisEvent.hh"
-#include "G4TouchableHistory.hh"
 #include "G4Step.hh"
 #include "G4SDManager.hh"
 
@@ -69,9 +67,9 @@ G4bool SensitiveDetector::ProcessHits(G4Step *aStep, G4TouchableHistory *)
   G4double edep = aStep->GetTotalEnergyDeposit();
   if (edep == 0.) return true;
 
-  G4TouchableHistory *aTouchable = (G4TouchableHistory *)(aStep->GetPreStepPoint()->GetTouchable());
+  std::size_t hitID = fScoringMap[aStep->GetPreStepPoint()->GetPhysicalVolume()->GetInstanceID()];
 
-  auto hit = RetrieveAndSetupHit(aTouchable);
+  auto hit = RetrieveAndSetupHit(hitID);
 
   // Add energy deposit from G4Step
   hit->AddEdep(edep);
@@ -89,9 +87,8 @@ G4bool SensitiveDetector::ProcessHits(G4Step *aStep, G4TouchableHistory *)
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-SimpleHit *SensitiveDetector::RetrieveAndSetupHit(G4TouchableHistory *aTouchable)
+SimpleHit *SensitiveDetector::RetrieveAndSetupHit(std::size_t hitID)
 {
-  std::size_t hitID = fScoringMap[aTouchable->GetHistory()->GetTopVolume()->GetInstanceID()];
 
   assert(hitID < fNumSensitive);
 
@@ -102,4 +99,4 @@ SimpleHit *SensitiveDetector::RetrieveAndSetupHit(G4TouchableHistory *aTouchable
   }
 
   return (*fHitsCollection)[hitID];
-}
+}

examples/IntegrationBenchmark/include/ActionInitialisation.hh

@@ -43,7 +43,8 @@ class DetectorConstruction;
 
 class ActionInitialisation : public G4VUserActionInitialization {
 public:
-  ActionInitialisation(G4String aOutputDirectory, G4String aOutputFilename, bool aDoBenchmark, bool aDoValidation);
+  ActionInitialisation(G4String aOutputDirectory, G4String aOutputFilename, bool aDoBenchmark, bool aDoValidation,
+                       bool aDoAccumulatdEvents);
   ~ActionInitialisation();
   /// Create all user actions.
   virtual void Build() const final;
@@ -57,6 +58,7 @@ private:
   G4String fOutputFilename;
   bool fDoBenchmark;
   bool fDoValidation;
+  bool fDoAccumulatedEvents;
 };
 
 #endif /* ACTIONINITIALISATION_HH */

examples/IntegrationBenchmark/include/RunAction.hh

@@ -48,7 +48,8 @@ class RunAction : public G4UserRunAction {
 public:
   /// Constructor. Defines the histograms.
   RunAction();
-  RunAction(G4String aOutputDirectory, G4String aOutputFilename, bool aDoBenchmark, bool aDoValidation);
+  RunAction(G4String aOutputDirectory, G4String aOutputFilename, bool aDoBenchmark, bool aDoValidation,
+            bool aDoAccumulatedEvents);
   virtual ~RunAction();
 
   /// Open the file for the analysis
@@ -58,10 +59,11 @@ public:
 
   G4Run *GenerateRun() override;
 
-  bool &GetDoBenchmark(){return fDoBenchmark;};
-  bool &GetDoValidation(){return fDoValidation;};
-  const G4String &GetOutputDirectory(){return fOutputDirectory;};
-  const G4String &GetOutputFilename(){return fOutputFilename;};
+  bool &GetDoBenchmark() { return fDoBenchmark; };
+  bool &GetDoValidation() { return fDoValidation; };
+  bool &GetDoAccumulatedEvents() { return fDoAccumulatedEvents; };
+  const G4String &GetOutputDirectory() { return fOutputDirectory; };
+  const G4String &GetOutputFilename() { return fOutputFilename; };
 
 private:
   /// Pointer to detector construction to retrieve the detector dimensions to
@@ -70,6 +72,7 @@ private:
   G4String fOutputFilename;
   bool fDoBenchmark;
   bool fDoValidation;
+  bool fDoAccumulatedEvents;
   vecgeom::Stopwatch fTimer;
   Run *fRun;
 };

examples/IntegrationBenchmark/include/SensitiveDetector.hh

@@ -58,7 +58,7 @@ public:
   /// Process energy deposit from the full simulation.
   virtual G4bool ProcessHits(G4Step *aStep, G4TouchableHistory *aROhist) final;
 
-  SimpleHit *RetrieveAndSetupHit(G4TouchableHistory *aTouchable);
+  SimpleHit *RetrieveAndSetupHit(std::size_t hitID);
 
   std::vector<G4LogicalVolume *> fSensitiveLogicalVolumes;
   /// Physical Volumes where we want to score

examples/IntegrationBenchmark/integrationBenchmark.cpp

@@ -24,6 +24,7 @@ int main(int argc, char **argv)
   G4String outputFilename   = "";
   bool doBenchmark          = false;
   bool doValidation         = false;
+  bool doAccumlatedEvents   = false; // whether the edep is accumulated across events in the validation csv file
   G4bool useInteractiveMode = true;
   G4bool useAdePT           = true;
   G4bool allSensitive = false; // If set, ignores the sensitive detector flags in the GDML and marks all volumes as
@@ -51,6 +52,8 @@ int main(int argc, char **argv)
       doBenchmark = true;
     } else if (argument == "--do_validation") {
       doValidation = true;
+    } else if (argument == "--accumulated_events") {
+      doAccumlatedEvents = true;
     } else if (argument == "--noadept") {
       useAdePT = false;
     } else if (argument == "--allsensitive") {
@@ -102,7 +105,7 @@ int main(int argc, char **argv)
   // UserAction classes
   //-------------------------------
   runManager->SetUserInitialization(
-      new ActionInitialisation(outputDirectory, outputFilename, doBenchmark, doValidation));
+      new ActionInitialisation(outputDirectory, outputFilename, doBenchmark, doValidation, doAccumlatedEvents));
 
   G4UImanager *UImanager = G4UImanager::GetUIpointer();
   G4String command       = "/control/execute ";

examples/IntegrationBenchmark/src/ActionInitialisation.cc

@@ -33,9 +33,9 @@
 #include "SteppingAction.hh"
 
 ActionInitialisation::ActionInitialisation(G4String aOutputDirectory, G4String aOutputFilename, bool aDoBenchmark,
-                                           bool aDoValidation)
+                                           bool aDoValidation, bool aDoAccumulatedEvents)
     : G4VUserActionInitialization(), fOutputDirectory(aOutputDirectory), fOutputFilename(aOutputFilename),
-      fDoBenchmark(aDoBenchmark), fDoValidation(aDoValidation)
+      fDoBenchmark(aDoBenchmark), fDoValidation(aDoValidation), fDoAccumulatedEvents(aDoAccumulatedEvents)
 {
 }
 
@@ -48,15 +48,16 @@ ActionInitialisation::~ActionInitialisation() {}
 void ActionInitialisation::BuildForMaster() const
 {
   new PrimaryGeneratorAction();
-  SetUserAction(new RunAction(fOutputDirectory, fOutputFilename, fDoBenchmark, fDoValidation));
+  SetUserAction(new RunAction(fOutputDirectory, fOutputFilename, fDoBenchmark, fDoValidation, fDoAccumulatedEvents));
 }
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 void ActionInitialisation::Build() const
 {
   SetUserAction(new PrimaryGeneratorAction());
-  RunAction *aRunAction = new RunAction(fOutputDirectory, fOutputFilename, fDoBenchmark, fDoValidation);
+  RunAction *aRunAction =
+      new RunAction(fOutputDirectory, fOutputFilename, fDoBenchmark, fDoValidation, fDoAccumulatedEvents);
   SetUserAction(aRunAction);
   SetUserAction(new EventAction(aRunAction));
   TrackingAction *aTrackingAction = new TrackingAction();

examples/IntegrationBenchmark/src/EventAction.cc

@@ -143,14 +143,24 @@ void EventAction::EndOfEventAction(const G4Event *aEvent)
       auto id = hit->GetPhysicalVolumeId() + Run::accumulators::NUM_ACCUMULATORS;
       // Reset the accumulator from the last event
       // aTestManager->setAccumulator(id, 0);
-      // Set the accumulator
-      aTestManager->setAccumulator(id, hit->GetEdep());
-      // aTestManager->addToAccumulator(id, hit->GetEdep());
+
+      if (!fRunAction->GetDoAccumulatedEvents()) {
+        // Set the accumulator
+        aTestManager->setAccumulator(id, hit->GetEdep());
+      } else {
+        // write all events to one accumulator:
+        aTestManager->addToAccumulator(id, hit->GetEdep());
+      }
     }
 
-    // Write data to output file. Validation data can take a lot of memory, and we don't need to aggregate
-    // the results of multiple events at runtime, so for better performance it's easier to write the output here
-    aTestManager->exportCSV(false);
+    if (!fRunAction->GetDoAccumulatedEvents()) {
+      // Write data to output file. Validation data can take a lot of memory, and we don't need to aggregate
+      // the results of multiple events at runtime, so for better performance it's easier to write the output here
+      // aTestManager->exportCSV(false);
+    } else {
+      // overwrite to have all validation data written into a single line for all events
+      aTestManager->exportCSV(true);
+    }
 
     // Store test manager
     // TestManagerStore<int>::GetInstance()->RecordState(aTestManager);

examples/IntegrationBenchmark/src/RunAction.cc

@@ -35,9 +35,10 @@
 
 RunAction::RunAction() : G4UserRunAction(), fOutputDirectory(""), fOutputFilename("") {}
 
-RunAction::RunAction(G4String aOutputDirectory, G4String aOutputFilename, bool aDoBenchmark, bool aDoValidation)
+RunAction::RunAction(G4String aOutputDirectory, G4String aOutputFilename, bool aDoBenchmark, bool aDoValidation,
+                     bool aDoAccumulatedEvents)
     : G4UserRunAction(), fOutputDirectory(aOutputDirectory), fOutputFilename(aOutputFilename),
-      fDoBenchmark(aDoBenchmark), fDoValidation(aDoValidation)
+      fDoBenchmark(aDoBenchmark), fDoValidation(aDoValidation), fDoAccumulatedEvents(aDoAccumulatedEvents)
 {
 }
 

examples/IntegrationBenchmark/src/SensitiveDetector.cc

@@ -26,10 +26,8 @@
 // ********************************************************************
 //
 #include "SensitiveDetector.hh"
-#include "SimpleHit.hh"
 
 #include "G4HCofThisEvent.hh"
-#include "G4TouchableHistory.hh"
 #include "G4Step.hh"
 #include "G4SDManager.hh"
 
@@ -70,9 +68,9 @@ G4bool SensitiveDetector::ProcessHits(G4Step *aStep, G4TouchableHistory *)
   G4double edep = aStep->GetTotalEnergyDeposit();
   if (edep == 0.) return true;
 
-  G4TouchableHistory *aTouchable = (G4TouchableHistory *)(aStep->GetPreStepPoint()->GetTouchable());
+  std::size_t hitID = fScoringMap[aStep->GetPreStepPoint()->GetPhysicalVolume()->GetInstanceID()];
 
-  auto hit = RetrieveAndSetupHit(aTouchable);
+  auto hit = RetrieveAndSetupHit(hitID);
 
   // Add energy deposit from G4Step
   hit->AddEdep(edep);
@@ -90,10 +88,8 @@ G4bool SensitiveDetector::ProcessHits(G4Step *aStep, G4TouchableHistory *)
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-SimpleHit *SensitiveDetector::RetrieveAndSetupHit(G4TouchableHistory *aTouchable)
+SimpleHit *SensitiveDetector::RetrieveAndSetupHit(std::size_t hitID)
 {
-  std::size_t hitID = fScoringMap[aTouchable->GetHistory()->GetTopVolume()->GetInstanceID()];
-
   assert(hitID < fNumSensitive);
 
   if (hitID >= fHitsCollection->entries()) {
@@ -103,4 +99,4 @@ SimpleHit *SensitiveDetector::RetrieveAndSetupHit(G4TouchableHistory *aTouchable
   }
 
   return (*fHitsCollection)[hitID];
-}
+}

src/AdePTGeant4Integration.cpp

@@ -45,9 +45,10 @@ struct ScoringObjects {
   G4Step *fG4Step = new (&stepStorage) G4Step;
 
   G4TouchableHandle *fPreG4TouchableHistoryHandle =
-      ::new (toucheableHandleStorage) G4TouchableHandle{::new (toucheableHistoryStorage) G4TouchableHistory};
+      ::new (&toucheableHandleStorage[0]) G4TouchableHandle{::new (&toucheableHistoryStorage[0]) G4TouchableHistory};
+
   G4TouchableHandle *fPostG4TouchableHistoryHandle =
-      ::new (toucheableHandleStorage) G4TouchableHandle{::new (toucheableHistoryStorage + 1) G4TouchableHistory};
+      ::new (&toucheableHandleStorage[1]) G4TouchableHandle{::new (&toucheableHistoryStorage[1]) G4TouchableHistory};
 
   // We need the dynamic particle associated to the track to have the correct particle definition, however this can
   // only be set at construction time. Similarly, we can only set the dynamic particle for a track when creating it
@@ -348,6 +349,8 @@ void AdePTGeant4Integration::FillG4NavigationHistory(vecgeom::NavigationState aN
     const auto item = fglobal_vecgeom_to_g4_map.find(aNavState.At(aLevel)->id());
     pnewvol         = item == fglobal_vecgeom_to_g4_map.end() ? nullptr : const_cast<G4VPhysicalVolume *>(item->second);
 
+    if (!pnewvol) throw std::runtime_error("VecGeom volume not found in G4 mapping!");
+
     if (aG4HistoryDepth && (aLevel <= aG4HistoryDepth)) {
       pvol = aG4NavigationHistory->GetVolume(aLevel);
       // If they match we do not need to update the history at this level