diff --git a/Common/SimConfig/CMakeLists.txt b/Common/SimConfig/CMakeLists.txt index 65d30935904ad..5737221789471 100644 --- a/Common/SimConfig/CMakeLists.txt +++ b/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) diff --git a/Common/SimConfig/include/SimConfig/FluenceWeightCalculator.h b/Common/SimConfig/include/SimConfig/FluenceWeightCalculator.h index 15d74ba27ab1b..0936264416e35 100644 --- a/Common/SimConfig/include/SimConfig/FluenceWeightCalculator.h +++ b/Common/SimConfig/include/SimConfig/FluenceWeightCalculator.h @@ -31,5 +31,6 @@ class FluenceWeightCalculator static std::unique_ptr neutronG; static std::unique_ptr protonG; static std::unique_ptr pionG; + static std::unique_ptr electronG; }; #endif diff --git a/Common/SimConfig/include/SimConfig/G4ScoringMerger.h b/Common/SimConfig/include/SimConfig/G4ScoringMerger.h new file mode 100644 index 0000000000000..2deb0792f2533 --- /dev/null +++ b/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 + +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 .worker.txt in a directory into .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 diff --git a/Common/SimConfig/src/FluenceWeightCalculator.cxx b/Common/SimConfig/src/FluenceWeightCalculator.cxx index 63828f71286e2..04a8950d88a7b 100644 --- a/Common/SimConfig/src/FluenceWeightCalculator.cxx +++ b/Common/SimConfig/src/FluenceWeightCalculator.cxx @@ -11,6 +11,7 @@ #include "SimConfig/FluenceWeightCalculator.h" #include +#include #include #include #include @@ -18,6 +19,20 @@ std::unique_ptr FluenceWeightCalculator::neutronG; std::unique_ptr FluenceWeightCalculator::protonG; std::unique_ptr FluenceWeightCalculator::pionG; +std::unique_ptr 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(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(); pionG->SetName("pionDW"); auto pioN = 0; + electronG = std::make_unique(); + 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(); } diff --git a/Common/SimConfig/src/G4ScoringMerger.cxx b/Common/SimConfig/src/G4ScoringMerger.cxx new file mode 100644 index 0000000000000..c59f65542f105 --- /dev/null +++ b/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 +#include +#include +#include +#include +#include +#include +#include + +namespace o2::conf +{ + +namespace +{ +// One scorer block of a Geant4 mesh dump: its header lines and the summed rows +struct ScorerBlock { + std::vector header; + std::vector keys; // "iZ,iPHI,iR" in file order + std::vector sum; + std::vector sum2; + std::vector entries; +}; + +// Read one mesh dump into scorer blocks; returns false on a format error +bool readDump(const std::string& fileName, std::vector& meshHeader, std::vector& 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 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> 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(files.size()) != expectedWorkers) { + LOG(error) << "Found " << files.size() << " Geant4 scoring dumps for mesh " << mesh << " but expected " << expectedWorkers; + return -1; + } + std::vector meshHeader; + std::vector total; + for (auto& file : files) { + std::vector header; + std::vector 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 diff --git a/Detectors/gconfig/g4Config.C b/Detectors/gconfig/g4Config.C index 77494c559ca51..1907c1aa9ebcd 100644 --- a/Detectors/gconfig/g4Config.C +++ b/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 #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(TVirtualMC::GetMC()); if (geant4) { + dumpScoringMeshesPerWorker(); // we need to call finish run for Geant4 ... Since we use ProcessEvent() interface; geant4->FinishRun(); } - terminated = true; } } diff --git a/run/CMakeLists.txt b/run/CMakeLists.txt index 0b88a6e6f68d5..abba055cc7cca 100644 --- a/run/CMakeLists.txt +++ b/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 diff --git a/run/O2SimDeviceRunner.cxx b/run/O2SimDeviceRunner.cxx index 524e609883926..609311809d5d9 100644 --- a/run/O2SimDeviceRunner.cxx +++ b/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); } diff --git a/run/o2sim_mergeg4scoring.cxx b/run/o2sim_mergeg4scoring.cxx new file mode 100644 index 0000000000000..46f751bef8334 --- /dev/null +++ b/run/o2sim_mergeg4scoring.cxx @@ -0,0 +1,28 @@ +// 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. + +// Sum the per-worker Geant4 scoring dumps of an o2-sim run into one file per mesh + +#include "SimConfig/G4ScoringMerger.h" +#include +#include + +int main(int argc, char* argv[]) +{ + const std::string directory = argc > 1 ? argv[1] : "."; + const int expectedWorkers = argc > 2 ? std::atoi(argv[2]) : 0; + const int merged = o2::conf::mergeG4ScoringDumps(directory, expectedWorkers); + if (merged < 0) { + return 1; + } + std::cout << "merged " << merged << " scoring mesh(es) in " << directory << "\n"; + return 0; +} diff --git a/run/o2sim_parallel.cxx b/run/o2sim_parallel.cxx index c060059ce4d6b..8a92a5f251cb0 100644 --- a/run/o2sim_parallel.cxx +++ b/run/o2sim_parallel.cxx @@ -11,6 +11,7 @@ /// @author Sandro Wenzel +#include "SimConfig/G4ScoringMerger.h" #include #include #include @@ -806,6 +807,18 @@ int main(int argc, char* argv[]) LOG(debug) << "ShmManager operation " << o2::utils::ShmManager::Instance().isOperational() << "\n"; + // forked sim workers can still be writing their scoring dumps after their parent exited + for (auto p : gChildProcesses) { + while (p != 0 && killpg(p, 0) == 0) { + usleep(100000); + } + } + + // sum the Geant4 scoring meshes written by the individual workers + if (!errored && o2::conf::mergeG4ScoringDumps(".", conf.getNSimWorkers()) < 0) { + errored = true; + } + // do a quick check to see if simulation produced something reasonable // (mainly useful for continuous integration / automated testing suite) auto returncode = errored ? 1 : checkresult();