Enhance DepositedChargeActor with uncertainty statistics

core/opengate_core/opengate_lib/GateDepositedChargeActor.cpp

@@ -16,13 +16,18 @@ GateDepositedChargeActor::GateDepositedChargeActor(py::dict &user_info)
   // Geant4 callbacks need by this actor.
   fActions.insert("StartSimulationAction");
   fActions.insert("BeginOfRunAction");
+  fActions.insert("BeginOfEventAction");
   fActions.insert("PreUserTrackingAction");
   fActions.insert("PostUserTrackingAction");
+  fActions.insert("EndOfEventAction");
   fActions.insert("EndOfSimulationWorkerAction");
 
-  // Initialize the total deposited charge to zero
+  // Initialize the merged accumulators to zero.
   fDepositedNominalCharge = 0.0;
   fDepositedDynamicCharge = 0.0;
+  fDepositedNominalChargeSquared = 0.0;
+  fDepositedDynamicChargeSquared = 0.0;
+  fNumberOfEvents = 0;
 }
 
 void GateDepositedChargeActor::InitializeUserInfo(py::dict &user_info) {
@@ -31,20 +36,35 @@ void GateDepositedChargeActor::InitializeUserInfo(py::dict &user_info) {
 }
 
 void GateDepositedChargeActor::StartSimulationAction() {
-  // Reset the total deposited charges at the start of the simulation
+  // Reset the merged accumulators at the start of the simulation
   fDepositedNominalCharge = 0.0;
   fDepositedDynamicCharge = 0.0;
+  fDepositedNominalChargeSquared = 0.0;
+  fDepositedDynamicChargeSquared = 0.0;
+  fNumberOfEvents = 0;
 }
 
 void GateDepositedChargeActor::BeginOfRunAction(const G4Run *run) {
-  // Reset the thread-local charge accumulators at the beginning of the first
-  // run.
+  // Reset the thread-local accumulators at the beginning of the first run.
   if (run->GetRunID() == 0) {
-    threadLocalData.Get().fNominalCharge = 0.0;
-    threadLocalData.Get().fDynamicCharge = 0.0;
+    auto &data = threadLocalData.Get();
+    data.fEventNominalCharge = 0.0;
+    data.fEventDynamicCharge = 0.0;
+    data.fSumNominalCharge = 0.0;
+    data.fSumNominalChargeSquared = 0.0;
+    data.fSumDynamicCharge = 0.0;
+    data.fSumDynamicChargeSquared = 0.0;
+    data.fNumberOfEvents = 0;
   }
 }
 
+void GateDepositedChargeActor::BeginOfEventAction(const G4Event * /*event*/) {
+  // Reset the per-event charge buffer.
+  auto &data = threadLocalData.Get();
+  data.fEventNominalCharge = 0.0;
+  data.fEventDynamicCharge = 0.0;
+}
+
 void GateDepositedChargeActor::PreUserTrackingAction(const G4Track *track) {
   // Net deposited charge =
   //     + charge of particles dead in the attached volume
@@ -71,8 +91,10 @@ void GateDepositedChargeActor::PreUserTrackingAction(const G4Track *track) {
   q_dynamic *= weight;
 
   auto &data = threadLocalData.Get();
-  data.fNominalCharge -= q_nominal; // being born -> subtract nominal charge
-  data.fDynamicCharge -= q_dynamic; // being born -> subtract dynamic charge
+  data.fEventNominalCharge -=
+      q_nominal; // being born -> subtract nominal charge
+  data.fEventDynamicCharge -=
+      q_dynamic; // being born -> subtract dynamic charge
 }
 
 void GateDepositedChargeActor::PostUserTrackingAction(const G4Track *track) {
@@ -97,16 +119,32 @@ void GateDepositedChargeActor::PostUserTrackingAction(const G4Track *track) {
   q_dynamic *= weight;
 
   auto &data = threadLocalData.Get();
-  data.fNominalCharge += q_nominal; // dying -> add nominal charge
-  data.fDynamicCharge += q_dynamic; // dying -> add dynamic charge
+  data.fEventNominalCharge += q_nominal; // dying -> add nominal charge
+  data.fEventDynamicCharge += q_dynamic; // dying -> add dynamic charge
+}
+
+void GateDepositedChargeActor::EndOfEventAction(const G4Event * /*event*/) {
+  // Fold the per-event net charge into the running first and second moments.
+  auto &data = threadLocalData.Get();
+
+  const double xn = data.fEventNominalCharge;
+  const double xd = data.fEventDynamicCharge;
+
+  data.fSumNominalCharge += xn;
+  data.fSumNominalChargeSquared += xn * xn;
+  data.fSumDynamicCharge += xd;
+  data.fSumDynamicChargeSquared += xd * xd;
+  data.fNumberOfEvents += 1;
 }
 
 void GateDepositedChargeActor::EndOfSimulationWorkerAction(
     const G4Run * /*lastRun*/) {
-  // Accumulate the thread-local charges into the total deposited charge
-  G4AutoLock mutex(
-      &GateDepositedChargeActorMutex); // Lock the mutex to protect access to
-                                       // the total deposited charge
-  fDepositedNominalCharge += threadLocalData.Get().fNominalCharge;
-  fDepositedDynamicCharge += threadLocalData.Get().fDynamicCharge;
+  // Accumulate the thread-local moments into the merged accumulators.
+  G4AutoLock mutex(&GateDepositedChargeActorMutex);
+  auto &data = threadLocalData.Get();
+  fDepositedNominalCharge += data.fSumNominalCharge;
+  fDepositedDynamicCharge += data.fSumDynamicCharge;
+  fDepositedNominalChargeSquared += data.fSumNominalChargeSquared;
+  fDepositedDynamicChargeSquared += data.fSumDynamicChargeSquared;
+  fNumberOfEvents += data.fNumberOfEvents;
 }

core/opengate_core/opengate_lib/GateDepositedChargeActor.h

@@ -25,25 +25,54 @@ class GateDepositedChargeActor : public GateVActor {
 
   void BeginOfRunAction(const G4Run *run) override;
 
+  void BeginOfEventAction(const G4Event *event) override;
+
   void PreUserTrackingAction(const G4Track *track) override;
 
   void PostUserTrackingAction(const G4Track *track) override;
 
+  void EndOfEventAction(const G4Event *event) override;
+
   void EndOfSimulationWorkerAction(const G4Run *lastRun) override;
 
+  // Net deposited charge, summed over all events
   double GetDepositedNominalCharge() const { return fDepositedNominalCharge; }
   double GetDepositedDynamicCharge() const { return fDepositedDynamicCharge; }
 
+  // Sum of the squared per-event net charge
+  double GetDepositedNominalChargeSquared() const {
+    return fDepositedNominalChargeSquared;
+  }
+  double GetDepositedDynamicChargeSquared() const {
+    return fDepositedDynamicChargeSquared;
+  }
+
+  // Number of primary events (histories) scored
+  long long GetNumberOfEvents() const { return fNumberOfEvents; }
+
 protected:
   struct threadLocal_t {
-    double fNominalCharge = 0.0; // nominal charge of the particle definition
-    double fDynamicCharge = 0.0; // effective charge due to ionisation
+    // Net charge accumulated during the current event only
+    double fEventNominalCharge = 0.0;
+    double fEventDynamicCharge = 0.0;
+
+    // Running first and second moments over the events scored by this worker.
+    double fSumNominalCharge = 0.0;
+    double fSumNominalChargeSquared = 0.0;
+    double fSumDynamicCharge = 0.0;
+    double fSumDynamicChargeSquared = 0.0;
+
+    // Number of events (histories) scored by this worker.
+    long long fNumberOfEvents = 0;
   };
   G4Cache<threadLocal_t> threadLocalData;
 
-  // Merged net deposited charge in units of eplus
-  double fDepositedNominalCharge;
-  double fDepositedDynamicCharge;
+  // Merged over all workers.
+  double fDepositedNominalCharge;        // Sum x
+  double fDepositedDynamicCharge;        // Sum x
+  double fDepositedNominalChargeSquared; // Sum x^2
+  double fDepositedDynamicChargeSquared; // Sum x^2
+  long long fNumberOfEvents;             // N
 };
 
 #endif // GateDepositedChargeActor_h

core/opengate_core/opengate_lib/pyGateDepositedChargeActor.cpp

@@ -16,5 +16,10 @@ void init_GateDepositedChargeActor(py::module &m) {
       .def("GetDepositedNominalCharge",
            &GateDepositedChargeActor::GetDepositedNominalCharge)
       .def("GetDepositedDynamicCharge",
-           &GateDepositedChargeActor::GetDepositedDynamicCharge);
+           &GateDepositedChargeActor::GetDepositedDynamicCharge)
+      .def("GetDepositedNominalChargeSquared",
+           &GateDepositedChargeActor::GetDepositedNominalChargeSquared)
+      .def("GetDepositedDynamicChargeSquared",
+           &GateDepositedChargeActor::GetDepositedDynamicChargeSquared)
+      .def("GetNumberOfEvents", &GateDepositedChargeActor::GetNumberOfEvents);
 }

docs/source/user_guide/user_guide_reference_actors.rst

@@ -387,6 +387,34 @@ for examples covering stopping electrons, stopping positrons, neutral beams, tra
 
    For repeated volume placements, the actor will score the total charge deposited across all instances of the volume.
 
+Statistical uncertainty
+~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The actor also provides a history-by-history estimate of the uncertainty on the deposited charge.
+
+After the simulation, the statistics are exposed as two dictionaries, ``nominal_charge_statistics`` and ``dynamic_charge_statistics``, one for each accumulated quantity. Each contains:
+
+  - ``mean``: mean net deposited charge per event, ``Sum x / N``.
+  - ``std``: sample standard deviation of the per-event charge (Bessel-corrected).
+  - ``sem``: standard error of the mean, ``std / sqrt(N)``.
+  - ``total``: total net charge, ``Sum x`` (identical to ``deposited_nominal_charge`` / ``deposited_dynamic_charge``).
+  - ``total_uncertainty``: absolute uncertainty on ``total``, i.e. the standard deviation of the sum, ``sqrt(N) * std = N * sem``.
+  - ``relative_uncertainty``: ``total_uncertainty / |total|``.
+
+.. code-block:: python
+
+    charge = sim.add_actor("DepositedChargeActor", "charge")
+    charge.attached_to = target.name
+
+    sim.run()
+
+    stats = charge.nominal_charge_statistics
+    print(f"total charge: {stats['total']} +/- {stats['total_uncertainty']} e")
+    print(f"relative uncertainty: {stats['relative_uncertainty']:.3%}")
+
+The test ``test099_deposited_charge_actor_uncertainty.py`` cross-checks it against an explicit batch method.
+
+
 Reference
 ~~~~~~~~~
 .. autoclass:: opengate.actors.miscactors.DepositedChargeActor

opengate/actors/miscactors.py

@@ -533,8 +533,13 @@ class DepositedChargeActor(ActorBase, g4.GateDepositedChargeActor):
 
     def __init__(self, *args, **kwargs):
         ActorBase.__init__(self, *args, **kwargs)
+        # Total net deposited charge (sum over events), in eplus.
         self.deposited_nominal_charge = 0.0
         self.deposited_dynamic_charge = 0.0
+        # History-by-history statistics, filled at EndSimulationAction.
+        self.deposited_nominal_charge_squared = 0.0
+        self.deposited_dynamic_charge_squared = 0.0
+        self.number_of_events = 0
         self.__initcpp__()
 
     def __initcpp__(self):
@@ -544,8 +549,10 @@ def __initcpp__(self):
                 "StartSimulationAction",
                 "EndSimulationAction",
                 "BeginOfRunAction",
+                "BeginOfEventAction",
                 "PreUserTrackingAction",
                 "PostUserTrackingAction",
+                "EndOfEventAction",
                 "EndOfSimulationWorkerAction",
             }
         )
