Draw transport random numbers from a per-thread generator

This adds o2::base::simRandom(), the generator used for seeding and drawing during transport, and moves the simulation to it.

- simRandom() returns a generator installed for the calling thread, and gRandom otherwise, so the sequential simulation is unchanged.
- Sub-event, per-track and stack seeding and the stepping draws in CPV, FDD, FT0, HMPID, PHOS, TPC, TRD and ZDC use it.
- VMCSeederService seeds Geant4 from it directly, with the recipe of TG4RunManager::SetRandomSeed, which reads gRandom.
- A unit test checks the default and that threads keep separate sequences.

This prepares the simulation for multithreaded Geant4, where gRandom would be shared by all workers.

https://its.cern.ch/jira/browse/O2-7305

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

Author: sawenzel

Detectors/Base/CMakeLists.txt

@@ -34,6 +34,7 @@ o2_add_library(DetectorsBase
                        src/GRPGeomHelper.cxx
                        src/Stack.cxx
                        src/VMCSeederService.cxx
+                       src/SimRandom.cxx
                        src/GlobalParams.cxx
                        src/O2Tessellated.cxx
                        src/TGeoGeometryUtils.cxx
@@ -89,6 +90,13 @@ o2_add_test(
   PUBLIC_LINK_LIBRARIES O2::DetectorsBase
   LABELS detectorsbase)
 
+o2_add_test(
+  SimRandom
+  SOURCES test/testSimRandom.cxx
+  COMPONENT_NAME DetectorsBase
+  PUBLIC_LINK_LIBRARIES O2::DetectorsBase
+  LABELS detectorsbase)
+
 o2_add_test(
   O2Tessellated
   SOURCES test/testO2Tessellated.cxx

Detectors/Base/include/DetectorsBase/SimRandom.h

@@ -0,0 +1,29 @@
+// 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 ALICEO2_BASE_SIMRANDOM_H
+#define ALICEO2_BASE_SIMRANDOM_H
+
+#include <TRandom.h>
+
+namespace o2::base
+{
+
+/// The random number generator used during particle transport, for seeding and for drawing.
+/// Each thread can install its own generator; a thread that has none uses gRandom.
+TRandom* simRandom();
+
+/// Install the transport generator of the calling thread (not owned); nullptr reverts to gRandom.
+void setSimRandom(TRandom* random);
+
+} // namespace o2::base
+
+#endif

Detectors/Base/include/DetectorsBase/VMCSeederService.h

@@ -13,7 +13,9 @@
 #define ALICEO2_VMC_SEEDERSERVICE_H
 
 class TVirtualMC;
+class TRandom;
 
+#include <atomic>
 #include <functional>
 
 namespace o2
@@ -33,19 +35,19 @@ class VMCSeederService
     return inst;
   }
 
-  void setSeed() const; // will propagate seed to the VMC engines
+  void setSeed() const; // will propagate the seed of o2::base::simRandom() to the VMC engine of this thread
 
   /// how often a seed was propagated; lets callers detect a silent no-op
   unsigned long long getSeedCount() const { return mSeedCount; }
 
-  typedef std::function<void()> SeederFcn;
+  typedef std::function<void(TRandom*)> SeederFcn;
 
  private:
   VMCSeederService();
   void initSeederFunction(TVirtualMC const*);
 
-  SeederFcn mSeederFcn;                     // the just-in-time compiled function talking to the VMC engines
-  mutable unsigned long long mSeedCount{0}; // number of setSeed() calls
+  SeederFcn mSeederFcn;                                  // the just-in-time compiled function talking to the VMC engines
+  mutable std::atomic<unsigned long long> mSeedCount{0}; // number of setSeed() calls
 };
 
 } // namespace base

Detectors/Base/src/SimRandom.cxx

@@ -0,0 +1,26 @@
+// 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 "DetectorsBase/SimRandom.h"
+
+namespace
+{
+thread_local TRandom* sThreadSimRandom = nullptr;
+}
+
+namespace o2::base
+{
+
+TRandom* simRandom() { return sThreadSimRandom ? sThreadSimRandom : gRandom; }
+
+void setSimRandom(TRandom* random) { sThreadSimRandom = random; }
+
+} // namespace o2::base

Detectors/Base/src/Stack.cxx

@@ -36,6 +36,7 @@
 #include <cassert>
 #include <cstddef> // for NULL
 #include <cmath>
+#include "DetectorsBase/SimRandom.h"
 
 using std::cout;
 using std::endl;
