Parallel Geant4 scoring and NIEL clamp

Common/SimConfig/CMakeLists.txt

@@ -21,6 +21,7 @@ o2_add_library(SimConfig
                        src/InteractionDiamondParam.cxx
                        src/GlobalProcessCutSimParam.cxx
 		       src/FluenceWeightCalculator.cxx
+                       src/G4ScoringMerger.cxx
                PUBLIC_LINK_LIBRARIES O2::CommonUtils
                                      O2::DetectorsCommonDataFormats O2::SimulationDataFormat
                                      FairRoot::Base Boost::program_options)

Common/SimConfig/include/SimConfig/FluenceWeightCalculator.h

@@ -31,5 +31,6 @@ class FluenceWeightCalculator
   static std::unique_ptr<TGraph> neutronG;
   static std::unique_ptr<TGraph> protonG;
   static std::unique_ptr<TGraph> pionG;
+  static std::unique_ptr<TGraph> electronG;
 };
 #endif

Common/SimConfig/include/SimConfig/G4ScoringMerger.h

@@ -0,0 +1,30 @@
+// Copyright 2019-2026 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+#ifndef O2_SIMCONFIG_G4SCORINGMERGER_H
+#define O2_SIMCONFIG_G4SCORINGMERGER_H
+
+#include <string>
+
+namespace o2::conf
+{
+
+/// Name of the Geant4 scoring dump written by one simulation worker
+std::string g4ScoringWorkerFileName(const std::string& meshName, int pid);
+
+/// Sum the per-worker Geant4 scoring dumps <mesh>.worker<pid>.txt in a directory into <mesh>.txt.
+/// If expectedWorkers > 0, each mesh must have exactly that many dumps.
+/// Returns the number of merged meshes, or -1 if the worker files are inconsistent.
+int mergeG4ScoringDumps(const std::string& directory, int expectedWorkers = 0);
+
+} // namespace o2::conf
+
+#endif

Common/SimConfig/src/FluenceWeightCalculator.cxx

@@ -11,13 +11,28 @@
 
 #include "SimConfig/FluenceWeightCalculator.h"
 #include <TFile.h>
+#include <algorithm>
 #include <fstream>
 #include <sstream>
 #include <iostream>
 
 std::unique_ptr<TGraph> FluenceWeightCalculator::neutronG;
 std::unique_ptr<TGraph> FluenceWeightCalculator::protonG;
 std::unique_ptr<TGraph> FluenceWeightCalculator::pionG;
+std::unique_ptr<TGraph> FluenceWeightCalculator::electronG;
+
+namespace
+{
+// Damage weight at an energy clamped to the tabulated range
+double evalClamped(const TGraph& g, double kineticEnergy)
+{
+  if (g.GetN() == 0) {
+    return 0.;
+  }
+  const double e = std::clamp(kineticEnergy, g.GetX()[0], g.GetX()[g.GetN() - 1]);
+  return g.Eval(e, nullptr, "S");
+}
+} // namespace
 
 double FluenceWeightCalculator::GetWeight(const int pdg, const double kineticEnergy)
 {
@@ -27,19 +42,22 @@ double FluenceWeightCalculator::GetWeight(const int pdg, const double kineticEne
     std::cerr << "FluenceWeightCalculator not initialized\n";
     return 0.;
   }
-  switch (std::abs(pdg)) {
-    case 2112: {
-      return neutronG->Eval(kineticEnergy, nullptr, "S");
-    }
-    case 2212: {
-      return ((kineticEnergy > 1e-3) ? protonG->Eval(kineticEnergy, nullptr, "S") : 0.);
-    }
-    case 211: {
-      return ((kineticEnergy > 10.) ? pionG->Eval(kineticEnergy, nullptr, "S") : 0.);
-    }
-    default:
-      return 0.0;
+  const int apdg = std::abs(pdg);
+  if (pdg == 2112) {
+    return evalClamped(*neutronG, kineticEnergy);
+  }
+  if (apdg == 11) {
+    return electronG ? evalClamped(*electronG, kineticEnergy) : 0.;
+  }
+  // other (anti)baryons use the proton weights
+  if (apdg >= 1000 && apdg < 10000) {
+    return ((kineticEnergy > 1e-3) ? evalClamped(*protonG, kineticEnergy) : 0.);
+  }
+  // mesons use the pion weights
+  if (apdg >= 100 && apdg < 1000) {
+    return ((kineticEnergy > 10.) ? evalClamped(*pionG, kineticEnergy) : 0.);
   }
+  return 0.;
 }
 
 void FluenceWeightCalculator::InitWeights(const std::string& filename)
@@ -74,6 +92,13 @@ void FluenceWeightCalculator::InitWeights(const std::string& filename)
     return;
   }
   pionG->SetBit(TGraph::kIsSortedX);
