Adding Example22

base/inc/AdePT/TestManager.h

@@ -11,8 +11,8 @@
 
 // The clock to use for measurements
 #define CLOCK std::chrono::system_clock
-// The precission with which we will provide the results
-#define PRECISSION std::chrono::microseconds
+// The precision with which we will provide the results
+#define PRECISION std::chrono::microseconds
 
 /**
  * @brief A type containing a timestamp, accumulated time, and whether the timer is running
@@ -75,8 +75,8 @@ class TestManager {
   /** @brief Returns the accumulated duration for the specified tag in seconds */
   double getDurationSeconds(TTag tag)
   {
-    return (double)(std::chrono::duration_cast<PRECISSION>(fTimers[tag].accumulatedDuration)).count() /
-           PRECISSION::period::den;
+    return (double)(std::chrono::duration_cast<PRECISION>(fTimers[tag].accumulatedDuration)).count() /
+           PRECISION::period::den;
   }
 
   /** @brief Sets the accumulator value for the specified tag */
@@ -96,8 +96,8 @@ class TestManager {
     }
   }
 
-  /** @brief Empties the timer map */
-  void reset() { fTimers.clear(); }
+  /** @brief Empties the maps */
+  void reset() { fTimers.clear(); fAccumulators.clear(); }
 
   /** @brief Removes a timer from the map */
   void removeTimer(TTag tag) { fTimers.erase(tag); }