@@ -361,7 +362,7 @@ TParticle* Stack::PopNextTrack(Int_t& iTrack)
         auto hash = getHash(mCurrentParticle);
         // LOG(info) << "SEEDING NEW TRACK USING HASH" << hash;
         // init seed per track
-        gRandom->SetSeed(hash);
+        o2::base::simRandom()->SetSeed(hash);
         o2::base::VMCSeederService::instance().setSeed();
         // NOTE: THE BETTER PLACE WOULD BE IN PRETRACK HOOK BUT THIS DOES NOT SEEM TO WORK
         // WORKS ONLY WITH G3 SINCE G4 DOES NOT CALL THIS FUNCTION

Detectors/Base/src/VMCSeederService.cxx

@@ -10,24 +10,31 @@
 // or submit itself to any jurisdiction.
 
 #include "DetectorsBase/VMCSeederService.h"
+#include "DetectorsBase/SimRandom.h"
 #include "TVirtualMC.h"
 #include <fairlogger/Logger.h>                    // for FairLogger
 #include <CommonUtils/ConfigurationMacroHelper.h> // for ROOT JIT helpers
+#include <TInterpreter.h>
+#include <TRandom.h>
 
 using namespace o2::base;
 
 void VMCSeederService::initSeederFunction(TVirtualMC const* vmc)
 {
   if (strcmp(vmc->GetName(), "TGeant3TGeo") == 0) {
     // Geant3 doesn't need anything special in our context
-    mSeederFcn = []() {};
+    mSeederFcn = [](TRandom*) {};
   } else if (strcmp(vmc->GetName(), "TGeant4") == 0) {
-    // dynamically get access to the Geant4_VMC seeding function (without this function linking against Geant4)
-    std::string G4func("std::function<void()> G4func() { gSystem->Load(\"libgeant4vmc\"); return [](){ ((TGeant4*)TVirtualMC::GetMC())->SetRandomSeed(); };}");
-    mSeederFcn = o2::conf::JITAndEvalFunction<SeederFcn>(G4func, "G4func()", "std::function<void()>", "VMCSEEDERFUNC123");
+    // Seed the Geant4 engine of the calling thread from the given generator, with the same recipe as
+    // TG4RunManager::SetRandomSeed (which always reads gRandom). JIT-compiled to avoid linking against Geant4.
+    if (!gInterpreter->Declare("#include <CLHEP/Random/Random.h>")) {
+      LOG(fatal) << "CLHEP/Random/Random.h not found on ROOT_INCLUDE_PATH; cannot seed Geant4";
+    }
+    std::string G4func("std::function<void(TRandom*)> G4func() { gSystem->Load(\"libgeant4vmc\"); return [](TRandom* r){ long seeds[10]; seeds[0] = r->GetSeed(); seeds[1] = r->GetSeed(); seeds[2] = 0; CLHEP::HepRandom::setTheSeeds(seeds); };}");
+    mSeederFcn = o2::conf::JITAndEvalFunction<SeederFcn>(G4func, "G4func()", "std::function<void(TRandom*)>", "VMCSEEDERFUNC123");
   } else {
     LOG(warn) << "Unknown VMC engine or unimplemented VMC seeding function";
-    mSeederFcn = []() {};
+    mSeederFcn = [](TRandom*) {};
   }
 }
 
@@ -45,10 +52,11 @@ VMCSeederService::VMCSeederService()
 
 void VMCSeederService::setSeed() const
 {
-  // Make sure gRandom is sufficiently well initialized
+  // Make sure the generator is sufficiently well initialized
   // by calling Rndm() once.
-  // This is ok since in any case gRandom->SetSeed(seed); gRandom->GetSeed() != seed;
-  gRandom->Rndm();
-  mSeederFcn();
+  // This is ok since in any case r->SetSeed(seed); r->GetSeed() != seed;
+  auto r = simRandom();
+  r->Rndm();
+  mSeederFcn(r);
   ++mSeedCount;
 }

Detectors/Base/test/testSimRandom.cxx

