Integrate geometry step killing

Common/SimConfig/include/SimConfig/SimParams.h

@@ -23,6 +23,7 @@ namespace conf
 // (mostly used in O2MCApplication stepping)
 struct SimCutParams : public o2::conf::ConfigurableParamHelper<SimCutParams> {
   bool stepFiltering = true; // if we activate the step filtering in O2BaseMCApplication
+  std::string stepFilteringMacro = "";                                                                        // ROOT macro providing keepStep(); empty = built-in z/R cut
   bool stepTrackRefHook = false;                                                                              // if we create track references during generic stepping
   std::string stepTrackRefHookFile = "${O2_ROOT}/share/Detectors/gconfig/StandardSteppingTrackRefHook.macro"; // the standard code holding the TrackRef callback
 

Detectors/gconfig/CMakeLists.txt

@@ -65,4 +65,5 @@ o2_add_test_root_macro(g3Config.C
 
 
 o2_data_file(COPY data  DESTINATION Detectors/gconfig/)
-install(FILES src/StandardSteppingTrackRefHook.macro src/FlukaRuntimeConfig.macro DESTINATION share/Detectors/gconfig/)
+install(FILES src/StandardSteppingTrackRefHook.macro src/FlukaRuntimeConfig.macro
+              src/KeepStepCylinders.macro DESTINATION share/Detectors/gconfig/)

Detectors/gconfig/include/SimSetup/MCReplayParam.h

@@ -29,6 +29,7 @@ struct MCReplayParam : public o2::conf::ConfigurableParamHelper<MCReplayParam> {
   std::string stepFilename = "MCStepLoggerOutput.root"; // filename where to find the stepTreename
   float energyCut = -1.;                                // minimum energy required for a step to continue tracking
   std::string cutFile = "";
+  bool allowStopTrack = false;
   O2ParamDef(MCReplayParam, "MCReplayParam");
 };
 } // end namespace o2

Detectors/gconfig/src/KeepStepCylinders.macro

@@ -0,0 +1,26 @@
+// Generated by run/SimExamples/Geometry_StepFiltering/makeKeepStepCylinders.macro
+// -- do not edit by hand.
+// Geometry: o2sim_geometry.root
+// Envelope: max radius of sensitive volumes and of material with rho > 0.01 g/cm3,
+//           excluding modules HALL,CAVE,
+//           sampled at 30000 z slices x 128 phi directions, margin 2% + 2 cm.
+// Steps outside the z range below are always kept.
+
+o2::steer::O2MCApplicationBase::KeepStepFcn keepStep()
+{
+  const float rSq[] = {50.0f * 50.0f, 176.0f * 176.0f, 437.0f * 437.0f, 51.0f * 51.0f, 420.0f * 420.0f, 48.0f * 48.0f, 378.0f * 378.0f, 576.0f * 576.0f, 330.0f * 330.0f, 899.0f * 899.0f, 97.0f * 97.0f, 50.0f * 50.0f};
+  const float edges[] = {-15000.0f, -1890.0f, -1740.0f, -1690.0f, -1640.0f, -1390.0f, -1330.0f, -1200.0f, -830.0f, -720.0f, 720.0f, 1200.0f, 15000.0f};
+  return [rSq, edges](TVirtualMC const* mc) {
+    float x, y, z;
+    mc->TrackPosition(x, y, z);
+    if (z < edges[0] || z >= edges[12]) {
+      return true; // beyond the traced region, e.g. the ZDC tunnel
+    }
+    for (auto i = 0U; i < 12; ++i) {
+      if (edges[i + 1] > z && z >= edges[i]) {
+        return (x * x + y * y) < rSq[i];
+      }
+    }
+    return true;
+  };
+}

Detectors/gconfig/src/MCReplayConfig.cxx

@@ -39,6 +39,7 @@ void Config()
   replay->SetCut("CUTALLE", params.energyCut);
   replay->cutsFromConfig(params.cutFile);
   replay->blockSetProcessesCuts();
+  replay->allowStopTrack(params.allowStopTrack);
 }
 
 void MCReplayConfig()

Steer/include/Steer/O2MCApplicationBase.h

