Trigger a fast simulation model on its envelope, not on its region

This makes a model cover a whole module, which the region mechanism alone
cannot express.

- A Geant4 region in O2 is always "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.
- G4FastTrack::GetEnvelopeSolid() therefore returns one absorber piece, and
  FastSimModel no longer uses it.
- G4.fastSimEnvelope names the volume a model stands in for. Containment and the
  exit distance are read from the track's own touchable, which already carries
  the full ancestry and the transform of every level.
- ModelTrigger now requires geometric containment in that envelope, which also
  excludes volumes that merely share a material with it, such as the absorber's
  support cradle.
- The regions needed for Geant4 to consult the model are derived by walking the
  envelope's subtree and collecting its media, so no list is maintained by hand.
  The walk happens in CreateUserPostDetConstruction, the one point after the
  geometry is built and before the media are turned into regions.
- G4.fastSimRegions still overrides the walk with an explicit list.

Author: sawenzel

Common/SimConfig/include/SimConfig/G4Params.h

@@ -58,7 +58,9 @@ struct G4Params : public o2::conf::ConfigurableParamHelper<G4Params> {
   // 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 fastSimRegions = "";  // tracking media the models apply to (wildcards allowed)
+  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");
 };

Detectors/FastSim/CMakeLists.txt

@@ -11,6 +11,7 @@
 
 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

Detectors/FastSim/include/FastSim/FastSimModel.h

@@ -39,7 +39,7 @@ struct FastSimInput {
   double kineticEnergy = 0.; ///< GeV
   double mass = 0.;          ///< GeV
   double time = 0.;          ///< ns
-  double exitDistance = 0.;  ///< cm from `position` to the envelope surface along `direction`
+  double exitDistance = 0.;  ///< cm from `position` to the ENVELOPE surface along `direction`
 };
 
 /// One particle leaving the envelope.
@@ -55,10 +55,23 @@ struct FastSimOutput {
 /// 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 measures against 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, double minEnergyGeV);
+  FastSimModel(const G4String& name, const G4String& envelopeVolume, double minEnergyGeV);
 
   G4bool IsApplicable(const G4ParticleDefinition& particle) override;
   G4bool ModelTrigger(const G4FastTrack& fastTrack) override;
@@ -74,7 +87,13 @@ class FastSimModel : public G4VFastSimulationModel
   virtual std::vector<FastSimOutput> 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 measures against
   double mMinEnergy = 0.; ///< internal Geant4 units; below this the detailed transport runs
+  mutable bool mWarned = false;
 };
 
 } // namespace o2::fastsim

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 <set>
+#include <string>
+#include <vector>
+
+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<std::string> mediaInSubtree(TGeoVolume* volume);
+
+/// Same, looked up by volume name in the current TGeo geometry. Empty if there
+/// is no such volume.
+std::set<std::string> 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<ModelRegions> models);
+  void Construct() override;
+
+ private:
+  std::vector<ModelRegions> mModels;
+};
+
+} // namespace o2::fastsim
+
+#endif // O2_FASTSIM_REGIONS_H_

Detectors/FastSim/include/FastSim/G4FastSimulation.h

@@ -18,17 +18,20 @@
 ///
 ///   o2-sim -n 10 -g pythia8pp -e TGeant4 -m PIPE ABSO
 ///          --configKeyValues "G4.fastSimModels=toyAbsorber;
-///                             G4.fastSimRegions=ABSO_AIR_ENVELOPE"
+///                             G4.fastSimEnvelope=AFaM"
 ///
-/// Regions are selected by TRACKING MEDIUM name (wildcards allowed), which
-/// Geant4-VMC maps to the medium's MATERIAL; every volume of that material joins
-/// the region. That is why the absorber mother volume AFaM carries a material of
-/// its own (see Detectors/Passive/src/Absorber.cxx). Selection by volume is not
-/// available: the VMC special cuts already root every logical volume in a
-/// per-material region, and Geant4 allows a volume in only one region.
+/// `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 <string>
 #include <vector>
@@ -40,12 +43,13 @@ namespace o2::fastsim
 class G4FastSimulation : public TG4VUserFastSimulation
 {
  public:
-  G4FastSimulation(std::vector<std::string> models, const std::string& regions,
+  G4FastSimulation(std::vector<std::string> models, const std::string& envelope,
                    double minEnergyGeV);
   void Construct() override;
 
  private:
   std::vector<std::string> mModels;
+  std::string mEnvelope;
   double mMinEnergy = 1.;
 };
 
@@ -56,6 +60,7 @@ class G4RunConfiguration : public TG4RunConfiguration
  public:
   using TG4RunConfiguration::TG4RunConfiguration;
   TG4VUserFastSimulation* CreateUserFastSimulation() override;
+  TG4VUserPostDetConstruction* CreateUserPostDetConstruction() override;
 };
 
 } // namespace o2::fastsim

Detectors/FastSim/src/FastSimModel.cxx

@@ -21,19 +21,51 @@
 #include <G4SystemOfUnits.hh>
 #include <G4ThreeVector.hh>
 #include <G4Track.hh>
+#include <G4NavigationHistory.hh>
+#include <G4TouchableHistory.hh>
+#include <G4VPhysicalVolume.hh>
 #include <G4VSolid.hh>
+#include <G4VTouchable.hh>
 
 #include <algorithm>
 
 namespace o2::fastsim
 {
 
 //_____________________________________________________________________________
-FastSimModel::FastSimModel(const G4String& name, double minEnergyGeV)
-  : G4VFastSimulationModel(name), mMinEnergy(minEnergyGeV * CLHEP::GeV)
+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&)
 {
@@ -45,9 +77,28 @@ G4bool FastSimModel::IsApplicable(const G4ParticleDefinition&)
 //_____________________________________________________________________________
 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.
-  return fastTrack.GetPrimaryTrack()->GetKineticEnergy() > mMinEnergy;
+  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;
 }
 
 //_____________________________________________________________________________
@@ -68,9 +119,17 @@ void FastSimModel::DoIt(const G4FastTrack& fastTrack, G4FastStep& fastStep)
   input.kineticEnergy = track->GetKineticEnergy() / CLHEP::GeV;
   input.mass = track->GetDefinition()->GetPDGMass() / CLHEP::GeV;
   input.time = track->GetGlobalTime() / CLHEP::ns;
-  input.exitDistance = fastTrack.GetEnvelopeSolid()->DistanceToOut(
-                         fastTrack.GetPrimaryTrackLocalPosition(),
-                         fastTrack.GetPrimaryTrackLocalDirection()) /
+  // Deliberately NOT GetEnvelopeSolid(): that is the region's root volume, i.e.
+  // one absorber piece. Measure against the envelope volume instead, using 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<FastSimOutput> outgoing = sample(input);