@@ -0,0 +1,63 @@
+// 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.
+
+/// \file testSimRandom.cxx
+/// \author Sandro Wenzel (CERN)
+/// \brief Checks that o2::base::simRandom() defaults to gRandom and is per thread
+
+#define BOOST_TEST_MODULE Test SimRandom
+#define BOOST_TEST_MAIN
+#define BOOST_TEST_DYN_LINK
+#include <boost/test/unit_test.hpp>
+
+#include "DetectorsBase/SimRandom.h"
+#include <TRandom.h>
+#include <TRandom3.h>
+#include <thread>
+#include <vector>
+
+BOOST_AUTO_TEST_CASE(defaults_to_gRandom)
+{
+  BOOST_CHECK(o2::base::simRandom() == gRandom);
+  TRandom3 own(1);
+  o2::base::setSimRandom(&own);
+  BOOST_CHECK(o2::base::simRandom() == &own);
+  o2::base::setSimRandom(nullptr);
+  BOOST_CHECK(o2::base::simRandom() == gRandom);
+}
+
+BOOST_AUTO_TEST_CASE(per_thread_generators)
+{
+  // each thread draws from its own generator; equal seeds give equal sequences, whatever the other threads do
+  constexpr int nthreads = 8;
+  constexpr int ndraws = 1000;
+  std::vector<std::vector<double>> draws(nthreads);
+  std::vector<std::thread> threads;
+  for (int t = 0; t < nthreads; ++t) {
+    threads.emplace_back([t, &draws]() {
+      TRandom3 own(t % 2 + 1);
+      o2::base::setSimRandom(&own);
+      for (int i = 0; i < ndraws; ++i) {
+        draws[t].push_back(o2::base::simRandom()->Rndm());
+      }
+      o2::base::setSimRandom(nullptr);
+    });
+  }
+  for (auto& th : threads) {
+    th.join();
+  }
+  for (int t = 2; t < nthreads; ++t) {
+    BOOST_CHECK(draws[t] == draws[t % 2]);
+  }
+  BOOST_CHECK(draws[0] != draws[1]);
+  // the main thread was never changed
+  BOOST_CHECK(o2::base::simRandom() == gRandom);
+}

Detectors/CPV/simulation/src/Detector.cxx

@@ -32,6 +32,7 @@
 
 #include <boost/algorithm/string/predicate.hpp>
 #include <boost/range/irange.hpp>
+#include "DetectorsBase/SimRandom.h"
 
 using namespace o2::cpv;
 
@@ -171,7 +172,7 @@ Bool_t Detector::ProcessHits(FairVolume* v)
   float dZY = pZ / pNorm * cpvparam.mCPVGasThickness;
   float dXY = pX / pNorm * cpvparam.mCPVGasThickness;
   float rnor1 = 0., rnor2 = 0.;
-  gRandom->Rannor(rnor1, rnor2);
+  o2::base::simRandom()->Rannor(rnor1, rnor2);
   float eloss = cpvparam.mdEdx * (1 + cpvparam.mDetR * rnor1) *
                 TMath::Sqrt((1. + (pow(dZY, 2) + pow(dXY, 2)) / pow(cpvparam.mCPVGasThickness, 2)));
   float zhit1 = hitZ + cpvparam.mPadSizeZ * cpvparam.mnCellZ / 2 - dZY / 2;

Detectors/FIT/FDD/simulation/src/Detector.cxx

@@ -33,6 +33,7 @@
 #include "FairRootManager.h"
 #include "FairVolume.h"
 #include "FairRootManager.h"
+#include "DetectorsBase/SimRandom.h"
 
 using namespace o2::fdd;
 using o2::fdd::Geometry;
@@ -149,7 +150,7 @@ Bool_t Detector::ProcessHits(FairVolume* vol)
 
   Float_t destep_ad = fMC->Edep();
   Int_t nPhotonsInStep_ad = Int_t(fFDDLightYield * destep_ad);
-  nPhotonsInStep_ad = gRandom->Poisson(nPhotonsInStep_ad);
+  nPhotonsInStep_ad = o2::base::simRandom()->Poisson(nPhotonsInStep_ad);
 
   mElossAD += destep_ad;
 

Detectors/FIT/FT0/simulation/src/Detector.cxx

@@ -37,6 +37,7 @@
 #include "FT0Base/Geometry.h"
 #include "FT0Simulation/Detector.h"
 #include "DetectorsBase/Stack.h"
+#include "DetectorsBase/SimRandom.h"
 
 using namespace o2::ft0;
 using o2::ft0::Geometry;
@@ -1373,7 +1374,7 @@ void Detector::FillOtherOptProperties()
 Bool_t Detector::RegisterPhotoE(float energy)
 {
   float eff = mPMTeff->Eval(energy);
-  float p = gRandom->Rndm();
+  float p = o2::base::simRandom()->Rndm();
   if (p > eff) {
     return kFALSE;
   }