@@ -35,9 +35,14 @@ namespace steer
 class O2MCApplicationBase : public FairMCApplication
 {
  public:
-  O2MCApplicationBase() : FairMCApplication(), mCutParams(o2::conf::SimCutParams::Instance()) { initTrackRefHook(); }
+  O2MCApplicationBase() : FairMCApplication(), mCutParams(o2::conf::SimCutParams::Instance())
+  {
+    initStepFilterHook();
+    initTrackRefHook();
+  }
   O2MCApplicationBase(const char* name, const char* title, TObjArray* ModList, const char* MatName) : FairMCApplication(name, title, ModList, MatName), mCutParams(o2::conf::SimCutParams::Instance())
   {
+    initStepFilterHook();
     initTrackRefHook();
   }
 
@@ -57,6 +62,7 @@ class O2MCApplicationBase : public FairMCApplication
   double TrackingZmax() const override { return mCutParams.maxAbsZTracking; }
 
   typedef std::function<void(TVirtualMC const*)> TrackRefFcn;
+  typedef std::function<bool(TVirtualMC const*)> KeepStepFcn;
 
   void fixTGeoRuntimeShapes();
 
@@ -77,6 +83,11 @@ class O2MCApplicationBase : public FairMCApplication
   void finishEventCommon();
   TrackRefFcn mTrackRefFcn; // a function hook that gets (optionally) called during Stepping
   void initTrackRefHook();
+  /// an optional extra per-step criterion, loaded from
+  /// SimCutParams.stepFilteringMacro; only consulted if mHasStepFilterMacro
+  KeepStepFcn mKeepStepFcn;
+  bool mHasStepFilterMacro = false;
+  void initStepFilterHook();
 
   ClassDefOverride(O2MCApplicationBase, 1);
 };

Steer/src/O2MCApplication.cxx

@@ -115,6 +115,14 @@ void O2MCApplicationBase::Stepping()
     }
   }
 
+  // an additional, user-provided criterion; only consulted when one is
+  // configured, so that SimCutParams.stepFilteringMacro being unset leaves the
+  // code above as the whole of the geometry cut
+  if (mHasStepFilterMacro && !mKeepStepFcn(fMC)) {
+    fMC->StopTrack();
+    return;
+  }
+
   if (mCutParams.stepTrackRefHook) {
     mTrackRefFcn(fMC);
   }
@@ -239,6 +247,28 @@ void O2MCApplicationBase::ConstructGeometry()
   }
 }
 