@@ -558,12 +565,89 @@ def initialize(self):
     def EndSimulationAction(self):
         self.deposited_nominal_charge = self.GetDepositedNominalCharge()
         self.deposited_dynamic_charge = self.GetDepositedDynamicCharge()
+        self.deposited_nominal_charge_squared = self.GetDepositedNominalChargeSquared()
+        self.deposited_dynamic_charge_squared = self.GetDepositedDynamicChargeSquared()
+        self.number_of_events = self.GetNumberOfEvents()
+
+    def recover_user_output(self, actor):
+        # The results are stored as plain scalar attributes (not in user_output),
+        # so when the simulation runs in a subprocess (start_new_process=True)
+        # they must be copied back explicitly from the subprocess actor.
+        super().recover_user_output(actor)
+        self.deposited_nominal_charge = actor.deposited_nominal_charge
+        self.deposited_dynamic_charge = actor.deposited_dynamic_charge
+        self.deposited_nominal_charge_squared = actor.deposited_nominal_charge_squared
+        self.deposited_dynamic_charge_squared = actor.deposited_dynamic_charge_squared
+        self.number_of_events = actor.number_of_events
+
+    @staticmethod
+    def _history_statistics(sum_x, sum_x2, n):
+        """History-by-history statistics from the first and second moments.
+
+        Returns a dict with, for the per-event charge distribution:
+            - ``mean``: mean net charge per event (sum_x / n)
+            - ``std``: sample standard deviation of the per-event charge (Bessel-corrected)
+            - ``sem``: standard error of the mean (std / sqrt(n))
+            - ``total``: total net charge (sum_x)
+            - ``total_uncertainty``: absolute uncertainty on ``total``
+            - ``relative_uncertainty``: total_uncertainty / |total|
+        """
+        stats = {
+            "mean": 0.0,
+            "std": 0.0,
+            "sem": 0.0,
+            "total": sum_x,
+            "total_uncertainty": 0.0,
+            "relative_uncertainty": 0.0,
+        }
+        if n < 2:
+            return stats
+        mean = sum_x / n
+        # Sum of squared deviations = sum_x2 - (sum_x)^2 / n, clamped to >= 0
+        sum_sq_dev = max(sum_x2 - sum_x * sum_x / n, 0.0)
+        variance = sum_sq_dev / (n - 1)
+        std = variance**0.5
+        sem = std / n**0.5
+        total_uncertainty = n * sem
+        stats.update(
+            mean=mean,
+            std=std,
+            sem=sem,
+            total_uncertainty=total_uncertainty,
+            relative_uncertainty=(
+                total_uncertainty / abs(sum_x) if sum_x != 0.0 else 0.0
+            ),
+        )
+        return stats
+
+    @property
+    def nominal_charge_statistics(self):
+        """History-by-history statistics for the nominal deposited charge."""
+        return self._history_statistics(
+            self.deposited_nominal_charge,
+            self.deposited_nominal_charge_squared,
+            self.number_of_events,
+        )
+
+    @property
+    def dynamic_charge_statistics(self):
+        """History-by-history statistics for the dynamic deposited charge."""
+        return self._history_statistics(
+            self.deposited_dynamic_charge,
+            self.deposited_dynamic_charge_squared,
+            self.number_of_events,
+        )
 
     def __str__(self):
+        nominal = self.nominal_charge_statistics
+        dynamic = self.dynamic_charge_statistics
         return (
             f"DepositedChargeActor {self.name}:"
-            f"  Nominal: {self.deposited_nominal_charge} e"
-            f"  Dynamic: {self.deposited_dynamic_charge} e"
+            f"  Nominal: {self.deposited_nominal_charge} "
+            f"+/- {nominal['total_uncertainty']} e"
+            f"  Dynamic: {self.deposited_dynamic_charge} "
+            f"+/- {dynamic['total_uncertainty']} e"
+            f"  (N events: {self.number_of_events})"
         )