Merge pull request #634 from ManuelHu/validation

tests: update nist and ge validation

Author: ManuelHu

docs/rmg-commands.md

@@ -11,14 +11,39 @@ electrons and compare them with data from the
 short description of the reported quantities is available
 [here](https://physics.nist.gov/PhysRefData/Star/Text/method.html)).
 
-Multiple Coulomb scattering is disabled through the
+The simulation physics is modified by the commands:
 
 ```
 /process/inactivate msc
+/RMG/Processes/DefaultProductionCut 1 km
+/RMG/Processes/SensitiveProductionCut 1 km
 ```
 
-command, since it prevents us to measure the true integrated particle path
-length in the simulation (unless a very short step size limit is set).
+- Multiple Coulomb scattering is disabled, since it prevents us to measure the
+  true integrated particle path length in the simulation (unless a very short
+  step size limit is set).
+- the high production cuts ensure that no secondaries are created. This is
+  required to correctly measure the bremsstrahlung component of the CSDA loss.
+
+An effect of the step length is visible. The larger the step length, the larger
+the difference to the expected integrated path length.
+
+The statistical uncertainties are dominated by energy-loss fluctuations. They
+are kept enabled here to show the effect of the user-chosen step length values.
+Especially the stopping power at low step sizes are the most affected.
+
+:::{note}
+
+The simulated stopping powers are derived from the energy deposition in the
+first step and its length. This leads to the stopping powers only matching with
+the ESTAR-provided values until a material-specific energy.
+
+Below that, the energy is fully deposited in one single small step, from which
+we cannot estimate the stopping power at the incident energy any more. For these
+low energies, we plot the incident energy over the full expected ranged from
+ESTAR data as the expectation value.
+
+:::
 
 ### Electron interactions in natural germanium
 

docs/validation/nist.md

@@ -44,6 +44,7 @@ class RMGSteppingAction : public G4UserSteppingAction {
     RMGRunAction* fRunAction;
 
     bool fSkipTracking = false;
+    bool fKillSecondaries = false;
     double fDaughterKillLifetime = -1;
 
     std::unique_ptr<G4GenericMessenger> fMessenger;

include/RMGSteppingAction.hh

@@ -32,7 +32,7 @@ RMGSteppingAction::RMGSteppingAction(RMGRunAction* run_action) : fRunAction(run_
 
 void RMGSteppingAction::UserSteppingAction(const G4Step* step) {
 
-  if (fSkipTracking) {
+  if (fSkipTracking || (fKillSecondaries && step->GetTrack()->GetParentID() > 0)) {
     step->GetTrack()->SetTrackStatus(fKillTrackAndSecondaries);
     return;
   }
@@ -114,6 +114,13 @@ void RMGSteppingAction::DefineCommands() {
       .SetParameterName("boolean", true)
       .SetDefaultValue("true")
       .SetStates(G4State_Idle);
+
+  fMessenger->DeclareProperty("KillSecondaries", fKillSecondaries)
+      .SetGuidance("Kill all secondary particles (tracks with parent ID > 0) at their first step.")
+      .SetGuidance(std::string("This is ") + (fKillSecondaries ? "enabled" : "disabled") + " by default")
+      .SetParameterName("boolean", true)
+      .SetDefaultValue("true")
+      .SetStates(G4State_Idle);
 }
 
 

src/RMGSteppingAction.cc

@@ -8,7 +8,16 @@
 # to the theoretical range
 /process/inactivate msc
 
+# disable energy loss fluctuations (straggling)
+/process/eLoss/fluct {FLUCT}
+
+# prevent any secondary production
+/RMG/Processes/DefaultProductionCut 1 km
+/RMG/Processes/SensitiveProductionCut 1 km
+
 /RMG/Output/Germanium/StoreTrackID true
+# store only post-steps so that we can correctly calculate step lengths
+/RMG/Output/Germanium/StepPositionMode PostStep
 
 /RMG/Generator/Confine Volume
 /RMG/Generator/Confinement/Geometrical/AddSolid Sphere

tests/nist/macros/run.mac

@@ -1,36 +1,69 @@
 from __future__ import annotations
 
 import os
+from multiprocessing import Pool
+from pathlib import Path
 
 import numpy as np
-from remage import remage_run_from_args
+from remage import remage_run
 
-n_events = 500 * int(os.environ["RMG_STATS_FACTOR"])
+Path("output/").mkdir(parents=True, exist_ok=True)
+
+n_proc = int(os.environ["RMG_STATS_FACTOR"])
+n_events = 500 * n_proc
+
+energies = np.geomspace(30, 50000, num=15).astype(int)  # in keV
+step_sizes = (1, 0.1, 0.01, 0.001, 0.0001)  # in m
+range_energies = np.array([300, 500, 1000, 1500, 2000, 3000, 4000])  # in keV
+range_step_sizes = (1,)  # in m
+
+sims = []
+
+for mat in ("ge", "ar", "cu"):
+    for step_size in step_sizes:  # in m
+        for energy in energies:  # in keV
+            sims.append(
+                {
+                    "macros": "macros/run.mac",
+                    "macro_substitutions": {
+                        "ENERGY": energy,
+                        "EVENTS": n_events,
+                        "MAX_STEP_SIZE_IN_M": step_size,
+                        "FLUCT": "true",
+                    },
+                    "gdml_files": f"gdml/geometry-{mat}.gdml",
+                    "output": f"output/electrons-{mat}-{energy}-keV-{step_size}-m.lh5",
+                    "overwrite_output": True,
+                    "flat_output": True,
+                    "log_level": "summary",
+                }
+            )
 
 for mat in ("ge", "ar", "cu"):
-    base_args = [
-        "--macro-substitutions",
-        "ENERGY={}",
-        f"EVENTS={n_events}",
-        "MAX_STEP_SIZE_IN_M=1",
-        "--gdml-files",
-        f"gdml/geometry-{mat}.gdml",
-        "--output-file",
-        f"electrons-{mat}-{{}}-keV.lh5",
-        "--flat-output",
-        "--overwrite",
-        "macros/run.mac",
-    ]
-
-    energies = np.concatenate(
-        [
-            np.arange(30, 100, step=20),
-            np.arange(100, 500, step=200),
-            np.arange(500, 5000, step=500),
-        ]
-    )
-
-    for energy in energies:  # in keV
-        args = [arg.format(energy) for arg in base_args]
-        print("remage", " ".join(args))
-        remage_run_from_args(args, raise_on_error=True)
+    for step_size in range_step_sizes:  # in m
+        for energy in range_energies:  # in keV
+            sims.append(
+                {
+                    "macros": "macros/run.mac",
+                    "macro_substitutions": {
+                        "ENERGY": energy,
+                        "EVENTS": n_events,
+                        "MAX_STEP_SIZE_IN_M": step_size,
+                        "FLUCT": "true",
+                    },
+                    "gdml_files": f"gdml/geometry-{mat}.gdml",
+                    "output": f"output/electrons-range-{mat}-{energy}-keV-{step_size}-m.lh5",
+                    "overwrite_output": True,
+                    "flat_output": True,
+                    "log_level": "summary",
+                }
+            )
+
+
+def _run(sim):
+    remage_run(**sim)
+
+
+if __name__ == "__main__":
+    with Pool(n_proc) as pool:
+        pool.map(_run, sims)