+  // electron weights are optional
+  tmp = nullptr;
+  inFile.GetObject("electronDW", tmp);
+  electronG.reset(tmp ? static_cast<TGraph*>(tmp->Clone()) : nullptr);
+  if (electronG) {
+    electronG->SetBit(TGraph::kIsSortedX);
+  }
 }
 
 void FluenceWeightCalculator::InitWeightsFromCSV(const std::string& filename)
@@ -89,6 +114,9 @@ void FluenceWeightCalculator::InitWeightsFromCSV(const std::string& filename)
   pionG = std::make_unique<TGraph>();
   pionG->SetName("pionDW");
   auto pioN = 0;
+  electronG = std::make_unique<TGraph>();
+  electronG->SetName("electronDW");
+  auto eleN = 0;
 
   std::ifstream in(filename);
   if (!in.is_open()) {
@@ -127,12 +155,21 @@ void FluenceWeightCalculator::InitWeightsFromCSV(const std::string& filename)
         pionG->SetPoint(pioN++, e, w);
         break;
       }
+      case 11: {
+        electronG->SetPoint(eleN++, e, w);
+        break;
+      }
       default:;
     }
   }
+  neutronG->Sort();
+  protonG->Sort();
+  pionG->Sort();
+  electronG->Sort();
   auto fout = new TFile("rd50_niel.root", "recreate");
   neutronG->Write();
   protonG->Write();
   pionG->Write();
+  electronG->Write();
   fout->Close();
 }

Common/SimConfig/src/G4ScoringMerger.cxx

@@ -0,0 +1,156 @@
+// Copyright 2019-2026 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+#include "SimConfig/G4ScoringMerger.h"
+#include <fairlogger/Logger.h>
+#include <filesystem>
+#include <fstream>
+#include <iomanip>
+#include <map>
+#include <regex>
+#include <sstream>
+#include <vector>
+
+namespace o2::conf
+{
+
+namespace
+{
+// One scorer block of a Geant4 mesh dump: its header lines and the summed rows
+struct ScorerBlock {
+  std::vector<std::string> header;
+  std::vector<std::string> keys; // "iZ,iPHI,iR" in file order
+  std::vector<double> sum;
+  std::vector<double> sum2;
+  std::vector<long> entries;
+};
+
+// Read one mesh dump into scorer blocks; returns false on a format error
+bool readDump(const std::string& fileName, std::vector<std::string>& meshHeader, std::vector<ScorerBlock>& blocks)
+{
+  std::ifstream in(fileName);
+  if (!in) {
+    return false;
+  }
+  std::string line;
+  ScorerBlock* current = nullptr;
+  while (std::getline(in, line)) {
+    if (line.rfind("# mesh name", 0) == 0) {
+      meshHeader.push_back(line);
+    } else if (line.rfind("# primitive scorer name", 0) == 0) {
+      blocks.emplace_back();
+      current = &blocks.back();
+      current->header.push_back(line);
+    } else if (line.rfind("#", 0) == 0) {
+      if (!current) {
+        return false;
+      }
+      current->header.push_back(line);
+    } else if (!line.empty()) {
+      if (!current) {
+        return false;
+      }
+      // iZ, iPHI, iR, total, total^2, entries
+      std::vector<std::string> fields;
+      std::stringstream ss(line);
+      std::string field;
+      while (std::getline(ss, field, ',')) {
+        fields.push_back(field);
+      }
+      if (fields.size() != 6) {
+        return false;
+      }
+      current->keys.push_back(fields[0] + "," + fields[1] + "," + fields[2]);
+      current->sum.push_back(std::stod(fields[3]));
+      current->sum2.push_back(std::stod(fields[4]));
+      current->entries.push_back(std::stol(fields[5]));
+    }
+  }
+  return !blocks.empty();
+}
+} // namespace
+
+std::string g4ScoringWorkerFileName(const std::string& meshName, int pid)
+{
+  return meshName + ".worker" + std::to_string(pid) + ".txt";
+}
+
+int mergeG4ScoringDumps(const std::string& directory, int expectedWorkers)
+{
+  namespace fs = std::filesystem;
+  const std::regex pattern(R"((.+)\.worker([0-9]+)\.txt)");
+  std::map<std::string, std::vector<fs::path>> filesPerMesh;
+  for (auto& entry : fs::directory_iterator(directory)) {
+    std::smatch match;
+    const auto name = entry.path().filename().string();
+    if (entry.is_regular_file() && std::regex_match(name, match, pattern)) {
+      filesPerMesh[match[1]].push_back(entry.path());
+    }
+  }
+
+  int merged = 0;
+  for (auto& [mesh, files] : filesPerMesh) {
+    if (expectedWorkers > 0 && static_cast<int>(files.size()) != expectedWorkers) {
+      LOG(error) << "Found " << files.size() << " Geant4 scoring dumps for mesh " << mesh << " but expected " << expectedWorkers;
+      return -1;
+    }
+    std::vector<std::string> meshHeader;
+    std::vector<ScorerBlock> total;
+    for (auto& file : files) {
+      std::vector<std::string> header;
+      std::vector<ScorerBlock> blocks;
+      if (!readDump(file.string(), header, blocks)) {
+        LOG(error) << "Cannot read Geant4 scoring dump " << file;
+        return -1;
+      }
+      if (total.empty()) {
+        meshHeader = header;
+        total = std::move(blocks);
+        continue;
+      }
+      if (blocks.size() != total.size()) {
+        LOG(error) << "Geant4 scoring dump " << file << " has a different set of scorers";
+        return -1;
+      }
+      for (size_t b = 0; b < blocks.size(); ++b) {
+        if (blocks[b].header != total[b].header || blocks[b].keys != total[b].keys) {
+          LOG(error) << "Geant4 scoring dump " << file << " does not match the mesh layout of the other workers";
+          return -1;
+        }
+        for (size_t i = 0; i < blocks[b].keys.size(); ++i) {
+          total[b].sum[i] += blocks[b].sum[i];
+          total[b].sum2[i] += blocks[b].sum2[i];
+          total[b].entries[i] += blocks[b].entries[i];
+        }
+      }
+    }
+
+    const auto outName = (fs::path(directory) / (mesh + ".txt")).string();
+    std::ofstream out(outName);
+    out << std::setprecision(16);
+    for (auto& line : meshHeader) {
+      out << line << "\n";
+    }
+    for (auto& block : total) {
+      for (auto& line : block.header) {
+        out << line << "\n";
+      }
+      for (size_t i = 0; i < block.keys.size(); ++i) {
+        out << block.keys[i] << "," << block.sum[i] << "," << block.sum2[i] << "," << block.entries[i] << "\n";
+      }
+    }
+    LOG(info) << "Merged " << files.size() << " Geant4 scoring dumps into " << outName;
+    ++merged;
+  }
+  return merged;
+}
+
+} // namespace o2::conf

