From e2da7786a665da9b87d8b4665f20a5a7d3f5901b Mon Sep 17 00:00:00 2001 From: swenzel Date: Wed, 22 Jul 2026 14:34:03 +0200 Subject: [PATCH 1/2] CAD->TGeo: support multiple external modules and sensitive detectors Improve the CAD->TGeo integration path for o2-sim with a configurable JSON-driven setup for externally provided CAD geometry. The configuration can describe multiple passive external modules and sensitive external detectors, both activated through the normal detector/module list. Add a --extGeomFile option carrying the external geometry configuration. Passive CAD modules no longer need to be hardcoded in build_geometry.C, and several configured CAD modules can be injected in one simulation. Introduce shared CAD geometry utilities to JIT-load O2_CADtoTGeo.py-style ROOT macros into unique namespaces, allowing several CAD modules to coexist in one Cling session despite identical exported builder symbols. The utilities also remap imported CAD media into the O2 MaterialManager. Add o2::ext::ExternalDetector, a generic sensitive detector implementation for externally supplied CAD geometry. Configured volumes or media can be marked sensitive, the detector is tied to a free DetID slot, and hits flow through the standard o2-sim forwarding and merging machinery. Add o2::ext::Hit as a detector-agnostic external hit payload storing entrance/exit position, momentum, energy, energy loss, time, length, PDG code and track status flags. All external detector instances share this wire format, so the hit merger only needs one external hit type. Support optional per-detector sensitive actions loaded from ROOT macros at runtime via GetFromMacro. When no custom action is configured, a built-in charged-track entrance/exit action is used. Register active external detectors in O2HitMerger so their hits are persisted in parallel mode, and key DetImpl hit collection buffers by detector instance to support multiple detectors sharing the same C++ type. Add a self-contained SimExample with two artificial sensitive external detectors, demonstrating configurable geometry injection, built-in and custom sensitive actions, and parallel-mode hit persistence. --- .clang-format | 4 + .../SimConfig/include/SimConfig/SimConfig.h | 4 +- Common/SimConfig/src/SimConfig.cxx | 4 +- Detectors/Base/CMakeLists.txt | 1 + .../include/DetectorsBase/CADGeometryUtils.h | 50 ++ .../Base/include/DetectorsBase/Detector.h | 11 +- Detectors/Base/src/CADGeometryUtils.cxx | 202 ++++++++ Detectors/CMakeLists.txt | 1 + Detectors/External/CMakeLists.txt | 22 + .../ExternalDetectors/ExternalDetector.h | 176 +++++++ .../External/include/ExternalDetectors/Hit.h | 140 ++++++ .../macro/sensitiveActionExample.macro | 68 +++ Detectors/External/src/ExternalDetector.cxx | 451 ++++++++++++++++++ .../External/src/ExternalDetectorsLinkDef.h | 23 + Detectors/Passive/CMakeLists.txt | 3 +- .../include/DetectorsPassive/ExternalModule.h | 15 +- Detectors/Passive/src/ExternalModule.cxx | 228 ++++----- macro/build_geometry.C | 36 +- run/CMakeLists.txt | 1 + run/O2HitMerger.h | 60 +++ .../External_Sensitive_Detectors/README.md | 54 +++ .../detectorlist.json | 6 + .../externalDetectors.json | 23 + .../geometry_innerCylinder.macro | 35 ++ .../geometry_outerDisk.macro | 35 ++ .../inspect_hits.macro | 50 ++ .../External_Sensitive_Detectors/run.sh | 32 ++ .../sensitive_action.macro | 44 ++ scripts/geometry/README.md | 236 ++++++++- scripts/geometry/TODO.md | 4 + 30 files changed, 1853 insertions(+), 166 deletions(-) create mode 100644 Detectors/Base/include/DetectorsBase/CADGeometryUtils.h create mode 100644 Detectors/Base/src/CADGeometryUtils.cxx create mode 100644 Detectors/External/CMakeLists.txt create mode 100644 Detectors/External/include/ExternalDetectors/ExternalDetector.h create mode 100644 Detectors/External/include/ExternalDetectors/Hit.h create mode 100644 Detectors/External/macro/sensitiveActionExample.macro create mode 100644 Detectors/External/src/ExternalDetector.cxx create mode 100644 Detectors/External/src/ExternalDetectorsLinkDef.h create mode 100644 run/SimExamples/External_Sensitive_Detectors/README.md create mode 100644 run/SimExamples/External_Sensitive_Detectors/detectorlist.json create mode 100644 run/SimExamples/External_Sensitive_Detectors/externalDetectors.json create mode 100644 run/SimExamples/External_Sensitive_Detectors/geometry_innerCylinder.macro create mode 100644 run/SimExamples/External_Sensitive_Detectors/geometry_outerDisk.macro create mode 100644 run/SimExamples/External_Sensitive_Detectors/inspect_hits.macro create mode 100755 run/SimExamples/External_Sensitive_Detectors/run.sh create mode 100644 run/SimExamples/External_Sensitive_Detectors/sensitive_action.macro create mode 100644 scripts/geometry/TODO.md diff --git a/.clang-format b/.clang-format index 93ba0f7b0c187..f0eb7d31df8e3 100644 --- a/.clang-format +++ b/.clang-format @@ -54,3 +54,7 @@ UseTab: Never # Do not format protobuf files Language: Proto DisableFormat: true +--- +# Do not format JSON configuration files +Language: Json +DisableFormat: true diff --git a/Common/SimConfig/include/SimConfig/SimConfig.h b/Common/SimConfig/include/SimConfig/SimConfig.h index be88d9fbd8c33..0be82ba1921b5 100644 --- a/Common/SimConfig/include/SimConfig/SimConfig.h +++ b/Common/SimConfig/include/SimConfig/SimConfig.h @@ -88,8 +88,9 @@ struct SimConfigData { bool mForwardKine = false; // true if tracks and event headers are to be published on a FairMQ channel (for reading by other consumers) bool mWriteToDisc = true; // whether we write simulation products (kine, hits) to disc VertexMode mVertexMode = VertexMode::kDiamondParam; // by default we should use die InteractionDiamond parameter + std::string mExtGeomFile = ""; // optional path to a JSON file describing external (CAD) geometry modules to inject - ClassDefNV(SimConfigData, 4); + ClassDefNV(SimConfigData, 5); }; // A singleton class which can be used @@ -178,6 +179,7 @@ class SimConfig bool forwardKine() const { return mConfigData.mForwardKine; } bool writeToDisc() const { return mConfigData.mWriteToDisc; } VertexMode getVertexMode() const { return mConfigData.mVertexMode; } + std::string getExtGeomFilename() const { return mConfigData.mExtGeomFile; } // returns the pair of collision context filename as well as event prefix encoded // in the mFromCollisionContext string. Returns empty string if information is not available or set. diff --git a/Common/SimConfig/src/SimConfig.cxx b/Common/SimConfig/src/SimConfig.cxx index 15879687872d5..14f99dd4ea580 100644 --- a/Common/SimConfig/src/SimConfig.cxx +++ b/Common/SimConfig/src/SimConfig.cxx @@ -75,7 +75,8 @@ void SimConfig::initOptions(boost::program_options::options_description& options "asservice", bpo::value()->default_value(false), "run in service/server mode")( "noGeant", bpo::bool_switch(), "prohibits any Geant transport/physics (by using tight cuts)")( "forwardKine", bpo::bool_switch(), "forward kinematics on a FairMQ channel")( - "noDiscOutput", bpo::bool_switch(), "switch off writing sim results to disc (useful in combination with forwardKine)"); + "noDiscOutput", bpo::bool_switch(), "switch off writing sim results to disc (useful in combination with forwardKine)")( + "extGeomFile", bpo::value()->default_value(""), "Path to a JSON file describing external (CAD) geometry modules to inject (see Detectors/Passive ExternalModule). Modules are added when their 'name' is part of the active module list."); options.add_options()("fromCollContext", bpo::value()->default_value(""), "Use a pregenerated collision context to infer number of events to simulate, how to embedd them, the vertex position etc. Takes precedence of other options such as \"--nEvents\". The format is COLLISIONCONTEXTFILE.root[:SIGNALNAME] where SIGNALNAME is the event part in the context which is relevant."); } @@ -354,6 +355,7 @@ bool SimConfig::resetFromParsedMap(boost::program_options::variables_map const& if (vm.count("noemptyevents")) { mConfigData.mFilterNoHitEvents = true; } + mConfigData.mExtGeomFile = vm["extGeomFile"].as(); mConfigData.mFromCollisionContext = vm["fromCollContext"].as(); auto collcontext_simprefix = getCollContextFilenameAndEventPrefix(); adjustFromCollContext(collcontext_simprefix.first, collcontext_simprefix.second); diff --git a/Detectors/Base/CMakeLists.txt b/Detectors/Base/CMakeLists.txt index 83a9193274e4f..b3162ba9caaf0 100644 --- a/Detectors/Base/CMakeLists.txt +++ b/Detectors/Base/CMakeLists.txt @@ -31,6 +31,7 @@ o2_add_library(DetectorsBase src/GlobalParams.cxx src/O2Tessellated.cxx src/TGeoGeometryUtils.cxx + src/CADGeometryUtils.cxx PUBLIC_LINK_LIBRARIES FairRoot::Base O2::CommonUtils O2::DetectorsCommonDataFormats diff --git a/Detectors/Base/include/DetectorsBase/CADGeometryUtils.h b/Detectors/Base/include/DetectorsBase/CADGeometryUtils.h new file mode 100644 index 0000000000000..12d626d0eb572 --- /dev/null +++ b/Detectors/Base/include/DetectorsBase/CADGeometryUtils.h @@ -0,0 +1,50 @@ +// 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 CADGeometryUtils.h +/// \brief Helpers to inject CAD-derived (TGeo) geometry into O2 simulation +/// +/// These utilities are shared between purely passive external modules +/// (o2::passive::ExternalModule) and sensitive external detectors +/// (o2::ext::ExternalDetector). They deal with the geometry produced by +/// scripts/geometry/O2_CADtoTGeo.py, which is emitted as a ROOT macro. + +#ifndef ALICEO2_BASE_CADGEOMETRYUTILS_H +#define ALICEO2_BASE_CADGEOMETRYUTILS_H + +#include + +class TGeoVolume; + +namespace o2::base +{ + +/// JIT-compile a CAD-derived ROOT geometry macro (as produced by O2_CADtoTGeo.py) +/// and execute it to obtain the top TGeoVolume of the described module. +/// +/// The macro body is wrapped into a unique namespace (derived from \a instanceTag) +/// so that several such macros — which all export identically named symbols +/// (build(), get_builder_hook_unchecked(), ...) — can coexist in the same Cling +/// session without colliding. Returns nullptr on failure. +/// +/// \param macroFile path to the geometry macro (shell variables are expanded) +/// \param instanceTag a short tag used to build a unique, human-readable namespace +TGeoVolume* buildCADVolumeFromMacro(const std::string& macroFile, const std::string& instanceTag); + +/// Re-register the TGeo media used in the volume tree rooted at \a top into the O2 +/// MaterialManager under ownership of \a modulename, rewiring the volumes to the +/// newly created media. This brings the CAD-imported media under O2's media/cut +/// handling (so that e.g. tracking cuts apply consistently). +void remapCADMedia(TGeoVolume* top, const char* modulename); + +} // namespace o2::base + +#endif diff --git a/Detectors/Base/include/DetectorsBase/Detector.h b/Detectors/Base/include/DetectorsBase/Detector.h index f1744086d6a05..54e264dd22666 100644 --- a/Detectors/Base/include/DetectorsBase/Detector.h +++ b/Detectors/Base/include/DetectorsBase/Detector.h @@ -530,9 +530,14 @@ class DetImpl : public o2::base::Detector { using Hit_t = typename std::remove_pointer(this)->Det::getHits(0))>::type; using Collector_t = tbb::concurrent_unordered_map>>>; - static Collector_t hitcollector; // note: we can't put this as member because - // decltype type deduction doesn't seem to work for class members; so we use a static member - // and will use some pointer member to communicate this data to other functions + // note: we can't put this as a member because decltype type deduction doesn't seem to work for + // class members; so we use a static and communicate it to other functions via a pointer member. + // The collector must be kept *per detector instance* (keyed by 'this'): for most detectors there + // is a single instance per C++ type, but several external detectors share the same type + // (o2::ext::ExternalDetector) and would otherwise clobber/double-free each other's buffers. + // tbb::concurrent_unordered_map is node-based, so the reference stays valid across insertions. + static tbb::concurrent_unordered_map hitcollectors; + auto& hitcollector = hitcollectors[this]; mHitCollectorBufferPtr = (char*)&hitcollector; int probe = 0; diff --git a/Detectors/Base/src/CADGeometryUtils.cxx b/Detectors/Base/src/CADGeometryUtils.cxx new file mode 100644 index 0000000000000..84c1890d844b5 --- /dev/null +++ b/Detectors/Base/src/CADGeometryUtils.cxx @@ -0,0 +1,202 @@ +// 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. + +#include "DetectorsBase/CADGeometryUtils.h" +#include "DetectorsBase/MaterialManager.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace o2::base +{ + +TGeoVolume* buildCADVolumeFromMacro(const std::string& macroFile, const std::string& instanceTag) +{ + if (macroFile.empty()) { + return nullptr; + } + auto expandedHookFileName = o2::utils::expandShellVarsInFileName(macroFile); + if (!std::filesystem::exists(expandedHookFileName)) { + LOG(error) << "External geometry macro " << expandedHookFileName << " does not exist"; + return nullptr; + } + + // We JIT the macro into a *unique* namespace per call. This is essential when several + // external geometries are present at the same time: every macro produced by + // O2_CADtoTGeo.py exports identically named symbols (build(), get_builder_hook_unchecked(), + // LoadFacets(), ...). Loading them all into the single global Cling scope would collide + // (the first definition wins and subsequent macros are silently ignored). By wrapping each + // macro body in its own namespace we keep the symbols separate. The preprocessor #include + // lines must stay at global scope, so we hoist them out of the namespace. + std::ifstream macroStream(expandedHookFileName, std::ios::in); + if (!macroStream.is_open()) { + LOG(error) << "Cannot open external geometry macro " << expandedHookFileName; + return nullptr; + } + std::string preamble; // #include (and other top-level preprocessor) lines -> global scope + std::string body; // everything else -> wrapped into a unique namespace + std::string line; + while (std::getline(macroStream, line)) { + auto firstNonSpace = line.find_first_not_of(" \t"); + if (firstNonSpace != std::string::npos && line[firstNonSpace] == '#') { + preamble += line + "\n"; + } else { + body += line + "\n"; + } + } + + // build a unique, valid C++ identifier for the namespace + static std::atomic instanceCounter{0}; + std::string ns = std::string("o2_cadgeom_") + instanceTag + "_" + std::to_string(instanceCounter++); + for (auto& c : ns) { + if (!std::isalnum(static_cast(c)) && c != '_') { + c = '_'; + } + } + + const std::string wrapped = preamble + "\nnamespace " + ns + " {\n" + body + "\n}\n"; + if (!gInterpreter->Declare(wrapped.c_str())) { + LOG(error) << "Failed to JIT external geometry macro " << expandedHookFileName; + return nullptr; + } + + // retrieve the builder hook from the unique namespace + const std::string globalName = "__" + ns + "_hook__"; + gROOT->ProcessLine(Form("std::function %s = %s::get_builder_hook_unchecked();", + globalName.c_str(), ns.c_str())); + auto global = gROOT->GetGlobal(globalName.c_str()); + if (!global) { + LOG(error) << "Could not retrieve geometry builder hook from macro " << expandedHookFileName; + return nullptr; + } + auto hook = *reinterpret_cast*>(global->GetAddress()); + LOG(info) << "CAD geometry hook initialized from file " << expandedHookFileName << " (namespace " << ns << ")"; + + auto top = hook(); + if (!top) { + LOG(error) << "CAD geometry macro " << expandedHookFileName << " did not return a top volume"; + } + return top; +} + +void remapCADMedia(TGeoVolume* top, const char* modulename) +{ + std::unordered_map medium_ptr_mapping; + std::unordered_set volumes_already_treated; + int counter = 1; + + // The transformer function + auto transform_media = [&](TGeoVolume* vol_) { + if (volumes_already_treated.find(vol_) != volumes_already_treated.end()) { + // this volume was already transformed + return; + } + volumes_already_treated.insert(vol_); + + if (dynamic_cast(vol_)) { + // do nothing for assemblies (they don't have a medium) + return; + } + + auto medium = vol_->GetMedium(); + if (!medium) { + return; + } + + auto iter = medium_ptr_mapping.find(medium); + if (iter != medium_ptr_mapping.end()) { + // This medium has already been transformed, so + // we just update the volume + vol_->SetMedium(iter->second); + return; + } else { + LOG(info) << "Transforming media with name " << medium->GetName() << " for volume " << vol_->GetName(); + + // we found a medium, not yet treated + auto curr_mat = medium->GetMaterial(); + auto& matmgr = o2::base::MaterialManager::Instance(); + + matmgr.Material(modulename, counter, curr_mat->GetName(), curr_mat->GetA(), curr_mat->GetZ(), curr_mat->GetDensity(), curr_mat->GetRadLen(), curr_mat->GetIntLen()); + // TGeo medium params are stored in a flat array with the following convention + // fParams[0] = isvol; + // fParams[1] = ifield; + // fParams[2] = fieldm; + // fParams[3] = tmaxfd; + // fParams[4] = stemax; + // fParams[5] = deemax; + // fParams[6] = epsil; + // fParams[7] = stmin; + const auto isvol = medium->GetParam(0); + const auto isxfld = medium->GetParam(1); + const auto sxmgmx = medium->GetParam(2); + const auto tmaxfd = medium->GetParam(3); + const auto stemax = medium->GetParam(4); + const auto deemax = medium->GetParam(5); + const auto epsil = medium->GetParam(6); + const auto stmin = medium->GetParam(7); + + matmgr.Medium(modulename, counter, medium->GetName(), counter, isvol, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); + + // there will be new Material and Medium objects; fetch them + auto new_med = matmgr.getTGeoMedium(modulename, counter); + + // insert into cache + medium_ptr_mapping[medium] = new_med; + vol_->SetMedium(new_med); + counter++; + } + }; // end transformer lambda + + // a generic volume walker + std::function visit_volume; + visit_volume = [&](TGeoVolume* vol) -> void { + if (!vol) { + return; + } + + // call the transformer + transform_media(vol); + + // Recurse into daughters + const int nd = vol->GetNdaughters(); + for (int i = 0; i < nd; ++i) { + TGeoNode* node = vol->GetNode(i); + if (!node) { + continue; + } + TGeoVolume* child = node->GetVolume(); + if (!child) { + continue; + } + + visit_volume(child); + } + }; + + visit_volume(top); +} + +} // namespace o2::base diff --git a/Detectors/CMakeLists.txt b/Detectors/CMakeLists.txt index eef692ff18ca7..9143761508998 100644 --- a/Detectors/CMakeLists.txt +++ b/Detectors/CMakeLists.txt @@ -25,6 +25,7 @@ add_subdirectory(TOF) add_subdirectory(ZDC) add_subdirectory(ITSMFT) +add_subdirectory(External) # sensitive external (CAD-derived) detectors; uses ITSMFT hit type add_subdirectory(TRD) add_subdirectory(MUON) diff --git a/Detectors/External/CMakeLists.txt b/Detectors/External/CMakeLists.txt new file mode 100644 index 0000000000000..ea5f0b53e2b8e --- /dev/null +++ b/Detectors/External/CMakeLists.txt @@ -0,0 +1,22 @@ +# 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. + +o2_add_library(ExternalDetectors + SOURCES src/ExternalDetector.cxx + PUBLIC_LINK_LIBRARIES O2::DetectorsBase + O2::SimulationDataFormat + O2::CommonUtils + RapidJSON::RapidJSON) + +o2_target_root_dictionary(ExternalDetectors + HEADERS include/ExternalDetectors/Hit.h + include/ExternalDetectors/ExternalDetector.h + LINKDEF src/ExternalDetectorsLinkDef.h) diff --git a/Detectors/External/include/ExternalDetectors/ExternalDetector.h b/Detectors/External/include/ExternalDetectors/ExternalDetector.h new file mode 100644 index 0000000000000..fbfc30bd5e148 --- /dev/null +++ b/Detectors/External/include/ExternalDetectors/ExternalDetector.h @@ -0,0 +1,176 @@ +// 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 ExternalDetector.h +/// \brief Sensitive detector built from an externally provided (CAD-derived) geometry +/// +/// ExternalDetector is the sensitive counterpart of o2::passive::ExternalModule. +/// It injects a CAD-derived TGeo geometry (produced by scripts/geometry/O2_CADtoTGeo.py) +/// and turns a configurable set of its volumes (selected by medium or volume name) into +/// sensitive volumes which produce hits. It derives from o2::base::DetImpl, so it +/// transparently participates in the full o2-sim hit forwarding/merging machinery +/// (FairMQ serialization, sub-event merging, ...). +/// +/// All instances share one generic hit type (o2::ext::Hit) so that an arbitrary number +/// of external detectors can coexist (each tied to a different o2::detectors::DetID) +/// without the hit merger needing to know more than this single wire format. +/// +/// The sensitive action itself is configurable: by default a generic entrance/exit hit +/// is produced, but a user can instead provide a ROOT macro (loaded at runtime via +/// o2::conf::GetFromMacro, the same mechanism used for generator/stepping hooks) whose +/// function receives the detector instance and may query the TVirtualMC singleton and +/// call addHit(...) to implement an arbitrary sensitive action -- without recompiling O2. + +#ifndef ALICEO2_EXT_EXTERNALDETECTOR_H +#define ALICEO2_EXT_EXTERNALDETECTOR_H + +#include "DetectorsBase/Detector.h" // for DetImpl +#include "DetectorsCommonDataFormats/DetID.h" // for DetID +#include "ExternalDetectors/Hit.h" // for the generic external hit type + +#include "Rtypes.h" +#include "TLorentzVector.h" + +#include +#include +#include +#include +#include + +class FairVolume; +class TGeoMatrix; +class TVector3; + +namespace o2::ext +{ + +/// Configuration of a single sensitive external detector. +struct ExternalDetectorOptions { + std::string root_macro_file; // ROOT macro describing the CAD geometry (O2_CADtoTGeo.py output) + std::string anchor_volume; // existing volume into which the geometry is hooked + TGeoMatrix const* placement = nullptr; // placement of the geometry inside the anchor (may be null) + std::vector sensitiveMedia; // media (substring match on the medium name) to be made sensitive + std::vector sensitiveVolumes; // volumes (substring match on the volume name) to be made sensitive + int detID = o2::detectors::DetID::ITS; // DetID this detector's hits are tied to (identity / output format) + std::string sensitiveMacro; // optional ROOT macro implementing the sensitive action + std::string sensitiveFunction; // global function in the macro returning the action (default "sensitiveAction()") +}; + +class ExternalDetector : public o2::base::DetImpl +{ + public: + /// Signature of a (JIT-able) sensitive action. The function is handed the detector + /// instance and is expected to query the TVirtualMC singleton (TVirtualMC::GetMC()) + /// for the current step and to call addHit(...) to produce hits. Returning true means + /// a hit-relevant step was processed (mirrors the ProcessHits return value). + using SensitiveFcn = std::function; + + ExternalDetector(const char* name, const char* title, ExternalDetectorOptions options); + ExternalDetector(); + ~ExternalDetector() override; + + /// Build a list of sensitive external detectors from a JSON description file. + /// The file must contain an "externalDetectors" array; each entry needs at least + /// "name", "macro", "anchor" and at least one of "sensitiveMedia" / "sensitiveVolumes" + /// (arrays of substrings matched against medium / volume names); an optional + /// "detID" (name, default "ITS") ties the hit output to an existing detector, and + /// an optional "placement" object may carry "translation"/"rotation_deg". + /// Ownership of the returned detectors is transferred to the caller. + static std::vector createFromJSON(const std::string& jsonfile); + + /// Build the CAD geometry, remap its media and register the sensitive volumes. + void ConstructGeometry() override; + + /// Resolve the Monte Carlo volume IDs of the sensitive volumes. + void InitializeO2Detector() override; + + /// Called for each tracking step; produces hits in the sensitive volumes. + Bool_t ProcessHits(FairVolume* v = nullptr) override; + + /// Register the hit collection with the FairRootManager. + void Register() override; + + /// Get the produced hit collection (probe interface used by DetImpl). + std::vector* getHits(Int_t iColl) const + { + if (iColl == 0) { + return mHits; + } + return nullptr; + } + + void Reset() override; + void EndOfEvent() override; + + void FinishPrimary() override {} + void BeginPrimary() override {} + void PostTrack() override {} + void PreTrack() override {} + + /// \name Helpers usable from a user-provided sensitive-action macro + /// These wrap the bookkeeping a sensitive action typically needs so that a macro can + /// stay focused on physics and the TVirtualMC queries. + ///@{ + /// Append a hit to the output collection and flag the MCTrack as having left a hit + /// in this detector. Returns a pointer to the stored hit. + o2::ext::Hit* addHit(int trackID, int sensorID, const TVector3& startPos, const TVector3& endPos, + const TVector3& startMom, double startE, double endTime, double eLoss, + unsigned char startStatus, unsigned char endStatus, int pdg = 0, float length = 0.f); + + /// Running sensor index of the volume currently being processed, or -1 if the current + /// volume is not one of the configured sensitive volumes. + int currentSensorID() const; + + /// MCTrack number of the track currently being stepped. + int currentTrackID() const; + ///@} + + protected: + /// the built-in sensitive action used when no macro is configured (generic entrance/exit hit) + Bool_t defaultProcessHits(); + + /// recursively collect names of volumes whose medium matches the configured sensitive media + void collectSensitiveVolumeNames(TGeoVolume* vol, std::set& visited); + + ExternalDetectorOptions mOptions; + + std::vector mSensitiveVolumeNames; //! names of the volumes to be made sensitive (filled at geometry build) + std::set mSensitiveVolIDs; //! MC volume IDs of the sensitive volumes + std::unordered_map mVolID2SensorID; //! dense sensor index per sensitive MC volume ID + + /// transient data about a track passing a sensor (mirrors the ITS approach) + struct TrackData { + bool mHitStarted; //! hit creation started + unsigned char mTrkStatusStart; //! track status flag at entrance + TLorentzVector mPositionStart; //! position at entrance + TLorentzVector mMomentumStart; //! momentum at entrance + double mEnergyLoss; //! accumulated energy loss + } mTrackData; //! + + std::vector* mHits = nullptr; //! container for produced hits + + SensitiveFcn mSensitiveAction; //! optional user-provided sensitive action (loaded from a macro) + FairVolume* mCurrentVolume = nullptr; //! volume currently passed to ProcessHits (for the action helpers) + + int mStepCount = 0; //! number of stepping calls inside our sensitive volumes this event (probe) + + private: + ExternalDetector(const ExternalDetector&); + ExternalDetector& operator=(const ExternalDetector&); + + template + friend class o2::base::DetImpl; + ClassDefOverride(ExternalDetector, 1); +}; + +} // namespace o2::ext + +#endif diff --git a/Detectors/External/include/ExternalDetectors/Hit.h b/Detectors/External/include/ExternalDetectors/Hit.h new file mode 100644 index 0000000000000..4a77f85050f06 --- /dev/null +++ b/Detectors/External/include/ExternalDetectors/Hit.h @@ -0,0 +1,140 @@ +// 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 Hit.h +/// \brief Generic hit type for externally injected (CAD-derived) sensitive detectors +/// +/// o2::ext::Hit is a deliberately rich, detector-agnostic hit. A single external +/// hit type lets an arbitrary number of o2::ext::ExternalDetector instances (each +/// tied to a different DetID) share one wire format, so that the o2-sim hit merger +/// only ever has to know how to (de)serialize this one type, independently of how +/// many external detectors are configured or what their sensitive action does. +/// +/// It stores entrance and exit position, the momentum, energy and energy loss, the +/// time, the track length in the volume, the PDG code and the MC status flags at +/// entrance/exit, so that most information an external sensitive action might want +/// to keep is available downstream. + +#ifndef ALICEO2_EXT_HIT_H +#define ALICEO2_EXT_HIT_H + +#include "SimulationDataFormat/BaseHits.h" // for BasicXYZEHit +#include "CommonUtils/ShmAllocator.h" +#include "Rtypes.h" +#include "TVector3.h" +#include + +namespace o2::ext +{ + +class Hit : public o2::BasicXYZEHit +{ + public: + enum HitStatus_t { + kTrackEntering = 0x1, + kTrackInside = 0x1 << 1, + kTrackExiting = 0x1 << 2, + kTrackOut = 0x1 << 3, + kTrackStopped = 0x1 << 4, + kTrackAlive = 0x1 << 5 + }; + + Hit() = default; + + /// \param trackID index of the MCTrack + /// \param sensorID index of the sensitive volume (per-detector running id) + /// \param startPos coordinates at entrance to the active volume [cm] + /// \param endPos coordinates at exit of the active volume [cm] + /// \param startMom momentum of the track at entrance [GeV] + /// \param startE total energy at entrance [GeV] + /// \param endTime time at exit [ns] + /// \param eLoss energy deposited in the volume [GeV] + /// \param startStatus MC status flags at entrance + /// \param endStatus MC status flags at exit + /// \param pdg PDG code of the track (optional) + /// \param length track length inside the volume [cm] (optional) + Hit(int trackID, unsigned short sensorID, const TVector3& startPos, const TVector3& endPos, + const TVector3& startMom, double startE, double endTime, double eLoss, + unsigned char startStatus, unsigned char endStatus, int pdg = 0, float length = 0.f) + : BasicXYZEHit(endPos.X(), endPos.Y(), endPos.Z(), endTime, eLoss, trackID, sensorID), + mMomentum(startMom.Px(), startMom.Py(), startMom.Pz()), + mPosStart(startPos.X(), startPos.Y(), startPos.Z()), + mE(startE), + mLength(length), + mPdg(pdg), + mTrackStatusEnd(endStatus), + mTrackStatusStart(startStatus) + { + } + + // entrance position + math_utils::Point3D GetPosStart() const { return mPosStart; } + float GetStartX() const { return mPosStart.X(); } + float GetStartY() const { return mPosStart.Y(); } + float GetStartZ() const { return mPosStart.Z(); } + void SetPosStart(const math_utils::Point3D& p) { mPosStart = p; } + + // momentum / energy + math_utils::Vector3D GetMomentum() const { return mMomentum; } + math_utils::Vector3D& GetMomentum() { return mMomentum; } + float GetPx() const { return mMomentum.X(); } + float GetPy() const { return mMomentum.Y(); } + float GetPz() const { return mMomentum.Z(); } + float GetE() const { return mE; } + float GetTotalEnergy() const { return mE; } + + // extra bookkeeping + float GetLength() const { return mLength; } + void SetLength(float l) { mLength = l; } + int GetPdg() const { return mPdg; } + void SetPdg(int pdg) { mPdg = pdg; } + + // status flags + unsigned char GetStatusStart() const { return mTrackStatusStart; } + unsigned char GetStatusEnd() const { return mTrackStatusEnd; } + bool IsEntering() const { return mTrackStatusEnd & kTrackEntering; } + bool IsInside() const { return mTrackStatusEnd & kTrackInside; } + bool IsExiting() const { return mTrackStatusEnd & kTrackExiting; } + bool IsOut() const { return mTrackStatusEnd & kTrackOut; } + bool IsStopped() const { return mTrackStatusEnd & kTrackStopped; } + bool IsAlive() const { return mTrackStatusEnd & kTrackAlive; } + + friend std::ostream& operator<<(std::ostream& of, const Hit& point) + { + of << "-I- o2::ext::Hit for track " << point.GetTrackID() << " in sensor " << point.GetDetectorID(); + return of; + } + + private: + math_utils::Vector3D mMomentum; ///< momentum at entrance + math_utils::Point3D mPosStart; ///< position at entrance (base mPos holds the exit position) + float mE; ///< total energy at entrance + float mLength; ///< track length inside the volume + int mPdg; ///< PDG code of the track + unsigned char mTrackStatusEnd; ///< MC status flag at exit + unsigned char mTrackStatusStart; ///< MC status flag at entrance + + ClassDefNV(Hit, 1); +}; + +} // namespace o2::ext + +#ifdef USESHM +namespace std +{ +template <> +class allocator : public o2::utils::ShmAllocator +{ +}; +} // namespace std +#endif + +#endif // ALICEO2_EXT_HIT_H diff --git a/Detectors/External/macro/sensitiveActionExample.macro b/Detectors/External/macro/sensitiveActionExample.macro new file mode 100644 index 0000000000000..0f496bfcece0a --- /dev/null +++ b/Detectors/External/macro/sensitiveActionExample.macro @@ -0,0 +1,68 @@ +// 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 sensitiveActionExample.macro +/// \brief Example sensitive action for an o2::ext::ExternalDetector +/// +/// Point this macro at an "externalDetectors" entry via the "sensitiveMacro" key +/// (an absolute path, or one using shell variables such as $O2_ROOT, is resolved at runtime): +/// "sensitiveMacro": ".../Detectors/External/macro/sensitiveActionExample.macro" +/// +/// The global function sensitiveAction() returns the callable that o2-sim invokes for +/// every tracking step inside one of the detector's sensitive volumes. It is loaded at +/// runtime with o2::conf::GetFromMacro (the same just-in-time mechanism used for +/// generator and stepping hooks), so the sensitive logic can be changed without +/// recompiling O2. +/// +/// Inside the action you have the full TVirtualMC singleton available (exactly like a +/// hand-written ProcessHits) plus a few convenience helpers on the detector: +/// - det->currentSensorID() : running index of the current sensitive volume (-1 if none) +/// - det->currentTrackID() : MCTrack number of the track being stepped +/// - det->addHit(...) : append an o2::ext::Hit and flag the MCTrack +/// +/// This trivial example records one hit each time a charged particle enters a sensitive +/// volume. + +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include "ExternalDetectors/ExternalDetector.h" +#include "ExternalDetectors/Hit.h" +#include +#include +#include +#endif + +// NOTE: the return type must be spelled exactly as the typedef name passed to +// GetFromMacro ("o2::ext::ExternalDetector::SensitiveFcn") -- the loader compares the +// function's return-type name textually, so std::function<...> would not match. +o2::ext::ExternalDetector::SensitiveFcn sensitiveAction() +{ + return [](o2::ext::ExternalDetector* det) -> bool { + auto vmc = TVirtualMC::GetMC(); + if (vmc->TrackCharge() == 0) { + return false; // ignore neutral particles + } + if (!vmc->IsTrackEntering()) { + return false; // record only the entrance point in this example + } + const int sensor = det->currentSensorID(); + if (sensor < 0) { + return false; // current volume is not one of our sensitive volumes + } + TLorentzVector pos, mom; + vmc->TrackPosition(pos); + vmc->TrackMomentum(mom); + // a point-like hit at the entrance (start == end position, no accumulated energy loss) + det->addHit(det->currentTrackID(), sensor, pos.Vect(), pos.Vect(), mom.Vect(), + mom.E(), pos.T(), 0. /*eLoss*/, o2::ext::Hit::kTrackEntering, o2::ext::Hit::kTrackEntering, + vmc->TrackPid(), vmc->TrackLength()); + return true; + }; +} diff --git a/Detectors/External/src/ExternalDetector.cxx b/Detectors/External/src/ExternalDetector.cxx new file mode 100644 index 0000000000000..f79d6ce6f2505 --- /dev/null +++ b/Detectors/External/src/ExternalDetector.cxx @@ -0,0 +1,451 @@ +// 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. + +#include "ExternalDetectors/ExternalDetector.h" +#include "DetectorsBase/CADGeometryUtils.h" +#include "DetectorsBase/Stack.h" +#include "CommonUtils/ConfigurationMacroHelper.h" +#include "CommonUtils/FileSystemUtils.h" +#include "CommonUtils/ShmManager.h" +#include "CommonUtils/ShmAllocator.h" + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +namespace o2::ext +{ + +ExternalDetector::ExternalDetector(const char* name, const char* title, ExternalDetectorOptions options) + : o2::base::DetImpl(name, true), + mOptions(options), + mTrackData(), + mHits(o2::utils::createSimVector()) +{ + (void)title; // the FairModule title is the second base ctor argument; kept for symmetry with other detectors + // Decouple the user-facing FairModule name (e.g. "IRIS") from the DetId: the base + // ctor derives fDetId from the name which is generally not a registered DetID, so we + // explicitly tie this detector to the configured (existing) DetID. This is what makes + // the hit output format / identity well defined, as discussed. + fDetId = mOptions.detID; +} + +ExternalDetector::ExternalDetector() + : o2::base::DetImpl("EXTDET", true), + mTrackData(), + mHits(o2::utils::createSimVector()) +{ +} + +ExternalDetector::ExternalDetector(const ExternalDetector& rhs) + : o2::base::DetImpl(rhs), + mOptions(rhs.mOptions), + mSensitiveVolumeNames(rhs.mSensitiveVolumeNames), + mSensitiveVolIDs(rhs.mSensitiveVolIDs), + mVolID2SensorID(rhs.mVolID2SensorID), + mTrackData(), + mHits(o2::utils::createSimVector()) +{ +} + +ExternalDetector::~ExternalDetector() +{ + if (mHits) { + o2::utils::freeSimVector(mHits); + } +} + +void ExternalDetector::collectSensitiveVolumeNames(TGeoVolume* vol, std::set& visited) +{ + if (!vol || visited.count(vol)) { + return; + } + visited.insert(vol); + + bool sensitive = false; + // match by volume name + const std::string volname = vol->GetName(); + for (const auto& token : mOptions.sensitiveVolumes) { + if (!token.empty() && volname.find(token) != std::string::npos) { + sensitive = true; + break; + } + } + // otherwise match by medium name + if (!sensitive) { + if (auto medium = vol->GetMedium()) { + const std::string medname = medium->GetName(); + for (const auto& token : mOptions.sensitiveMedia) { + if (!token.empty() && medname.find(token) != std::string::npos) { + sensitive = true; + break; + } + } + } + } + if (sensitive) { + mSensitiveVolumeNames.emplace_back(volname); + } + + const int nd = vol->GetNdaughters(); + for (int i = 0; i < nd; ++i) { + if (auto node = vol->GetNode(i)) { + collectSensitiveVolumeNames(node->GetVolume(), visited); + } + } +} + +void ExternalDetector::ConstructGeometry() +{ + // build the CAD geometry and obtain its top volume + auto module_top = o2::base::buildCADVolumeFromMacro(mOptions.root_macro_file, GetName()); + if (!module_top) { + LOG(error) << "No geometry could be built for external detector " << GetName(); + return; + } + + // bring the CAD media under O2's MaterialManager + o2::base::remapCADMedia(module_top, GetName()); + + // determine which volumes should become sensitive (selected by medium name) + mSensitiveVolumeNames.clear(); + std::set visited; + collectSensitiveVolumeNames(module_top, visited); + if (mSensitiveVolumeNames.empty()) { + LOG(warning) << "External detector " << GetName() << ": no volume matched the configured sensitive media; " + << "no hits will be produced"; + } else { + LOG(info) << "External detector " << GetName() << ": " << mSensitiveVolumeNames.size() + << " sensitive volume(s) selected"; + } + + // place it into the provided anchor volume (needs to exist) + auto anchor = gGeoManager->FindVolumeFast(mOptions.anchor_volume.c_str()); + if (!anchor) { + LOG(error) << "Anchor volume " << mOptions.anchor_volume << " not found. Aborting"; + return; + } + anchor->AddNode(module_top, 1, const_cast(mOptions.placement)); +} + +void ExternalDetector::InitializeO2Detector() +{ + // resolve the MC volume IDs of the sensitive volumes and register them with FairRoot + mSensitiveVolIDs.clear(); + mVolID2SensorID.clear(); + int sensorID = 0; + for (const auto& name : mSensitiveVolumeNames) { + const int volID = registerSensitiveVolumeAndGetVolID(name); + if (volID <= 0) { + continue; + } + mSensitiveVolIDs.insert(volID); + mVolID2SensorID[volID] = sensorID++; + LOG(info) << "External detector " << GetName() << ": registered sensitive volume '" << name + << "' (MC volID " << volID << ", sensor " << mVolID2SensorID[volID] << ")"; + } + + // optionally load a user-provided sensitive action from a ROOT macro (same mechanism as + // generator/stepping hooks). When given, it fully replaces the built-in action. + if (!mOptions.sensitiveMacro.empty()) { + const auto file = o2::utils::expandShellVarsInFileName(mOptions.sensitiveMacro); + const auto func = mOptions.sensitiveFunction.empty() ? std::string("sensitiveAction()") : mOptions.sensitiveFunction; + const auto unique = std::string("o2ext_sensitive_action_") + GetName(); + mSensitiveAction = o2::conf::GetFromMacro(file, func, "o2::ext::ExternalDetector::SensitiveFcn", unique); + if (mSensitiveAction) { + LOG(info) << "External detector " << GetName() << ": using sensitive action '" << func + << "' from macro '" << file << "'"; + } else { + LOG(fatal) << "External detector " << GetName() << ": could not load sensitive action '" << func + << "' from macro '" << file << "'"; + } + } +} + +Bool_t ExternalDetector::ProcessHits(FairVolume* vol) +{ + // This method is called from the MC stepping for the registered sensitive volumes. + // Remember the current volume so the action helpers (currentSensorID()) can resolve it, + // then either run the user-provided action or the built-in one. + mCurrentVolume = vol; + ++mStepCount; // probe: count stepping calls inside our sensitive volumes + if (mSensitiveAction) { + return mSensitiveAction(this) ? kTRUE : kFALSE; + } + return defaultProcessHits(); +} + +Bool_t ExternalDetector::defaultProcessHits() +{ + if (!(fMC->TrackCharge())) { + return kFALSE; + } + + const int sensorID = currentSensorID(); + if (sensorID < 0) { + return kFALSE; // not one of our sensitive volumes + } + + bool startHit = false, stopHit = false; + unsigned char status = 0; + if (fMC->IsTrackEntering()) { + status |= o2::ext::Hit::kTrackEntering; + } + if (fMC->IsTrackInside()) { + status |= o2::ext::Hit::kTrackInside; + } + if (fMC->IsTrackExiting()) { + status |= o2::ext::Hit::kTrackExiting; + } + if (fMC->IsTrackOut()) { + status |= o2::ext::Hit::kTrackOut; + } + if (fMC->IsTrackStop()) { + status |= o2::ext::Hit::kTrackStopped; + } + if (fMC->IsTrackAlive()) { + status |= o2::ext::Hit::kTrackAlive; + } + + // track is entering or created in the volume + if ((status & o2::ext::Hit::kTrackEntering) || (status & o2::ext::Hit::kTrackInside && !mTrackData.mHitStarted)) { + startHit = true; + } else if ((status & (o2::ext::Hit::kTrackExiting | o2::ext::Hit::kTrackOut | o2::ext::Hit::kTrackStopped))) { + stopHit = true; + } + + // increment energy loss at all steps except entrance + if (!startHit) { + mTrackData.mEnergyLoss += fMC->Edep(); + } + if (!(startHit | stopHit)) { + return kFALSE; // do nothing + } + + if (startHit) { + mTrackData.mEnergyLoss = 0.; + fMC->TrackMomentum(mTrackData.mMomentumStart); + fMC->TrackPosition(mTrackData.mPositionStart); + mTrackData.mTrkStatusStart = status; + mTrackData.mHitStarted = true; + } + if (stopHit) { + TLorentzVector positionStop; + fMC->TrackPosition(positionStop); + addHit(currentTrackID(), sensorID, mTrackData.mPositionStart.Vect(), positionStop.Vect(), + mTrackData.mMomentumStart.Vect(), mTrackData.mMomentumStart.E(), positionStop.T(), + mTrackData.mEnergyLoss, mTrackData.mTrkStatusStart, status, fMC->TrackPid(), fMC->TrackLength()); + mTrackData.mHitStarted = false; + } + return kTRUE; +} + +int ExternalDetector::currentSensorID() const +{ + const int volID = mCurrentVolume ? mCurrentVolume->getMCid() : -1; + auto it = mVolID2SensorID.find(volID); + return it == mVolID2SensorID.end() ? -1 : it->second; +} + +int ExternalDetector::currentTrackID() const +{ + return static_cast(fMC->GetStack())->GetCurrentTrackNumber(); +} + +o2::ext::Hit* ExternalDetector::addHit(int trackID, int sensorID, const TVector3& startPos, const TVector3& endPos, + const TVector3& startMom, double startE, double endTime, double eLoss, + unsigned char startStatus, unsigned char endStatus, int pdg, float length) +{ + mHits->emplace_back(trackID, sensorID, startPos, endPos, startMom, startE, endTime, eLoss, + startStatus, endStatus, pdg, length); + // register that this track left a hit in our detector (sets the hit bit on the MCTrack) + static_cast(fMC->GetStack())->addHit(GetDetId()); + return &(mHits->back()); +} + +void ExternalDetector::Register() +{ + // Create a branch (named "Hit") holding the produced hits. + if (FairRootManager::Instance()) { + FairRootManager::Instance()->RegisterAny(addNameTo("Hit").data(), mHits, kTRUE); + } +} + +void ExternalDetector::Reset() +{ + if (!o2::utils::ShmManager::Instance().isOperational()) { + mHits->clear(); + } +} + +void ExternalDetector::EndOfEvent() +{ + // probe: report how often our sensitive volumes were stepped through and how many hits resulted + LOG(info) << "External detector " << GetName() << " EndOfEvent: " << mStepCount + << " sensitive step(s) -> " << (mHits ? mHits->size() : 0) << " hit(s)"; + mStepCount = 0; + Reset(); +} + +namespace +{ +// Build a TGeoCombiTrans from an optional JSON "placement" object carrying +// "translation":[x,y,z] (cm) and/or "rotation_deg":[rx,ry,rz] (deg, applied X,Y,Z). +TGeoMatrix* makePlacementFromJSON(const rapidjson::Value& placement) +{ + auto combi = new TGeoCombiTrans(); + if (placement.HasMember("rotation_deg") && placement["rotation_deg"].IsArray()) { + const auto& r = placement["rotation_deg"]; + if (r.Size() == 3) { + combi->RotateX(r[0].GetDouble()); + combi->RotateY(r[1].GetDouble()); + combi->RotateZ(r[2].GetDouble()); + } else { + LOG(warning) << "ExternalDetector placement 'rotation_deg' must have 3 entries; ignoring"; + } + } + if (placement.HasMember("translation") && placement["translation"].IsArray()) { + const auto& t = placement["translation"]; + if (t.Size() == 3) { + combi->SetDx(t[0].GetDouble()); + combi->SetDy(t[1].GetDouble()); + combi->SetDz(t[2].GetDouble()); + } else { + LOG(warning) << "ExternalDetector placement 'translation' must have 3 entries; ignoring"; + } + } + return combi; +} +} // namespace + +std::vector ExternalDetector::createFromJSON(const std::string& jsonfile) +{ + std::vector result; + + auto expanded = o2::utils::expandShellVarsInFileName(jsonfile); + std::ifstream fileStream(expanded, std::ios::in); + if (!fileStream.is_open()) { + LOG(error) << "Cannot open external geometry config file '" << expanded << "'"; + return result; + } + + rapidjson::IStreamWrapper isw(fileStream); + rapidjson::Document doc; + doc.ParseStream(isw); + if (doc.HasParseError()) { + LOG(error) << "Error parsing external geometry JSON '" << expanded << "': " + << rapidjson::GetParseError_En(doc.GetParseError()) + << " (offset " << doc.GetErrorOffset() << ")"; + return result; + } + // the array of sensitive external detectors is optional (the same file may only + // configure passive external modules) + if (!doc.HasMember("externalDetectors")) { + return result; + } + if (!doc["externalDetectors"].IsArray()) { + LOG(error) << "External geometry JSON '" << expanded << "': 'externalDetectors' must be an array"; + return result; + } + + auto getString = [](const rapidjson::Value& v, const char* key) -> std::string { + if (v.HasMember(key) && v[key].IsString()) { + return v[key].GetString(); + } + return std::string(); + }; + + for (const auto& entry : doc["externalDetectors"].GetArray()) { + if (!entry.IsObject()) { + LOG(error) << "Skipping non-object entry in 'externalDetectors'"; + continue; + } + const auto name = getString(entry, "name"); + if (name.empty()) { + LOG(error) << "Skipping external detector entry without 'name'"; + continue; + } + ExternalDetectorOptions options; + options.root_macro_file = getString(entry, "macro"); + options.anchor_volume = getString(entry, "anchor"); + if (options.root_macro_file.empty() || options.anchor_volume.empty()) { + LOG(error) << "External detector '" << name << "' requires both 'macro' and 'anchor'; skipping"; + continue; + } + + if (entry.HasMember("sensitiveMedia") && entry["sensitiveMedia"].IsArray()) { + for (const auto& m : entry["sensitiveMedia"].GetArray()) { + if (m.IsString()) { + options.sensitiveMedia.emplace_back(m.GetString()); + } + } + } + if (entry.HasMember("sensitiveVolumes") && entry["sensitiveVolumes"].IsArray()) { + for (const auto& v : entry["sensitiveVolumes"].GetArray()) { + if (v.IsString()) { + options.sensitiveVolumes.emplace_back(v.GetString()); + } + } + } + if (options.sensitiveMedia.empty() && options.sensitiveVolumes.empty()) { + LOG(error) << "External detector '" << name + << "' requires a non-empty 'sensitiveMedia' or 'sensitiveVolumes' array; skipping"; + continue; + } + + const auto detIDName = getString(entry, "detID"); + if (!detIDName.empty()) { + const auto did = o2::detectors::DetID::nameToID(detIDName.c_str()); + if (did < 0 || did >= o2::detectors::DetID::nDetectors) { + LOG(error) << "External detector '" << name << "': unknown detID '" << detIDName << "'; skipping"; + continue; + } + options.detID = did; + } + + if (entry.HasMember("placement") && entry["placement"].IsObject()) { + options.placement = makePlacementFromJSON(entry["placement"]); + } + + // optional user-provided sensitive action (a ROOT macro). When absent, the built-in + // generic entrance/exit hit action is used. + options.sensitiveMacro = getString(entry, "sensitiveMacro"); + options.sensitiveFunction = getString(entry, "sensitiveFunction"); + + auto title = getString(entry, "title"); + if (title.empty()) { + title = name; + } + LOG(info) << "Configured external detector '" << name << "' from macro '" << options.root_macro_file + << "' anchored to '" << options.anchor_volume << "', tied to DetID '" + << o2::detectors::DetID::getName(options.detID) << "'"; + result.push_back(new ExternalDetector(name.c_str(), title.c_str(), options)); + } + return result; +} + +} // namespace o2::ext + +ClassImp(o2::ext::ExternalDetector); diff --git a/Detectors/External/src/ExternalDetectorsLinkDef.h b/Detectors/External/src/ExternalDetectorsLinkDef.h new file mode 100644 index 0000000000000..f6e13b70200fc --- /dev/null +++ b/Detectors/External/src/ExternalDetectorsLinkDef.h @@ -0,0 +1,23 @@ +// 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. + +#ifdef __CLING__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class o2::ext::Hit + ; +#pragma link C++ class std::vector < o2::ext::Hit> + ; +#pragma link C++ class o2::ext::ExternalDetector + ; +#pragma link C++ class o2::base::DetImpl < o2::ext::ExternalDetector> + ; + +#endif diff --git a/Detectors/Passive/CMakeLists.txt b/Detectors/Passive/CMakeLists.txt index a24954ad10539..e3bc45f8e8153 100644 --- a/Detectors/Passive/CMakeLists.txt +++ b/Detectors/Passive/CMakeLists.txt @@ -24,7 +24,8 @@ o2_add_library(DetectorsPassive src/HallSimParam.cxx src/PassiveBase.cxx src/ExternalModule.cxx - PUBLIC_LINK_LIBRARIES O2::Field O2::DetectorsBase O2::SimConfig) + PUBLIC_LINK_LIBRARIES O2::Field O2::DetectorsBase O2::SimConfig + RapidJSON::RapidJSON) o2_target_root_dictionary(DetectorsPassive HEADERS include/DetectorsPassive/Absorber.h diff --git a/Detectors/Passive/include/DetectorsPassive/ExternalModule.h b/Detectors/Passive/include/DetectorsPassive/ExternalModule.h index 155870ae42a6d..01c6dfea16947 100644 --- a/Detectors/Passive/include/DetectorsPassive/ExternalModule.h +++ b/Detectors/Passive/include/DetectorsPassive/ExternalModule.h @@ -14,6 +14,8 @@ #include "DetectorsPassive/PassiveBase.h" // base class of passive modules #include "Rtypes.h" // for Pipe::Class, ClassDef, Pipe::Streamer +#include +#include class TGeoVolume; class TGeoTransformation; @@ -41,22 +43,23 @@ class ExternalModule : public PassiveBase ~ExternalModule() override = default; void ConstructGeometry() override; + /// Build a list of external (passive) modules from a JSON description file. + /// The file must contain an "externalModules" array; each entry needs at least + /// "name", "macro" and "anchor"; an optional "placement" object may carry + /// "translation":[x,y,z] (cm) and "rotation_deg":[rx,ry,rz] (degrees). + /// Ownership of the returned modules is transferred to the caller. + static std::vector createFromJSON(const std::string& jsonfile); + /// Clone this object (used in MT mode only) FairModule* CloneModule() const override { return nullptr; } - typedef std::function GeomBuilderFcn; // function hook for external geometry builder - private: // void createMaterials(); ExternalModule(const ExternalModule& orig); ExternalModule& operator=(const ExternalModule&); - GeomBuilderFcn mGeomHook; ExternalModuleOptions mOptions; - bool initGeomBuilderHook(); // function to load/JIT Geometry builder hook - void remapMedia(TGeoVolume* vol); // performs a remapping of materials/media IDs after registration with VMC - // ClassDefOverride(ExternalModule, 0); }; } // namespace passive diff --git a/Detectors/Passive/src/ExternalModule.cxx b/Detectors/Passive/src/ExternalModule.cxx index fc6bd6953b82d..ebfa405d877de 100644 --- a/Detectors/Passive/src/ExternalModule.cxx +++ b/Detectors/Passive/src/ExternalModule.cxx @@ -12,16 +12,15 @@ // Sandro Wenzel (CERN), 2026 #include -#include -#include +#include +#include #include #include #include -#include -#include -#include -#include -#include +#include +#include +#include +#include // ClassImp(o2::passive::ExternalModule) @@ -32,121 +31,17 @@ ExternalModule::ExternalModule(const char* name, const char* long_title, Externa { } -void ExternalModule::remapMedia(TGeoVolume* top_volume) -{ - std::unordered_map medium_ptr_mapping; - std::unordered_set volumes_already_treated; - int counter = 1; - - auto modulename = GetName(); - - // The transformer function - auto transform_media = [&](TGeoVolume* vol_) { - if (volumes_already_treated.find(vol_) != volumes_already_treated.end()) { - // this volume was already transformed - return; - } - volumes_already_treated.insert(vol_); - - if (dynamic_cast(vol_)) { - // do nothing for assemblies (they don't have a medium) - return; - } - - auto medium = vol_->GetMedium(); - if (!medium) { - return; - } - - auto iter = medium_ptr_mapping.find(medium); - if (iter != medium_ptr_mapping.end()) { - // This medium has already been transformed, so - // we just update the volume - vol_->SetMedium(iter->second); - return; - } else { - std::cout << "Transforming media with name " << medium->GetName() << " for volume " << vol_->GetName() << "\n"; - - // we found a medium, not yet treated - auto curr_mat = medium->GetMaterial(); - auto& matmgr = o2::base::MaterialManager::Instance(); - - matmgr.Material(modulename, counter, curr_mat->GetName(), curr_mat->GetA(), curr_mat->GetZ(), curr_mat->GetDensity(), curr_mat->GetRadLen(), curr_mat->GetIntLen()); - // TGeo medium params are stored in a flat array with the following convention - // fParams[0] = isvol; - // fParams[1] = ifield; - // fParams[2] = fieldm; - // fParams[3] = tmaxfd; - // fParams[4] = stemax; - // fParams[5] = deemax; - // fParams[6] = epsil; - // fParams[7] = stmin; - const auto isvol = medium->GetParam(0); - const auto isxfld = medium->GetParam(1); - const auto sxmgmx = medium->GetParam(2); - const auto tmaxfd = medium->GetParam(3); - const auto stemax = medium->GetParam(4); - const auto deemax = medium->GetParam(5); - const auto epsil = medium->GetParam(6); - const auto stmin = medium->GetParam(7); - - matmgr.Medium(modulename, counter, medium->GetName(), counter, isvol, isxfld, sxmgmx, tmaxfd, stemax, deemax, epsil, stmin); - - // there will be new Material and Medium objects; fetch them - auto new_med = matmgr.getTGeoMedium(modulename, counter); - - // insert into cache - medium_ptr_mapping[medium] = new_med; - vol_->SetMedium(new_med); - counter++; - } - }; // end transformer lambda - - // a generic volume walker - std::function visit_volume; - visit_volume = [&](TGeoVolume* vol) -> void { - if (!vol) { - return; - } - - // call the transformer - transform_media(vol); - - // Recurse into daughters - const int nd = vol->GetNdaughters(); - for (int i = 0; i < nd; ++i) { - TGeoNode* node = vol->GetNode(i); - if (!node) { - continue; - } - TGeoVolume* child = node->GetVolume(); - if (!child) { - continue; - } - - visit_volume(child); - } - }; - - visit_volume(top_volume); -} - void ExternalModule::ConstructGeometry() { - // JIT the geom builder hook - if (!initGeomBuilderHook()) { - LOG(error) << " Could not load geometry builder hook"; - return; - } - - // otherwise execute it and obtain pointer to top most module volume - auto module_top = mGeomHook(); + // JIT the geom builder macro and obtain the top most module volume + auto module_top = o2::base::buildCADVolumeFromMacro(mOptions.root_macro_file, GetName()); if (!module_top) { - LOG(error) << "No module found\n"; + LOG(error) << "No module geometry could be built from " << mOptions.root_macro_file; return; } - remapMedia(const_cast(module_top)); + // bring the CAD media under O2's MaterialManager + o2::base::remapCADMedia(module_top, GetName()); // place it into the provided anchor volume (needs to exist) auto anchor = gGeoManager->FindVolumeFast(mOptions.anchor_volume.c_str()); @@ -154,22 +49,101 @@ void ExternalModule::ConstructGeometry() LOG(error) << "Anchor volume " << mOptions.anchor_volume << " not found. Aborting"; return; } - anchor->AddNode(const_cast(module_top), 1, const_cast(mOptions.placement)); + anchor->AddNode(module_top, 1, const_cast(mOptions.placement)); +} + +namespace +{ +// Build a TGeoCombiTrans from an optional JSON "placement" object carrying +// "translation":[x,y,z] (cm) and/or "rotation_deg":[rx,ry,rz] (deg, applied X,Y,Z). +TGeoMatrix* makePlacementFromJSON(const rapidjson::Value& placement) +{ + auto combi = new TGeoCombiTrans(); + if (placement.HasMember("rotation_deg") && placement["rotation_deg"].IsArray()) { + const auto& r = placement["rotation_deg"]; + if (r.Size() == 3) { + combi->RotateX(r[0].GetDouble()); + combi->RotateY(r[1].GetDouble()); + combi->RotateZ(r[2].GetDouble()); + } else { + LOG(warning) << "ExternalModule placement 'rotation_deg' must have 3 entries; ignoring"; + } + } + if (placement.HasMember("translation") && placement["translation"].IsArray()) { + const auto& t = placement["translation"]; + if (t.Size() == 3) { + combi->SetDx(t[0].GetDouble()); + combi->SetDy(t[1].GetDouble()); + combi->SetDz(t[2].GetDouble()); + } else { + LOG(warning) << "ExternalModule placement 'translation' must have 3 entries; ignoring"; + } + } + return combi; } +} // namespace -bool ExternalModule::initGeomBuilderHook() +std::vector ExternalModule::createFromJSON(const std::string& jsonfile) { - if (mOptions.root_macro_file.size() > 0) { - LOG(info) << "Initializing the hook for geometry module building"; - auto expandedHookFileName = o2::utils::expandShellVarsInFileName(mOptions.root_macro_file); - if (std::filesystem::exists(expandedHookFileName)) { - // if this file exists we will compile the hook on the fly (the last one is an identifier --> maybe make it dependent on this class) - mGeomHook = o2::conf::GetFromMacro(mOptions.root_macro_file, "get_builder_hook_unchecked()", "function", "o2_passive_extmodule_builder"); - LOG(info) << "Hook initialized from file " << expandedHookFileName; - return true; + std::vector result; + + auto expanded = o2::utils::expandShellVarsInFileName(jsonfile); + std::ifstream fileStream(expanded, std::ios::in); + if (!fileStream.is_open()) { + LOG(error) << "Cannot open external geometry config file '" << expanded << "'"; + return result; + } + + rapidjson::IStreamWrapper isw(fileStream); + rapidjson::Document doc; + doc.ParseStream(isw); + if (doc.HasParseError()) { + LOG(error) << "Error parsing external geometry JSON '" << expanded << "': " + << rapidjson::GetParseError_En(doc.GetParseError()) + << " (offset " << doc.GetErrorOffset() << ")"; + return result; + } + if (!doc.HasMember("externalModules") || !doc["externalModules"].IsArray()) { + LOG(error) << "External geometry JSON '" << expanded << "' must contain an 'externalModules' array"; + return result; + } + + auto getString = [](const rapidjson::Value& v, const char* key) -> std::string { + if (v.HasMember(key) && v[key].IsString()) { + return v[key].GetString(); + } + return std::string(); + }; + + for (const auto& entry : doc["externalModules"].GetArray()) { + if (!entry.IsObject()) { + LOG(error) << "Skipping non-object entry in 'externalModules'"; + continue; + } + const auto name = getString(entry, "name"); + if (name.empty()) { + LOG(error) << "Skipping external module entry without 'name'"; + continue; + } + ExternalModuleOptions options; + options.root_macro_file = getString(entry, "macro"); + options.anchor_volume = getString(entry, "anchor"); + if (options.root_macro_file.empty() || options.anchor_volume.empty()) { + LOG(error) << "External module '" << name << "' requires both 'macro' and 'anchor'; skipping"; + continue; + } + if (entry.HasMember("placement") && entry["placement"].IsObject()) { + options.placement = makePlacementFromJSON(entry["placement"]); + } + auto title = getString(entry, "title"); + if (title.empty()) { + title = name; } + LOG(info) << "Configured external module '" << name << "' from macro '" << options.root_macro_file + << "' anchored to volume '" << options.anchor_volume << "'"; + result.push_back(new ExternalModule(name.c_str(), title.c_str(), options)); } - return false; + return result; } } // namespace o2::passive \ No newline at end of file diff --git a/macro/build_geometry.C b/macro/build_geometry.C index e46b45254d308..d9d69f82363dc 100644 --- a/macro/build_geometry.C +++ b/macro/build_geometry.C @@ -64,6 +64,7 @@ #endif #include +#include using Return = o2::base::Detector*; @@ -184,16 +185,18 @@ void build_geometry(FairRunSim* run = nullptr) } #endif - if (isActivated("EXT")) { - // EXAMPLE!! how to pick geometry generated from external (CAD) module via `O2_CADtoTGeo.py` - o2::passive::ExternalModuleOptions options; - options.root_macro_file = "PATH_TO_EXTERNAL_GEOM_MODULE/geom.C"; - options.anchor_volume = "barrel"; // hook this into barrel - auto rot = new TGeoCombiTrans(); - rot->RotateX(90); - rot->SetDy(30); // we need to compensate for a shift of barrel with respect to zero - options.placement = rot; - run->AddModule(new o2::passive::ExternalModule("FOO", "BAR", options)); + // external (e.g. CAD-derived) geometry modules are injected from the outside via a JSON + // description file given with `--extGeomFile` (geometry generated via `O2_CADtoTGeo.py`). + // Each module is added when its 'name' is part of the active module list (so it can be + // switched on/off via the detector list, like any other module). + if (auto extGeomFile = confref.getExtGeomFilename(); !extGeomFile.empty()) { + for (auto* extmod : o2::passive::ExternalModule::createFromJSON(extGeomFile)) { + if (isActivated(extmod->GetName())) { + run->AddModule(extmod); + } else { + delete extmod; // not requested in the active module list + } + } } // the absorber @@ -226,6 +229,19 @@ void build_geometry(FairRunSim* run = nullptr) } }; + // sensitive external (CAD-derived) detectors, injected from the same JSON used for passive + // external modules (entries under "externalDetectors"). These derive from o2::base::Detector, + // so they produce hits and participate in the regular hit forwarding/merging machinery. + if (auto extGeomFile = confref.getExtGeomFilename(); !extGeomFile.empty()) { + for (auto* extdet : o2::ext::ExternalDetector::createFromJSON(extGeomFile)) { + if (isActivated(extdet->GetName())) { + addReadoutDetector(extdet); + } else { + delete extdet; // not requested in the active module list + } + } + } + if (isActivated("TOF")) { // TOF addReadoutDetector(new o2::tof::Detector(isReadout("TOF"))); diff --git a/run/CMakeLists.txt b/run/CMakeLists.txt index 474aa7e41eb7c..973e17ebe2a8c 100644 --- a/run/CMakeLists.txt +++ b/run/CMakeLists.txt @@ -22,6 +22,7 @@ target_link_libraries(allsim FairMQ::FairMQ O2::CPVSimulation O2::DetectorsPassive + O2::ExternalDetectors O2::EMCALSimulation O2::FDDSimulation O2::Field diff --git a/run/O2HitMerger.h b/run/O2HitMerger.h index 2f54a8b5d202a..a0a79ac4fef96 100644 --- a/run/O2HitMerger.h +++ b/run/O2HitMerger.h @@ -59,6 +59,7 @@ #include #include #include +#include #include "CommonUtils/ShmManager.h" #include @@ -890,6 +891,7 @@ class O2HitMerger : public fair::mq::Device int mPipeToDriver = -1; std::vector> mDetectorInstances; //! + std::vector mExternalDetIDs; //! DetID slots occupied by external (CAD) detectors // output folder configuration std::string mInitialOutputDir; // initial output folder of the process (initialized during construction) @@ -900,6 +902,7 @@ class O2HitMerger : public fair::mq::Device // init detector instances void initDetInstances(); + void initExternalDetInstances(); void initHitFiles(std::string prefix); }; @@ -921,6 +924,12 @@ void O2HitMerger::initHitFiles(std::string prefix) // init the detector specific output files initHitTreeAndOutFile(prefix, i); } + + // external (CAD) detectors are not part of the readout-detector list (their module names + // are not DetID names); their slots were determined in initDetInstances() + for (auto detID : mExternalDetIDs) { + initHitTreeAndOutFile(prefix, detID); + } } // init detector instances used to write hit data to a TTree @@ -1052,6 +1061,57 @@ void O2HitMerger::initDetInstances() if (counter != DetID::nDetectors) { LOG(warning) << " O2HitMerger: Some Detectors are potentially missing in this initialization "; } + + // also register external (CAD-derived) sensitive detectors so their hits are persisted + // in parallel (multi-worker) mode + initExternalDetInstances(); +} + +// init detector instances for external (CAD-derived) sensitive detectors. +// These are not part of the hard-coded DetID switch above: they are described in the +// external geometry JSON (the same file used by build_geometry.C on the worker side) and +// tied to an existing (free) DetID. The merger only needs an instance able to interpret the +// generic o2::ext::Hit wire format and write the "Hit" branch; no geometry is built here. +void O2HitMerger::initExternalDetInstances() +{ + using o2::detectors::DetID; + + auto& simConfig = o2::conf::SimConfig::Instance(); + const auto extGeomFile = simConfig.getExtGeomFilename(); + if (extGeomFile.empty()) { + return; + } + + // mirror the worker-side activation: an external detector participates when its module + // name is part of the active module list + auto const& activeModules = simConfig.getActiveModules(); + auto isActivated = [&activeModules](std::string const& s) -> bool { + return std::find(activeModules.begin(), activeModules.end(), s) != activeModules.end(); + }; + + for (auto* extdet : o2::ext::ExternalDetector::createFromJSON(extGeomFile)) { + const std::string name = extdet->GetName(); + if (!isActivated(name)) { + delete extdet; // not requested in the active module list + continue; + } + const int detID = extdet->GetDetId(); + if (detID < DetID::First || detID > DetID::Last) { + LOG(error) << "O2HitMerger: external detector " << name << " has invalid DetID " << detID << "; skipping"; + delete extdet; + continue; + } + if (mDetectorInstances[detID]) { + LOG(error) << "O2HitMerger: DetID " << DetID::getName(detID) << " requested by external detector " << name + << " is already occupied; its hits will not be persisted. Assign a free DetID."; + delete extdet; + continue; + } + mDetectorInstances[detID].reset(extdet); + mExternalDetIDs.emplace_back(detID); + LOG(info) << "O2HitMerger: registered external detector " << name << " on DetID " << DetID::getName(detID) + << " (branch " << name << "Hit)"; + } } } // namespace devices diff --git a/run/SimExamples/External_Sensitive_Detectors/README.md b/run/SimExamples/External_Sensitive_Detectors/README.md new file mode 100644 index 0000000000000..d36a6efa5e297 --- /dev/null +++ b/run/SimExamples/External_Sensitive_Detectors/README.md @@ -0,0 +1,54 @@ +# External sensitive detectors + +A minimal example showing how to add **sensitive** detectors to `o2-sim` without compiling +anything into O2. Two artificial detectors are injected purely from data: + +| name | geometry (runtime macro) | DetID slot | sensitive action | +|---------|-------------------------------|------------|---------------------------------| +| `ACYL` | `geometry_innerCylinder.macro`| `ITS` | built-in entrance/exit action | +| `BDISK` | `geometry_outerDisk.macro` | `TST` | custom action `sensitive_action.macro` | + +Both are instances of the single compiled `o2::ext::ExternalDetector` class and produce the +single generic hit type `o2::ext::Hit`. They differ only in their data: geometry, the DetID +slot they occupy, and (optionally) a sensitive-action macro. + +## Run + +```bash +./run.sh +``` + +This transports a few `boxgen` events and prints the per-detector hit counts, e.g. + +``` +External sensitive detector hits: + ACYLHit : ~430 hits over 5 events, mean radius 20.0 cm [o2sim_HitsITS.root] + BDISKHit : ~120 hits over 5 events, mean radius ... cm [o2sim_HitsTST.root] +``` + +## How it works + +* **Geometry** — each `"macro"` is a ROOT macro exporting `get_builder_hook_unchecked()`, the + same symbol `O2_CADtoTGeo.py` emits when converting CAD (STEP) files to TGeo. The macros here + are hand-written so the example is self-contained (no CAD binaries). The volumes whose names + match `"sensitiveVolumes"` are registered as sensitive. + +* **Sensitive action** — `ACYL` has none, so it uses the built-in action that records an + entrance/exit hit per track. `BDISK` points `"sensitiveMacro"` at a macro whose + `sensitiveAction()` returns the per-step callable; it is JIT-compiled at runtime via + `o2::conf::GetFromMacro` (the same mechanism as the generator / stepping hooks — no + recompilation, no ACLIC). The callable has the full `TVirtualMC` singleton and a few helpers + (`currentSensorID()`, `currentTrackID()`, `addHit()`). + +* **Persistence** — each detector is tied to a free **DetID** slot (`ITS`, `TST` here, chosen + because no real detector of that name is active). The DetID is the scarce resource: it fixes + the hit-file name (`o2sim_Hits.root`) and lets the hit merger instantiate a matching + receiver, so hits are written in parallel multi-worker mode just like for any built-in + detector. The branch keeps the detector's own name (`ACYLHit`, `BDISKHit`). + +## Adding more detectors + +Append entries to `externalDetectors.json` (each on a different free DetID slot) and list their +names in `detectorlist.json`. The mechanism is fully data-driven; nothing needs to be rebuilt. + +See also `Detectors/External` and `scripts/geometry/O2_CADtoTGeo.py`. diff --git a/run/SimExamples/External_Sensitive_Detectors/detectorlist.json b/run/SimExamples/External_Sensitive_Detectors/detectorlist.json new file mode 100644 index 0000000000000..2284e67f19748 --- /dev/null +++ b/run/SimExamples/External_Sensitive_Detectors/detectorlist.json @@ -0,0 +1,6 @@ +{ + "EXTEXAMPLE": [ + "ACYL", + "BDISK" + ] +} diff --git a/run/SimExamples/External_Sensitive_Detectors/externalDetectors.json b/run/SimExamples/External_Sensitive_Detectors/externalDetectors.json new file mode 100644 index 0000000000000..abccc209c4819 --- /dev/null +++ b/run/SimExamples/External_Sensitive_Detectors/externalDetectors.json @@ -0,0 +1,23 @@ +{ + "externalDetectors": [ + { + "name": "ACYL", + "title": "artificial barrel cylinder - built-in entrance/exit action", + "macro": "geometry_innerCylinder.macro", + "anchor": "barrel", + "detID": "ITS", + "sensitiveVolumes": ["ACYL_SENS"], + "placement": { "translation": [0, 0, 0] } + }, + { + "name": "BDISK", + "title": "artificial endcap disk - custom JITed sensitive action", + "macro": "geometry_outerDisk.macro", + "anchor": "barrel", + "detID": "TST", + "sensitiveVolumes": ["BDISK_SENS"], + "sensitiveMacro": "sensitive_action.macro", + "placement": { "translation": [0, 0, 0] } + } + ] +} diff --git a/run/SimExamples/External_Sensitive_Detectors/geometry_innerCylinder.macro b/run/SimExamples/External_Sensitive_Detectors/geometry_innerCylinder.macro new file mode 100644 index 0000000000000..78bec83f42213 --- /dev/null +++ b/run/SimExamples/External_Sensitive_Detectors/geometry_innerCylinder.macro @@ -0,0 +1,35 @@ +// Geometry builder for the "ACYL" artificial detector: a thin silicon cylindrical shell. +// +// This is a hand-written stand-in for the macros that O2_CADtoTGeo.py produces from CAD +// files. It exports the same builder-hook symbol (get_builder_hook_unchecked) that +// o2::base::buildCADVolumeFromMacro expects, so it is injected exactly like a CAD module +// through the external-geometry JSON. Media are placeholders here; remapCADMedia() +// re-registers them under O2's MaterialManager at construction time. +#include +#include +#include +#include +#include +#include +#include + +TGeoVolume* build(bool /*check*/ = true) +{ + if (!gGeoManager) { + throw std::runtime_error("gGeoManager is null when building ACYL"); + } + auto* mat = new TGeoMaterial("SiACYL", 28.0855, 14., 2.33, 9.37, 45.5); + double params[8] = {1., 0., 0., 0., 0., 0., 0., 0.}; // isvol=1, no field, defaults + auto* med = new TGeoMedium("SiACYL", 1, mat, params); + + // 2 mm thick shell, R = 20 cm, half-length 40 cm + auto* tube = new TGeoTube("ACYL_SENS", 20.0, 20.2, 40.0); + auto* sens = new TGeoVolume("ACYL_SENS", tube, med); + + auto* top = new TGeoVolumeAssembly("ACYL"); + top->AddNode(sens, 1, nullptr); + return top; +} + +std::function get_builder_hook_checked() { return []() { return build(true); }; } +std::function get_builder_hook_unchecked() { return []() { return build(false); }; } diff --git a/run/SimExamples/External_Sensitive_Detectors/geometry_outerDisk.macro b/run/SimExamples/External_Sensitive_Detectors/geometry_outerDisk.macro new file mode 100644 index 0000000000000..13907ebb94a5b --- /dev/null +++ b/run/SimExamples/External_Sensitive_Detectors/geometry_outerDisk.macro @@ -0,0 +1,35 @@ +// Geometry builder for the "BDISK" artificial detector: a thin silicon endcap disk. +// +// Hand-written stand-in for a CAD-exported macro (see geometry_innerCylinder.macro). Exports the +// same builder-hook symbol consumed by o2::base::buildCADVolumeFromMacro. Together with ACYL +// this gives two geometrically distinct artificial detectors (a barrel shell and an endcap +// disk) sharing one generic o2::ext::Hit type but writing into separate DetID slots. +#include +#include +#include +#include +#include +#include +#include +#include + +TGeoVolume* build(bool /*check*/ = true) +{ + if (!gGeoManager) { + throw std::runtime_error("gGeoManager is null when building BDISK"); + } + auto* mat = new TGeoMaterial("SiBDISK", 28.0855, 14., 2.33, 9.37, 45.5); + double params[8] = {1., 0., 0., 0., 0., 0., 0., 0.}; // isvol=1, no field, defaults + auto* med = new TGeoMedium("SiBDISK", 1, mat, params); + + // a 2 mm thick disk: inner r = 5 cm, outer r = 40 cm, placed at z = +70 cm (endcap) + auto* disk = new TGeoTube("BDISK_SENS", 5.0, 40.0, 0.1); + auto* sens = new TGeoVolume("BDISK_SENS", disk, med); + + auto* top = new TGeoVolumeAssembly("BDISK"); + top->AddNode(sens, 1, new TGeoTranslation(0., 0., 70.0)); + return top; +} + +std::function get_builder_hook_checked() { return []() { return build(true); }; } +std::function get_builder_hook_unchecked() { return []() { return build(false); }; } diff --git a/run/SimExamples/External_Sensitive_Detectors/inspect_hits.macro b/run/SimExamples/External_Sensitive_Detectors/inspect_hits.macro new file mode 100644 index 0000000000000..caebe1496b851 --- /dev/null +++ b/run/SimExamples/External_Sensitive_Detectors/inspect_hits.macro @@ -0,0 +1,50 @@ +// Inspect the hits produced by the two artificial external detectors of this example. +// +// In parallel mode (o2-sim) the hit merger writes one file per detector, named after the +// DetID slot the detector borrows: ACYL -> o2sim_HitsITS.root, BDISK -> o2sim_HitsTST.root. +// The branch inside keeps the detector's own name: "ACYLHit" / "BDISKHit". +// +// Run: root -l -b -q inspect_hits.macro +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include "ExternalDetectors/Hit.h" +#include +#include +#include +#include +#endif + +void inspectOne(const char* fname, const char* branch) +{ + auto* f = TFile::Open(fname); + if (!f || f->IsZombie()) { + printf(" %-22s : file not found (%s)\n", branch, fname); + return; + } + auto* t = (TTree*)f->Get("o2sim"); + if (!t || !t->GetBranch(branch)) { + printf(" %-22s : no '%s' branch in %s\n", branch, branch, fname); + delete f; + return; + } + std::vector* hits = nullptr; + t->SetBranchAddress(branch, &hits); + long total = 0; + double rsum = 0.; + for (Long64_t e = 0; e < t->GetEntries(); ++e) { + t->GetEntry(e); + for (auto& h : *hits) { + ++total; + rsum += std::hypot(h.GetX(), h.GetY()); + } + } + printf(" %-22s : %ld hits over %lld events, mean radius %.1f cm [%s]\n", + branch, total, t->GetEntries(), total ? rsum / total : 0., fname); + delete f; +} + +void inspect_hits() +{ + printf("External sensitive detector hits:\n"); + inspectOne("o2sim_HitsITS.root", "ACYLHit"); // barrel cylinder, built-in action + inspectOne("o2sim_HitsTST.root", "BDISKHit"); // endcap disk, custom JITed action +} diff --git a/run/SimExamples/External_Sensitive_Detectors/run.sh b/run/SimExamples/External_Sensitive_Detectors/run.sh new file mode 100755 index 0000000000000..2cd94aa49c8c9 --- /dev/null +++ b/run/SimExamples/External_Sensitive_Detectors/run.sh @@ -0,0 +1,32 @@ +#!/usr/bin/env bash +# +# Minimal example injecting two artificial *sensitive* external detectors into o2-sim. +# +# Neither detector is compiled into O2: both are described purely by data. The geometry of +# each is built at runtime from a ROOT macro (here hand-written stand-ins for what +# O2_CADtoTGeo.py produces from CAD files), and each is tied to a free DetID slot so its hits +# are written like those of any built-in detector -- including in parallel (multi-worker) mode, +# where the hit merger instantiates a matching receiver per detector. +# +# ACYL : thin silicon barrel cylinder, DetID slot ITS, built-in entrance/exit action +# BDISK : thin silicon endcap disk, DetID slot TST, custom JITed sensitive action +# +# Both share the single generic hit type o2::ext::Hit. Multiplicity is data-driven: add more +# entries (on more free DetID slots) to externalDetectors.json and detectorlist.json. + +set -x + +# run from this example's directory so the relative macro/JSON paths resolve +cd "$(dirname "$0")" + +NWORKERS=2 +EVENTS=5 + +o2-sim -j ${NWORKERS} -n ${EVENTS} -g boxgen \ + --detectorList EXTEXAMPLE:detectorlist.json \ + --extGeomFile externalDetectors.json \ + --configKeyValues 'BoxGun.number=50' \ + > sim.log 2>&1 + +# count and locate the produced hits (one file per detector / DetID slot) +root -l -b -q inspect_hits.macro diff --git a/run/SimExamples/External_Sensitive_Detectors/sensitive_action.macro b/run/SimExamples/External_Sensitive_Detectors/sensitive_action.macro new file mode 100644 index 0000000000000..29734a8aea9f9 --- /dev/null +++ b/run/SimExamples/External_Sensitive_Detectors/sensitive_action.macro @@ -0,0 +1,44 @@ +// Custom sensitive action for the "BDISK" artificial detector. +// +// The action is JIT-compiled at runtime via o2::conf::GetFromMacro (the same mechanism used +// for generator and stepping hooks) and fully replaces ExternalDetector's built-in +// entrance/exit action. It runs for every MC step inside one of the detector's sensitive +// volumes and has the full TVirtualMC singleton available, plus a few convenience helpers on +// the detector (currentSensorID(), currentTrackID(), addHit()). +// +// IMPORTANT: the return type must be spelled exactly as the typedef name passed to +// GetFromMacro ("o2::ext::ExternalDetector::SensitiveFcn"); the loader compares the function's +// return-type name textually, so std::function<...> would not match. +// +// This example records a single point-like hit whenever a charged particle enters a sensitive +// volume. +#if !defined(__CLING__) || defined(__ROOTCLING__) +#include "ExternalDetectors/ExternalDetector.h" +#include "ExternalDetectors/Hit.h" +#include +#include +#endif + +o2::ext::ExternalDetector::SensitiveFcn sensitiveAction() +{ + return [](o2::ext::ExternalDetector* det) -> bool { + auto vmc = TVirtualMC::GetMC(); + if (vmc->TrackCharge() == 0) { + return false; // ignore neutral particles + } + if (!vmc->IsTrackEntering()) { + return false; // record only the entrance point + } + const int sensor = det->currentSensorID(); + if (sensor < 0) { + return false; // not one of our sensitive volumes + } + TLorentzVector pos, mom; + vmc->TrackPosition(pos); + vmc->TrackMomentum(mom); + det->addHit(det->currentTrackID(), sensor, pos.Vect(), pos.Vect(), mom.Vect(), + mom.E(), pos.T(), 0. /*eLoss*/, o2::ext::Hit::kTrackEntering, o2::ext::Hit::kTrackEntering, + vmc->TrackPid(), vmc->TrackLength()); + return true; + }; +} diff --git a/scripts/geometry/README.md b/scripts/geometry/README.md index 4fb2d1ec610d4..8c02cadef51a4 100644 --- a/scripts/geometry/README.md +++ b/scripts/geometry/README.md @@ -1,27 +1,229 @@ -This is the tool O2_CADtoTGeo.py which translates from geometries in STEP format (CAD export) to -TGeo. +# CAD-to-TGeo geometry import -To use the tool, setup a conda environment with python-occ core installed. -The following should work on standard linux x86: +`O2_CADtoTGeo.py` converts CAD geometries exported as STEP files into ROOT TGeo geometry. +The converter emits a small ROOT macro plus compact binary facet payloads. The generated +macro can be loaded directly in ROOT, or injected into `o2-sim` as an external passive +module or as a sensitive external detector. +The current integration path is data-driven: the CAD geometry is converted once, then a +JSON file tells `o2-sim` which generated macro to load, where to anchor it in the existing +geometry, and, for sensitive detectors, which volumes or media should produce hits. + +## Software setup + +The preferred setup is the normal ALICE software environment. The `pythonOCC` package pulls +in OpenCascade and the Python bindings needed by the converter: + +```bash +alienv enter O2sim/latest,pythonOCC/latest +python PATH_TO_ALICEO2_SOURCES/scripts/geometry/O2_CADtoTGeo.py --help ``` -# -) download miniconda into $HOME/miniconda (if not already done) -if [ ! -d $HOME/miniconda ]; then - curl -fsSL https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-aarch64.sh -o miniconda.sh - bash miniconda.sh -b -p $HOME/miniconda -fi -# -) source conda into the environment (in every shell you want to use this) -source $HOME/miniconda/etc/profile.d/conda.sh +If you are working from a local O2 checkout, `PATH_TO_ALICEO2_SOURCES` is the directory that +contains this `scripts/geometry` folder. -# -) Create an OCC environment (for OpenCacade) +For standalone studies outside the ALICE software stack, a conda environment with +`pythonocc-core` can also be used: + +```bash conda create -n occ python=3.10 -y conda activate occ - -# 3) Install OpenCascade Python bindings conda install -c conda-forge pythonocc-core -y -# 4) Run the tool, e.g. -conda activate occ python PATH_TO_ALICEO2_SOURCES/scripts/geometry/O2_CADtoTGeo.py --help -``` \ No newline at end of file +``` + +## Convert a STEP file + +For a quick, robust geometry preview, convert leaves to bounding boxes: + +```bash +mkdir -p cad_out/mydet +python PATH_TO_ALICEO2_SOURCES/scripts/geometry/O2_CADtoTGeo.py \ + my_detector.step \ + --output-folder cad_out/mydet \ + -o geom.C \ + --step-unit auto +``` + +For a more detailed faceted representation, enable meshing: + +```bash +python PATH_TO_ALICEO2_SOURCES/scripts/geometry/O2_CADtoTGeo.py \ + my_detector.step \ + --output-folder cad_out/mydet \ + -o geom.C \ + --mesh \ + --mesh-prec 0.05 \ + --step-unit auto +``` + +The output folder contains: + +- `geom.C`, a ROOT macro exporting `get_builder_hook_unchecked()` +- `facets_*.bin`, one compact triangle payload per leaf logical volume + +The generated macro is the file referenced from the external-geometry JSON examples below. +The macro and its facet binaries should stay together, because the macro loads the facet files +relative to its own location. + +## Optional material mapping + +The converter can use a BOM CSV to assign materials and, when part masses and CAD volumes are +available, derive effective densities. A Geant4 NIST material JSON dump enables richer material +and tracking-length information: + +```bash +python PATH_TO_ALICEO2_SOURCES/scripts/geometry/O2_CADtoTGeo.py \ + my_detector.step \ + --output-folder cad_out/mydet \ + -o geom.C \ + --mesh \ + --materials-csv detector_bom.csv \ + --bom-mass-unit kg \ + --g4-nist-json g4_nist_materials.json +``` + +The expected BOM rows are mechanical part rows of the form: + +```text +CAD,Mechanical/Part,,,,,,... +``` + +Material names are matched to the NIST database when possible. If a match is ambiguous or not +available, the generated macro falls back to a simple material and leaves comments in `geom.C` +for follow-up. + +## Inject passive CAD geometry into `o2-sim` + +Passive external modules are configured under an `externalModules` array. Each entry needs a +module name, the generated macro, and an anchor volume already present in the O2 geometry. An +optional placement can translate and rotate the imported geometry inside the anchor volume: + +```json +{ + "externalModules": [ + { + "name": "IRIS", + "title": "IRIS support from CAD", + "macro": "cad_out/iris/geom.C", + "anchor": "barrel", + "placement": { + "translation": [0.0, 0.0, 0.0], + "rotation_deg": [0.0, 0.0, 15.0] + } + } + ] +} +``` + +The module is only added when its `name` is present in the active detector/module list. One +way to make such a list is a detector-list JSON file: + +```json +{ + "EXTCAD": ["IRIS"] +} +``` + +Run `o2-sim` with both files: + +```bash +o2-sim -n 1 -g boxgen \ + --detectorList EXTCAD:detectorlist.json \ + --extGeomFile externalGeometry.json +``` + +Multiple passive modules can be listed in the same file. The CAD macro loader JIT-compiles each +macro into a unique namespace, so several `O2_CADtoTGeo.py` outputs can coexist even though they +export the same builder-hook symbol names. + +## Inject sensitive external CAD detectors + +Sensitive external detectors are configured under an `externalDetectors` array. They use the +same generated geometry macro, but additionally select sensitive volumes or media and bind the +detector to a free O2 `DetID` slot. All such detectors are instances of +`o2::ext::ExternalDetector` and write the generic `o2::ext::Hit` format. + +```json +{ + "externalDetectors": [ + { + "name": "ECYL", + "title": "External silicon cylinder", + "macro": "cad_out/ecyl/geom.C", + "anchor": "barrel", + "detID": "ITS", + "sensitiveVolumes": ["ECYL_SENSOR"], + "placement": { "translation": [0.0, 0.0, 0.0] } + }, + { + "name": "EDISK", + "title": "External endcap disk with custom action", + "macro": "cad_out/edisk/geom.C", + "anchor": "barrel", + "detID": "TST", + "sensitiveMedia": ["Silicon"], + "sensitiveMacro": "sensitive_action.macro", + "sensitiveFunction": "sensitiveAction()" + } + ] +} +``` + +Selection rules: + +- `sensitiveVolumes` matches substrings of TGeo volume names. +- `sensitiveMedia` matches substrings of TGeo medium names. +- At least one of the two arrays must be non-empty. + +The `detID` determines the hit-file identity, for example `o2sim_HitsITS.root` or +`o2sim_HitsTST.root`. Choose a DetID that is not already occupied by an active built-in detector. +The branch name keeps the external detector name, for example `ECYLHit` or `EDISKHit`. + +If no `sensitiveMacro` is provided, the built-in action records a charged-track entrance/exit hit. +With a custom action, the macro is JIT-compiled at runtime and must return an +`o2::ext::ExternalDetector::SensitiveFcn`. The action can query `TVirtualMC::GetMC()` and use +helpers such as `currentSensorID()`, `currentTrackID()`, and `addHit()`. + +Run the detectors just like passive modules, with their names in the detector-list JSON: + +```json +{ + "EXTCAD": ["ECYL", "EDISK"] +} +``` + +```bash +o2-sim -j 2 -n 5 -g boxgen \ + --detectorList EXTCAD:detectorlist.json \ + --extGeomFile externalGeometry.json \ + --configKeyValues 'BoxGun.number=50' +``` + +In parallel mode, the hit merger reads the same `--extGeomFile`, registers the configured active +external detectors, and persists their generic external hits like built-in detector hits. + +## Complete runnable example + +A self-contained example is available in: + +```bash +run/SimExamples/External_Sensitive_Detectors +``` + +It defines two artificial sensitive detectors entirely from data: + +- `ACYL`, a silicon barrel cylinder using the built-in entrance/exit action +- `BDISK`, a silicon endcap disk using a custom JITed sensitive action + +The example uses hand-written geometry macros that mimic `O2_CADtoTGeo.py` output, so it does not +require CAD input files. Run it from its directory: + +```bash +cd run/SimExamples/External_Sensitive_Detectors +./run.sh +``` + +The script transports a few box-generator events and prints the hit counts for the produced +external-detector branches. \ No newline at end of file diff --git a/scripts/geometry/TODO.md b/scripts/geometry/TODO.md new file mode 100644 index 0000000000000..8add645e06357 --- /dev/null +++ b/scripts/geometry/TODO.md @@ -0,0 +1,4 @@ +- implement a BVHSurface solid as an exact representation of a CAD solid +- complete geometry configurable as JSON --> even the world volume ? + + \ No newline at end of file From 1ef44171e01f6fa9e8edca6f321e36e89dda4df5 Mon Sep 17 00:00:00 2001 From: swenzel Date: Thu, 23 Jul 2026 11:02:08 +0200 Subject: [PATCH 2/2] Add CAD clip-box and name-based part selection to O2_CADtoTGeo Support culling/cutting the CAD model against an axis-aligned --clip-box (with --clip-deduplicate {none,intact}) and whitelisting/blacklisting parts via --include-name/--exclude-name regexes before TGeo conversion. Aided by GPT-5.5 (Copilot CERN pilot licence) --- scripts/geometry/O2_CADtoTGeo.py | 421 ++++++++++++++++++++++++++++--- scripts/geometry/README.md | 54 ++++ 2 files changed, 443 insertions(+), 32 deletions(-) diff --git a/scripts/geometry/O2_CADtoTGeo.py b/scripts/geometry/O2_CADtoTGeo.py index 3de2fd75973df..836488444c50e 100644 --- a/scripts/geometry/O2_CADtoTGeo.py +++ b/scripts/geometry/O2_CADtoTGeo.py @@ -50,11 +50,14 @@ import struct from dataclasses import dataclass from pathlib import Path as _Path -from typing import Dict, List, Optional, Tuple +from typing import Dict, List, Optional, Pattern, Tuple from OCC.Core.Bnd import Bnd_Box +from OCC.Core.BRepAlgoAPI import BRepAlgoAPI_Common from OCC.Core.BRepBndLib import brepbndlib +from OCC.Core.BRepBuilderAPI import BRepBuilderAPI_Transform from OCC.Core.BRepMesh import BRepMesh_IncrementalMesh +from OCC.Core.BRepPrimAPI import BRepPrimAPI_MakeBox from OCC.Core.BRep import BRep_Tool from OCC.Core.TopLoc import TopLoc_Location from OCC.Core.TopAbs import TopAbs_REVERSED @@ -67,7 +70,7 @@ from OCC.Core.TDF import TDF_Label, TDF_LabelSequence, TDF_Tool from OCC.Core.TCollection import TCollection_AsciiString -from OCC.Core.gp import gp_Trsf +from OCC.Core.gp import gp_Pnt, gp_Trsf # volume properties for density calcs (may not be present in older pythonOCC builds) try: @@ -161,6 +164,61 @@ def detect_step_length_unit(step_path: str) -> str: return "mm" +@dataclass(frozen=True) +class ClipBox: + xmin: float + ymin: float + zmin: float + xmax: float + ymax: float + zmax: float + + @classmethod + def from_values(cls, values: List[float]) -> "ClipBox": + if len(values) != 6: + raise ValueError("--clip-box expects 6 values: xmin ymin zmin xmax ymax zmax") + xmin, ymin, zmin, xmax, ymax, zmax = (float(v) for v in values) + if not (xmin < xmax and ymin < ymax and zmin < zmax): + raise ValueError("--clip-box requires xmin Tuple[float, float, float, float, float, float]: + return (self.xmin, self.ymin, self.zmin, self.xmax, self.ymax, self.zmax) + + +@dataclass(frozen=True) +class NameFilter: + include: Tuple[Pattern[str], ...] + exclude: Tuple[Pattern[str], ...] + + @classmethod + def from_patterns(cls, include: List[str], exclude: List[str], case_sensitive: bool = False) -> "NameFilter": + flags = 0 if case_sensitive else re.IGNORECASE + return cls( + tuple(re.compile(pattern, flags) for pattern in include), + tuple(re.compile(pattern, flags) for pattern in exclude), + ) + + @property + def active(self) -> bool: + return bool(self.include or self.exclude) + + @property + def has_include(self) -> bool: + return bool(self.include) + + def _text(self, lid: str, name: str) -> str: + return f"{name} {lid}".strip() + + def matches_include(self, lid: str, name: str) -> bool: + text = self._text(lid, name) + return any(pattern.search(text) for pattern in self.include) + + def matches_exclude(self, lid: str, name: str) -> bool: + text = self._text(lid, name) + return any(pattern.search(text) for pattern in self.exclude) + + # ------------------------------- # Triangulation helpers # ------------------------------- @@ -921,6 +979,113 @@ def emit_assembly_cpp(lid: str, asm_display_name: str) -> str: return f' TGeoVolumeAssembly *asm_{safe} = new TGeoVolumeAssembly("{name}");' +# ------------------------------- +# CAD clipping helpers +# ------------------------------- + +def make_clip_box_shape(clip_box: ClipBox): + return BRepPrimAPI_MakeBox( + gp_Pnt(clip_box.xmin, clip_box.ymin, clip_box.zmin), + gp_Pnt(clip_box.xmax, clip_box.ymax, clip_box.zmax), + ).Shape() + + +def _compose_trsf(parent_to_world: gp_Trsf, local_to_parent: gp_Trsf) -> gp_Trsf: + return parent_to_world.Multiplied(local_to_parent) + + +def _shape_is_empty(shape) -> bool: + if shape is None: + return True + try: + if shape.IsNull(): + return True + except Exception: + pass + try: + for _ in TopologyExplorer(shape).faces(): + return False + return True + except Exception: + return False + + +def _transformed_bbox(shape, trsf: gp_Trsf) -> Optional[Tuple[float, float, float, float, float, float]]: + box = Bnd_Box() + brepbndlib.Add(shape, box) + try: + xmin, ymin, zmin, xmax, ymax, zmax = box.Get() + except Exception: + return None + + points = [] + for x in (xmin, xmax): + for y in (ymin, ymax): + for z in (zmin, zmax): + p = gp_Pnt(x, y, z) + p.Transform(trsf) + points.append((p.X(), p.Y(), p.Z())) + + return ( + min(p[0] for p in points), + min(p[1] for p in points), + min(p[2] for p in points), + max(p[0] for p in points), + max(p[1] for p in points), + max(p[2] for p in points), + ) + + +def _bbox_outside_clip_box(bbox: Tuple[float, float, float, float, float, float], clip_box: ClipBox) -> bool: + xmin, ymin, zmin, xmax, ymax, zmax = bbox + return ( + xmax < clip_box.xmin or xmin > clip_box.xmax or + ymax < clip_box.ymin or ymin > clip_box.ymax or + zmax < clip_box.zmin or zmin > clip_box.zmax + ) + + +def _bbox_inside_clip_box(bbox: Tuple[float, float, float, float, float, float], clip_box: ClipBox) -> bool: + xmin, ymin, zmin, xmax, ymax, zmax = bbox + return ( + xmin >= clip_box.xmin and xmax <= clip_box.xmax and + ymin >= clip_box.ymin and ymax <= clip_box.ymax and + zmin >= clip_box.zmin and zmax <= clip_box.zmax + ) + + +def _classify_shape_against_clip_box(shape, clip_box: ClipBox, local_to_world: gp_Trsf) -> Optional[str]: + world_bbox = _transformed_bbox(shape, local_to_world) + if world_bbox is None: + return None + if _bbox_outside_clip_box(world_bbox, clip_box): + return "outside" + if _bbox_inside_clip_box(world_bbox, clip_box): + return "inside" + return "overlap" + + +def clip_shape_to_box(shape, clip_box: ClipBox, clip_box_shape, local_to_world: gp_Trsf, lid: str): + clip_state = _classify_shape_against_clip_box(shape, clip_box, local_to_world) + if clip_state is None: + return None + if clip_state == "outside": + return None + if clip_state == "inside": + return shape + + local_clip = BRepBuilderAPI_Transform(clip_box_shape, local_to_world.Inverted(), True).Shape() + common = BRepAlgoAPI_Common(shape, local_clip) + common.Build() + if not common.IsDone(): + raise RuntimeError(f"Failed to clip CAD shape {lid} against --clip-box") + + clipped = common.Shape() + if _shape_is_empty(clipped): + return None + return clipped + + # ------------------------------- # Definition graph extraction # ------------------------------- @@ -939,61 +1104,188 @@ def cpp_var_for_def(lid: str) -> str: return f"asm_{safe}" if lid in assemblies else f"vol_{safe}" -def expand_definition(def_label: TDF_Label, shape_tool, meshparam=None, scale_to_cm: float = 1.0): - def_lid = label_id(def_label) - if def_lid in visited_defs: - return - visited_defs.add(def_lid) +def expand_definition( + def_label: TDF_Label, + shape_tool, + meshparam=None, + scale_to_cm: float = 1.0, + clip_box: Optional[ClipBox] = None, + clip_box_shape=None, + clip_deduplicate: str = "intact", + name_filter: Optional[NameFilter] = None, + include_subtree: bool = False, + world_trsf: Optional[gp_Trsf] = None, + occ_path: str = "r1", +) -> Optional[str]: + clip_enabled = clip_box_shape is not None + if world_trsf is None: + world_trsf = gp_Trsf() + def_lid = label_id(def_label) nm = label_name(def_label) - if nm and def_lid not in def_names: - def_names[def_lid] = nm - elif def_lid not in def_names: - def_names[def_lid] = "" + + subtree_included = include_subtree + if name_filter is not None: + if name_filter.matches_exclude(def_lid, nm): + return None + if name_filter.has_include and name_filter.matches_include(def_lid, nm): + subtree_included = True + + if clip_enabled and clip_box is not None: + try: + shape_for_clip = shape_tool.GetShape(def_label) + except Exception: + shape_for_clip = None + if shape_for_clip is not None: + clip_state = _classify_shape_against_clip_box(shape_for_clip, clip_box, world_trsf) + if clip_state == "outside": + return None + if clip_state == "inside" and clip_deduplicate == "intact": + return expand_definition( + def_label, + shape_tool, + meshparam=meshparam, + scale_to_cm=scale_to_cm, + clip_box=None, + clip_box_shape=None, + clip_deduplicate=clip_deduplicate, + name_filter=name_filter, + include_subtree=subtree_included, + ) + + def_key = f"{def_lid}@{occ_path}" if clip_enabled else def_lid + if not clip_enabled and def_lid in visited_defs: + return def_lid + if not clip_enabled: + visited_defs.add(def_lid) + + if nm and def_key not in def_names: + def_names[def_key] = nm + elif def_key not in def_names: + def_names[def_key] = "" children = TDF_LabelSequence() shape_tool.GetComponents(def_label, children) has_children = children.Length() > 0 if has_children or shape_tool.IsAssembly(def_label): - assemblies.add(def_lid) + assemblies.add(def_key) + kept_children = 0 for i in range(children.Length()): child = children.Value(i + 1) + child_occ_path = f"{occ_path}_{i + 1}" if shape_tool.IsReference(child): referred = TDF_Label() shape_tool.GetReferredShape(child, referred) - child_def_lid = label_id(referred) loc = shape_tool.GetLocation(child) trsf = loc.Transformation() - placements.append((def_lid, child_def_lid, trsf)) - - expand_definition(referred, shape_tool, meshparam=meshparam, scale_to_cm=scale_to_cm) + if clip_enabled: + child_key = expand_definition( + referred, + shape_tool, + meshparam=meshparam, + scale_to_cm=scale_to_cm, + clip_box=clip_box, + clip_box_shape=clip_box_shape, + clip_deduplicate=clip_deduplicate, + name_filter=name_filter, + include_subtree=subtree_included, + world_trsf=_compose_trsf(world_trsf, trsf), + occ_path=child_occ_path, + ) + if child_key is None: + continue + placements.append((def_key, child_key, trsf)) + else: + child_key = expand_definition( + referred, + shape_tool, + meshparam=meshparam, + scale_to_cm=scale_to_cm, + clip_deduplicate=clip_deduplicate, + name_filter=name_filter, + include_subtree=subtree_included, + ) + if child_key is None: + continue + placements.append((def_key, child_key, trsf)) else: - child_def_lid = label_id(child) - placements.append((def_lid, child_def_lid, gp_Trsf())) - expand_definition(child, shape_tool, meshparam=meshparam, scale_to_cm=scale_to_cm) - return + trsf = gp_Trsf() + if clip_enabled: + child_key = expand_definition( + child, + shape_tool, + meshparam=meshparam, + scale_to_cm=scale_to_cm, + clip_box=clip_box, + clip_box_shape=clip_box_shape, + clip_deduplicate=clip_deduplicate, + name_filter=name_filter, + include_subtree=subtree_included, + world_trsf=world_trsf, + occ_path=child_occ_path, + ) + if child_key is None: + continue + placements.append((def_key, child_key, trsf)) + else: + child_key = expand_definition( + child, + shape_tool, + meshparam=meshparam, + scale_to_cm=scale_to_cm, + clip_deduplicate=clip_deduplicate, + name_filter=name_filter, + include_subtree=subtree_included, + ) + if child_key is None: + continue + placements.append((def_key, child_key, trsf)) + kept_children += 1 + + if (clip_enabled or (name_filter is not None and name_filter.has_include)) and kept_children == 0: + assemblies.discard(def_key) + return None + return def_key if shape_tool.IsSimpleShape(def_label): - if def_lid not in logical_volumes: + if name_filter is not None and name_filter.has_include and not subtree_included: + return None + + if def_key not in logical_volumes: shape = shape_tool.GetShape(def_label) # store volume (for density estimation) try: - def_volumes_cm3[def_lid] = volume_cm3_of_shape(shape, scale_to_cm=scale_to_cm) + volume_cm3 = volume_cm3_of_shape(shape, scale_to_cm=scale_to_cm) except Exception: - def_volumes_cm3[def_lid] = 0.0 + volume_cm3 = 0.0 + + if clip_enabled: + shape = clip_shape_to_box(shape, clip_box, clip_box_shape, world_trsf, def_lid) + if shape is None: + return None + + def_volumes_cm3[def_key] = volume_cm3 do_meshing = (meshparam is not None) and meshparam.get("do_meshing", None) is True - logical_volumes[def_lid] = triangulate_CAD_solid(shape, meshparam=meshparam, scale_to_cm=scale_to_cm) if do_meshing else triangulate_asbbox(shape, scale_to_cm=scale_to_cm) - return + logical_volumes[def_key] = triangulate_CAD_solid(shape, meshparam=meshparam, scale_to_cm=scale_to_cm) if do_meshing else triangulate_asbbox(shape, scale_to_cm=scale_to_cm) + return def_key - assemblies.add(def_lid) + assemblies.add(def_key) + return def_key -def extract_graph(step_path: str, meshparam=None, scale_to_cm: float = 1.0): +def extract_graph( + step_path: str, + meshparam=None, + scale_to_cm: float = 1.0, + clip_box: Optional[ClipBox] = None, + clip_deduplicate: str = "intact", + name_filter: Optional[NameFilter] = None, +): global logical_volumes, def_names, def_volumes_cm3, assemblies, placements, top_defs, visited_defs logical_volumes = {} def_names = {} @@ -1004,20 +1296,42 @@ def extract_graph(step_path: str, meshparam=None, scale_to_cm: float = 1.0): visited_defs = set() doc, shape_tool = load_step_with_xcaf(step_path) + clip_box_shape = make_clip_box_shape(clip_box) if clip_box is not None else None roots = TDF_LabelSequence() shape_tool.GetFreeShapes(roots) for i in range(roots.Length()): root = roots.Value(i + 1) + root_occ_path = f"r{i + 1}" if shape_tool.IsReference(root): ref = TDF_Label() shape_tool.GetReferredShape(root, ref) - top_defs.add(label_id(ref)) - expand_definition(ref, shape_tool, meshparam=meshparam, scale_to_cm=scale_to_cm) + top = expand_definition( + ref, + shape_tool, + meshparam=meshparam, + scale_to_cm=scale_to_cm, + clip_box=clip_box, + clip_box_shape=clip_box_shape, + clip_deduplicate=clip_deduplicate, + name_filter=name_filter, + occ_path=root_occ_path, + ) else: - top_defs.add(label_id(root)) - expand_definition(root, shape_tool, meshparam=meshparam, scale_to_cm=scale_to_cm) + top = expand_definition( + root, + shape_tool, + meshparam=meshparam, + scale_to_cm=scale_to_cm, + clip_box=clip_box, + clip_box_shape=clip_box_shape, + clip_deduplicate=clip_deduplicate, + name_filter=name_filter, + occ_path=root_occ_path, + ) + if top is not None: + top_defs.add(top) return doc, shape_tool @@ -1074,6 +1388,9 @@ def emit_root_macro( out_folder: _Path, meshparam=None, step_unit: str = "auto", + clip_box: Optional[ClipBox] = None, + clip_deduplicate: str = "intact", + name_filter: Optional[NameFilter] = None, materials_csv: Optional[str] = None, bom_mass_unit: str = "kg", g4_nist_json: Optional[str] = None, @@ -1087,7 +1404,21 @@ def emit_root_macro( scale_to_cm = step_unit_scale_to_cm(step_unit) print(f"Using overridden STEP length unit: {step_unit} (scale to cm = {scale_to_cm})") - extract_graph(step_path, meshparam=meshparam, scale_to_cm=scale_to_cm) + if clip_box is not None: + print(f"Clipping CAD geometry to STEP-coordinate bounding box: {clip_box.as_tuple()}") + print(f"Clip deduplication mode: {clip_deduplicate}") + + if name_filter is not None and name_filter.active: + print(f"CAD name filters: {len(name_filter.include)} include regex(es), {len(name_filter.exclude)} exclude regex(es)") + + extract_graph( + step_path, + meshparam=meshparam, + scale_to_cm=scale_to_cm, + clip_box=clip_box, + clip_deduplicate=clip_deduplicate, + name_filter=name_filter, + ) out_folder = out_folder.expanduser().resolve() out_folder.mkdir(parents=True, exist_ok=True) @@ -1278,6 +1609,11 @@ def main(): ap.add_argument("--print-tree", action="store_true", help="Just prints the geometry tree") ap.add_argument("--mesh-prec", default=0.1, help="meshing precision. lower --> slower") ap.add_argument("--step-unit", default="auto", choices=["auto", "mm", "cm", "m", "in", "ft"], help="STEP length unit override (default: auto-detect); TGeo expects cm") + ap.add_argument("--clip-box", nargs=6, type=float, metavar=("XMIN", "YMIN", "ZMIN", "XMAX", "YMAX", "ZMAX"), default=None, help="Clip CAD geometry to this axis-aligned bounding box before meshing (coordinates in STEP file units, before conversion to cm)") + ap.add_argument("--clip-deduplicate", default="intact", choices=["none", "intact"], help="When clipping, reuse original logical definitions for subtrees fully inside the clip box (default: intact); use 'none' for one volume per surviving occurrence") + ap.add_argument("--include-name", action="append", default=[], help="Only convert CAD labels whose XCAF name or label entry matches this regex; may be repeated. Matching an assembly includes its subtree.") + ap.add_argument("--exclude-name", action="append", default=[], help="Skip CAD labels/subtrees whose XCAF name or label entry matches this regex; may be repeated.") + ap.add_argument("--name-filter-case-sensitive", action="store_true", help="Make --include-name/--exclude-name matching case-sensitive (default: case-insensitive)") # NEW: BOM / material support ap.add_argument("--materials-csv", default=None, help="BOM CSV file providing material + mass per part (optional)") @@ -1304,6 +1640,24 @@ def main(): if args.out_path is not None: out_folder = _Path(args.out_path) + clip_box = None + if args.clip_box is not None: + try: + clip_box = ClipBox.from_values(args.clip_box) + except ValueError as exc: + ap.error(str(exc)) + + name_filter = None + if args.include_name or args.exclude_name: + try: + name_filter = NameFilter.from_patterns( + args.include_name, + args.exclude_name, + case_sensitive=args.name_filter_case_sensitive, + ) + except re.error as exc: + ap.error(f"Invalid CAD name filter regex: {exc}") + meshparam = {"do_meshing": args.mesh, "lin_defl": args.mesh_prec, "ang_defl": args.mesh_prec} @@ -1325,6 +1679,9 @@ def main(): out_folder, meshparam=meshparam, step_unit=args.step_unit, + clip_box=clip_box, + clip_deduplicate=args.clip_deduplicate, + name_filter=name_filter, materials_csv=args.materials_csv, bom_mass_unit=args.bom_mass_unit, g4_nist_json=args.g4_nist_json, diff --git a/scripts/geometry/README.md b/scripts/geometry/README.md index 8c02cadef51a4..a59da13c12716 100644 --- a/scripts/geometry/README.md +++ b/scripts/geometry/README.md @@ -67,6 +67,60 @@ The generated macro is the file referenced from the external-geometry JSON examp The macro and its facet binaries should stay together, because the macro loads the facet files relative to its own location. +## Restrict the converted region with a clip box + +Large CAD assemblies often contain far more than the region of interest. The `--clip-box` +option restricts the conversion to an axis-aligned bounding box, so only the geometry inside +(or overlapping) that box is written out: + +```bash +python PATH_TO_ALICEO2_SOURCES/scripts/geometry/O2_CADtoTGeo.py \ + my_detector.step \ + --output-folder cad_out/mydet \ + -o geom.C \ + --mesh \ + --clip-box XMIN YMIN ZMIN XMAX YMAX ZMAX +``` + +Notes on the coordinates: + +- The six values are `xmin ymin zmin xmax ymax zmax` and must satisfy `xmin < xmax`, + `ymin < ymax`, and `zmin < zmax`. +- Coordinates are given in the **STEP file units** (before conversion to cm), and are applied + in the global/world coordinate system of the assembly. + +Each solid is classified against the box before meshing: + +- Solids fully outside the box are dropped. +- Solids fully inside the box are kept unchanged. +- Solids straddling the box boundary are cut against it (a boolean intersection), so only the + part inside the box is meshed. + +Assemblies that end up with no surviving children are removed from the output tree. + +### Deduplication mode + +The `--clip-deduplicate` option controls how subtrees that fall entirely inside the box are +emitted: + +- `intact` (default): subtrees fully inside the box reuse their original shared logical + definitions, keeping the output compact. +- `none`: every surviving occurrence becomes its own volume, which is useful when you need a + flat, per-instance representation. + +```bash +python PATH_TO_ALICEO2_SOURCES/scripts/geometry/O2_CADtoTGeo.py \ + my_detector.step \ + --output-folder cad_out/mydet \ + -o geom.C \ + --mesh \ + --clip-box -50 -50 -20 50 50 20 \ + --clip-deduplicate none +``` + +Clipping can be combined with the name-based selection options (`--include-name` / +`--exclude-name`) to further narrow down which parts are converted. + ## Optional material mapping The converter can use a BOM CSV to assign materials and, when part masses and CAD volumes are