@@ -143,7 +143,7 @@ class TestManager {
 
   /** @brief Export a CSV file with the timer names as labels and the accumulated time for each */
   // If the file is not empty, write only the times
-  void exportCSV()
+  void exportCSV(bool overwrite = true)
   { 
     std::string aOutputDir;
     std::string aOutputFilename;
@@ -158,10 +158,10 @@ class TestManager {
     bool first_write = !std::filesystem::exists(path);
 
     std::ofstream output_file;
-    output_file.open(path, std::ofstream::app);
+    output_file.open(path, overwrite ? std::ofstream::trunc : std::ofstream::app);
 
     // Write the header only the first time
-    if (first_write) {
+    if (first_write || overwrite) {
       // Timers
       for (auto iter = fTimers.begin(); iter != fTimers.end(); ++iter) {
         output_file << iter->first;

examples/CMakeLists.txt

@@ -23,6 +23,7 @@ add_subdirectory(Example18)
 add_subdirectory(Example19)
 add_subdirectory(Example20)
 add_subdirectory(Example21)
+add_subdirectory(Example22)
 add_subdirectory(TestEm3)
 add_subdirectory(TestEm3MT)
 add_subdirectory(TestEm3Compact)

examples/Example21/AdeptIntegration.cu

@@ -336,7 +336,7 @@ void AdeptIntegration::ShowerGPU(int event, TrackBuffer &buffer) // const &buffe
   int numLeaked         = 0;
   int num_compact       = 0;
   int loopingNo         = 0;
-  int previousElectrons = -1, previousPositrons = -1;
+  int previousElectrons = -1, previousPositrons = -1, previousGammas = -1;
   LeakedTracks leakedTracks = {.leakedElectrons = electrons.leakedTracks,
                             .leakedPositrons = positrons.leakedTracks,
                             .leakedGammas    = gammas.leakedTracks};
@@ -429,11 +429,12 @@ void AdeptIntegration::ShowerGPU(int event, TrackBuffer &buffer) // const &buffe
     if (fDebugLevel > 1) {
       printf("iter %d: elec %d, pos %d, gam %d, leak %d\n", niter++, numElectrons, numPositrons, numGammas, numLeaked);
     }
-    if (numElectrons == previousElectrons && numPositrons == previousPositrons && numGammas == 0) {
+    if (numElectrons == previousElectrons && numPositrons == previousPositrons && numGammas == previousGammas) {
       loopingNo++;
     } else {
       previousElectrons = numElectrons;
       previousPositrons = numPositrons;
+      previousGammas = numGammas;
       loopingNo         = 0;
     }
 

examples/Example21/electrons.cuh

@@ -108,7 +108,30 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
     }
 
     // Call G4HepEm to compute the physics step limit.
-    G4HepEmElectronManager::HowFar(&g4HepEmData, &g4HepEmPars, &elTrack, &rnge);
+    G4HepEmElectronManager::HowFarToDiscreteInteraction(&g4HepEmData, &g4HepEmPars, &elTrack);
+
+    bool restrictedPhysicalStepLength = false;
+    if (BzFieldValue != 0) {
+      const double momentumMag = sqrt(energy * (energy + 2.0 * Mass));
+      // Distance along the track direction to reach the maximum allowed error
+      const double safeLength = fieldPropagatorBz.ComputeSafeLength(momentumMag, Charge, dir);
+
+      constexpr int MaxSafeLength = 10;
+      double limit                = MaxSafeLength * safeLength;
+      limit                       = safety > limit ? safety : limit;
+
+      double physicalStepLength = elTrack.GetPStepLength();
+      if (physicalStepLength > limit) {
+        physicalStepLength           = limit;
+        restrictedPhysicalStepLength = true;
+        elTrack.SetPStepLength(physicalStepLength);
+        // Note: We are limiting the true step length, which is converted to
+        // a shorter geometry step length in HowFarToMSC. In that sense, the
+        // limit is an over-approximation, but that is fine for our purpose.
+      }
+    }
+
+    G4HepEmElectronManager::HowFarToMSC(&g4HepEmData, &g4HepEmPars, &elTrack, &rnge);
 
     // Remember MSC values for the next step(s).
     currentTrack.initialRange       = mscData->fInitialRange;
@@ -257,7 +280,7 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
         }
       }
       continue;
-    } else if (!propagated) {
+    } else if (!propagated || restrictedPhysicalStepLength) {
       // Did not yet reach the interaction point due to error in the magnetic
       // field propagation. Try again next time.
       survive();

examples/Example21/include/Run.hh

@@ -59,6 +59,7 @@ public:
     ECAL_SQ,
     EVENT_SUM,
     EVENT_SQ,
+    NUM_PARTICLES,
     NUM_ACCUMULATORS
   };
 

examples/Example21/include/SteppingAction.hh

@@ -15,14 +15,18 @@ class TrackingAction;
 class SteppingAction : public G4UserSteppingAction {
 
 public:
-  SteppingAction(DetectorConstruction *aDetector, RunAction *aRunAction, TrackingAction *aTrackingAction);
+  SteppingAction(DetectorConstruction *aDetector, RunAction *aRunAction, TrackingAction *aTrackingAction,
+                 bool aDoBenchmark);
   ~SteppingAction() override;
   void UserSteppingAction(const G4Step *step) override;
+  void SetNumSteps(int aNumSteps) { fNumSteps = aNumSteps; }
 
 private:
   DetectorConstruction *fDetector;
   RunAction *fRunAction;
   TrackingAction *fTrackingAction;
+  int fNumSteps{0};
+  bool fDoBenchmark{false};
 };
 
 #endif

examples/Example21/include/TrackingAction.hh

@@ -35,6 +35,7 @@
 #define TRACKINGACTION_HH
 
 #include "G4UserTrackingAction.hh"
+#include "SteppingAction.hh"
 #include "G4Region.hh"
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -58,13 +59,15 @@ public:
   inline void setCurrentRegion(G4Region* aCurrentRegion){fCurrentRegion = aCurrentRegion;}
   inline G4VPhysicalVolume* getCurrentVolume(){return fCurrentVolume;}
   inline void setCurrentVolume(G4VPhysicalVolume* aCurrentVolume){fCurrentVolume = aCurrentVolume;}
+  inline void setSteppingAction(SteppingAction* aSteppingAction){fSteppingAction = aSteppingAction;}
 
 private:
   DetectorConstruction* fDetector;
   bool fInsideEcal;
   G4Region* fCurrentRegion;
   G4VPhysicalVolume* fCurrentVolume;
   G4Region* fGPURegion;
+  SteppingAction* fSteppingAction;
 };
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

examples/Example21/macros/example21.mac.in

@@ -15,7 +15,8 @@
 ##
 /example21/setSeed 1
 
-/example21/detector/filename @CMS2018_GDML@
+#/example21/detector/filename @CMS2018_GDML@
+/example21/detector/filename cms2018_sd.gdml
 /example21/adept/verbose 0
 ## Threshold for buffering tracks before sending to GPU
 /example21/adept/threshold 200
@@ -144,5 +145,5 @@
 
 # run events with parametrised simulation
 # by default all created models are active
-/run/beamOn 1
+/run/beamOn 10
 

examples/Example21/macros/example21_ttbar.mac.in

@@ -7,15 +7,16 @@
 ## Geant4 macro for modelling simplified sampling calorimeters
 ## =============================================================================
 ##
-/run/numberOfThreads 16
+/run/numberOfThreads 8
 /control/verbose 0
 /run/verbose 0
 /process/verbose 0
 /tracking/verbose 0
 ##
 /example21/setSeed 1
 
-/example21/detector/filename @CMS2018_GDML@
+#/example21/detector/filename @CMS2018_GDML@
+/example21/detector/filename cms2018_sd.gdml
 /example21/detector/regionname EcalRegion
 /example21/adept/verbose 0
 ## Threshold for buffering tracks before sending to GPU
@@ -144,4 +145,4 @@
 /example21/setSeed 1
 
 # run events with parametrised simulation
-/run/beamOn 128
+/run/beamOn 16

examples/Example21/src/ActionInitialisation.cc

@@ -64,12 +64,8 @@ void ActionInitialisation::Build() const
   SetUserAction(aRunAction);
   TrackingAction *aTrackingAction = new TrackingAction(fDetector);
   SetUserAction(aTrackingAction);
+  SteppingAction *aSteppingAction = new SteppingAction(fDetector, aRunAction, aTrackingAction, fDoBenchmark);
+  SetUserAction(aSteppingAction);
+  aTrackingAction->setSteppingAction(aSteppingAction);
 
-// Do not register this if the TestManager is not active or if benchmark is not selected
-#if defined TEST
-  if(fDoBenchmark)
-  {
-    SetUserAction(new SteppingAction(fDetector, aRunAction, aTrackingAction));
-  }
-#endif
 }

examples/Example21/src/EventAction.cc

@@ -91,6 +91,7 @@ void EventAction::EndOfEventAction(const G4Event *aEvent)
     {
       aTestManager->timerStop(Run::timers::EVENT);
     }
+    aTestManager->addToAccumulator(Run::accumulators::NUM_PARTICLES, aEvent->GetPrimaryVertex()->GetNumberOfParticle());
   #endif
 
   // Get hits collection ID (only once)

examples/Example21/src/Run.cc

@@ -37,6 +37,8 @@ void Run::Merge(const G4Run *run)
     fTestManager->addToAccumulator(accumulators::ECAL_SUM,
                                         aTestManager->getAccumulator(accumulators::ECAL_SUM));
     fTestManager->addToAccumulator(accumulators::ECAL_SQ, aTestManager->getAccumulator(accumulators::ECAL_SQ));
+    fTestManager->addToAccumulator(accumulators::NUM_PARTICLES, 
+                                        aTestManager->getAccumulator(accumulators::NUM_PARTICLES));
   }
 
   G4Run::Merge(run);