Detectors/gconfig/g4Config.C

@@ -61,6 +61,10 @@ R__LOAD_LIBRARY(libgeant4vmc)
 #include "TG4RunConfiguration.h"
 #include "SimConfig/G4Params.h"
 #include "SimConfig/FluenceWeightCalculator.h"
+#include "SimConfig/G4ScoringMerger.h"
+#include "G4ScoringManager.hh"
+#include "G4VScoringMesh.hh"
+#include <unistd.h>
 #include "FastSim/G4FastSimulation.h"
 #endif
 #include "commonConfig.C"
@@ -159,16 +163,30 @@ void Config()
   std::cout << "g4Config.C finished" << std::endl;
 }
 
+// Write each Geant4 scoring mesh to a file named after this process, so that parallel workers do not overwrite each other
+void dumpScoringMeshesPerWorker()
+{
+  auto scoringManager = G4ScoringManager::GetScoringManagerIfExist();
+  if (!scoringManager) {
+    return;
+  }
+  for (size_t i = 0; i < scoringManager->GetNumberOfMesh(); ++i) {
+    const auto meshName = scoringManager->GetMesh(i)->GetWorldName();
+    scoringManager->DumpAllQuantitiesToFile(meshName, o2::conf::g4ScoringWorkerFileName(meshName, getpid()));
+  }
+}
+
 void Terminate()
 {
   static bool terminated = false;
   if (!terminated) {
+    terminated = true;
     std::cout << "Executing G4 terminate\n";
     TGeant4* geant4 = dynamic_cast<TGeant4*>(TVirtualMC::GetMC());
     if (geant4) {
+      dumpScoringMeshesPerWorker();
       // we need to call finish run for Geant4 ... Since we use ProcessEvent() interface;
       geant4->FinishRun();
     }
-    terminated = true;
   }
 }

run/CMakeLists.txt

@@ -62,6 +62,11 @@ o2_add_executable(serial
                   COMPONENT_NAME sim
                   SOURCES o2sim.cxx
                   PUBLIC_LINK_LIBRARIES internal::allsim)
+o2_add_executable(merge-g4scoring
+                  COMPONENT_NAME sim
+                  SOURCES o2sim_mergeg4scoring.cxx
+                  PUBLIC_LINK_LIBRARIES O2::SimConfig)
+
 o2_add_executable(evalmat
                   COMPONENT_NAME sim
                   SOURCES o2sim_evalmat.cxx

run/O2SimDeviceRunner.cxx

@@ -59,6 +59,8 @@ void sigaction_handler(int signal, siginfo_t* signal_info, void*)
     // signal was sent from driver process --> not error
     // or it was a standard SIGTERM
 
+    // shut down before waiting, so that the master worker finalises (e.g. writes its scoring dumps) before the driver's kill timer
+    o2::SimSetup::shutdown();
     // need to wait for potential children before exiting itself
     // ... in order to have correct resource accounting
     int status, cpid;
@@ -67,7 +69,6 @@ void sigaction_handler(int signal, siginfo_t* signal_info, void*)
         break;
       }
     }
-    o2::SimSetup::shutdown();
     _exit(0);
   }