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,13 @@ 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) {
+    //if (numElectrons == previousElectrons && numPositrons == previousPositrons && numGammas == 0) {
       loopingNo++;
     } else {
       previousElectrons = numElectrons;
       previousPositrons = numPositrons;
+      previousGammas = numGammas;
       loopingNo         = 0;
     }
 

examples/Example21/electrons.cuh

@@ -108,7 +108,32 @@ static __device__ __forceinline__ void TransportElectrons(adept::TrackManager<Tr
     }
 
     // Call G4HepEm to compute the physics step limit.
-    G4HepEmElectronManager::HowFar(&g4HepEmData, &g4HepEmPars, &elTrack, &rnge);
+    //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 +282,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

@@ -18,11 +18,13 @@ public:
   SteppingAction(DetectorConstruction *aDetector, RunAction *aRunAction, TrackingAction *aTrackingAction);
   ~SteppingAction() override;
   void UserSteppingAction(const G4Step *step) override;
+  void SetNumSteps(int aNumSteps){fNumSteps = aNumSteps;}
 
 private:
   DetectorConstruction *fDetector;
   RunAction *fRunAction;
   TrackingAction *fTrackingAction;
+  int fNumSteps{0};
 };
 
 #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

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
@@ -106,7 +107,7 @@
 ## -----------------------------------------------------------------------------
 ## Optionally, set a constant magnetic filed:
 ## -----------------------------------------------------------------------------
-/example21/detector/setField 0 0 0 tesla
+/example21/detector/setField 0 0 3.8 tesla
 #/example21/detector/setField 0 0 3.8 tesla
 
 ##
@@ -144,4 +145,4 @@
 /example21/setSeed 1
 
 # run events with parametrised simulation
-/run/beamOn 128
+/run/beamOn 16

examples/Example21/src/ActionInitialisation.cc

@@ -66,10 +66,13 @@ void ActionInitialisation::Build() const
   SetUserAction(aTrackingAction);
 
 // 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
+//#if defined TEST
+  //if(fDoBenchmark)
+  //{
+    SteppingAction *aSteppingAction = new SteppingAction(fDetector, aRunAction, aTrackingAction);
+    SetUserAction(aSteppingAction);
+    aTrackingAction->setSteppingAction(aSteppingAction);
+  //}
+//#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

@@ -26,6 +26,14 @@ SteppingAction::~SteppingAction() {}
 
 void SteppingAction::UserSteppingAction(const G4Step *theStep)
 {
+  // 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);
+  }
+
   // Check if we moved to a new volume
   if (theStep->IsLastStepInVolume()) {
     G4VPhysicalVolume *nextVolume = theStep->GetTrack()->GetNextVolume();

examples/Example21/src/TrackingAction.cc

@@ -60,6 +60,8 @@ TrackingAction::TrackingAction(DetectorConstruction* aDetector) : G4UserTracking
 
 void TrackingAction::PreUserTrackingAction(const G4Track *aTrack) 
 {
+  // Reset step counter
+  fSteppingAction->SetNumSteps(0);
   #if defined TEST
     //For leptons, get the Run object associated to this thread and start the timer for this track, only if it is outside 
     //the GPU region
@@ -93,6 +95,8 @@ void TrackingAction::PreUserTrackingAction(const G4Track *aTrack)
 
 void TrackingAction::PostUserTrackingAction(const G4Track *aTrack)
 {
+  // Reset step counter
+  fSteppingAction->SetNumSteps(0);
   #if defined TEST
     //Get the Run object associated to this thread and end the timer for this track
     Run* currentRun = static_cast< Run* > ( G4RunManager::GetRunManager()->GetNonConstCurrentRun() );