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/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 4fb2d1ec610d4..a59da13c12716 100644 --- a/scripts/geometry/README.md +++ b/scripts/geometry/README.md @@ -1,27 +1,283 @@ -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. + +For standalone studies outside the ALICE software stack, a conda environment with +`pythonocc-core` can also be used: -# -) Create an OCC environment (for OpenCacade) +```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. + +## 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 +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