Draw transport random numbers from a per-thread generator

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,27 @@ 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;
+
+  /// installs the function that seeds Geant4 from a generator; called by the Geant4 configuration, which links Geant4
+  static void setGeant4Seeder(SeederFcn fcn) { geant4Seeder() = std::move(fcn); }
 
  private:
   VMCSeederService();
+  static SeederFcn& geant4Seeder()
+  {
+    static SeederFcn fcn;
+    return fcn;
+  }
   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 function talking to the VMC engine
+  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,27 @@
 // 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 <fairlogger/Logger.h> // for FairLogger
+#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");
+    // installed by the Geant4 configuration, so that this library does not link against Geant4
+    mSeederFcn = geant4Seeder();
+    if (!mSeederFcn) {
+      LOG(fatal) << "No Geant4 seeding function installed (VMCSeederService::setGeant4Seeder)";
+    }
   } else {
     LOG(warn) << "Unknown VMC engine or unimplemented VMC seeding function";
-    mSeederFcn = []() {};
+    mSeederFcn = [](TRandom*) {};
   }
 }
 
@@ -45,10 +48,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

@@ -33,6 +33,7 @@
 
 #include <boost/algorithm/string/predicate.hpp>
 #include <boost/range/irange.hpp>
+#include "DetectorsBase/SimRandom.h"
 
 using namespace o2::cpv;
 
@@ -172,7 +173,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);
 
   static Float_t eloss_ad = 0.;
   static Float_t tlength_ad = 0.;

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;
   }

Detectors/HMPID/simulation/src/Detector.cxx

@@ -33,6 +33,7 @@
 #include "DetectorsBase/MaterialManager.h"
 
 #include "Framework/Logger.h"
+#include "DetectorsBase/SimRandom.h"
 
 namespace o2
 {
@@ -187,22 +188,22 @@ void Detector::GenFee(Float_t qtot)
   // eloss=0 means photon so only pulse height distribution is to be analysed.
   TLorentzVector x4;
   fMC->TrackPosition(x4);
-  Int_t iNphotons = fMC->GetRandom()->Poisson(0.02 * qtot); // # of feedback photons is proportional to the charge of hit
+  Int_t iNphotons = o2::base::simRandom()->Poisson(0.02 * qtot); // # of feedback photons is proportional to the charge of hit
   // AliDebug(1,Form("N photons=%i",iNphotons));
   Int_t j;
   Float_t cthf, phif, enfp = 0, sthf, e1[3], e2[3], e3[3], vmod, uswop, dir[3], phi, pol[3], mom[4];
   // Generate photons
   for (Int_t i = 0; i < iNphotons; i++) { // feedbacks loop
     Double_t ranf[2];
-    fMC->GetRandom()->RndmArray(2, ranf); // Sample direction
+    o2::base::simRandom()->RndmArray(2, ranf); // Sample direction
     cthf = ranf[0] * 2 - 1.0;
     if (cthf < 0) {
       continue;
     }
     sthf = TMath::Sqrt((1. - cthf) * (1. + cthf));
     phif = ranf[1] * 2 * TMath::Pi();
 
-    if (Double_t randomNumber = fMC->GetRandom()->Rndm() <= 0.57) {
+    if (Double_t randomNumber = o2::base::simRandom()->Rndm() <= 0.57) {
       enfp = 7.5e-9;
     } else if (randomNumber <= 0.7) {
       enfp = 6.4e-9;
@@ -270,7 +271,7 @@ void Detector::GenFee(Float_t qtot)
       e2[j] *= vmod;
     }
 
-    phi = fMC->GetRandom()->Rndm() * 2 * TMath::Pi();
+    phi = o2::base::simRandom()->Rndm() * 2 * TMath::Pi();
     for (j = 0; j < 3; j++) {
       pol[j] = e1[j] * TMath::Sin(phi) + e2[j] * TMath::Cos(phi);
     }
@@ -295,7 +296,7 @@ Bool_t Detector::IsLostByFresnel()
   Double_t localTc = localMom[0] * localMom[0] + localMom[2] * localMom[2];
   Double_t localTheta = TMath::ATan2(TMath::Sqrt(localTc), localMom[1]);
   Double_t cotheta = TMath::Abs(TMath::Cos(localTheta));
-  if (fMC->GetRandom()->Rndm() < Fresnel(p4.E() * 1e9, cotheta, 1)) {
+  if (o2::base::simRandom()->Rndm() < Fresnel(p4.E() * 1e9, cotheta, 1)) {
     // AliDebug(1,"Photon lost");
     return kTRUE;
   } else {

Detectors/PHOS/simulation/src/Detector.cxx

@@ -33,6 +33,7 @@
 
 #include <boost/algorithm/string/predicate.hpp>
 #include <boost/range/irange.hpp>
+#include "DetectorsBase/SimRandom.h"
 
 using namespace o2::phos;
 
@@ -108,7 +109,7 @@ void Detector::FinishEvent()
   first = mHits->begin();
   last = mHits->end();
   while (first != last) {
-    float light = gRandom->Poisson(first->GetEnergyLoss() * o2::phos::PHOSSimParams::Instance().mLightYieldPerGeV);
+    float light = o2::base::simRandom()->Poisson(first->GetEnergyLoss() * o2::phos::PHOSSimParams::Instance().mLightYieldPerGeV);
     first->SetEnergyLoss(light / o2::phos::PHOSSimParams::Instance().mLightYieldPerGeV);
     first++;
   }

Detectors/Passive/src/Cave.cxx

@@ -33,6 +33,7 @@
 #include "TGeoPgon.h"
 #include "TGeoTube.h"
 #include "TGeoCompositeShape.h"
+#include "DetectorsBase/SimRandom.h"
 using namespace o2::passive;
 
 void Cave::createMaterials()
@@ -156,7 +157,7 @@ void Cave::BeginPrimary()
       if (primcounter % chunks == 0) {
         static int counter = 1;
         auto seed = counter + 10;
-        gRandom->SetSeed(seed);
+        o2::base::simRandom()->SetSeed(seed);
         counter++;
       }
     }