+void O2MCApplicationBase::initStepFilterHook()
+{
+  if (mCutParams.stepFilteringMacro.empty()) {
+    return;
+  }
+  const auto macro = o2::utils::expandShellVarsInFileName(mCutParams.stepFilteringMacro);
+  if (!std::filesystem::exists(macro)) {
+    LOG(error) << "Macro for step filtering does not exist at " << macro << "; ignoring it";
+    return;
+  }
+  LOG(info) << "Initializing step filtering from macro " << macro;
+  mKeepStepFcn = o2::conf::GetFromMacro<KeepStepFcn>(macro, "keepStep()",
+                                                     "o2::steer::O2MCApplicationBase::KeepStepFcn",
+                                                     "o2mc_stepping_keep_step");
+  if (!mKeepStepFcn) {
+    LOG(error) << "Could not set up keepStep() from " << macro << "; ignoring it";
+    return;
+  }
+  mHasStepFilterMacro = true;
+  LOG(info) << "Step filtering initialized from macro " << macro;
+}
+
 void O2MCApplicationBase::InitGeometry()
 {
   // load special cuts which might be given from the outside first.

run/SimExamples/Geometry_StepFiltering/README.md

@@ -0,0 +1,76 @@
+# Geometry step filtering
+
+Stops transport of a track once it steps outside a hull wrapped around the
+detector, and derives that hull from the geometry.
+
+`run.sh` runs the whole chain: a reference simulation with the MCStepLogger, the
+derivation of a hull from its geometry, a check that every recorded hit lies
+inside that hull, and two MCReplay runs — one without cuts, one with the hull —
+whose step and hit counts can be compared.
+
+| file | |
+|---|---|
+| `makeKeepStepCylinders.macro` | derives the hull and writes a `keepStep()` macro |
+| `checkKeepStepCylinders.macro` | verifies that no recorded hit lies outside a hull |
+| `countHits.macro` | hits per detector of an `o2-sim` output |
+
+The hull is applied through `SimCutParams.stepFilteringMacro`, in addition to
+the built-in z/R cut, so it can only remove steps. A ready-made one is installed
+at `$O2_ROOT/share/Detectors/gconfig/KeepStepCylinders.macro`; regenerate it
+whenever the geometry changes.
+
+`MCReplayParam.allowStopTrack=true` is what lets a replay act on the
+`StopTrack()` calls the hull makes.
+
+Measured on 10 pp minimum-bias events, Pythia8 and Geant4, with the hull
+generated from the geometry of the run itself:
+
+| module list | steps removed | hits lost |
+|---|---|---|
+| default minus ZDC | about 5 % | 0.18 % |
+| default | about 2 % | 0.20 % |
+
+Hits exclude FT0, which draws a random number while creating hits and does not
+reproduce under replay.
+
+**Read the replay's step accounting with care.** `MCReplayParam.allowStopTrack`
+makes a replay honour *every* `StopTrack()` of a replayed step, not only the
+ones the hull causes. FT0's photocathode efficiency and HMPID's Fresnel loss
+stop tracks on a random draw, and those draws do not repeat on replay, so a
+replay with `allowStopTrack=true` and a `keepStep()` that returns `true`
+already reports about 8 % of steps skipped with no cut in play; running it with
+`--skipModules ZDC FT0 HMP` reports exactly zero. Use the replay only as a
+difference against such a no-op macro, and expect that difference to still read
+a little high, because a secondary is dropped whenever its parent was skipped
+before it was born, which approximates rather than reproduces what Geant4 does.
+
+The number to quote comes from a real Geant4 pair with `SimCutParams.trackSeed`
+on, summing the per-event `This event/chunk did N steps`. Without per-track
+seeding the two runs diverge into different physics and the comparison is
+worthless: on one sample it returned 0.19 % where the truth was 4.4 %.
+
+Do not expect the step reduction to become a CPU reduction. Without ZDC,
+measured as user time over sequential runs on an idle machine with the same
+generator seed:
+
+| run | Geant4 steps | CPU |
+|---|---|---|
+| no macro | 9 855 230 | 88.7 s |
+| macro returning `true` | 9 855 230 | 91.1 s |
+| the hull | 9 380 571 | 90.4 s |
+
+Two things follow. Reaching a cling-compiled `keepStep()` through a
+`std::function` costs about 235 ns per step on its own -- the middle row removes
+no steps at all -- while the cylinder scan itself is negligible next to it. And
+the steps the hull removes are cheap ones: the ~5 % of steps it removes are
+worth only about 0.7 % of the runtime, so even at zero hook overhead the gain
+here would be small. The macro is compiled once at start-up, not per step.
+
+Evaluating a cylinder set natively rather than through a macro would remove the
+overhead; whether a geometry hull is worth it without ZDC is a separate
+question.
+
+Background: A. Swain, *Geometric Hyperparameter Optimisation of ALICE Monte
+Carlo Transport Simulations*, CERN-STUDENTS-Note-2023-164, and B. Völkel,
+*Geometry cuts in MC transport*, WP12/13 meeting, 11.10.2023
+(https://indico.cern.ch/event/1334852/contributions/5620122/).

run/SimExamples/Geometry_StepFiltering/checkKeepStepCylinders.macro

@@ -0,0 +1,164 @@
+// Copyright 2019-2020 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 checkKeepStepCylinders.macro
+/// \brief Verify that a generated hull contains every recorded hit
+///
+/// A hull that excludes a position where a hit was produced will lose that hit
+/// in a real run. This reads the hull straight out of a generated
+/// KeepStepCylinders.macro and tests it against the hit positions of an
+/// o2-sim output, per detector. It is a necessary condition, cheap to run, and
+/// independent of the replay machinery.
+///
+///   root -l -b -q 'checkKeepStepCylinders.macro("KeepStepCylinders.macro", "o2sim")'
+
+#include <TFile.h>
+#include <TTree.h>
+#include <TSystem.h>
+#include <TSystemDirectory.h>
+#include <TString.h>
+#include <TH1.h>
+#include <TObjArray.h>
+#include <TBranch.h>
+
+#include <cstdio>
+#include <fstream>
+#include <regex>
+#include <string>
+#include <vector>
+
+namespace
+{
+/// pull the rSq[] and edges[] initialisers out of a generated macro
+bool readHull(const char* macro, std::vector<double>& rad, std::vector<double>& edges)
+{
+  std::ifstream in(macro);
+  if (!in) {
+    printf("[check] cannot open %s\n", macro);
+    return false;
+  }
+  std::string all((std::istreambuf_iterator<char>(in)), std::istreambuf_iterator<char>());
+  auto grab = [&all](const char* name, std::vector<double>& out, bool sqrtIt) {
+    std::smatch m;
+    std::regex e(std::string(name) + R"(\[\]\s*=\s*\{([^}]*)\})");
+    if (!std::regex_search(all, m, e)) {
+      return false;
+    }
+    std::string body = m[1].str();
+    std::regex num(R"(-?[0-9]+(\.[0-9]+)?)");
+    for (auto it = std::sregex_iterator(body.begin(), body.end(), num); it != std::sregex_iterator(); ++it) {
+      out.push_back(std::stod(it->str()));
+    }
+    if (sqrtIt) { // rSq is written as "R.0f * R.0f", so every value appears twice
+      std::vector<double> half;
+      for (size_t i = 0; i < out.size(); i += 2) {
+        half.push_back(out[i]);
+      }
+      out.swap(half);
+    }
+    return true;
+  };
+  return grab("rSq", rad, true) && grab("edges", edges, false);
+}
+} // namespace
+
+void checkKeepStepCylinders(const char* macro = "KeepStepCylinders.macro",
+                            const char* prefix = "o2sim")
+{
+  std::vector<double> rad, edges;
+  if (!readHull(macro, rad, edges) || edges.size() != rad.size() + 1) {
+    printf("[check] could not parse a hull out of %s\n", macro);
+    return;
+  }
+  printf("[check] hull from %s:\n", macro);
+  for (size_t k = 0; k < rad.size(); ++k) {
+    printf("   z in [%8.1f, %8.1f)  R < %6.1f\n", edges[k], edges[k + 1], rad[k]);
+  }
+
+  auto inside = [&](double x, double y, double z) {
+    if (z < edges.front() || z >= edges.back()) {
+      return true;
+    }
+    for (size_t k = 0; k < rad.size(); ++k) {
+      if (z >= edges[k] && z < edges[k + 1]) {
+        return (x * x + y * y) < rad[k] * rad[k];
+      }
+    }
+    return true;
+  };
+
+  TString dir = gSystem->DirName(prefix);
+  TString base = gSystem->BaseName(prefix);
+  TSystemDirectory sdir(dir, dir);
+  TList* files = sdir.GetListOfFiles();
+  if (!files) {
+    printf("[check] no such directory %s\n", dir.Data());
+    return;
+  }
+  long long nTot = 0, nOut = 0;
+  TIter next(files);
+  while (auto* f = (TSystemFile*)next()) {
+    TString name = f->GetName();
+    if (!name.BeginsWith(base + "_Hits") || !name.EndsWith(".root")) {
+      continue;
+    }
+    TString det = name;
+    det.ReplaceAll(base + "_Hits", "");
+    det.ReplaceAll(".root", "");
+    TFile* file = TFile::Open(dir + "/" + name);
+    TTree* t = file ? (TTree*)file->Get("o2sim") : nullptr;
+    if (!t) {
+      if (file) {
+        file->Close();
+      }
+      continue;
+    }
+    // BasicXYZxHit and friends store the position as a ROOT::Math::PositionVector3D
+    TString br;
+    for (int i = 0; i < t->GetListOfBranches()->GetEntries(); ++i) {
+      TString b = t->GetListOfBranches()->At(i)->GetName();
+      if (b.Contains("Hit")) {
+        br = b;
+        break;
+      }
+    }
+    long long n = 0, bad = 0;
+    const TString pos = br + ".mPos.fCoordinates.f";
+    if (!br.IsNull() && t->GetLeaf(pos + "X")) {
+      const TString expr = Form("%sX:%sY:%sZ", pos.Data(), pos.Data(), pos.Data());
+      // one hit is one value, not one entry: the value buffer has to be sized
+      // to the number of selected values or GetV1..3 return truncated arrays
+      Long64_t sel = t->Draw(expr, "", "goff");
+      if (sel > t->GetEstimate()) {
+        t->SetEstimate(sel + 1);
+        sel = t->Draw(expr, "", "goff");
+      }
+      const double *vx = t->GetV1(), *vy = t->GetV2(), *vz = t->GetV3();
+      for (Long64_t i = 0; i < sel; ++i) {
+        ++n;
+        if (!inside(vx[i], vy[i], vz[i])) {
+          ++bad;
+        }
+      }
+    }
+    if (n > 0) {
+      printf("[check] %-5s hits %10lld   outside hull %8lld  (%.3f%%)%s\n",
+             det.Data(), n, bad, 100. * bad / n, bad ? "   <-- HULL TOO SMALL" : "");
+      nTot += n;
+      nOut += bad;
+    } else {
+      printf("[check] %-5s no position leaf (%s), skipped\n", det.Data(), br.Data());
+    }
+    file->Close();
+  }
+  printf("[check] TOTAL hits %lld, outside hull %lld (%.4f%%)\n", nTot, nOut,
+         nTot ? 100. * nOut / nTot : 0.);
+}