@@ -65,12 +67,6 @@ void Run::EndOfRunSummary(G4String aOutputDirectory, G4String aOutputFilename,
     G4cout << "BENCHMARK: Run time: " << runTime << "\n";
     G4cout << "BENCHMARK: Mean Event time: " << eventMean << "\n";
     G4cout << "BENCHMARK: Event Standard Deviation: " << eventStdev << "\n";
-    G4cout << "BENCHMARK: Mean Non EM time: " << nonEMMean << "\n";
-    G4cout << "BENCHMARK: Non EM Standard Deviation: " << nonEMStdev << "\n";
-    G4cout << "BENCHMARK: Mean ECAL e-, e+ and gammas time: " << ecalMean << "\n";
-    G4cout << "BENCHMARK: ECAL e-, e+ and gammas Standard Deviation: " << ecalStdev << "\n";
-    G4cout << "BENCHMARK: Mean proportion of time spent simulating e-, e+ and gammas in ECAL: "
-          << 100 * ecalMean / eventMean << "%\n";
   }
 
   // Export the results per event
@@ -105,10 +101,19 @@ void Run::EndOfRunSummary(G4String aOutputDirectory, G4String aOutputFilename,
 
       aOutputTestManager.setOutputDirectory(aOutputDirectory);
       aOutputTestManager.setOutputFilename(aOutputFilename);
-      aOutputTestManager.exportCSV();
+      aOutputTestManager.exportCSV(false);
 
       aOutputTestManager.reset();
     }
   }
   TestManagerStore<int>::GetInstance()->Reset();
