diff --git a/Common/SimConfig/include/SimConfig/G4Params.h b/Common/SimConfig/include/SimConfig/G4Params.h index 2a333a39e4242..97c3cc4ddb412 100644 --- a/Common/SimConfig/include/SimConfig/G4Params.h +++ b/Common/SimConfig/include/SimConfig/G4Params.h @@ -54,6 +54,14 @@ struct G4Params : public o2::conf::ConfigurableParamHelper { bool g4scoring = false; bool g4fluenceweight = false; + + // Fast simulation. Empty fastSimModels (the default) disables the feature + // entirely; see Detectors/gconfig/include/SimSetup/G4FastSimulation.h. + std::string fastSimModels = ""; // comma-separated model names to activate + std::string fastSimEnvelope = ""; // volume a model stands in for, e.g. AFaM; the media of its + // subtree are collected automatically + std::string fastSimRegions = ""; // optional explicit media, overriding the subtree walk + float fastSimMinEnergy = 1.f; // GeV; below this the detailed transport runs O2ParamDef(G4Params, "G4"); }; diff --git a/Detectors/CMakeLists.txt b/Detectors/CMakeLists.txt index 9143761508998..09e784b0338ee 100644 --- a/Detectors/CMakeLists.txt +++ b/Detectors/CMakeLists.txt @@ -53,6 +53,7 @@ add_subdirectory(ForwardAlign) if(BUILD_SIMULATION) + add_subdirectory(FastSim) add_subdirectory(gconfig) o2_data_file(COPY gconfig DESTINATION Detectors) endif() diff --git a/Detectors/FastSim/CMakeLists.txt b/Detectors/FastSim/CMakeLists.txt new file mode 100644 index 0000000000000..ae4de342e391b --- /dev/null +++ b/Detectors/FastSim/CMakeLists.txt @@ -0,0 +1,18 @@ +# 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. + +o2_add_library(FastSim + SOURCES src/FastSimModel.cxx + src/FastSimRegions.cxx + src/ToyAbsorberFastSim.cxx + src/G4FastSimulation.cxx + PUBLIC_LINK_LIBRARIES MC::Geant4VMC MC::Geant4 O2::SimConfig +) diff --git a/Detectors/FastSim/include/FastSim/FastSimModel.h b/Detectors/FastSim/include/FastSim/FastSimModel.h new file mode 100644 index 0000000000000..3bcd5cbe52fd8 --- /dev/null +++ b/Detectors/FastSim/include/FastSim/FastSimModel.h @@ -0,0 +1,105 @@ +// 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_FASTSIM_MODEL_H_ +#define O2_FASTSIM_MODEL_H_ + +/// Base class for fast simulation models. +/// +/// A fast simulation model replaces the detailed transport through a region of +/// the geometry by a function from the particle that enters it to the particles +/// that leave it. `DoIt()` is Geant4's own entry point +/// (`G4VFastSimulationModel::DoIt`); the implementation here wraps the logic +/// common to every model and delegates the physics to `sample()`. A model that +/// needs a different shape can still override `DoIt()`. + +#include "G4VFastSimulationModel.hh" + +#include + +class G4FastStep; +class G4FastTrack; +class G4ParticleDefinition; + +namespace o2::fastsim +{ + +/// The particle entering the envelope, plus the geometric context a model needs. +/// Units are the O2/VMC ones: cm, GeV, ns. +struct FastSimInput { + int pdg = 0; + double position[3] = {}; ///< global, on the envelope surface + double direction[3] = {}; ///< unit vector + double kineticEnergy = 0.; ///< GeV + double mass = 0.; ///< GeV + double time = 0.; ///< ns + double exitDistance = 0.; ///< cm from `position` to the ENVELOPE surface along `direction` +}; + +/// One particle leaving the envelope. +struct FastSimOutput { + int pdg = 0; + double position[3] = {}; ///< global; put it outside the envelope surface + double momentum[3] = {}; ///< GeV/c + double time = 0.; ///< ns +}; + +/// A secondary created exactly on the envelope surface is located by the +/// navigator in whichever daughter owns that point, which costs two extra +/// zero-length steps before it gets out. Models should emit just beyond it. +constexpr double kSurfaceEpsilonCm = 1e-5; + +/// Base class for fast simulation models. +/// +/// The model is attached to regions (see G4FastSimulation.h) purely so that +/// Geant4 consults it; what it encloses is the ENVELOPE VOLUME named below, +/// which is normally the mother volume of a whole module. The two are +/// deliberately separate, because a Geant4 region in O2 can only ever be "every +/// volume of a given material" -- the VMC special cuts make every logical volume +/// a root of its own material's region, and Geant4 stops propagating a region at +/// any such daughter. So `G4FastTrack::GetEnvelopeSolid()` would hand back one +/// absorber piece rather than the absorber, and this class does not use it. +/// +/// Containment and the exit distance are taken from the track's own touchable, +/// which already carries the full ancestry and the transform of every level. +class FastSimModel : public G4VFastSimulationModel +{ + public: + FastSimModel(const G4String& name, const G4String& envelopeVolume, double minEnergyGeV); + + G4bool IsApplicable(const G4ParticleDefinition& particle) override; + G4bool ModelTrigger(const G4FastTrack& fastTrack) override; + + /// Wraps the common logic of a fast simulation action and delegates the + /// physics to `sample()`: it measures the distance to the envelope surface, + /// kills the incident particle, stacks what `sample()` returned and books the + /// energy difference as a deposit. Override it for a model that does not fit + /// that shape. + void DoIt(const G4FastTrack& fastTrack, G4FastStep& fastStep) override; + + protected: + /// Given the particle that entered, return everything that leaves. This is + /// the function a trained model implements. + virtual std::vector sample(const FastSimInput& input) const = 0; + + private: + /// The track's ancestry level at which the envelope volume sits, or -1 when + /// the track is not inside it at all. + int envelopeDepth(const G4Track* track) const; + + G4String mEnvelope; ///< logical volume the model encloses + double mMinEnergy = 0.; ///< internal Geant4 units; below this the detailed transport runs + mutable bool mWarned = false; +}; + +} // namespace o2::fastsim + +#endif // O2_FASTSIM_MODEL_H_ diff --git a/Detectors/FastSim/include/FastSim/FastSimRegions.h b/Detectors/FastSim/include/FastSim/FastSimRegions.h new file mode 100644 index 0000000000000..a47a9319a18ba --- /dev/null +++ b/Detectors/FastSim/include/FastSim/FastSimRegions.h @@ -0,0 +1,65 @@ +// 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_FASTSIM_REGIONS_H_ +#define O2_FASTSIM_REGIONS_H_ + +/// Deriving a model's regions from its envelope volume. +/// +/// Geant4-VMC attaches a fast simulation model to regions, and in O2 a region +/// can only be "every volume of a given material" (see FastSimModel.h). To have +/// the model consulted everywhere inside a module, it must therefore be attached +/// to every material that module is built from -- which is a list nobody should +/// maintain by hand, because it changes whenever the geometry does. +/// +/// So walk the envelope's subtree and collect the media as they actually are. + +#include "TG4VUserPostDetConstruction.h" + +#include +#include +#include + +class TGeoVolume; + +namespace o2::fastsim +{ + +/// Every tracking medium used by `volume` or anything below it. +/// Takes a non-const pointer because TGeo's accessors are not const. +std::set mediaInSubtree(TGeoVolume* volume); + +/// Same, looked up by volume name in the current TGeo geometry. Empty if there +/// is no such volume. +std::set mediaInSubtree(const std::string& volumeName); + +/// Sets each model's regions from its envelope, in the one window where that is +/// possible: after the geometry is built and before Geant4-VMC turns the media +/// into regions. +class FastSimRegionConstruction : public TG4VUserPostDetConstruction +{ + public: + struct ModelRegions { + std::string model; + std::string envelope; ///< volume whose subtree supplies the media + std::string regions; ///< explicit media, used instead of the walk if given + }; + + explicit FastSimRegionConstruction(std::vector models); + void Construct() override; + + private: + std::vector mModels; +}; + +} // namespace o2::fastsim + +#endif // O2_FASTSIM_REGIONS_H_ diff --git a/Detectors/FastSim/include/FastSim/G4FastSimulation.h b/Detectors/FastSim/include/FastSim/G4FastSimulation.h new file mode 100644 index 0000000000000..900e5e2d01ef7 --- /dev/null +++ b/Detectors/FastSim/include/FastSim/G4FastSimulation.h @@ -0,0 +1,68 @@ +// 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_FASTSIM_G4_FAST_SIMULATION_H_ +#define O2_FASTSIM_G4_FAST_SIMULATION_H_ + +/// Wiring of the fast simulation models into the Geant4 engine. +/// +/// The feature is OFF unless `G4.fastSimModels` names a model. +/// +/// o2-sim -n 10 -g pythia8pp -e TGeant4 -m PIPE ABSO +/// --configKeyValues "G4.fastSimModels=toyAbsorber; +/// G4.fastSimEnvelope=AFaM" +/// +/// `G4.fastSimEnvelope` names the VOLUME the model stands in for. The regions +/// Geant4 needs in order to consult the model are derived from it by walking its +/// subtree and collecting the media (FastSimRegions.h) -- a region in O2 can only +/// be "every volume of a given material", so covering a module means naming all +/// of its materials, and that list should not be maintained by hand. +/// +/// `G4.fastSimRegions` overrides the walk with an explicit space-separated list +/// of media, for when a model should see less than a whole subtree. + +#include "TG4RunConfiguration.h" +#include "TG4VUserFastSimulation.h" +#include "TG4VUserPostDetConstruction.h" + +#include +#include + +namespace o2::fastsim +{ + +/// Creates and registers the models named in `G4.fastSimModels`. +class G4FastSimulation : public TG4VUserFastSimulation +{ + public: + G4FastSimulation(std::vector models, const std::string& envelope, + double minEnergyGeV); + void Construct() override; + + private: + std::vector mModels; + std::string mEnvelope; + double mMinEnergy = 1.; +}; + +/// Supplies Geant4-VMC with the fast simulation models and their regions. +/// Returns nullptr when no model is configured, so nothing is set up. +class G4RunConfiguration : public TG4RunConfiguration +{ + public: + using TG4RunConfiguration::TG4RunConfiguration; + TG4VUserFastSimulation* CreateUserFastSimulation() override; + TG4VUserPostDetConstruction* CreateUserPostDetConstruction() override; +}; + +} // namespace o2::fastsim + +#endif // O2_FASTSIM_G4_FAST_SIMULATION_H_ diff --git a/Detectors/FastSim/include/FastSim/ToyAbsorberFastSim.h b/Detectors/FastSim/include/FastSim/ToyAbsorberFastSim.h new file mode 100644 index 0000000000000..95ed180bb9527 --- /dev/null +++ b/Detectors/FastSim/include/FastSim/ToyAbsorberFastSim.h @@ -0,0 +1,34 @@ +// 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_FASTSIM_TOY_ABSORBER_H_ +#define O2_FASTSIM_TOY_ABSORBER_H_ + +#include "FastSim/FastSimModel.h" + +namespace o2::fastsim +{ + +/// A toy fast simulation model for the absorber: one particle out, continuing +/// along the incident direction with the energy exponentially attenuated over +/// the path through the envelope. +class ToyAbsorberFastSim : public FastSimModel +{ + public: + using FastSimModel::FastSimModel; + + protected: + std::vector sample(const FastSimInput& input) const override; +}; + +} // namespace o2::fastsim + +#endif // O2_FASTSIM_TOY_ABSORBER_H_ diff --git a/Detectors/FastSim/src/FastSimModel.cxx b/Detectors/FastSim/src/FastSimModel.cxx new file mode 100644 index 0000000000000..2237d0cd5e69c --- /dev/null +++ b/Detectors/FastSim/src/FastSimModel.cxx @@ -0,0 +1,174 @@ +// 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 "FastSim/FastSimModel.h" + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace o2::fastsim +{ + +//_____________________________________________________________________________ +FastSimModel::FastSimModel(const G4String& name, const G4String& envelopeVolume, + double minEnergyGeV) + : G4VFastSimulationModel(name), + mEnvelope(envelopeVolume), + mMinEnergy(minEnergyGeV * CLHEP::GeV) +{ +} + +//_____________________________________________________________________________ +int FastSimModel::envelopeDepth(const G4Track* track) const +{ + /// Where the envelope volume sits in the track's ancestry, or -1 if the track + /// is not inside it. + /// + /// The touchable is the cheapest exact answer to "is this track inside the + /// module": it is the navigator's own record of the volume and every one of + /// its ancestors, so no geometry lookup, no cached transform and no name list + /// is needed, and it stays correct if the envelope is ever placed more than + /// once. + const G4VTouchable* touchable = track->GetTouchable(); + if (touchable == nullptr) { + return -1; + } + const G4int depth = touchable->GetHistoryDepth(); + for (G4int level = 0; level <= depth; ++level) { + const G4VPhysicalVolume* volume = touchable->GetVolume(level); + if (volume != nullptr && volume->GetLogicalVolume()->GetName() == mEnvelope) { + return level; + } + } + return -1; +} + +//_____________________________________________________________________________ +G4bool FastSimModel::IsApplicable(const G4ParticleDefinition&) +{ + // Which particles a model sees is decided by the `setParticles` selection, + // not here. + return true; +} + +//_____________________________________________________________________________ +G4bool FastSimModel::ModelTrigger(const G4FastTrack& fastTrack) +{ + const G4Track* track = fastTrack.GetPrimaryTrack(); + + // Below the threshold the detailed transport is cheap and a surrogate would + // be extrapolating. + if (track->GetKineticEnergy() <= mMinEnergy) { + return false; + } + + // Geometric containment rather than a name list. This is what excludes, for + // instance, the absorber's steel support cradle: it shares its material with + // parts of the absorber, so no selection by material can separate them, but + // it sits outside the envelope and so fails here. + if (envelopeDepth(track) < 0) { + if (!mWarned) { + mWarned = true; + LOG(warn) << "fast simulation: model " << GetName() << " was consulted for a track " + << "outside its envelope '" << mEnvelope << "'; the region selection is " + << "wider than the envelope, which is allowed but wasteful"; + } + return false; + } + return true; +} + +//_____________________________________________________________________________ +void FastSimModel::DoIt(const G4FastTrack& fastTrack, G4FastStep& fastStep) +{ + const G4Track* track = fastTrack.GetPrimaryTrack(); + const G4ThreeVector& position = track->GetPosition(); + const G4ThreeVector& direction = track->GetMomentumDirection(); + + FastSimInput input; + input.pdg = track->GetDefinition()->GetPDGEncoding(); + input.position[0] = position.x() / CLHEP::cm; + input.position[1] = position.y() / CLHEP::cm; + input.position[2] = position.z() / CLHEP::cm; + input.direction[0] = direction.x(); + input.direction[1] = direction.y(); + input.direction[2] = direction.z(); + input.kineticEnergy = track->GetKineticEnergy() / CLHEP::GeV; + input.mass = track->GetDefinition()->GetPDGMass() / CLHEP::GeV; + input.time = track->GetGlobalTime() / CLHEP::ns; + // Deliberately NOT GetEnvelopeSolid(): that is the region's root volume, i.e. + // one absorber piece. Use the envelope volume instead, with the transform the + // touchable already holds for that level. + const G4int level = envelopeDepth(track); + const G4VTouchable* touchable = track->GetTouchable(); + const G4AffineTransform& toLocal = + touchable->GetHistory()->GetTransform(touchable->GetHistoryDepth() - level); + const G4VSolid* envelopeSolid = touchable->GetVolume(level)->GetLogicalVolume()->GetSolid(); + + input.exitDistance = envelopeSolid->DistanceToOut(toLocal.TransformPoint(position), + toLocal.TransformAxis(direction)) / + CLHEP::cm; + + const std::vector outgoing = sample(input); + + fastStep.KillPrimaryTrack(); + fastStep.ProposePrimaryTrackPathLength(input.exitDistance * CLHEP::cm); + + // NOTE: a fast step defaults to AvoidHitInvocation, so Geant4 does not call + // the sensitive detector and TVirtualMCApplication::Stepping() is not invoked + // for it. For a passive envelope that is what we want -- there are no hits to + // lose, and the steps disappearing from the step log is the saving. A model + // covering a region that scores would add + // fastStep.ProposeSteppingControl(NormalCondition); + // here. + + double outgoingKineticEnergy = 0.; + fastStep.SetNumberOfSecondaryTracks(outgoing.size()); + for (const auto& out : outgoing) { + const G4ParticleDefinition* definition = + G4ParticleTable::GetParticleTable()->FindParticle(out.pdg); + if (definition == nullptr) { + LOG(error) << "fast simulation: model " << GetName() << " returned unknown pdg " << out.pdg + << "; particle dropped"; + continue; + } + const G4ThreeVector momentum(out.momentum[0] * CLHEP::GeV, out.momentum[1] * CLHEP::GeV, + out.momentum[2] * CLHEP::GeV); + G4DynamicParticle particle(definition, momentum); + outgoingKineticEnergy += particle.GetKineticEnergy(); + fastStep.CreateSecondaryTrack(particle, + G4ThreeVector(out.position[0] * CLHEP::cm, + out.position[1] * CLHEP::cm, + out.position[2] * CLHEP::cm), + out.time * CLHEP::ns, /*localCoordinates=*/false); + } + + // Whatever did not come out stayed in. + fastStep.ProposeTotalEnergyDeposited( + std::max(0., track->GetKineticEnergy() - outgoingKineticEnergy)); +} + +} // namespace o2::fastsim diff --git a/Detectors/FastSim/src/FastSimRegions.cxx b/Detectors/FastSim/src/FastSimRegions.cxx new file mode 100644 index 0000000000000..5179fd458c619 --- /dev/null +++ b/Detectors/FastSim/src/FastSimRegions.cxx @@ -0,0 +1,117 @@ +// 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 "FastSim/FastSimRegions.h" + +#include "TG4GeometryManager.h" +#include "TG4ModelConfigurationManager.h" + +#include +#include +#include +#include + +#include + +namespace o2::fastsim +{ + +namespace +{ +void collect(TGeoVolume* volume, std::set& media, std::set& seen) +{ + if (volume == nullptr || !seen.insert(volume).second) { + return; // a volume can be placed many times; visit it once + } + if (const TGeoMedium* medium = volume->GetMedium()) { + media.insert(medium->GetName()); + } + TObjArray* nodes = volume->GetNodes(); + if (nodes == nullptr) { + return; + } + for (int i = 0; i < nodes->GetEntriesFast(); ++i) { + auto* node = static_cast(nodes->UncheckedAt(i)); + if (node != nullptr) { + collect(node->GetVolume(), media, seen); + } + } +} +} // namespace + +//_____________________________________________________________________________ +std::set mediaInSubtree(TGeoVolume* volume) +{ + std::set media; + std::set seen; + collect(volume, media, seen); + return media; +} + +//_____________________________________________________________________________ +std::set mediaInSubtree(const std::string& volumeName) +{ + if (gGeoManager == nullptr) { + LOG(error) << "fast simulation: no TGeo geometry when resolving '" << volumeName << "'"; + return {}; + } + auto* volume = gGeoManager->GetVolume(volumeName.c_str()); + if (volume == nullptr) { + LOG(error) << "fast simulation: no volume named '" << volumeName << "' in the geometry"; + return {}; + } + return mediaInSubtree(volume); +} + +//_____________________________________________________________________________ +FastSimRegionConstruction::FastSimRegionConstruction(std::vector models) + : mModels(std::move(models)) +{ +} + +//_____________________________________________________________________________ +void FastSimRegionConstruction::Construct() +{ + /// Called by Geant4-VMC after the geometry is built and immediately before + /// the media are turned into regions, which is the only moment at which this + /// can be done: the geometry does not exist when the model is created, and + /// the regions are fixed once they are made. + auto* manager = TG4GeometryManager::Instance()->GetFastModelsManager(); + + for (const auto& model : mModels) { + std::string regions = model.regions; + if (regions.empty()) { + const auto media = mediaInSubtree(model.envelope); + for (const auto& medium : media) { + // The setter tokenizes on whitespace, so a medium whose name contains a + // space cannot go through it. O2 composes medium names as + // _ and none contain spaces, but say so if that changes + // rather than silently selecting the wrong thing. + if (medium.find(' ') != std::string::npos) { + LOG(warn) << "fast simulation: medium '" << medium << "' contains a space and cannot " + << "be selected; it is skipped"; + continue; + } + regions += (regions.empty() ? "" : " ") + medium; + } + LOG(info) << "fast simulation: model " << model.model << " covers " << media.size() + << " media found under '" << model.envelope << "'"; + } + if (regions.empty()) { + LOG(error) << "fast simulation: model " << model.model << " ended up with no regions"; + continue; + } + LOG(debug) << "fast simulation: regions for " << model.model << ": " << regions; + manager->SetModelRegions(model.model, regions); + } +} + +} // namespace o2::fastsim diff --git a/Detectors/FastSim/src/G4FastSimulation.cxx b/Detectors/FastSim/src/G4FastSimulation.cxx new file mode 100644 index 0000000000000..6349084fdf04d --- /dev/null +++ b/Detectors/FastSim/src/G4FastSimulation.cxx @@ -0,0 +1,97 @@ +// 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 "FastSim/G4FastSimulation.h" +#include "FastSim/FastSimRegions.h" +#include "FastSim/ToyAbsorberFastSim.h" +#include "SimConfig/G4Params.h" + +#include + +#include + +namespace o2::fastsim +{ + +namespace +{ +std::vector split(const std::string& value, char sep) +{ + std::vector out; + std::stringstream stream(value); + std::string token; + while (std::getline(stream, token, sep)) { + if (!token.empty()) { + out.push_back(token); + } + } + return out; +} +} // namespace + +//_____________________________________________________________________________ +G4FastSimulation::G4FastSimulation(std::vector models, + const std::string& envelope, double minEnergyGeV) + : TG4VUserFastSimulation(), mModels(std::move(models)), mEnvelope(envelope), mMinEnergy(minEnergyGeV) +{ + // Only the model itself can be declared here: this constructor runs before the + // geometry exists, so the regions cannot be derived yet. They are set later by + // FastSimRegionConstruction, in the window between geometry construction and + // region creation. + for (const auto& model : mModels) { + SetModel(model); + SetModelParticles(model, "all"); + } +} + +//_____________________________________________________________________________ +void G4FastSimulation::Construct() +{ + for (const auto& model : mModels) { + if (model == "toyAbsorber") { + LOG(info) << "fast simulation: registering model " << model << " on envelope '" << mEnvelope + << "' above " << mMinEnergy << " GeV"; + Register(new ToyAbsorberFastSim(model, mEnvelope, mMinEnergy)); + } else { + LOG(error) << "fast simulation: unknown model " << model << "; ignored"; + } + } +} + +//_____________________________________________________________________________ +TG4VUserFastSimulation* G4RunConfiguration::CreateUserFastSimulation() +{ + const auto& params = o2::conf::G4Params::Instance(); + auto models = split(params.fastSimModels, ','); + if (models.empty()) { + return nullptr; // the default: no fast simulation, unchanged behaviour + } + LOG(info) << "fast simulation is ENABLED on envelope '" << params.fastSimEnvelope << "'"; + return new G4FastSimulation(std::move(models), params.fastSimEnvelope, params.fastSimMinEnergy); +} + +//_____________________________________________________________________________ +TG4VUserPostDetConstruction* G4RunConfiguration::CreateUserPostDetConstruction() +{ + const auto& params = o2::conf::G4Params::Instance(); + auto models = split(params.fastSimModels, ','); + if (models.empty()) { + return TG4RunConfiguration::CreateUserPostDetConstruction(); + } + std::vector wanted; + wanted.reserve(models.size()); + for (auto& model : models) { + wanted.push_back({model, params.fastSimEnvelope, params.fastSimRegions}); + } + return new FastSimRegionConstruction(std::move(wanted)); +} + +} // namespace o2::fastsim diff --git a/Detectors/FastSim/src/ToyAbsorberFastSim.cxx b/Detectors/FastSim/src/ToyAbsorberFastSim.cxx new file mode 100644 index 0000000000000..f91b85cc141f2 --- /dev/null +++ b/Detectors/FastSim/src/ToyAbsorberFastSim.cxx @@ -0,0 +1,47 @@ +// 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 "FastSim/ToyAbsorberFastSim.h" + +#include + +namespace o2::fastsim +{ + +namespace +{ +/// Attenuation length of the toy transformation. +constexpr double kAbsorptionLengthCm = 60.; +} // namespace + +//_____________________________________________________________________________ +std::vector ToyAbsorberFastSim::sample(const FastSimInput& input) const +{ + // A toy transformation, not a physics model: the incident particle carries on + // in its direction with the energy attenuated over the path through the + // envelope. A trained model returns a shower here instead. + // Always positive: ModelTrigger only calls a model above its threshold, and + // an exponential of a finite path cannot reach zero. + const double kinetic = input.kineticEnergy * std::exp(-input.exitDistance / kAbsorptionLengthCm); + const double momentum = std::sqrt(kinetic * (kinetic + 2. * input.mass)); + + FastSimOutput out; + out.pdg = input.pdg; + out.time = input.time; + for (int i = 0; i < 3; ++i) { + out.position[i] = + input.position[i] + (input.exitDistance + kSurfaceEpsilonCm) * input.direction[i]; + out.momentum[i] = momentum * input.direction[i]; + } + return {out}; +} + +} // namespace o2::fastsim diff --git a/Detectors/Passive/src/Absorber.cxx b/Detectors/Passive/src/Absorber.cxx index 97b091f5965a6..4798734cbf043 100644 --- a/Detectors/Passive/src/Absorber.cxx +++ b/Detectors/Passive/src/Absorber.cxx @@ -171,6 +171,21 @@ void Absorber::createMaterials() matmgr.Medium("ABSO", 35, "AIR_C1", 35, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); matmgr.Medium("ABSO", 55, "AIR_C2", 55, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); + // + // Air of the absorber envelope + // + // Chemically identical to AIR0$ above. It exists so that the absorber + // mother volume AFaM has a material no other volume shares: Geant4-VMC + // selects fast-simulation regions by MATERIAL name and adds every volume of + // that material to the region, so a volume can only be addressed as a + // region of its own if its material is its own. Naming ABSO_AIR_ENVELOPE + // therefore means exactly "the front absorber", with all of its daughters + // inside it. Cuts and processes are the global defaults, the same ones + // ABSO_AIR_C0 gets, so the physics is unchanged. + matmgr.Mixture("ABSO", 20, "AIR_ENVELOPE0$", aAir, zAir, dAir, 4, wAir); + matmgr.Medium("ABSO", 20, "AIR_ENVELOPE", 20, 0, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, + stmin); + // // Vacuum matmgr.Mixture("ABSO", 16, "VACUUM0$", aAir, zAir, dAir1, 4, wAir); @@ -249,6 +264,8 @@ void Absorber::ConstructGeometry() auto kMedSteelSh = matmgr.getTGeoMedium("ABSO_ST_C3"); // auto kMedAir = matmgr.getTGeoMedium("ABSO_AIR_C0"); + // Air again, under a name only AFaM uses -- see the comment at its definition. + auto kMedAirEnvelope = matmgr.getTGeoMedium("ABSO_AIR_ENVELOPE"); // auto kMedPb = matmgr.getTGeoMedium("ABSO_PB_C0"); auto kMedPbSh = matmgr.getTGeoMedium("ABSO_PB_C2"); @@ -867,7 +884,9 @@ void Absorber::ConstructGeometry() shFaM->DefineSection(14, z, rInFaCH2Cone2 - dz * angle10, rOuSteelEnvelopeR2); z += dzSteelEnvelopeR / 2.; shFaM->DefineSection(15, z, rInFaCH2Cone2, rOuSteelEnvelopeR2); - TGeoVolume* voFaM = new TGeoVolume("AFaM", shFaM, kMedAir); + // AFaM is the mother of the whole absorber, and its dedicated medium is what + // makes "the absorber" addressable as one fast-simulation region. + TGeoVolume* voFaM = new TGeoVolume("AFaM", shFaM, kMedAirEnvelope); voFaM->SetVisibility(0); // diff --git a/Detectors/gconfig/CMakeLists.txt b/Detectors/gconfig/CMakeLists.txt index be50c9bcfe72f..282fa4c9d124e 100644 --- a/Detectors/gconfig/CMakeLists.txt +++ b/Detectors/gconfig/CMakeLists.txt @@ -16,7 +16,7 @@ o2_add_library(G3Setup o2_add_library(G4Setup SOURCES src/G4Config.cxx - PUBLIC_LINK_LIBRARIES MC::Geant4VMC MC::Geant4 FairRoot::Base O2::SimulationDataFormat O2::Generators O2::SimSetup + PUBLIC_LINK_LIBRARIES MC::Geant4VMC MC::Geant4 FairRoot::Base O2::SimulationDataFormat O2::Generators O2::SimSetup O2::FastSim ) o2_add_library(FLUKASetup diff --git a/Detectors/gconfig/g4Config.C b/Detectors/gconfig/g4Config.C index 16374cf9fd4a3..77494c559ca51 100644 --- a/Detectors/gconfig/g4Config.C +++ b/Detectors/gconfig/g4Config.C @@ -61,6 +61,7 @@ R__LOAD_LIBRARY(libgeant4vmc) #include "TG4RunConfiguration.h" #include "SimConfig/G4Params.h" #include "SimConfig/FluenceWeightCalculator.h" +#include "FastSim/G4FastSimulation.h" #endif #include "commonConfig.C" @@ -114,8 +115,12 @@ void Config() LOG(fatal) << "Unsupported geometry navigation mode"; } - auto runConfiguration = new TG4RunConfiguration(geomNavStr, physicsSetup, "stepLimiter+specialCuts", - specialStacking, mtMode); + // o2::fastsim::G4RunConfiguration differs from TG4RunConfiguration only in + // providing the fast-simulation hook; with G4.fastSimModels empty it behaves + // identically. + auto runConfiguration = new o2::fastsim::G4RunConfiguration(geomNavStr, physicsSetup, + "stepLimiter+specialCuts", + specialStacking, mtMode); if (g4Params.g4scoring) { runConfiguration->SetUseOfG4Scoring(); if (g4Params.g4fluenceweight) { diff --git a/run/SimExamples/FastSim_Absorber/README.md b/run/SimExamples/FastSim_Absorber/README.md new file mode 100644 index 0000000000000..3a7c51a717e8c --- /dev/null +++ b/run/SimExamples/FastSim_Absorber/README.md @@ -0,0 +1,106 @@ +# Fast simulation of the front absorber + +Replaces the detailed transport through the ALICE front absorber by a model, and +compares the result against a full simulation of the same events. + +`run.sh` runs both: a reference `o2-sim` with PIPE and ABSO only, and the same +setup with the fast simulation switched on. Both use the same seed, so the two +simulations see the same primaries and the comparison is paired. + +| file | | +|---|---| +| `run.sh` | the two simulations and the comparison | +| `countTracks.macro` | tracks per event of an `o2-sim` output | + +## Switching it on + +The feature does nothing unless a model is named: + +``` +--configKeyValues "G4.fastSimModels=toyAbsorber;G4.fastSimEnvelope=AFaM" +``` + +`G4.fastSimEnvelope` is the volume the model stands in for — here `AFaM`, the +mother of the whole absorber. `G4.fastSimMinEnergy` (GeV, default 1) is the +threshold below which the detailed transport still runs. + +## What the model does + +`toyAbsorber` is a placeholder: it returns the incident particle continuing in +its direction with the energy attenuated exponentially over its path through the +envelope. A real absorber turns one incident hadron into a shower, so these are +not physics numbers. + +A model implements one function, `sample()`, which maps the particle that +entered the region to the particles that leave it +(`Detectors/FastSim/include/FastSim/FastSimModel.h`). The surrounding work — +measuring the distance to the envelope surface, killing the incident particle, +stacking what comes back, booking the energy difference as a deposit — is shared. + +## How the envelope and the regions relate + +Geant4 consults a fast simulation model through *regions*, and in O2 a region can +only ever be "every volume of a given material": the VMC special cuts make every +logical volume a root of its own material's region, and Geant4 stops propagating +a region at any such daughter. A region named after `AFaM` would therefore +contain `AFaM` alone, which tracks skip entirely because its daughters touch its +surface. + +So the two jobs are separated. The regions are derived by walking the envelope's +subtree and collecting its media, which for `AFaM` finds fourteen: + +``` +fast simulation: model toyAbsorber covers 14 media found under 'AFaM' +Adding fast simulation model toyAbsorber to regions ABSO_AIR0$ ABSO_AIR_ENVELOPE0$ + ABSO_CONCRETE2$ ABSO_POLYETHYLEN2$ ABSO_CARBON0$ ABSO_CARBON2$ ABSO_MAGNESIUM$ + ABSO_Ni-W-Cu0$ ABSO_Ni-W-Cu2$ ABSO_LEAD0$ ABSO_LEAD2$ ABSO_STAINLESS STEEL0$ + ABSO_STAINLESS STEEL2$ +``` + +What the model measures against is the envelope itself, read from the track's own +touchable. `ModelTrigger` also requires geometric containment in it, which is +what excludes the steel support cradle — it shares its material with the end +plate, so no selection by material could separate them. + +## Measured + +Five pp minimum-bias events, Pythia8 and Geant4, `-m PIPE ABSO`, same seed, one +EPN node against `O2PDPSuite/daily-20260819-0000-1`: + +| | full | fast | +|---|---|---| +| tracks per event | 1547 | 1144 | +| transport real time | 14.9 s | 14.2 s | + +A 20 GeV muon fired into the muon-arm acceptance with `/tracking/verbose 1` +shows what the model does to a single track: + +``` + 10 54.4 -0.137 -900 1.99e+04 0.0749 287 902 AFaMgRing Transportation + 11 54.4 -0.137 -900 0 1.99e+04 4.1e+03 5e+03 AFaMgRing G4FastSimulationManagerProcess +``` + +One step of 4.1 m across the whole absorber, in place of roughly twenty `muIoni` +steps through graphite, concrete and steel. + +## What these numbers are not + +**Not a performance result.** A quarter fewer tracks buys only five percent of +wall clock, because the tracks the absorber's shower contributes are cheap +low-energy ones and most of the CPU in this setup is spent elsewhere. And with +minimum-bias pp only the forward cone reaches the absorber at all. A CPU number +worth quoting needs a workload where the absorber is on the critical path — a +forward-biased generator, or the full detector where the muon arm is the point. + +**Not a physics validation.** The toy returns one particle where a real absorber +returns a shower, so the track counts above say more about the placeholder than +about the absorber. Comparing the outgoing multiplicity and spectrum against a +full simulation is what a trained model has to pass, and that is the measurement +this example is scaffolding for. + +Two further caveats worth knowing when reading any output of this: a fast step +does not call the sensitive detector, so its steps disappear from MCStepLogger +(harmless for a passive envelope, which has no hits); and secondaries the model +creates carry `TMCProcess` `kPNull`, the code geant4_vmc gives to everything it +has no VMC equivalent for, so they cannot be told apart from other tracks by +process alone. diff --git a/run/SimExamples/FastSim_Absorber/countTracks.macro b/run/SimExamples/FastSim_Absorber/countTracks.macro new file mode 100644 index 0000000000000..0c5e10106aec4 --- /dev/null +++ b/run/SimExamples/FastSim_Absorber/countTracks.macro @@ -0,0 +1,20 @@ +// Tracks per event of an o2-sim output, as a cheap proxy for how much transport +// happened. Usage: root -l -b -q 'countTracks.macro("o2sim")' +void countTracks(const char* prefix = "o2sim") +{ + TFile file(Form("%s_Kine.root", prefix)); + auto* tree = (TTree*)file.Get("o2sim"); + if (!tree) { + printf("no o2sim tree in %s_Kine.root\n", prefix); + return; + } + std::vector* tracks = nullptr; + tree->SetBranchAddress("MCTrack", &tracks); + Long64_t total = 0; + for (Long64_t entry = 0; entry < tree->GetEntries(); ++entry) { + tree->GetEntry(entry); + total += tracks->size(); + } + printf("%s: %lld events, %lld tracks, %.1f per event\n", prefix, tree->GetEntries(), total, + tree->GetEntries() ? double(total) / tree->GetEntries() : 0.); +} diff --git a/run/SimExamples/FastSim_Absorber/run.sh b/run/SimExamples/FastSim_Absorber/run.sh new file mode 100755 index 0000000000000..057b649fbc391 --- /dev/null +++ b/run/SimExamples/FastSim_Absorber/run.sh @@ -0,0 +1,62 @@ +#!/usr/bin/env bash +# +# This example exercises the Geant4 fast simulation route on the front absorber. +# +# A fast simulation model replaces the detailed transport through a region of +# the geometry by a function from the particle that enters it to the particles +# that leave it. The model shipped here, `toyAbsorber`, is a placeholder: it +# returns the incident particle carrying on in its direction with an +# exponentially attenuated energy. What the example demonstrates is the +# machinery, not the physics. +# +# `G4.fastSimEnvelope` names the volume the model stands in for: AFaM, the +# mother of the whole absorber. The regions Geant4 needs in order to consult the +# model are derived from it by walking its subtree and collecting the media, +# because a region can only ever be "every volume of a given material". +# +# The setup is PIPE and ABSO only, which keeps the run short and puts the +# absorber in the path of everything. + +set -x + +EVENTS=5 +MODULES="-m PIPE ABSO" +GEN="-g pythia8pp" +# The same seed in both runs, so the two simulations see the SAME primaries and +# the comparison is paired. Without it o2-sim picks a seed per run and the +# difference between the two is mostly different events. +SEED="--seed 12345" +# Alignment is irrelevant here and switching it off keeps the example from +# needing a CCDB connection and an alien token. +COMMON="align-geom.mDetectors=none" + +# --------------------------------------------------------------- 1. reference +# Detailed transport, for comparison. +mkdir -p full && cd full +o2-sim-serial -n ${EVENTS} ${GEN} ${SEED} -e TGeant4 ${MODULES} -o full \ + --configKeyValues "${COMMON}" > logfull 2>&1 +cd .. + +# -------------------------------------------------------------------- 2. fast +# G4.fastSimModels is what switches the feature on; with it empty (the default) +# nothing about the simulation changes. +mkdir -p fast && cd fast +o2-sim-serial -n ${EVENTS} ${GEN} ${SEED} -e TGeant4 ${MODULES} -o fast \ + --configKeyValues "${COMMON};\ +G4.fastSimModels=toyAbsorber;\ +G4.fastSimEnvelope=AFaM;\ +G4.fastSimMinEnergy=1.0" > logfast 2>&1 +cd .. + +# ----------------------------------------------------------------- 3. compare +# The model prints once at setup; if this line is missing the region name did +# not resolve to a medium and nothing was applied. +grep -h "fast simulation" fast/logfast + +# Tracks written per event. The absorber normally turns one incident hadron into +# a shower, and the toy model returns a single particle instead, so the fast run +# has to produce far fewer. +for d in full fast; do + echo "=== ${d}" + root -l -b -q "$(dirname "$0")/countTracks.macro(\"${d}/${d}\")" +done