From 2137c445cbfa26aa0d5e0cfb806673c7af5d5fa7 Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Mon, 17 Aug 2026 12:35:31 +0200 Subject: [PATCH 1/2] Integrate geometry step killing This commit integrates the geometry-step killing mechanism originally studied by Benedikt Voelkel and Anthony Swain. A step filter is loaded from a ROOT macro named by SimCutParams.stepFilteringMacro and evaluated for every step in addition to the built-in z/R cut, and MCReplayParam.allowStopTrack lets a replay act on it. Both are unset by default. The code is taken from https://github.com/AliceO2Group/AliceO2/pull/12119 https://github.com/benedikt-voelkel/AliceO2/tree/geom-cuts https://github.com/benedikt-voelkel/AliceO2/tree/envelope https://github.com/AnthonySwain/AliceO2/tree/ParticleSpecificMaps and the studies behind it are presented in B. Voelkel, Geometry cuts in MC transport, WP12/13, 11.10.2023 https://indico.cern.ch/event/1334852/contributions/5620122/ A. Swain, Geometric hyperparameter optimisation of ALICE MC transport https://indico.cern.ch/event/1323900/contributions/5571214/ CERN-STUDENTS-Note-2023-164 In addition, we provide run/SimExamples/Geometry_StepFiltering, which runs the whole chain from a reference simulation to a hull-filtered replay, and makeKeepStepCylinders.macro, which is able to rederive the geometry hull from a closed geometry to some extent. Thanks to @galickazuzanna for help in compiling this commit. Co-authored-by: Benedikt Volkel Co-authored-by: Anthony Swain Co-authored-by: Zuzanna Galicka Co-authored-by: Claude Opus 5 --- .../SimConfig/include/SimConfig/SimParams.h | 1 + Detectors/gconfig/CMakeLists.txt | 3 +- .../gconfig/include/SimSetup/MCReplayParam.h | 1 + Detectors/gconfig/src/KeepStepCylinders.macro | 26 ++ Detectors/gconfig/src/MCReplayConfig.cxx | 1 + Steer/include/Steer/O2MCApplicationBase.h | 13 +- Steer/src/O2MCApplication.cxx | 30 ++ .../Geometry_StepFiltering/README.md | 76 ++++ .../checkKeepStepCylinders.macro | 164 +++++++ .../Geometry_StepFiltering/countHits.macro | 75 +++ .../makeKeepStepCylinders.macro | 430 ++++++++++++++++++ run/SimExamples/Geometry_StepFiltering/run.sh | 86 ++++ 12 files changed, 904 insertions(+), 2 deletions(-) create mode 100644 Detectors/gconfig/src/KeepStepCylinders.macro create mode 100644 run/SimExamples/Geometry_StepFiltering/README.md create mode 100644 run/SimExamples/Geometry_StepFiltering/checkKeepStepCylinders.macro create mode 100644 run/SimExamples/Geometry_StepFiltering/countHits.macro create mode 100644 run/SimExamples/Geometry_StepFiltering/makeKeepStepCylinders.macro create mode 100755 run/SimExamples/Geometry_StepFiltering/run.sh diff --git a/Common/SimConfig/include/SimConfig/SimParams.h b/Common/SimConfig/include/SimConfig/SimParams.h index b5f975d1b0c6e..40d9883c576a2 100644 --- a/Common/SimConfig/include/SimConfig/SimParams.h +++ b/Common/SimConfig/include/SimConfig/SimParams.h @@ -23,6 +23,7 @@ namespace conf // (mostly used in O2MCApplication stepping) struct SimCutParams : public o2::conf::ConfigurableParamHelper { 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 diff --git a/Detectors/gconfig/CMakeLists.txt b/Detectors/gconfig/CMakeLists.txt index 444f125b7cbb7..be50c9bcfe72f 100644 --- a/Detectors/gconfig/CMakeLists.txt +++ b/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/) diff --git a/Detectors/gconfig/include/SimSetup/MCReplayParam.h b/Detectors/gconfig/include/SimSetup/MCReplayParam.h index dd39c87ab312c..31d51c62e62f2 100644 --- a/Detectors/gconfig/include/SimSetup/MCReplayParam.h +++ b/Detectors/gconfig/include/SimSetup/MCReplayParam.h @@ -29,6 +29,7 @@ struct MCReplayParam : public o2::conf::ConfigurableParamHelper { 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 diff --git a/Detectors/gconfig/src/KeepStepCylinders.macro b/Detectors/gconfig/src/KeepStepCylinders.macro new file mode 100644 index 0000000000000..67f319bf5eb89 --- /dev/null +++ b/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; + }; +} diff --git a/Detectors/gconfig/src/MCReplayConfig.cxx b/Detectors/gconfig/src/MCReplayConfig.cxx index 764c43088fb90..3554ddde3197d 100644 --- a/Detectors/gconfig/src/MCReplayConfig.cxx +++ b/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() diff --git a/Steer/include/Steer/O2MCApplicationBase.h b/Steer/include/Steer/O2MCApplicationBase.h index bd730c0f2fcb2..cdbf99b3b4b01 100644 --- a/Steer/include/Steer/O2MCApplicationBase.h +++ b/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 TrackRefFcn; + typedef std::function 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); }; diff --git a/Steer/src/O2MCApplication.cxx b/Steer/src/O2MCApplication.cxx index 7f37fd4bccf03..dba61328c2d9c 100644 --- a/Steer/src/O2MCApplication.cxx +++ b/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(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. diff --git a/run/SimExamples/Geometry_StepFiltering/README.md b/run/SimExamples/Geometry_StepFiltering/README.md new file mode 100644 index 0000000000000..3f93977042198 --- /dev/null +++ b/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/). diff --git a/run/SimExamples/Geometry_StepFiltering/checkKeepStepCylinders.macro b/run/SimExamples/Geometry_StepFiltering/checkKeepStepCylinders.macro new file mode 100644 index 0000000000000..264dd84671d38 --- /dev/null +++ b/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 +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace +{ +/// pull the rSq[] and edges[] initialisers out of a generated macro +bool readHull(const char* macro, std::vector& rad, std::vector& edges) +{ + std::ifstream in(macro); + if (!in) { + printf("[check] cannot open %s\n", macro); + return false; + } + std::string all((std::istreambuf_iterator(in)), std::istreambuf_iterator()); + auto grab = [&all](const char* name, std::vector& 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 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 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.); +} diff --git a/run/SimExamples/Geometry_StepFiltering/countHits.macro b/run/SimExamples/Geometry_StepFiltering/countHits.macro new file mode 100644 index 0000000000000..544f26d83e23f --- /dev/null +++ b/run/SimExamples/Geometry_StepFiltering/countHits.macro @@ -0,0 +1,75 @@ +// Count hits per detector for an o2-sim output prefix. +// Usage: root -b -q 'countHits.macro("dir/prefix")' +// Prints: and a TOTAL line. +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +void countHits(const char* prefix) +{ + TString p(prefix); + TString dir = gSystem->DirName(p); + TString base = gSystem->BaseName(p); + + TSystemDirectory sdir(dir, dir); + TList* files = sdir.GetListOfFiles(); + if (!files) { + std::cerr << "no dir " << dir << "\n"; + return; + } + std::map counts; + long long total = 0; + TIter next(files); + TSystemFile* f; + while ((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); + if (!file || file->IsZombie()) { + continue; + } + TTree* t = (TTree*)file->Get("o2sim"); + if (!t) { + file->Close(); + continue; + } + long long n = 0; + TObjArray* brs = t->GetListOfBranches(); + for (int i = 0; i < brs->GetEntries(); ++i) { + TBranch* b = (TBranch*)brs->At(i); + TString bn = b->GetName(); + if (!bn.Contains("Hit")) { + continue; + } + // number of elements of the vector branch, summed over entries + Long64_t nsel = t->Draw(Form("@%s.size()", bn.Data()), "", "goff"); + if (nsel <= 0) { + continue; + } + double* v = t->GetV1(); + for (Long64_t e = 0; e < nsel; ++e) { + n += (long long)v[e]; + } + } + counts[det.Data()] += n; + total += n; + file->Close(); + } + for (auto& kv : counts) { + std::cout << "HITS " << kv.first << " " << kv.second << "\n"; + } + std::cout << "HITS TOTAL " << total << "\n"; +} diff --git a/run/SimExamples/Geometry_StepFiltering/makeKeepStepCylinders.macro b/run/SimExamples/Geometry_StepFiltering/makeKeepStepCylinders.macro new file mode 100644 index 0000000000000..bfb4c871d550c --- /dev/null +++ b/run/SimExamples/Geometry_StepFiltering/makeKeepStepCylinders.macro @@ -0,0 +1,430 @@ +// 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 makeKeepStepCylinders.macro +/// \brief Derive the geometry hull used by SimCutParams.stepFilteringMacro +/// +/// Ray-traces the closed TGeo geometry to obtain the outermost radius at which +/// ALICE material sits, R(z), then approximates R(z) by a small number of +/// z-slabs (cylinders) and writes a ready-to-use keepStep() macro. +/// +/// The hull produced here is an *over*-approximation of R(z) by construction: +/// each slab takes the maximum of the sampled envelope over its z range, and a +/// safety margin is added on top. A step outside the hull can therefore not be +/// a step inside ALICE material. +/// +/// Method. For every z slice and every azimuthal direction phi, a ray is +/// started far out in the cavern and walked *inwards* boundary by boundary with +/// TGeoManager::FindNextBoundaryAndStep() until it enters the first qualifying +/// volume. The radius of that entry point is the envelope for that (z, phi). +/// Walking inwards means each ray stops at the outermost material instead of +/// crossing the whole detector, which is what makes a 1 cm z scan affordable. +/// Because the radius comes from a geometry boundary and not from a binned +/// scan, R(z) is exact in r; the only discretisation is in z and phi. +/// +/// A volume qualifies if its module is not excluded (default: HALL and CAVE, +/// i.e. the experimental cavern and the world/mother envelopes) and either +/// - its material density is above rhoMin, i.e. it is neither vacuum nor air, or +/// - it is a sensitive volume. Sensitive volumes must be included regardless +/// of density: several of them are gas (MID's RPCs, the TPC drift volume) +/// and would otherwise fall outside the hull that is supposed to contain +/// every hit. +/// The volume-to-module map and the sensitive-volume list are the +/// MCStepLoggerVolMap.dat and MCStepLoggerSenVol.dat that o2-sim writes next to +/// the geometry. Without the map, the module of a volume is taken to be the +/// top-level node the ray is currently under. +/// +/// Usage: +/// root -l -b -q 'makeKeepStepCylinders.macro("o2sim_geometry.root")' +/// +/// \author derived from the procedure presented by B. Voelkel, WP12/13, 11.10.2023 + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace +{ + +/// volume name -> module name, as written by O2MCApplicationBase::ConstructGeometry +std::map readVolumeToModule(const char* fname) +{ + std::map m; + std::ifstream in(fname); + if (!in) { + printf("[hull] no volume/module map at %s, falling back to top-node names\n", fname); + return m; + } + std::string line; + while (std::getline(in, line)) { + auto pos = line.rfind(':'); + if (pos == std::string::npos) { + continue; + } + m[line.substr(0, pos)] = line.substr(pos + 1); + } + printf("[hull] volume/module map: %zu volumes from %s\n", m.size(), fname); + return m; +} + +/// sensitive volume names, as written by O2MCApplicationBase::ConstructGeometry +/// (format ":") +std::set readSensitiveVolumes(const char* fname) +{ + std::set s; + std::ifstream in(fname); + if (!in) { + printf("[hull] no sensitive-volume list at %s; gas-only sensitive volumes may be clipped\n", fname); + return s; + } + std::string line; + while (std::getline(in, line)) { + auto pos = line.find(':'); + if (pos == std::string::npos) { + continue; + } + s.insert(line.substr(pos + 1)); + } + printf("[hull] sensitive volumes: %zu from %s\n", s.size(), fname); + return s; +} + +std::set splitCSV(const char* s) +{ + std::set out; + std::stringstream ss(s); + std::string item; + while (std::getline(ss, item, ',')) { + if (!item.empty()) { + out.insert(item); + } + } + return out; +} + +/// Best piecewise-constant over-approximation of the envelope with exactly +/// nSlab slabs. Minimises sum_i Rslab(i)^2 * len(i), i.e. the transported +/// volume that the hull encloses, by exact dynamic programming over the fine +/// z bins. Rslab is the maximum of R over the slab, so the result always +/// contains the fine envelope. +void segment(const std::vector& r, int nSlab, std::vector& firstBin, std::vector& radius) +{ + const int n = r.size(); + const double kInf = std::numeric_limits::max() / 4; + // dp[k][j]: best cost of covering bins [0, j) with k slabs + std::vector> dp(nSlab + 1, std::vector(n + 1, kInf)); + std::vector> back(nSlab + 1, std::vector(n + 1, -1)); + dp[0][0] = 0.; + for (int k = 1; k <= nSlab; ++k) { + for (int j = 1; j <= n; ++j) { + double runMax = 0.; + for (int i = j; i >= 1; --i) { // slab covers bins [i-1, j) + runMax = std::max(runMax, r[i - 1]); + if (dp[k - 1][i - 1] >= kInf) { + continue; + } + const double cost = dp[k - 1][i - 1] + runMax * runMax * (j - i + 1); + if (cost < dp[k][j]) { + dp[k][j] = cost; + back[k][j] = i - 1; + } + } + } + } + firstBin.assign(nSlab, 0); + radius.assign(nSlab, 0.); + int j = n; + for (int k = nSlab; k >= 1; --k) { + const int i = back[k][j]; + firstBin[k - 1] = i; + radius[k - 1] = *std::max_element(r.begin() + i, r.begin() + j); + j = i; + } +} + +} // namespace + +void makeKeepStepCylinders(const char* geomfile = "o2sim_geometry.root", + const char* volmapfile = "MCStepLoggerVolMap.dat", + const char* senvolfile = "MCStepLoggerSenVol.dat", + const char* outmacro = "KeepStepCylinders.macro", + int nSlab = 12, + double zmin = 0., double zmax = 0., int nz = 0, + int nphi = 128, double rscan = 2000., + double marginAbs = 2., double marginRel = 0.02, + double rhoMin = 0.01, + double rFloor = 50., double rFloorBeyondZ = 1900., + const char* excludeModules = "HALL,CAVE", + const char* outroot = "KeepStepCylinders.root", + const char* method = "radial", + int maxSegmentBins = 3000) +{ + if (!gGeoManager) { + TGeoManager::Import(geomfile); + } + if (!gGeoManager) { + printf("[hull] could not load %s\n", geomfile); + return; + } + const auto vol2mod = readVolumeToModule(volmapfile); + const auto sensitive = readSensitiveVolumes(senvolfile); + const auto excluded = splitCSV(excludeModules); + + // cache the qualification decision per TGeoVolume, the ray walk revisits them constantly + std::map cache; + auto qualifies = [&](const TGeoNode* node) { + if (!node) { + return false; + } + const TGeoVolume* v = node->GetVolume(); + auto it = cache.find(v); + if (it != cache.end()) { + return it->second; + } + const bool isSensitive = sensitive.count(v->GetName()) > 0; + bool ok = true; + const TGeoMaterial* mat = v->GetMaterial(); + if (!isSensitive && (!mat || mat->GetDensity() < rhoMin)) { + ok = false; + } + if (ok) { + std::string mod; + auto mit = vol2mod.find(v->GetName()); + if (mit != vol2mod.end()) { + mod = mit->second; + } else if (gGeoManager->GetLevel() > 1) { + // no map available: approximate the module by the top-level node + mod = gGeoManager->GetMother(gGeoManager->GetLevel() - 1)->GetVolume()->GetName(); + } + if (excluded.count(mod)) { + ok = false; + } + } + cache[v] = ok; + return ok; + }; + + // start the rays just inside the world volume + const auto* worldBox = dynamic_cast(gGeoManager->GetTopVolume()->GetShape()); + const double rstart = worldBox ? std::min(rscan, 0.98 * std::min(worldBox->GetDX(), worldBox->GetDY())) : rscan; + + // Cover the whole world in z unless told otherwise. Outside the traced range + // the hull keeps every step, so a shorter range simply leaves that region to + // the built-in z/R cut of O2MCApplicationBase. + if (zmin == 0. && zmax == 0.) { + const double dzWorld = worldBox ? worldBox->GetDZ() : 3000.; + zmin = -dzWorld; + zmax = dzWorld; + } + if (nz <= 0) { + nz = static_cast(zmax - zmin); // 1 cm slices + } + printf("[hull] z range [%.0f, %.0f] cm in %d slices\n", zmin, zmax, nz); + + const double dz = (zmax - zmin) / nz; + std::vector renv(nz, 0.), rprojXZ(nz, 0.), rprojYZ(nz, 0.); + printf("[hull] ray-tracing %d z slices x %d phi directions, inwards from r = %.0f cm ...\n", nz, nphi, rstart); + + for (int iz = 0; iz < nz; ++iz) { + const double z = zmin + (iz + 0.5) * dz; + for (int ip = 0; ip < nphi; ++ip) { + const double phi = TMath::TwoPi() * ip / nphi; + const double dx = std::cos(phi), dy = std::sin(phi); + gGeoManager->InitTrack(rstart * dx, rstart * dy, z, -dx, -dy, 0.); + double rhit = 0.; + if (qualifies(gGeoManager->GetCurrentNode())) { + rhit = rstart; // should not happen unless the cavern itself qualifies + } + while (rhit == 0. && !gGeoManager->IsOutside()) { + gGeoManager->FindNextBoundaryAndStep(); + if (gGeoManager->GetStep() <= 0.) { + break; + } + const double* p = gGeoManager->GetCurrentPoint(); + const double r = std::sqrt(p[0] * p[0] + p[1] * p[1]); + if (qualifies(gGeoManager->GetCurrentNode())) { + rhit = r; + } + if (r < 1e-3) { + break; // reached the beam line without finding material + } + } + if (rhit > renv[iz]) { + renv[iz] = rhit; + } + const double ax = std::fabs(rhit * dx), ay = std::fabs(rhit * dy); + if (ax > rprojXZ[iz]) { + rprojXZ[iz] = ax; + } + if (ay > rprojYZ[iz]) { + rprojYZ[iz] = ay; + } + } + if (nz >= 20 && iz % (nz / 20) == 0) { + printf("[hull] z = %8.1f cm R = %7.1f cm\n", z, renv[iz]); + } + } + + // "radial" is the maximal radius of material and is a true over-approximation. + // "projection" reproduces the 2023 procedure, max|x| and max|y| taken from the + // XZ and YZ projections; it is smaller than the radial envelope wherever the + // cross section is not a disc -- notably the octagonal L3 magnet -- and is + // therefore not conservative. Kept for comparison. + const bool useProjection = TString(method).EqualTo("projection", TString::kIgnoreCase); + if (useProjection) { + for (int iz = 0; iz < nz; ++iz) { + renv[iz] = std::max(rprojXZ[iz], rprojYZ[iz]); + } + } + printf("[hull] envelope method: %s\n", useProjection ? "projection (XZ/YZ, 2023)" : "radial (max r)"); + + // safety margin: covers the finite phi sampling and the finite z slicing + std::vector rsafe(nz); + for (int iz = 0; iz < nz; ++iz) { + rsafe[iz] = renv[iz] > 0. ? renv[iz] * (1. + marginRel) + marginAbs : 0.; + } + + // Beyond rFloorBeyondZ the hull must not be tighter than the built-in ZDC + // tunnel cut (SimCutParams.maxRTrackingZDC, 50 cm). There is no material to + // trace along the beam line there, but particles travel to the ZDC through + // the vacuum inside the pipe and are bent onto it by the beam optics, so an + // envelope taken from material alone does not contain their trajectories. + if (rFloor > 0.) { + for (int iz = 0; iz < nz; ++iz) { + const double z = zmin + (iz + 0.5) * dz; + if (std::fabs(z) > rFloorBeyondZ) { + rsafe[iz] = std::max(rsafe[iz], rFloor); + } + } + printf("[hull] floor of %.0f cm applied for |z| > %.0f cm\n", rFloor, rFloorBeyondZ); + } + + // The segmentation is O(nSlab * n^2), so run it on a coarser grid when the + // scan is fine. Coarsening takes the maximum inside each group, which keeps + // the result an over-approximation; only the z resolution of the slab edges + // is reduced. + const int group = std::max(1, (nz + maxSegmentBins - 1) / maxSegmentBins); + std::vector rseg((nz + group - 1) / group, 0.); + for (int iz = 0; iz < nz; ++iz) { + rseg[iz / group] = std::max(rseg[iz / group], rsafe[iz]); + } + if (group > 1) { + printf("[hull] segmenting on %zu bins of %.1f cm\n", rseg.size(), group * dz); + } + + std::vector firstBin; + std::vector radius; + segment(rseg, nSlab, firstBin, radius); + + // round up to whole cm; a larger hull can only be safer + std::vector edges(nSlab + 1); + for (int k = 0; k < nSlab; ++k) { + firstBin[k] *= group; + edges[k] = zmin + firstBin[k] * dz; + radius[k] = std::ceil(radius[k]); + } + edges[nSlab] = zmax; + + printf("\n[hull] %d-cylinder approximation of R(z):\n", nSlab); + double kept = 0., env = 0.; + for (int k = 0; k < nSlab; ++k) { + const int i0 = firstBin[k]; + const int i1 = (k + 1 < nSlab) ? firstBin[k + 1] : nz; + for (int i = i0; i < i1; ++i) { + kept += radius[k] * radius[k] * dz; + env += rsafe[i] * rsafe[i] * dz; + } + printf(" z in [%8.1f, %8.1f) R < %6.0f cm\n", edges[k], edges[k + 1], radius[k]); + } + printf("[hull] hull volume / exact-envelope volume = %.3f\n", env > 0 ? kept / env : 0.); + + // ---------------------------------------------------------------- the macro + std::ofstream out(outmacro); + out << "// Generated by run/SimExamples/Geometry_StepFiltering/makeKeepStepCylinders.macro\n" + << "// -- do not edit by hand.\n" + << "// Geometry: " << geomfile << "\n" + << "// Envelope: max radius of sensitive volumes and of material with rho > " << rhoMin + << " g/cm3,\n// excluding modules " << excludeModules << ",\n" + << "// sampled at " << nz << " z slices x " << nphi + << " phi directions, margin " << marginRel * 100 << "% + " << marginAbs << " cm.\n" + << "// Steps outside the z range below are always kept.\n\n" + << "o2::steer::O2MCApplicationBase::KeepStepFcn keepStep()\n" + << "{\n" + << " const float rSq[] = {"; + for (int k = 0; k < nSlab; ++k) { + out << (k ? ", " : "") << Form("%.1ff * %.1ff", radius[k], radius[k]); + } + out << "};\n" + << " const float edges[] = {"; + for (int k = 0; k <= nSlab; ++k) { + out << (k ? ", " : "") << Form("%.1ff", edges[k]); + } + out << "};\n" + << " return [rSq, edges](TVirtualMC const* mc) {\n" + << " float x, y, z;\n" + << " mc->TrackPosition(x, y, z);\n" + << " if (z < edges[0] || z >= edges[" << nSlab << "]) {\n" + << " return true; // beyond the traced region, e.g. the ZDC tunnel\n" + << " }\n" + << " for (auto i = 0U; i < " << nSlab << "; ++i) {\n" + << " if (edges[i + 1] > z && z >= edges[i]) {\n" + << " return (x * x + y * y) < rSq[i];\n" + << " }\n" + << " }\n" + << " return true;\n" + << " };\n" + << "}\n"; + out.close(); + printf("[hull] wrote %s\n", outmacro); + + // ------------------------------------------------------- envelope for plots + if (outroot && outroot[0]) { + TFile f(outroot, "RECREATE"); + TGraph gEnv(nz), gSafe(nz), gProjX(nz), gProjY(nz); + for (int iz = 0; iz < nz; ++iz) { + const double z = zmin + (iz + 0.5) * dz; + gEnv.SetPoint(iz, z, renv[iz]); + gSafe.SetPoint(iz, z, rsafe[iz]); + gProjX.SetPoint(iz, z, rprojXZ[iz]); + gProjY.SetPoint(iz, z, rprojYZ[iz]); + } + TGraph gHull(2 * nSlab); + for (int k = 0; k < nSlab; ++k) { + gHull.SetPoint(2 * k, edges[k], radius[k]); + gHull.SetPoint(2 * k + 1, edges[k + 1], radius[k]); + } + gEnv.Write("envelope"); + gSafe.Write("envelopeWithMargin"); + gProjX.Write("projectionXZ"); + gProjY.Write("projectionYZ"); + gHull.Write("hull"); + f.Close(); + printf("[hull] wrote %s\n", outroot); + } +} diff --git a/run/SimExamples/Geometry_StepFiltering/run.sh b/run/SimExamples/Geometry_StepFiltering/run.sh new file mode 100755 index 0000000000000..f19c040665cb2 --- /dev/null +++ b/run/SimExamples/Geometry_StepFiltering/run.sh @@ -0,0 +1,86 @@ +#!/usr/bin/env bash +# +# This is a simulation example showing how to reduce the number of transport +# steps with a geometry hull: a set of cylinders wrapped around the detector, +# outside of which transport of a track is stopped. +# +# The hull is plugged into the per-step decision of O2MCApplicationBase via +# ``` +# [SimCutParams] +# stepFilteringMacro=$O2_ROOT/share/Detectors/gconfig/KeepStepCylinders.macro +# ``` +# The macro provides a `keepStep()` function returning the predicate that is +# evaluated for every step. It is applied in addition to the built-in z/R cut, +# so it can only remove steps, never add them. +# +# The chain below +# 1. records a reference simulation with the MCStepLogger, +# 2. derives a hull from the geometry of that simulation, +# 3. checks that every recorded hit lies inside the hull, +# 4. replays the reference without cuts and with the hull, and compares +# the number of steps and the number of hits per detector. +# +# Requires the step logger in the environment: +# alienv enter O2/latest MCStepLogger/latest + +set -x + +MODULES="" # empty: default detector list +EVENTS=10 +HERE=$(cd "$(dirname "$0")" && pwd) +TOP=$PWD + +# Alignment is irrelevant to step and hit counting and switching it off keeps +# the example from needing a CCDB connection and an alien token. +CONFIG="align-geom.mDetectors=none" + +# ---------------------------------------------------------------- 1. reference +mkdir -p ref && cd ref +MCSTEPLOG_TTREE=1 \ +LD_PRELOAD=$MCSTEPLOGGER_ROOT/lib/libMCStepLoggerInterceptSteps.so \ +o2-sim-serial -n ${EVENTS} -g pythia8pp -e TGeant4 ${MODULES} \ + --configKeyValues "${CONFIG}" > logref 2>&1 +cd .. + +# ------------------------------------------------------------------- 2. derive +# The defaults scan the full world in 1 cm slices and 128 azimuthal directions +# and approximate the result with 12 cylinders; this takes about half a minute. +( cd ref && root -l -b -q "${HERE}/makeKeepStepCylinders.macro(\ +\"o2sim_geometry.root\", \"MCStepLoggerVolMap.dat\", \"MCStepLoggerSenVol.dat\", \ +\"${TOP}/KeepStepCylinders.macro\")" ) + +# -------------------------------------------------------------------- 3. check +# Necessary condition: no hit may sit outside the hull. +( cd ref && root -l -b -q "${HERE}/checkKeepStepCylinders.macro(\ +\"${TOP}/KeepStepCylinders.macro\", \"o2sim\")" ) + +# ------------------------------------------------------------------- 4. replay +# MCReplayParam.allowStopTrack lets the replay act on the StopTrack() calls the +# hull makes; without it a replay reproduces its reference and nothing is cut. +mkdir -p baseline && cd baseline +o2-sim-serial -n ${EVENTS} -e MCReplay -g extkinO2 ${MODULES} \ + --extKinFile ../ref/o2sim_Kine.root -o baseline \ + --configKeyValues "${CONFIG};MCReplayParam.stepFilename=../ref/MCStepLoggerOutput.root" \ + > logbaseline 2>&1 +cd .. + +mkdir -p hull && cd hull +o2-sim-serial -n ${EVENTS} -e MCReplay -g extkinO2 ${MODULES} \ + --extKinFile ../ref/o2sim_Kine.root -o hull \ + --configKeyValues "${CONFIG};MCReplayParam.stepFilename=../ref/MCStepLoggerOutput.root;\ +MCReplayParam.allowStopTrack=true;\ +SimCutParams.stepFilteringMacro=${TOP}/KeepStepCylinders.macro" \ + > loghull 2>&1 +cd .. + +# ------------------------------------------------------------------ 5. compare +# MCReplay prints "Original number, skipped, kept, ..." once per event. +grep -h "Original number" baseline/logbaseline hull/loghull + +# Hit counts per detector. FT0 draws a random number while creating hits and +# does not reproduce under replay, so compare it against the baseline replay +# rather than against the reference simulation, and expect it to be noisy. +for d in baseline hull; do + echo "=== $d" + root -l -b -q "${HERE}/countHits.macro(\"${d}/${d}\")" +done From 50437c2e43de724aa76cf3105318aa856e5a7c8a Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Mon, 17 Aug 2026 11:53:31 +0000 Subject: [PATCH 2/2] Please consider the following formatting changes --- Common/SimConfig/include/SimConfig/SimParams.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Common/SimConfig/include/SimConfig/SimParams.h b/Common/SimConfig/include/SimConfig/SimParams.h index 40d9883c576a2..c56e46c91c2ae 100644 --- a/Common/SimConfig/include/SimConfig/SimParams.h +++ b/Common/SimConfig/include/SimConfig/SimParams.h @@ -23,7 +23,7 @@ namespace conf // (mostly used in O2MCApplication stepping) struct SimCutParams : public o2::conf::ConfigurableParamHelper { bool stepFiltering = true; // if we activate the step filtering in O2BaseMCApplication - std::string stepFilteringMacro = ""; // ROOT macro providing keepStep(); empty = built-in z/R cut + 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