+
+  // Export global results
+
+  aOutputTestManager.setAccumulator("Totaltime", fTestManager->getDurationSeconds(timers::TOTAL));
+  aOutputTestManager.setAccumulator("NumParticles", fTestManager->getAccumulator(accumulators::NUM_PARTICLES));
+
+  aOutputTestManager.setOutputDirectory(aOutputDirectory);
+  aOutputTestManager.setOutputFilename(aOutputFilename + "_global");
+  aOutputTestManager.exportCSV();
 }

examples/Example21/src/SteppingAction.cc

@@ -17,57 +17,70 @@
 #include "Run.hh"
 #include "TestManager.h"
 
-SteppingAction::SteppingAction(DetectorConstruction *aDetector, RunAction *aRunAction, TrackingAction *aTrackingAction)
-    : G4UserSteppingAction(), fDetector(aDetector), fRunAction(aRunAction), fTrackingAction(aTrackingAction)
+SteppingAction::SteppingAction(DetectorConstruction *aDetector, RunAction *aRunAction, TrackingAction *aTrackingAction,
+                               bool aDoBenchmark)
+    : G4UserSteppingAction(), fDetector(aDetector), fRunAction(aRunAction), fTrackingAction(aTrackingAction),
+      fDoBenchmark(aDoBenchmark)
 {
 }
 
 SteppingAction::~SteppingAction() {}
 
 void SteppingAction::UserSteppingAction(const G4Step *theStep)
 {
-  // Check if we moved to a new volume
-  if (theStep->IsLastStepInVolume()) {
-    G4VPhysicalVolume *nextVolume = theStep->GetTrack()->GetNextVolume();
-    if (nextVolume != nullptr) {
-      if (!fTrackingAction->getInsideEcal()) {
-        // Check if the new volume is in the EM calorimeter region
-        if (nextVolume->GetLogicalVolume()->GetRegion() == fTrackingAction->getGPURegion()) {
-          // If it is, stop the timer for this track and store the result
-          G4Track *aTrack = theStep->GetTrack();
+  // Kill the particle if it has done over 1000 steps
+  fNumSteps++;
+  if (fNumSteps > 10000) {
+    G4cout << "Warning: Killing track over 10000 steps" << G4endl;
+    theStep->GetTrack()->SetTrackStatus(fStopAndKill);
+  }
+
+#if defined TEST
+  if (fDoBenchmark) {
+    // Check if we moved to a new volume
+    if (theStep->IsLastStepInVolume()) {
+      G4VPhysicalVolume *nextVolume = theStep->GetTrack()->GetNextVolume();
+      if (nextVolume != nullptr) {
+        if (!fTrackingAction->getInsideEcal()) {
+          // Check if the new volume is in the EM calorimeter region
+          if (nextVolume->GetLogicalVolume()->GetRegion() == fTrackingAction->getGPURegion()) {
+            // If it is, stop the timer for this track and store the result
+            G4Track *aTrack = theStep->GetTrack();
 
-          // Make sure this will only run on the first step when we enter the ECAL
-          fTrackingAction->setInsideEcal(true);
+            // Make sure this will only run on the first step when we enter the ECAL
+            fTrackingAction->setInsideEcal(true);
 
-          // We are only interested in the processing of e-, e+ and gammas in the EM calorimeter, for other
-          // particles the time keeps running
-          if (aTrack->GetDefinition() == G4Gamma::Gamma() || aTrack->GetDefinition() == G4Electron::Electron() ||
-              aTrack->GetDefinition() == G4Positron::Positron()) {
-            // Get the Run object associated to this thread and end the timer for this track
-            Run *currentRun        = static_cast<Run *>(G4RunManager::GetRunManager()->GetNonConstCurrentRun());
-            auto aTestManager = currentRun->GetTestManager();
+            // We are only interested in the processing of e-, e+ and gammas in the EM calorimeter, for other
+            // particles the time keeps running
+            if (aTrack->GetDefinition() == G4Gamma::Gamma() || aTrack->GetDefinition() == G4Electron::Electron() ||
+                aTrack->GetDefinition() == G4Positron::Positron()) {
+              // Get the Run object associated to this thread and end the timer for this track
+              Run *currentRun   = static_cast<Run *>(G4RunManager::GetRunManager()->GetNonConstCurrentRun());
+              auto aTestManager = currentRun->GetTestManager();
 
-            aTestManager->timerStop(Run::timers::NONEM);
-            aTestManager->addToAccumulator(Run::accumulators::NONEM_EVT,
-                                                aTestManager->getDurationSeconds(Run::timers::NONEM));
-            aTestManager->removeTimer(Run::timers::NONEM);
+              aTestManager->timerStop(Run::timers::NONEM);
+              aTestManager->addToAccumulator(Run::accumulators::NONEM_EVT,
+                                             aTestManager->getDurationSeconds(Run::timers::NONEM));
+              aTestManager->removeTimer(Run::timers::NONEM);
+            }
           }
-        }
-      } else {
-        // In case this track is exiting the EM calorimeter, start the timer
-        if (nextVolume->GetLogicalVolume()->GetRegion() != fTrackingAction->getGPURegion()) {
-          G4Track *aTrack = theStep->GetTrack();
+        } else {
+          // In case this track is exiting the EM calorimeter, start the timer
+          if (nextVolume->GetLogicalVolume()->GetRegion() != fTrackingAction->getGPURegion()) {
+            G4Track *aTrack = theStep->GetTrack();
 
-          fTrackingAction->setInsideEcal(false);
+            fTrackingAction->setInsideEcal(false);
 
-          if (aTrack->GetDefinition() == G4Gamma::Gamma() || aTrack->GetDefinition() == G4Electron::Electron() ||
-              aTrack->GetDefinition() == G4Positron::Positron()) {
-            // Get the Run object associated to this thread and start the timer for this track
-            Run *currentRun = static_cast<Run *>(G4RunManager::GetRunManager()->GetNonConstCurrentRun());
-            currentRun->GetTestManager()->timerStart(Run::timers::NONEM);
+            if (aTrack->GetDefinition() == G4Gamma::Gamma() || aTrack->GetDefinition() == G4Electron::Electron() ||
+                aTrack->GetDefinition() == G4Positron::Positron()) {
+              // Get the Run object associated to this thread and start the timer for this track
+              Run *currentRun = static_cast<Run *>(G4RunManager::GetRunManager()->GetNonConstCurrentRun());
+              currentRun->GetTestManager()->timerStart(Run::timers::NONEM);
+            }
           }
         }
       }
     }
   }
+#endif
 }