diff --git a/Detectors/Base/include/DetectorsBase/VMCSeederService.h b/Detectors/Base/include/DetectorsBase/VMCSeederService.h index 1669c73b39620..5f8f70f48840f 100644 --- a/Detectors/Base/include/DetectorsBase/VMCSeederService.h +++ b/Detectors/Base/include/DetectorsBase/VMCSeederService.h @@ -35,13 +35,17 @@ class VMCSeederService void setSeed() const; // will propagate seed to the VMC engines + /// how often a seed was propagated; lets callers detect a silent no-op + unsigned long long getSeedCount() const { return mSeedCount; } + typedef std::function SeederFcn; private: VMCSeederService(); void initSeederFunction(TVirtualMC const*); - SeederFcn mSeederFcn; // the just-in-time compiled function talking to the VMC engines + SeederFcn mSeederFcn; // the just-in-time compiled function talking to the VMC engines + mutable unsigned long long mSeedCount{0}; // number of setSeed() calls }; } // namespace base diff --git a/Detectors/Base/src/VMCSeederService.cxx b/Detectors/Base/src/VMCSeederService.cxx index 5bf4e1ed5641b..8fc36d9074fab 100644 --- a/Detectors/Base/src/VMCSeederService.cxx +++ b/Detectors/Base/src/VMCSeederService.cxx @@ -50,4 +50,5 @@ void VMCSeederService::setSeed() const // This is ok since in any case gRandom->SetSeed(seed); gRandom->GetSeed() != seed; gRandom->Rndm(); mSeederFcn(); + ++mSeedCount; } diff --git a/Steer/include/Steer/O2MCApplicationBase.h b/Steer/include/Steer/O2MCApplicationBase.h index d61199baba0ae..bd730c0f2fcb2 100644 --- a/Steer/include/Steer/O2MCApplicationBase.h +++ b/Steer/include/Steer/O2MCApplicationBase.h @@ -68,6 +68,11 @@ class O2MCApplicationBase : public FairMCApplication // keeping track of volumeIds and volume names double mLongestTrackTime = 0; + bool mTrackSeedWarned{false}; // whether we already complained that seeding never fired + + /// whether this engine needs per-track seeding in PreTrack (Geant3 seeds at + /// stack-pop time instead, see O2MCApplicationBase::seedsInPreTrack) + bool seedsInPreTrack() const; /// some common parts of finishEvent void finishEventCommon(); TrackRefFcn mTrackRefFcn; // a function hook that gets (optionally) called during Stepping diff --git a/Steer/src/O2MCApplication.cxx b/Steer/src/O2MCApplication.cxx index 1e3f925042d01..3a39b7d15be8b 100644 --- a/Steer/src/O2MCApplication.cxx +++ b/Steer/src/O2MCApplication.cxx @@ -35,6 +35,9 @@ #include #include #include "SimConfig/GlobalProcessCutSimParam.h" +#include +#include // full type: FairField derives from TVirtualMagField +#include #include "DetectorsBase/GeometryManagerParam.h" #include #include @@ -43,6 +46,10 @@ #include #include #include "SimConfig/G4Params.h" +#include "DetectorsBase/VMCSeederService.h" // per-track seeding of the engine +#include +#include +#include namespace o2 { @@ -116,14 +123,98 @@ void O2MCApplicationBase::Stepping() FairMCApplication::Stepping(); } +namespace +{ +// Hash of a track's initial state (vertex, global time, momentum, PDG). Used as +// the random seed for that track, so that a track's random stream depends only +// on the track itself and not on how many randoms earlier tracks happened to +// consume. +// +// The values are read from the transport engine, not from +// o2::data::Stack::GetCurrentTrack(): under Geant4 the stack's "current track" +// is only meaningful for primaries -- Stack::SetCurrentTrack() falls back to +// mCurrentParticle0 (the last particle *pushed*) for anything beyond the +// primary array, so every secondary would hash the wrong particle. Both engines +// have the track's initial state loaded by the time PreTrack is called (Geant4 +// sets the step to kVertex first; Geant3 calls GLTRAC before GUTRAK). +ULong_t hashCurrentTrack(TVirtualMC* vmc) +{ + auto asLong = [](double x) { + ULong_t l; + std::memcpy(&l, &x, sizeof(l)); + return l; + }; + + TLorentzVector pos, mom; + vmc->TrackPosition(pos); + vmc->TrackMomentum(mom); + + ULong_t hash = asLong(pos.X()); + hash ^= asLong(pos.Y()); + hash ^= asLong(pos.Z()); + hash ^= asLong(pos.T()); + hash ^= asLong(mom.Px()); + hash ^= asLong(mom.Py()); + hash ^= asLong(mom.Pz()); + hash += (ULong_t)vmc->TrackPid(); + return hash; +} +} // namespace + +bool O2MCApplicationBase::seedsInPreTrack() const +{ + // Geant3 seeds at stack-pop time, in o2::data::Stack::PopNextTrack(). That is + // strictly earlier than its PreTrack hook (gutrak) and measurably stronger: + // with the TOF module removed from an otherwise identical setup, pop-time + // seeding keeps all 603 ITS hits bit-identical, PreTrack seeding only 68 %. + // Do not seed Geant3 here as well -- it is already covered, and reseeding a + // second time mid-track would undo the first. + static const bool inPreTrack = [this]() { + const char* name = (fMC != nullptr) ? fMC->GetName() : ""; + return strncmp(name, "TGeant3", 7) != 0; + }(); + return inPreTrack; +} + void O2MCApplicationBase::PreTrack() { - // dispatch first to function in FairRoot + if (mCutParams.trackSeed && seedsInPreTrack()) { + // Per-track seeding for engines that do not go through + // o2::data::Stack::PopNextTrack(). Geant4 is one: it takes primaries via + // PopPrimaryForTracking and keeps secondaries internally, so the stack hook + // never fires and this is the only per-track hook available. It is called + // for primaries and secondaries alike + // (TG4TrackingAction::PreUserTrackingAction), and only on a track's first + // step, so a suspended track is not reseeded mid-flight. + auto hash = hashCurrentTrack(fMC); + // TRandom::SetSeed(0) means "seed from the clock" -- never let that happen. + gRandom->SetSeed(hash == 0 ? 1 : hash); + o2::base::VMCSeederService::instance().setSeed(); + } + + // dispatch now to function in FairRoot FairMCApplication::PreTrack(); } void O2MCApplicationBase::ConstructGeometry() { + // The transport engine constructs the geometry from inside its own + // constructor, long before FairMCApplication::InitMC() attaches the magnetic + // field to it. The media built below read the field through + // Detector::initFieldTrackingParams(), so without this they all silently fall + // back to hardcoded defaults. The run has known the field since + // build_geometry.C, which runs before Init() -- hand it over now. + if (auto* vmc = TVirtualMC::GetMC(); vmc != nullptr && vmc->GetMagField() == nullptr) { + auto* run = FairRunSim::Instance(); + if (run != nullptr && run->GetField() != nullptr) { + vmc->SetMagField(run->GetField()); + LOG(info) << "Magnetic field attached to the engine before media creation"; + } else { + LOG(warn) << "No magnetic field available at geometry construction; media " + "will be initialised with default tracking parameters"; + } + } + // fill the mapping mModIdToName.clear(); o2::detectors::DetID::mask_t dmask{}; @@ -289,6 +380,17 @@ void O2MCApplicationBase::finishEventCommon() header->setDetId2HitBitLUT(o2::base::Detector::getDetId2HitBitIndex()); static_cast(GetStack())->updateEventStats(); + + // Per-track seeding used to be wired to a stack callback that one of the two + // engines never invoked, and it failed silently. Never again: if it was asked + // for and nothing was seeded, say so. + if (mCutParams.trackSeed && o2::base::VMCSeederService::instance().getSeedCount() == 0 && + !mTrackSeedWarned) { + mTrackSeedWarned = true; + LOG(warn) << "Per-track seeding (SimCutParams.trackSeed) was requested but not a single track " + "was seeded -- neither the stack nor the PreTrack hook fired for this engine. " + "Seeding is NOT active."; + } } void O2MCApplicationBase::FinishEvent()