From 1f38c8bca94c8c97232d165f6b66cd8a9b6380f7 Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Sun, 8 Feb 2026 19:17:34 +0100 Subject: [PATCH 1/2] Support to plug-and-play external (CAD) geometry This commit adds support for - conversion from CAD geometries (STEP) file to (meshed) TGeo geometry via the tool O2_CADtoTGeo.py - ability to setup and include simulations modules from external geometries via the ExternalModule mechanism - ability to navigate complex TGeoTessellated solids efficiently (will stay here until fully integrated in the official ROOT repo) It includes all the (basic) functionality to pick a detector layout from a CAD design and actually simulate it with o2-sim. Instructions are provided in a dedicated documentation markdown. The commit is related to epic https://its.cern.ch/jira/browse/O2-6616 --- Detectors/Base/CMakeLists.txt | 7 +- .../include/DetectorsBase/O2Tessellated.h | 142 ++ .../include/DetectorsBase/TGeoGeometryUtils.h | 38 + Detectors/Base/src/DetectorsBaseLinkDef.h | 2 + Detectors/Base/src/O2Tessellated.cxx | 1487 +++++++++++++++++ Detectors/Base/src/TGeoGeometryUtils.cxx | 142 ++ Detectors/Base/src/bvh2_extra_kernels.h | 66 + Detectors/Base/src/bvh2_third_party.h | 36 + Detectors/Passive/CMakeLists.txt | 2 + .../include/DetectorsPassive/ExternalModule.h | 64 + Detectors/Passive/src/ExternalModule.cxx | 162 ++ Steer/include/Steer/O2MCApplicationBase.h | 2 + Steer/src/O2MCApplication.cxx | 60 + macro/build_geometry.C | 14 + scripts/geometry/O2_CADtoTGeo.py | 602 +++++++ scripts/geometry/README.md | 27 + scripts/geometry/simulating_CAD_modules.md | 72 + 17 files changed, 2924 insertions(+), 1 deletion(-) create mode 100644 Detectors/Base/include/DetectorsBase/O2Tessellated.h create mode 100644 Detectors/Base/include/DetectorsBase/TGeoGeometryUtils.h create mode 100644 Detectors/Base/src/O2Tessellated.cxx create mode 100644 Detectors/Base/src/TGeoGeometryUtils.cxx create mode 100644 Detectors/Base/src/bvh2_extra_kernels.h create mode 100644 Detectors/Base/src/bvh2_third_party.h create mode 100644 Detectors/Passive/include/DetectorsPassive/ExternalModule.h create mode 100644 Detectors/Passive/src/ExternalModule.cxx create mode 100644 scripts/geometry/O2_CADtoTGeo.py create mode 100644 scripts/geometry/README.md create mode 100644 scripts/geometry/simulating_CAD_modules.md diff --git a/Detectors/Base/CMakeLists.txt b/Detectors/Base/CMakeLists.txt index 30ab4c4fe8a40..83a9193274e4f 100644 --- a/Detectors/Base/CMakeLists.txt +++ b/Detectors/Base/CMakeLists.txt @@ -29,6 +29,8 @@ o2_add_library(DetectorsBase src/Stack.cxx src/VMCSeederService.cxx src/GlobalParams.cxx + src/O2Tessellated.cxx + src/TGeoGeometryUtils.cxx PUBLIC_LINK_LIBRARIES FairRoot::Base O2::CommonUtils O2::DetectorsCommonDataFormats @@ -46,6 +48,7 @@ o2_add_library(DetectorsBase O2::GPUDataTypes MC::VMC TBB::tbb + ROOT::Gdml ) o2_target_root_dictionary(DetectorsBase @@ -62,7 +65,9 @@ o2_target_root_dictionary(DetectorsBase include/DetectorsBase/Aligner.h include/DetectorsBase/Stack.h include/DetectorsBase/SimFieldUtils.h - include/DetectorsBase/GlobalParams.h) + include/DetectorsBase/GlobalParams.h + include/DetectorsBase/O2Tessellated.h + ) if(BUILD_SIMULATION) if (NOT APPLE) diff --git a/Detectors/Base/include/DetectorsBase/O2Tessellated.h b/Detectors/Base/include/DetectorsBase/O2Tessellated.h new file mode 100644 index 0000000000000..0a1cee8b3e01f --- /dev/null +++ b/Detectors/Base/include/DetectorsBase/O2Tessellated.h @@ -0,0 +1,142 @@ +// 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. + +#ifndef ALICEO2_BASE_O2TESSELLATED_ +#define ALICEO2_BASE_O2TESSELLATED_ + +#include "TGeoShape.h" +#include "TGeoBBox.h" +#include "TGeoVector3.h" +#include "TGeoTypedefs.h" +#include "TGeoTessellated.h" + +namespace o2 +{ +namespace base +{ + +class O2Tessellated : public TGeoBBox +{ + + public: + using Vertex_t = Tessellated::Vertex_t; + + private: + int fNfacets = 0; // Number of facets + int fNvert = 0; // Number of vertices + int fNseg = 0; // Number of segments + bool fDefined = false; //! Shape fully defined + bool fClosedBody = false; // The faces are making a closed body + + // for now separate vectors but might be better to group per face + std::vector fVertices; // List of vertices + std::vector fFacets; // List of facets + std::vector fOutwardNormals; // Vector of outward-facing normals (to be streamed !) + + std::multimap fVerticesMap; //! Temporary map used to deduplicate vertices + bool fIsClosed = false; //! to know if shape still needs closure/initialization + void* fBVH = nullptr; //! BVH acceleration structure for safety and navigation + + O2Tessellated(const O2Tessellated&) = delete; + O2Tessellated& operator=(const O2Tessellated&) = delete; + + // bvh helper functions + void BuildBVH(); + void CalculateNormals(); + + public: + // constructors + O2Tessellated() {} + O2Tessellated(const char* name, int nfacets = 0); + O2Tessellated(const char* name, const std::vector& vertices); + // from a TGeoTessellated + O2Tessellated(TGeoTessellated const&, bool check = false); + + // destructor + ~O2Tessellated() override {} + + void ComputeBBox() override; + void CloseShape(bool check = true, bool fixFlipped = true, bool verbose = true); + + bool AddFacet(const Vertex_t& pt0, const Vertex_t& pt1, const Vertex_t& pt2); + bool AddFacet(const Vertex_t& pt0, const Vertex_t& pt1, const Vertex_t& pt2, const Vertex_t& pt3); + bool AddFacet(int i1, int i2, int i3); + bool AddFacet(int i1, int i2, int i3, int i4); + int AddVertex(const Vertex_t& vert); + + bool FacetCheck(int ifacet) const; + Vertex_t FacetComputeNormal(int ifacet, bool& degenerated) const; + + int GetNfacets() const { return fFacets.size(); } + int GetNsegments() const { return fNseg; } + int GetNvertices() const { return fNvert; } + bool IsClosedBody() const { return fClosedBody; } + bool IsDefined() const { return fDefined; } + + const TGeoFacet& GetFacet(int i) const { return fFacets[i]; } + const Vertex_t& GetVertex(int i) const { return fVertices[i]; } + + int DistancetoPrimitive(int, int) override { return 99999; } + const TBuffer3D& GetBuffer3D(int reqSections, Bool_t localFrame) const override; + void GetMeshNumbers(int& nvert, int& nsegs, int& npols) const override; + int GetNmeshVertices() const override { return fNvert; } + void InspectShape() const override {} + TBuffer3D* MakeBuffer3D() const override; + void Print(Option_t* option = "") const override; + void SavePrimitive(std::ostream&, Option_t*) override {} + void SetPoints(double* points) const override; + void SetPoints(Float_t* points) const override; + void SetSegsAndPols(TBuffer3D& buff) const override; + void Sizeof3D() const override {} + + /// Resize and center the shape in a box of size maxsize + void ResizeCenter(double maxsize); + + /// Flip all facets + void FlipFacets() + { + for (auto facet : fFacets) + facet.Flip(); + } + + bool CheckClosure(bool fixFlipped = true, bool verbose = true); + + /// Reader from .obj format + static O2Tessellated* ImportFromObjFormat(const char* objfile, bool check = false, bool verbose = false); + + // navigation functions used by TGeoNavigator (attention: only the iact == 3 cases implemented for now) + Double_t DistFromOutside(const Double_t* point, const Double_t* dir, Int_t iact = 1, + Double_t step = TGeoShape::Big(), Double_t* safe = nullptr) const override; + Double_t DistFromInside(const Double_t* point, const Double_t* dir, Int_t iact = 1, Double_t step = TGeoShape::Big(), + Double_t* safe = nullptr) const override; + bool Contains(const Double_t* point) const override; + Double_t Safety(const Double_t* point, Bool_t in = kTRUE) const override; + void ComputeNormal(const Double_t* point, const Double_t* dir, Double_t* norm) const override; + + // these are trivial implementations, just for debugging + Double_t DistFromInside_Loop(const Double_t* point, const Double_t* dir) const; + Double_t DistFromOutside_Loop(const Double_t* point, const Double_t* dir) const; + bool Contains_Loop(const Double_t* point) const; + + Double_t Capacity() const override; + + private: + // a safety kernel used in multiple implementations + template + Double_t SafetyKernel(const Double_t* point, bool in, int* closest_facet_id = nullptr) const; + + ClassDefOverride(O2Tessellated, 1) // tessellated shape class +}; + +} // namespace base +} // namespace o2 + +#endif diff --git a/Detectors/Base/include/DetectorsBase/TGeoGeometryUtils.h b/Detectors/Base/include/DetectorsBase/TGeoGeometryUtils.h new file mode 100644 index 0000000000000..5ec85f1c14702 --- /dev/null +++ b/Detectors/Base/include/DetectorsBase/TGeoGeometryUtils.h @@ -0,0 +1,38 @@ +// 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 TGeoGeometryUtils.h +/// \author Sandro Wenzel (CERN) +/// \brief Collection of utility functions for TGeo + +#ifndef ALICEO2_BASE_TGEOGEOMETRYUTILS_H_ +#define ALICEO2_BASE_TGEOGEOMETRYUTILS_H_ + +class TGeoShape; +class TGeoTessellated; + +namespace o2 +{ +namespace base +{ + +/// A few utility functions to operate on TGeo geometries (transformations, printing, ...) +class TGeoGeometryUtils +{ + public: + ///< Transform any (primitive) TGeoShape to a tessellated representation + static TGeoTessellated* TGeoShapeToTGeoTessellated(TGeoShape const*); +}; + +} // namespace base +} // namespace o2 + +#endif diff --git a/Detectors/Base/src/DetectorsBaseLinkDef.h b/Detectors/Base/src/DetectorsBaseLinkDef.h index bd76e9bfbe2e4..8255c143ebb4a 100644 --- a/Detectors/Base/src/DetectorsBaseLinkDef.h +++ b/Detectors/Base/src/DetectorsBaseLinkDef.h @@ -42,4 +42,6 @@ #pragma link C++ class o2::data::Stack + ; +#pragma link C++ class o2::base::O2Tessellated - ; + #endif diff --git a/Detectors/Base/src/O2Tessellated.cxx b/Detectors/Base/src/O2Tessellated.cxx new file mode 100644 index 0000000000000..1bd209cfeab33 --- /dev/null +++ b/Detectors/Base/src/O2Tessellated.cxx @@ -0,0 +1,1487 @@ +#include +#include + +#include "TGeoManager.h" +#include "TGeoMatrix.h" +#include "TGeoVolume.h" +#include "TVirtualGeoPainter.h" +#include "DetectorsBase/O2Tessellated.h" +#include "TBuffer3D.h" +#include "TBuffer3DTypes.h" +#include "TMath.h" +#include "TBuffer.h" + +#include +#include + +// include the Third-party BVH headers +#include "bvh2_third_party.h" +// some kernels on top of BVH +#include "bvh2_extra_kernels.h" + +#include +#include + +using namespace o2::base; +ClassImp(O2Tessellated); + +using Vertex_t = Tessellated::Vertex_t; + +//////////////////////////////////////////////////////////////////////////////// +/// Compact consecutive equal vertices + +int TGeoFacet::CompactFacet(Vertex_t* vert, int nvertices) +{ + // Compact the common vertices and return new facet + if (nvertices < 2) + return nvertices; + int nvert = nvertices; + int i = 0; + while (i < nvert) { + if (vert[(i + 1) % nvert] == vert[i]) { + // shift last vertices left by one element + for (int j = i + 2; j < nvert; ++j) + vert[j - 1] = vert[j]; + nvert--; + } + i++; + } + return nvert; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Check if a connected neighbour facet has compatible normal + +bool TGeoFacet::IsNeighbour(const TGeoFacet& other, bool& flip) const +{ + + // Find a connecting segment + bool neighbour = false; + int line1[2], line2[2]; + int npoints = 0; + for (int i = 0; i < fNvert; ++i) { + auto ivert = fIvert[i]; + // Check if the other facet has the same vertex + for (int j = 0; j < other.GetNvert(); ++j) { + if (ivert == other[j]) { + line1[npoints] = i; + line2[npoints] = j; + if (++npoints == 2) { + neighbour = true; + bool order1 = line1[1] == line1[0] + 1; + bool order2 = line2[1] == (line2[0] + 1) % other.GetNvert(); + flip = (order1 == order2); + return neighbour; + } + } + } + } + return neighbour; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Constructor. In case nfacets is zero, it is user's responsibility to +/// call CloseShape once all faces are defined. + +O2Tessellated::O2Tessellated(const char* name, int nfacets) : TGeoBBox(name, 0, 0, 0) +{ + fNfacets = nfacets; + if (nfacets) + fFacets.reserve(nfacets); +} + +//////////////////////////////////////////////////////////////////////////////// +/// Constructor providing directly the array of vertices. Facets have to be added +/// providing vertex indices rather than coordinates. + +O2Tessellated::O2Tessellated(const char* name, const std::vector& vertices) : TGeoBBox(name, 0, 0, 0) +{ + fVertices = vertices; + fNvert = fVertices.size(); +} + +//////////////////////////////////////////////////////////////////////////////// +/// Construct from TGeoTessellated + +O2Tessellated::O2Tessellated(TGeoTessellated const& tsl, bool check) : TGeoBBox(tsl.GetName(), 0, 0, 0) +{ + fNfacets = tsl.GetNfacets(); + fNvert = tsl.GetNvertices(); + fNseg = tsl.GetNsegments(); + + // copy facet and vertex done + fVertices.reserve(fNvert); + fFacets.reserve(fNfacets); + for (int i = 0; i < fNfacets; ++i) { + fFacets.push_back(tsl.GetFacet(i)); + } + for (int i = 0; i < fNvert; ++i) { + fVertices.push_back(tsl.GetVertex(i)); + } + // finish remaining structures + CloseShape(check); +} + +//////////////////////////////////////////////////////////////////////////////// +/// Add a vertex checking for duplicates, returning the vertex index + +int O2Tessellated::AddVertex(Vertex_t const& vert) +{ + constexpr double tolerance = 1.e-10; + auto vertexHash = [&](Vertex_t const& vertex) { + // Compute hash for the vertex + long hash = 0; + // helper function to generate hash from integer numbers + auto hash_combine = [](long seed, const long value) { + return seed ^ (std::hash{}(value) + 0x9e3779b9 + (seed << 6) + (seed >> 2)); + }; + for (int i = 0; i < 3; i++) { + // use tolerance to generate int with the desired precision from a real number for hashing + hash = hash_combine(hash, std::roundl(vertex[i] / tolerance)); + } + return hash; + }; + + auto hash = vertexHash(vert); + bool isAdded = false; + int ivert = -1; + // Get the compatible vertices + auto range = fVerticesMap.equal_range(hash); + for (auto it = range.first; it != range.second; ++it) { + ivert = it->second; + if (fVertices[ivert] == vert) { + isAdded = true; + break; + } + } + if (!isAdded) { + ivert = fVertices.size(); + fVertices.push_back(vert); + fVerticesMap.insert(std::make_pair(hash, ivert)); + } + return ivert; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Adding a triangular facet from vertex positions in absolute coordinates + +bool O2Tessellated::AddFacet(const Vertex_t& pt0, const Vertex_t& pt1, const Vertex_t& pt2) +{ + if (fDefined) { + Error("AddFacet", "Shape %s already fully defined. Not adding", GetName()); + return false; + } + + Vertex_t vert[3]; + vert[0] = pt0; + vert[1] = pt1; + vert[2] = pt2; + int nvert = TGeoFacet::CompactFacet(vert, 3); + if (nvert < 3) { + Error("AddFacet", "Triangular facet at index %d degenerated. Not adding.", GetNfacets()); + return false; + } + int ind[3]; + for (auto i = 0; i < 3; ++i) + ind[i] = AddVertex(vert[i]); + fNseg += 3; + fFacets.emplace_back(ind[0], ind[1], ind[2]); + + return true; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Adding a triangular facet from indices of vertices + +bool O2Tessellated::AddFacet(int i0, int i1, int i2) +{ + if (fDefined) { + Error("AddFacet", "Shape %s already fully defined. Not adding", GetName()); + return false; + } + if (fVertices.empty()) { + Error("AddFacet", "Shape %s Cannot add facets by indices without vertices. Not adding", GetName()); + return false; + } + + fNseg += 3; + fFacets.emplace_back(i0, i1, i2); + return true; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Adding a quadrilateral facet from vertex positions in absolute coordinates + +bool O2Tessellated::AddFacet(const Vertex_t& pt0, const Vertex_t& pt1, const Vertex_t& pt2, const Vertex_t& pt3) +{ + if (fDefined) { + Error("AddFacet", "Shape %s already fully defined. Not adding", GetName()); + return false; + } + Vertex_t vert[4]; + vert[0] = pt0; + vert[1] = pt1; + vert[2] = pt2; + vert[3] = pt3; + int nvert = TGeoFacet::CompactFacet(vert, 4); + if (nvert < 3) { + Error("AddFacet", "Quadrilateral facet at index %d degenerated. Not adding.", GetNfacets()); + return false; + } + + int ind[4]; + for (auto i = 0; i < nvert; ++i) + ind[i] = AddVertex(vert[i]); + fNseg += nvert; + if (nvert == 3) + fFacets.emplace_back(ind[0], ind[1], ind[2]); + else + fFacets.emplace_back(ind[0], ind[1], ind[2], ind[3]); + + if (fNfacets > 0 && GetNfacets() == fNfacets) + CloseShape(false); + return true; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Adding a quadrilateral facet from indices of vertices + +bool O2Tessellated::AddFacet(int i0, int i1, int i2, int i3) +{ + if (fDefined) { + Error("AddFacet", "Shape %s already fully defined. Not adding", GetName()); + return false; + } + if (fVertices.empty()) { + Error("AddFacet", "Shape %s Cannot add facets by indices without vertices. Not adding", GetName()); + return false; + } + + fNseg += 4; + fFacets.emplace_back(i0, i1, i2, i3); + return true; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Compute normal for a given facet + +Vertex_t O2Tessellated::FacetComputeNormal(int ifacet, bool& degenerated) const +{ + // Compute normal using non-zero segments + constexpr double kTolerance = 1.e-20; + auto const& facet = fFacets[ifacet]; + int nvert = facet.GetNvert(); + degenerated = true; + Vertex_t normal; + for (int i = 0; i < nvert - 1; ++i) { + Vertex_t e1 = fVertices[facet[i + 1]] - fVertices[facet[i]]; + if (e1.Mag2() < kTolerance) + continue; + for (int j = i + 1; j < nvert; ++j) { + Vertex_t e2 = fVertices[facet[(j + 1) % nvert]] - fVertices[facet[j]]; + if (e2.Mag2() < kTolerance) + continue; + normal = Vertex_t::Cross(e1, e2); + // e1 and e2 may be colinear + if (normal.Mag2() < kTolerance) + continue; + normal.Normalize(); + degenerated = false; + break; + } + if (!degenerated) + break; + } + return normal; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Check validity of facet + +bool O2Tessellated::FacetCheck(int ifacet) const +{ + constexpr double kTolerance = 1.e-10; + auto const& facet = fFacets[ifacet]; + int nvert = facet.GetNvert(); + bool degenerated = true; + FacetComputeNormal(ifacet, degenerated); + if (degenerated) { + std::cout << "Facet: " << ifacet << " is degenerated\n"; + return false; + } + + // Compute surface area + double surfaceArea = 0.; + for (int i = 1; i < nvert - 1; ++i) { + Vertex_t e1 = fVertices[facet[i]] - fVertices[facet[0]]; + Vertex_t e2 = fVertices[facet[i + 1]] - fVertices[facet[0]]; + surfaceArea += 0.5 * Vertex_t::Cross(e1, e2).Mag(); + } + if (surfaceArea < kTolerance) { + std::cout << "Facet: " << ifacet << " has zero surface area\n"; + return false; + } + + return true; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Close the shape: calculate bounding box and compact vertices + +void O2Tessellated::CloseShape(bool check, bool fixFlipped, bool verbose) +{ + if (fIsClosed && fBVH) { + return; + } + // Compute bounding box + fDefined = true; + fNvert = fVertices.size(); + fNfacets = fFacets.size(); + ComputeBBox(); + + BuildBVH(); + if (fOutwardNormals.size() == 0) { + CalculateNormals(); + } else { + // short check if the normal container is of correct size + if (!fOutwardNormals.size() == fFacets.size()) { + std::cerr << "Inconsistency in normal container"; + } + } + fIsClosed = true; + + // Cleanup the vertex map + std::multimap().swap(fVerticesMap); + + if (fVertices.size() > 0) { + if (!check) + return; + + // Check facets + for (auto i = 0; i < fNfacets; ++i) + FacetCheck(i); + + fClosedBody = CheckClosure(fixFlipped, verbose); + } +} + +//////////////////////////////////////////////////////////////////////////////// +/// Check closure of the solid and check/fix flipped normals + +bool O2Tessellated::CheckClosure(bool fixFlipped, bool verbose) +{ + int* nn = new int[fNfacets]; + bool* flipped = new bool[fNfacets]; + bool hasorphans = false; + bool hasflipped = false; + for (int i = 0; i < fNfacets; ++i) { + nn[i] = 0; + flipped[i] = false; + } + + for (int icrt = 0; icrt < fNfacets; ++icrt) { + // all neighbours checked? + if (nn[icrt] >= fFacets[icrt].GetNvert()) + continue; + for (int i = icrt + 1; i < fNfacets; ++i) { + bool isneighbour = fFacets[icrt].IsNeighbour(fFacets[i], flipped[i]); + if (isneighbour) { + if (flipped[icrt]) + flipped[i] = !flipped[i]; + if (flipped[i]) + hasflipped = true; + nn[icrt]++; + nn[i]++; + if (nn[icrt] == fFacets[icrt].GetNvert()) + break; + } + } + if (nn[icrt] < fFacets[icrt].GetNvert()) + hasorphans = true; + } + + if (hasorphans && verbose) { + Error("Check", "Tessellated solid %s has following not fully connected facets:", GetName()); + for (int icrt = 0; icrt < fNfacets; ++icrt) { + if (nn[icrt] < fFacets[icrt].GetNvert()) + std::cout << icrt << " (" << fFacets[icrt].GetNvert() << " edges, " << nn[icrt] << " neighbours)\n"; + } + } + fClosedBody = !hasorphans; + int nfixed = 0; + if (hasflipped) { + if (verbose) + Warning("Check", "Tessellated solid %s has following facets with flipped normals:", GetName()); + for (int icrt = 0; icrt < fNfacets; ++icrt) { + if (flipped[icrt]) { + if (verbose) + std::cout << icrt << "\n"; + if (fixFlipped) { + fFacets[icrt].Flip(); + nfixed++; + } + } + } + if (nfixed && verbose) + Info("Check", "Automatically flipped %d facets to match first defined facet", nfixed); + } + delete[] nn; + delete[] flipped; + + return !hasorphans; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Compute bounding box + +void O2Tessellated::ComputeBBox() +{ + const double kBig = TGeoShape::Big(); + double vmin[3] = {kBig, kBig, kBig}; + double vmax[3] = {-kBig, -kBig, -kBig}; + for (const auto& facet : fFacets) { + for (int i = 0; i < facet.GetNvert(); ++i) { + for (int j = 0; j < 3; ++j) { + vmin[j] = TMath::Min(vmin[j], fVertices[facet[i]].operator[](j)); + vmax[j] = TMath::Max(vmax[j], fVertices[facet[i]].operator[](j)); + } + } + } + fDX = 0.5 * (vmax[0] - vmin[0]); + fDY = 0.5 * (vmax[1] - vmin[1]); + fDZ = 0.5 * (vmax[2] - vmin[2]); + for (int i = 0; i < 3; ++i) + fOrigin[i] = 0.5 * (vmax[i] + vmin[i]); +} + +//////////////////////////////////////////////////////////////////////////////// +/// Returns numbers of vertices, segments and polygons composing the shape mesh. + +void O2Tessellated::GetMeshNumbers(int& nvert, int& nsegs, int& npols) const +{ + nvert = fNvert; + nsegs = fNseg; + npols = GetNfacets(); +} + +//////////////////////////////////////////////////////////////////////////////// +/// Creates a TBuffer3D describing *this* shape. +/// Coordinates are in local reference frame. + +TBuffer3D* O2Tessellated::MakeBuffer3D() const +{ + const int nvert = fNvert; + const int nsegs = fNseg; + const int npols = GetNfacets(); + auto buff = new TBuffer3D(TBuffer3DTypes::kGeneric, nvert, 3 * nvert, nsegs, 3 * nsegs, npols, 6 * npols); + if (buff) { + SetPoints(buff->fPnts); + SetSegsAndPols(*buff); + } + return buff; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Prints basic info + +void O2Tessellated::Print(Option_t*) const +{ + std::cout << "=== Tessellated shape " << GetName() << " having " << GetNvertices() << " vertices and " + << GetNfacets() << " facets\n"; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Fills TBuffer3D structure for segments and polygons. + +void O2Tessellated::SetSegsAndPols(TBuffer3D& buff) const +{ + const int c = GetBasicColor(); + int* segs = buff.fSegs; + int* pols = buff.fPols; + + int indseg = 0; // segment internal data index + int indpol = 0; // polygon internal data index + int sind = 0; // segment index + for (const auto& facet : fFacets) { + auto nvert = facet.GetNvert(); + pols[indpol++] = c; + pols[indpol++] = nvert; + for (auto j = 0; j < nvert; ++j) { + int k = (j + 1) % nvert; + // segment made by next consecutive points + segs[indseg++] = c; + segs[indseg++] = facet[j]; + segs[indseg++] = facet[k]; + // add segment to current polygon and increment segment index + pols[indpol + nvert - j - 1] = sind++; + } + indpol += nvert; + } +} + +//////////////////////////////////////////////////////////////////////////////// +/// Fill tessellated points to an array. + +void O2Tessellated::SetPoints(double* points) const +{ + int ind = 0; + for (const auto& vertex : fVertices) { + vertex.CopyTo(&points[ind]); + ind += 3; + } +} + +//////////////////////////////////////////////////////////////////////////////// +/// Fill tessellated points in float. + +void O2Tessellated::SetPoints(Float_t* points) const +{ + int ind = 0; + for (const auto& vertex : fVertices) { + points[ind++] = vertex.x(); + points[ind++] = vertex.y(); + points[ind++] = vertex.z(); + } +} + +//////////////////////////////////////////////////////////////////////////////// +/// Resize the shape by scaling vertices within maxsize and center to origin + +void O2Tessellated::ResizeCenter(double maxsize) +{ + using Vector3_t = Vertex_t; + + if (!fDefined) { + Error("ResizeCenter", "Not all faces are defined"); + return; + } + Vector3_t origin(fOrigin[0], fOrigin[1], fOrigin[2]); + double maxedge = TMath::Max(TMath::Max(fDX, fDY), fDZ); + double scale = maxsize / maxedge; + for (size_t i = 0; i < fVertices.size(); ++i) { + fVertices[i] = scale * (fVertices[i] - origin); + } + fOrigin[0] = fOrigin[1] = fOrigin[2] = 0; + fDX *= scale; + fDY *= scale; + fDZ *= scale; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Fills a static 3D buffer and returns a reference. + +const TBuffer3D& O2Tessellated::GetBuffer3D(int reqSections, Bool_t localFrame) const +{ + static TBuffer3D buffer(TBuffer3DTypes::kGeneric); + + FillBuffer3D(buffer, reqSections, localFrame); + + const int nvert = fNvert; + const int nsegs = fNseg; + const int npols = GetNfacets(); + + if (reqSections & TBuffer3D::kRawSizes) { + if (buffer.SetRawSizes(nvert, 3 * nvert, nsegs, 3 * nsegs, npols, 6 * npols)) { + buffer.SetSectionsValid(TBuffer3D::kRawSizes); + } + } + if ((reqSections & TBuffer3D::kRaw) && buffer.SectionsValid(TBuffer3D::kRawSizes)) { + SetPoints(buffer.fPnts); + if (!buffer.fLocalFrame) { + TransformPoints(buffer.fPnts, buffer.NbPnts()); + } + + SetSegsAndPols(buffer); + buffer.SetSectionsValid(TBuffer3D::kRaw); + } + + return buffer; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Reads a single tessellated solid from an .obj file. + +O2Tessellated* O2Tessellated::ImportFromObjFormat(const char* objfile, bool check, bool verbose) +{ + using std::vector, std::string, std::ifstream, std::stringstream, std::endl; + + vector vertices; + vector sfacets; + + struct FacetInd_t { + int i0 = -1; + int i1 = -1; + int i2 = -1; + int i3 = -1; + int nvert = 0; + FacetInd_t(int a, int b, int c) + { + i0 = a; + i1 = b; + i2 = c; + nvert = 3; + }; + FacetInd_t(int a, int b, int c, int d) + { + i0 = a; + i1 = b; + i2 = c; + i3 = d; + nvert = 4; + }; + }; + + vector facets; + // List of geometric vertices, with (x, y, z [,w]) coordinates, w is optional and defaults to 1.0. + // struct vtx_t { double x = 0; double y = 0; double z = 0; double w = 1; }; + + // Texture coordinates in u, [,v ,w]) coordinates, these will vary between 0 and 1. v, w are optional and default to + // 0. + // struct tex_t { double u; double v; double w; }; + + // List of vertex normals in (x,y,z) form; normals might not be unit vectors. + // struct vn_t { double x; double y; double z; }; + + // Parameter space vertices in ( u [,v] [,w] ) form; free form geometry statement + // struct vp_t { double u; double v; double w; }; + + // Faces are defined using lists of vertex, texture and normal indices which start at 1. + // Polygons such as quadrilaterals can be defined by using more than three vertex/texture/normal indices. + // f v1//vn1 v2//vn2 v3//vn3 ... + + // Records starting with the letter "l" specify the order of the vertices which build a polyline. + // l v1 v2 v3 v4 v5 v6 ... + + string line; + int ind[4] = {0}; + ifstream file(objfile); + if (!file.is_open()) { + ::Error("O2Tessellated::ImportFromObjFormat", "Unable to open %s", objfile); + return nullptr; + } + + while (getline(file, line)) { + stringstream ss(line); + string tag; + + // We ignore everything which is not a vertex or a face + if (line.rfind('v', 0) == 0 && line.rfind("vt", 0) != 0 && line.rfind("vn", 0) != 0 && line.rfind("vn", 0) != 0) { + // Decode the vertex + double pos[4] = {0, 0, 0, 1}; + ss >> tag >> pos[0] >> pos[1] >> pos[2] >> pos[3]; + vertices.emplace_back(pos[0] * pos[3], pos[1] * pos[3], pos[2] * pos[3]); + } + + else if (line.rfind('f', 0) == 0) { + // Decode the face + ss >> tag; + string word; + sfacets.clear(); + while (ss >> word) + sfacets.push_back(word); + if (sfacets.size() > 4 || sfacets.size() < 3) { + ::Error("O2Tessellated::ImportFromObjFormat", "Detected face having unsupported %zu vertices", + sfacets.size()); + return nullptr; + } + int nvert = 0; + for (auto& sword : sfacets) { + stringstream ssword(sword); + string token; + getline(ssword, token, '/'); // just need the vertex index, which is the first token + // Convert string token to integer + + ind[nvert++] = stoi(token) - 1; + if (ind[nvert - 1] < 0) { + ::Error("O2Tessellated::ImportFromObjFormat", "Unsupported relative vertex index definition in %s", + objfile); + return nullptr; + } + } + if (nvert == 3) + facets.emplace_back(ind[0], ind[1], ind[2]); + else + facets.emplace_back(ind[0], ind[1], ind[2], ind[3]); + } + } + + int nvertices = (int)vertices.size(); + int nfacets = (int)facets.size(); + if (nfacets < 3) { + ::Error("O2Tessellated::ImportFromObjFormat", "Not enough faces detected in %s", objfile); + return nullptr; + } + + string sobjfile(objfile); + if (verbose) + std::cout << "Read " << nvertices << " vertices and " << nfacets << " facets from " << sobjfile << endl; + + auto tsl = new O2Tessellated(sobjfile.erase(sobjfile.find_last_of('.')).c_str(), vertices); + + for (int i = 0; i < nfacets; ++i) { + auto facet = facets[i]; + if (facet.nvert == 3) + tsl->AddFacet(facet.i0, facet.i1, facet.i2); + else + tsl->AddFacet(facet.i0, facet.i1, facet.i2, facet.i3); + } + tsl->CloseShape(check, true, verbose); + tsl->Print(); + return tsl; +} + +// implementation of some geometry helper functions in anonymous namespace +namespace +{ + +using Vertex_t = Tessellated::Vertex_t; +// The classic Moeller-Trumbore ray triangle-intersection kernel: +// - Compute triangle edges e1, e2 +// - Compute determinant det +// - Reject parallel rays +// - Compute barycentric coordinates u, v +// - Compute ray parameter t +double rayTriangle(const Vertex_t& orig, const Vertex_t& dir, const Vertex_t& v0, const Vertex_t& v1, + const Vertex_t& v2, double rayEPS = 1e-8) +{ + constexpr double EPS = 1e-8; + const double INF = std::numeric_limits::infinity(); + Vertex_t e1{v1[0] - v0[0], v1[1] - v0[1], v1[2] - v0[2]}; + Vertex_t e2{v2[0] - v0[0], v2[1] - v0[1], v2[2] - v0[2]}; + auto p = Vertex_t::Cross(dir, e2); + auto det = e1.Dot(p); + if (std::abs(det) <= EPS) { + return INF; + } + + Vertex_t tvec{orig[0] - v0[0], orig[1] - v0[1], orig[2] - v0[2]}; + auto invDet = 1.0 / det; + auto u = tvec.Dot(p) * invDet; + if (u < 0.0 || u > 1.0) { + return INF; + } + auto q = Vertex_t::Cross(tvec, e1); + auto v = dir.Dot(q) * invDet; + if (v < 0.0 || u + v > 1.0) { + return INF; + } + auto t = e2.Dot(q) * invDet; + return (t > rayEPS) ? t : INF; +} + +template +struct Vec3f { + T x, y, z; +}; + +template +inline Vec3f operator-(const Vec3f& a, const Vec3f& b) +{ + return {a.x - b.x, a.y - b.y, a.z - b.z}; +} + +template +inline Vec3f cross(const Vec3f& a, const Vec3f& b) +{ + return {a.y * b.z - a.z * b.y, a.z * b.x - a.x * b.z, a.x * b.y - a.y * b.x}; +} + +template +inline T dot(const Vec3f& a, const Vec3f& b) +{ + return a.x * b.x + a.y * b.y + a.z * b.z; +} + +// Kernel to get closest/shortest distance between a point and a triangl (a,b,c). +// Performed by default in float since Safety can be approximate. +// Project point onto triangle plane +// If projection lies inside → distance to plane +// Otherwise compute min distance to the three edges +// Return squared distance +template +T pointTriangleDistSq(const Vec3f& p, const Vec3f& a, const Vec3f& b, const Vec3f& c) +{ + // Edges + Vec3f ab = b - a; + Vec3f ac = c - a; + Vec3f ap = p - a; + + auto d1 = dot(ab, ap); + auto d2 = dot(ac, ap); + if (d1 <= T(0.0) && d2 <= T(0.0)) { + return dot(ap, ap); // barycentric (1,0,0) + } + + Vec3f bp = p - b; + auto d3 = dot(ab, bp); + auto d4 = dot(ac, bp); + if (d3 >= T(0.0) && d4 <= d3) { + return dot(bp, bp); // (0,1,0) + } + + T vc = d1 * d4 - d3 * d2; + if (vc <= 0.0f && d1 >= 0.0f && d3 <= 0.0f) { + T v = d1 / (d1 - d3); + Vec3f proj = {a.x + v * ab.x, a.y + v * ab.y, a.z + v * ab.z}; + Vec3f d = p - proj; + return dot(d, d); // edge AB + } + + Vec3f cp = p - c; + T d5 = dot(ab, cp); + T d6 = dot(ac, cp); + if (d6 >= T(0.0f) && d5 <= d6) { + return dot(cp, cp); // (0,0,1) + } + + T vb = d5 * d2 - d1 * d6; + if (vb <= 0.0f && d2 >= 0.0f && d6 <= 0.0f) { + T w = d2 / (d2 - d6); + Vec3f proj = {a.x + w * ac.x, a.y + w * ac.y, a.z + w * ac.z}; + Vec3f d = p - proj; + return dot(d, d); // edge AC + } + + T va = d3 * d6 - d5 * d4; + if (va <= 0.0f && (d4 - d3) >= 0.0f && (d5 - d6) >= 0.0f) { + T w = (d4 - d3) / ((d4 - d3) + (d5 - d6)); + Vec3f proj = {b.x + w * (c.x - b.x), b.y + w * (c.y - b.y), b.z + w * (c.z - b.z)}; + Vec3f d = p - proj; + return dot(d, d); // edge BC + } + + // Inside face region + T denom = T(1.0f) / (va + vb + vc); + T v = vb * denom; + T w = vc * denom; + + Vec3f proj = {a.x + ab.x * v + ac.x * w, a.y + ab.y * v + ac.y * w, a.z + ab.z * v + ac.z * w}; + + Vec3f d = p - proj; + return dot(d, d); +} + +template +inline Vec3f normalize(const Vec3f& v) +{ + T len2 = dot(v, v); + if (len2 == T(0.0f)) { + std::cerr << "Degnerate triangle. Cannot determine normal"; + return {0, 0, 0}; + } + T invLen = T(1.0f) / std::sqrt(len2); + return {v.x * invLen, v.y * invLen, v.z * invLen}; +} + +template +inline Vec3f triangleNormal(const Vec3f& a, const Vec3f& b, const Vec3f& c) +{ + const Vec3f e1 = b - a; + const Vec3f e2 = c - a; + return normalize(cross(e1, e2)); +} + +} // end anonymous namespace + +//////////////////////////////////////////////////////////////////////////////// +/// DistFromOutside + +Double_t O2Tessellated::DistFromOutside(const Double_t* point, const Double_t* dir, Int_t /*iact*/, Double_t stepmax, + Double_t* /*safe*/) const +{ + // use the BVH intersector in combination with leaf ray-triangle testing + double local_step = Big(); // we need this otherwise the lambda get's confused + + using Scalar = float; + using Vec3 = bvh::v2::Vec; + using Node = bvh::v2::Node; + using Bvh = bvh::v2::Bvh; + using Ray = bvh::v2::Ray; + + // let's fetch the bvh + auto mybvh = (Bvh*)fBVH; + if (!mybvh) { + assert(false); + return -1.; + } + + auto truncate_roundup = [](double orig) { + float epsilon = std::numeric_limits::epsilon() * std::fabs(orig); + // Add the bias to x before assigning it to y + return static_cast(orig + epsilon); + }; + + // let's do very quick checks against the top node + const auto topnode_bbox = mybvh->get_root().get_bbox(); + if ((-point[0] + topnode_bbox.min[0]) > stepmax) { + return Big(); + } + if ((-point[1] + topnode_bbox.min[1]) > stepmax) { + return Big(); + } + if ((-point[2] + topnode_bbox.min[2]) > stepmax) { + return Big(); + } + if ((point[0] - topnode_bbox.max[0]) > stepmax) { + return Big(); + } + if ((point[1] - topnode_bbox.max[1]) > stepmax) { + return Big(); + } + if ((point[2] - topnode_bbox.max[2]) > stepmax) { + return Big(); + } + + // the ray used for bvh interaction + Ray ray(Vec3(point[0], point[1], point[2]), // origin + Vec3(dir[0], dir[1], dir[2]), // direction + 0.0f, // minimum distance (could give stepmax ?) + truncate_roundup(local_step)); + + static constexpr bool use_robust_traversal = true; + + Vertex_t dir_v{dir[0], dir[1], dir[2]}; + // Traverse the BVH and apply concrete object intersection in BVH leafs + bvh::v2::GrowingStack stack; + mybvh->intersect(ray, mybvh->get_root().index, stack, [&](size_t begin, size_t end) { + for (size_t prim_id = begin; prim_id < end; ++prim_id) { + auto objectid = mybvh->prim_ids[prim_id]; + const auto& facet = fFacets[objectid]; + const auto& n = fOutwardNormals[objectid]; + + // quick normal test. Coming from outside, the dot product must be negative + if (n.Dot(dir_v) > 0.) { + continue; + } + + auto thisdist = rayTriangle(Vertex_t(point[0], point[1], point[2]), dir_v, + fVertices[facet[0]], fVertices[facet[1]], fVertices[facet[2]], 0.); + + if (thisdist < local_step) { + local_step = thisdist; + } + } + return false; // go on after this + }); + + return local_step; +} + +//////////////////////////////////////////////////////////////////////////////// +/// DistFromOutside + +Double_t O2Tessellated::DistFromInside(const Double_t* point, const Double_t* dir, Int_t /*iact*/, Double_t /*stepmax*/, + Double_t* /*safe*/) const +{ + // use the BVH intersector in combination with leaf ray-triangle testing + double local_step = Big(); // we need this otherwise the lambda get's confused + + using Scalar = float; + using Vec3 = bvh::v2::Vec; + using Node = bvh::v2::Node; + using Bvh = bvh::v2::Bvh; + using Ray = bvh::v2::Ray; + + // let's fetch the bvh + auto mybvh = (Bvh*)fBVH; + if (!mybvh) { + assert(false); + return -1.; + } + + auto truncate_roundup = [](double orig) { + float epsilon = std::numeric_limits::epsilon() * std::fabs(orig); + // Add the bias to x before assigning it to y + return static_cast(orig + epsilon); + }; + + // the ray used for bvh interaction + Ray ray(Vec3(point[0], point[1], point[2]), // origin + Vec3(dir[0], dir[1], dir[2]), // direction + 0., // minimum distance (could give stepmax ?) + truncate_roundup(local_step)); + + static constexpr bool use_robust_traversal = true; + + Vertex_t dir_v{dir[0], dir[1], dir[2]}; + // Traverse the BVH and apply concrete object intersection in BVH leafs + bvh::v2::GrowingStack stack; + mybvh->intersect(ray, mybvh->get_root().index, stack, [&](size_t begin, size_t end) { + for (size_t prim_id = begin; prim_id < end; ++prim_id) { + auto objectid = mybvh->prim_ids[prim_id]; + auto facet = fFacets[objectid]; + const auto& n = fOutwardNormals[objectid]; + + // Only exiting surfaces are relevant (from inside--> dot product must be positive) + if (n.Dot(dir_v) <= 0.) { + continue; + } + + const auto& v0 = fVertices[facet[0]]; + const auto& v1 = fVertices[facet[1]]; + const auto& v2 = fVertices[facet[2]]; + + const double t = + rayTriangle(Vertex_t{point[0], point[1], point[2]}, dir_v, v0, v1, v2, 0.); + if (t < local_step) { + local_step = t; + } + } + return false; // go on after this + }); + + return local_step; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Capacity + +Double_t O2Tessellated::Capacity() const +{ + // For explanation of the following algorithm see: + // https://en.wikipedia.org/wiki/Polyhedron#Volume + // http://wwwf.imperial.ac.uk/~rn/centroid.pdf + + double vol = 0.0; + for (size_t i = 0; i < fFacets.size(); ++i) { + auto& facet = fFacets[i]; + auto a = fVertices[facet[0]]; + auto b = fVertices[facet[1]]; + auto c = fVertices[facet[2]]; + vol += + a[0] * (b[1] * c[2] - b[2] * c[1]) + b[0] * (c[1] * a[2] - c[2] * a[1]) + c[0] * (a[1] * b[2] - a[2] * b[1]); + } + return vol / 6.0; +} + +//////////////////////////////////////////////////////////////////////////////// +/// BuildBVH + +void O2Tessellated::BuildBVH() +{ + using Scalar = float; + using BBox = bvh::v2::BBox; + using Vec3 = bvh::v2::Vec; + using Node = bvh::v2::Node; + using Bvh = bvh::v2::Bvh; + + // helper determining axis aligned bounding box from a facet; + auto GetBoundingBox = [this](TGeoFacet const& facet) { +#ifndef NDEBUG + const auto nvertices = facet.GetNvert(); + assert(nvertices == 3); // for now only triangles +#endif + const auto& v1 = fVertices[facet[0]]; + const auto& v2 = fVertices[facet[1]]; + const auto& v3 = fVertices[facet[2]]; + BBox bbox; + bbox.min[0] = std::min(std::min(v1[0], v2[0]), v3[0]) - 0.001f; + bbox.min[1] = std::min(std::min(v1[1], v2[1]), v3[1]) - 0.001f; + bbox.min[2] = std::min(std::min(v1[2], v2[2]), v3[2]) - 0.001f; + bbox.max[0] = std::max(std::max(v1[0], v2[0]), v3[0]) + 0.001f; + bbox.max[1] = std::max(std::max(v1[1], v2[1]), v3[1]) + 0.001f; + bbox.max[2] = std::max(std::max(v1[2], v2[2]), v3[2]) + 0.001f; + return bbox; + }; + + // we need bounding boxes enclosing the primitives and centers of primitives + // (replaced here by centers of bounding boxes) to build the bvh + std::vector bboxes; + std::vector centers; + + // loop over all the triangles/Facets; + int nd = fFacets.size(); + for (int i = 0; i < nd; ++i) { + auto& facet = fFacets[i]; + + // fetch the bounding box of this node and add to the vector of bounding boxes + (bboxes).push_back(GetBoundingBox(facet)); + centers.emplace_back((bboxes).back().get_center()); + } + + // check if some previous object is registered and delete if necessary + if (fBVH) { + delete (Bvh*)fBVH; + fBVH = nullptr; + } + + // create the bvh + typename bvh::v2::DefaultBuilder::Config config; + config.quality = bvh::v2::DefaultBuilder::Quality::High; + auto bvh = bvh::v2::DefaultBuilder::build(bboxes, centers, config); + auto bvhptr = new Bvh; + *bvhptr = std::move(bvh); // copy structure + fBVH = (void*)(bvhptr); + + return; +} + +//////////////////////////////////////////////////////////////////////////////// +/// Contains + +bool O2Tessellated::Contains(Double_t const* point) const +{ + // we do the parity test + using Scalar = float; + using Vec3 = bvh::v2::Vec; + using Node = bvh::v2::Node; + using Bvh = bvh::v2::Bvh; + using Ray = bvh::v2::Ray; + + // let's fetch the bvh + auto mybvh = (Bvh*)fBVH; + if (!mybvh) { + assert(false); + return -1.; + } + + auto truncate_roundup = [](double orig) { + float epsilon = std::numeric_limits::epsilon() * std::fabs(orig); + // Add the bias to x before assigning it to y + return static_cast(orig + epsilon); + }; + + // let's do very quick checks against the top node + if (!TGeoBBox::Contains(point)) { + return false; + } + + // An arbitrary test direction. + // Doesn't need to be normalized and probes all normals. Also ensuring to be skewed somewhat + // without evident symmetries. + Vertex_t test_dir{1.0, 1.41421356237, 1.73205080757}; + + double local_step = Big(); + // the ray used for bvh interaction + Ray ray(Vec3(point[0], point[1], point[2]), // origin + Vec3(test_dir[0], test_dir[1], test_dir[2]), // direction + 0.0f, // minimum distance (could give stepmax ?) + truncate_roundup(local_step)); + + static constexpr bool use_robust_traversal = true; + + // Traverse the BVH and apply concrete object intersection in BVH leafs + bvh::v2::GrowingStack stack; + size_t crossings = 0; + mybvh->intersect(ray, mybvh->get_root().index, stack, [&](size_t begin, size_t end) { + for (size_t prim_id = begin; prim_id < end; ++prim_id) { + auto objectid = mybvh->prim_ids[prim_id]; + auto& facet = fFacets[objectid]; + + // for the parity test, we probe all crossing surfaces + const auto& v0 = fVertices[facet[0]]; + const auto& v1 = fVertices[facet[1]]; + const auto& v2 = fVertices[facet[2]]; + + const double t = rayTriangle(Vertex_t(point[0], point[1], point[2]), + test_dir, v0, v1, v2, 0.); + + if (t != std::numeric_limits::infinity()) { + ++crossings; + } + } + return false; + }); + + return crossings & 1; +} + +namespace +{ + +// Helper classes/structs used for priority queue - BVH traversal +// structure keeping cost (value) for a BVH index +struct BVHPrioElement { + size_t bvh_node_id; + float value; +}; + +// A priority queue for BVHPrioElement with an additional clear method +// for quick reset. We intentionally derive from std::priority_queue here to expose a +// clear() convenience method via access to the protected container `c`. +// This is internal, non-polymorphic code and relies on standard-library +// implementation details that are stable across supported platforms. +template +class BVHPrioQueue : public std::priority_queue, Comparator> +{ + public: + using std::priority_queue, + Comparator>::priority_queue; // constructor inclusion + + // convenience method to quickly clear/reset the queue (instead of having to pop one by one) + void clear() { this->c.clear(); } +}; + +} // namespace + +/// a reusable safety kernel, which optionally returns the closest face +template +inline Double_t O2Tessellated::SafetyKernel(const Double_t* point, bool in, int* closest_facet_id) const +{ + // This is the classic traversal/pruning of a BVH based on priority queue search + + float smallest_safety_sq = TGeoShape::Big(); + + using Scalar = float; + using Vec3 = bvh::v2::Vec; + using Node = bvh::v2::Node; + using Bvh = bvh::v2::Bvh; + + // let's fetch the bvh + auto mybvh = (Bvh*)fBVH; + + // testpoint object in float for quick BVH interaction + Vec3 testpoint(point[0], point[1], point[2]); + + auto currnode = mybvh->nodes[0]; // we start from the top BVH node + // we do a quick check on the top node (in case we are outside shape) + bool outside_top = false; + if (!in) { + outside_top = !bvh::v2::extra::contains(currnode.get_bbox(), testpoint); + if (outside_top) { + const auto safety_sq_to_top = bvh::v2::extra::SafetySqToNode(currnode.get_bbox(), testpoint); + // we simply return safety to the outer bounding box as an estimate + return std::sqrt(safety_sq_to_top); + } + } + + // comparator bringing out "smallest" value on top + auto cmp = [](BVHPrioElement a, BVHPrioElement b) { return a.value > b.value; }; + static thread_local BVHPrioQueue queue(cmp); + queue.clear(); + + // algorithm is based on standard iterative tree traversal with priority queues + float current_safety_to_node_sq = 0.f; + + if (returnFace) { + *closest_facet_id = -1; + } + + do { + if (currnode.is_leaf()) { + // we are in a leaf node and actually talk to a face/triangular primitive + const auto begin_prim_id = currnode.index.first_id(); + const auto end_prim_id = begin_prim_id + currnode.index.prim_count(); + + for (auto p_id = begin_prim_id; p_id < end_prim_id; p_id++) { + const auto object_id = mybvh->prim_ids[p_id]; + + const auto& facet = fFacets[object_id]; + const auto& v1 = fVertices[facet[0]]; + const auto& v2 = fVertices[facet[1]]; + const auto& v3 = fVertices[facet[2]]; + + auto thissafetySQ = pointTriangleDistSq(Vec3f{point[0], point[1], point[2]}, Vec3f{v1[0], v1[1], v1[2]}, + Vec3f{v2[0], v2[1], v2[2]}, Vec3f{v3[0], v3[1], v3[2]}); + + if (thissafetySQ < smallest_safety_sq) { + smallest_safety_sq = thissafetySQ; + if (returnFace) { + *closest_facet_id = object_id; + } + } + } + } else { + // not a leave node ... for further traversal, + // we inject the children into priority queue based on distance to it's bounding box + const auto leftchild_id = currnode.index.first_id(); + const auto rightchild_id = leftchild_id + 1; + + for (size_t childid : {leftchild_id, rightchild_id}) { + if (childid >= mybvh->nodes.size()) { + continue; + } + + const auto& node = mybvh->nodes[childid]; + const auto inside = bvh::v2::extra::contains(node.get_bbox(), testpoint); + + if (inside) { + // this must be further considered because we are inside the bounding box + queue.push(BVHPrioElement{childid, -1.}); + } else { + auto safety_to_node_square = bvh::v2::extra::SafetySqToNode(node.get_bbox(), testpoint); + if (safety_to_node_square <= smallest_safety_sq) { + // this should be further considered + queue.push(BVHPrioElement{childid, safety_to_node_square}); + } + } + } + } + + if (queue.size() > 0) { + auto currElement = queue.top(); + currnode = mybvh->nodes[currElement.bvh_node_id]; + current_safety_to_node_sq = currElement.value; + queue.pop(); + } else { + break; + } + } while (current_safety_to_node_sq <= smallest_safety_sq); + + return std::nextafter(std::sqrt(smallest_safety_sq), 0.0f); +} + +//////////////////////////////////////////////////////////////////////////////// +/// Safety + +Double_t O2Tessellated::Safety(const Double_t* point, Bool_t in) const +{ + // we could use some caching here (in future) since queries to the solid will likely + // be made with some locality + + // fall-back to precise safety kernel + return SafetyKernel(point, in); +} + +//////////////////////////////////////////////////////////////////////////////// +/// ComputeNormal interface + +void O2Tessellated::ComputeNormal(const Double_t* point, const Double_t* dir, Double_t* norm) const +{ + // We take the approach to identify closest facet to the point via safety + // and returning the normal from this face. + + // TODO: Before doing that we could check for cached points from other queries + + // use safety kernel + int closest_face_id = -1; + SafetyKernel(point, true, &closest_face_id); + + if (closest_face_id < 0) { + norm[0] = 1.; + norm[1] = 0.; + norm[2] = 0.; + return; + } + + const auto& n = fOutwardNormals[closest_face_id]; + norm[0] = n[0]; + norm[1] = n[1]; + norm[2] = n[2]; + + // change sign depending on dir + if (norm[0] * dir[0] + norm[1] * dir[1] + norm[2] * dir[2] < 0) { + norm[0] = -norm[0]; + norm[1] = -norm[1]; + norm[2] = -norm[2]; + } + return; +} + +//////////////////////////////////////////////////////////////////////////////// +/// trivial (non-BVH) DistFromInside function + +Double_t O2Tessellated::DistFromInside_Loop(const Double_t* point, const Double_t* dir) const +{ + Vertex_t p(point[0], point[1], point[2]); + Vertex_t d(dir[0], dir[1], dir[2]); + + double dist = Big(); + for (size_t i = 0; i < fFacets.size(); ++i) { + const auto& facet = fFacets[i]; + const auto& n = fOutwardNormals[i]; + + // Only exiting surfaces are relevant (from inside--> dot product must be positive) + if (n.Dot(d) <= 0.0) { + continue; + } + + const auto& v0 = fVertices[facet[0]]; + const auto& v1 = fVertices[facet[1]]; + const auto& v2 = fVertices[facet[2]]; + + const double t = rayTriangle(p, d, v0, v1, v2, 0.); + + if (t < dist) { + dist = t; + } + } + return dist; +} + +//////////////////////////////////////////////////////////////////////////////// +/// trivial (non-BVH) DistFromOutside function + +Double_t O2Tessellated::DistFromOutside_Loop(const Double_t* point, const Double_t* dir) const +{ + Vertex_t p(point[0], point[1], point[2]); + Vertex_t d(dir[0], dir[1], dir[2]); + + double dist = Big(); + for (size_t i = 0; i < fFacets.size(); ++i) { + const auto& facet = fFacets[i]; + const auto& n = fOutwardNormals[i]; + + // Only exiting surfaces are relevant (from outside, the dot product must be negative) + if (n.Dot(d) > 0.0) { + continue; + } + + const auto& v0 = fVertices[facet[0]]; + const auto& v1 = fVertices[facet[1]]; + const auto& v2 = fVertices[facet[2]]; + + const double t = rayTriangle(p, d, v0, v1, v2, 0.); + + if (t < dist) { + dist = t; + } + } + return dist; +} + +//////////////////////////////////////////////////////////////////////////////// +/// trivial (non-BVH) Contains + +bool O2Tessellated::Contains_Loop(const Double_t* point) const +{ + // Fixed ray direction + const Vertex_t test_dir{1.0, 1.41421356237, 1.73205080757}; + + Vertex_t p(point[0], point[1], point[2]); + + int crossings = 0; + for (size_t i = 0; i < fFacets.size(); ++i) { + const auto& facet = fFacets[i]; + + const auto& v0 = fVertices[facet[0]]; + const auto& v1 = fVertices[facet[1]]; + const auto& v2 = fVertices[facet[2]]; + + const double t = rayTriangle(p, test_dir, v0, v1, v2, 0.); + if (t != std::numeric_limits::infinity()) { + ++crossings; + } + } + return (crossings & 1); +} + +//////////////////////////////////////////////////////////////////////////////// +/// Custom streamer which performs Closing on read. +/// Recalculation of BVH and normals is fast + +void O2Tessellated::Streamer(TBuffer& b) +{ + if (b.IsReading()) { + b.ReadClassBuffer(O2Tessellated::Class(), this); + CloseShape(false); // close shape but do not re-perform checks + } else { + b.WriteClassBuffer(O2Tessellated::Class(), this); + } +} + +//////////////////////////////////////////////////////////////////////////////// +/// Calculate the normals + +void O2Tessellated::CalculateNormals() +{ + fOutwardNormals.clear(); + for (auto& facet : fFacets) { + auto& v1 = fVertices[facet[0]]; + auto& v2 = fVertices[facet[1]]; + auto& v3 = fVertices[facet[2]]; + using Vec3d = Vec3f; + auto norm = triangleNormal(Vec3d{v1[0], v1[1], v1[2]}, Vec3d{v2[0], v2[1], v2[2]}, Vec3d{v3[0], v3[1], v3[2]}); + fOutwardNormals.emplace_back(Vertex_t{norm.x, norm.y, norm.z}); + } +} diff --git a/Detectors/Base/src/TGeoGeometryUtils.cxx b/Detectors/Base/src/TGeoGeometryUtils.cxx new file mode 100644 index 0000000000000..13d142cfb29be --- /dev/null +++ b/Detectors/Base/src/TGeoGeometryUtils.cxx @@ -0,0 +1,142 @@ +// 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 TGeoGeometryUtils.cxx +/// \author Sandro Wenzel (CERN) +/// \brief Collection of utility functions for TGeo + +#include +#include +#include +#include +#include + +namespace o2 +{ +namespace base +{ + +namespace +{ +// some helpers to interpret TGeo TBuffer3D output +// and convert it to surface triangles (reengineered from TGeo code) + +std::vector BuildVertexLoop(const TBuffer3D& buf, + const std::vector& segs) +{ + // adjacency list + std::unordered_map> adj; + + for (int s : segs) { + int a = buf.fSegs[3 * s + 1]; + int b = buf.fSegs[3 * s + 2]; + adj[a].push_back(b); + adj[b].push_back(a); + } + + // start from any vertex + int start = adj.begin()->first; + int prev = -1; + int curr = start; + + std::vector loop; + + while (true) { + loop.push_back(curr); + + const auto& nbrs = adj[curr]; + int next = -1; + + for (int n : nbrs) { + if (n != prev) { + next = n; + break; + } + } + + if (next == -1 || next == start) + break; + + prev = curr; + curr = next; + } + return loop; +} + +std::vector> ExtractPolygons(const TBuffer3D& buf) +{ + std::vector> polys; + Int_t idx = 0; + + for (Int_t ip = 0; ip < buf.NbPols(); ++ip) { + + idx++; // color + Int_t nseg = buf.fPols[idx++]; + + std::vector segs(nseg); + for (Int_t i = 0; i < nseg; ++i) { + segs[i] = buf.fPols[idx++]; + } + + auto verts = BuildVertexLoop(buf, segs); + if (verts.size() >= 3) { + polys.push_back(std::move(verts)); + } + } + + return polys; +} + +std::vector> + Triangulate(const std::vector>& polys) +{ + std::vector> tris; + for (const auto& poly : polys) { + int nv = poly.size(); + if (nv < 3) + continue; + + int v0 = poly[0]; + for (int i = 1; i < nv - 1; ++i) { + tris.push_back({{v0, poly[i], poly[i + 1]}}); + } + } + return tris; +} + +TGeoTessellated* MakeTessellated(const TBuffer3D& buf) +{ + auto polys = ExtractPolygons(buf); + auto tris = Triangulate(polys); + int i = 0; + auto* tess = new TGeoTessellated("tess"); + const Double_t* p = buf.fPnts; + for (auto& t : tris) { + tess->AddFacet( + TGeoTessellated::Vertex_t{p[3 * t[0]], p[3 * t[0] + 1], p[3 * t[0] + 2]}, + TGeoTessellated::Vertex_t{p[3 * t[1]], p[3 * t[1] + 1], p[3 * t[1] + 2]}, + TGeoTessellated::Vertex_t{p[3 * t[2]], p[3 * t[2] + 1], p[3 * t[2] + 2]}); + } + tess->CloseShape(); + return tess; +} +} // end anonymous namespace + +///< Transform any (primitive) TGeoShape to a TGeoTessellated +TGeoTessellated* TGeoGeometryUtils::TGeoShapeToTGeoTessellated(TGeoShape const* shape) +{ + auto& buf = shape->GetBuffer3D(TBuffer3D::kRawSizes | TBuffer3D::kRaw | TBuffer3D::kCore, false); + auto tes = MakeTessellated(buf); + return tes; +} + +} // namespace base +} // namespace o2 diff --git a/Detectors/Base/src/bvh2_extra_kernels.h b/Detectors/Base/src/bvh2_extra_kernels.h new file mode 100644 index 0000000000000..e466c58c9071c --- /dev/null +++ b/Detectors/Base/src/bvh2_extra_kernels.h @@ -0,0 +1,66 @@ +#ifndef ROOT_GEOM_BVH2_EXTRA + +namespace bvh::v2::extra +{ + +// reusable geometry kernels used in multiple places +// for interaction with BVH2 structures + +// determines if a point is inside the bounding box +template +bool contains(bvh::v2::BBox const& box, bvh::v2::Vec const& p) +{ + auto min = box.min; + auto max = box.max; + return (p[0] >= min[0] && p[0] <= max[0]) && (p[1] >= min[1] && p[1] <= max[1]) && + (p[2] >= min[2] && p[2] <= max[2]); +} + +// determines the largest squared distance of point to any of the bounding box corners +template +auto RmaxSqToNode(bvh::v2::BBox const& box, bvh::v2::Vec const& p) +{ + // construct the 8 corners to get the maximal distance + const auto minCorner = box.min; + const auto maxCorner = box.max; + using Vec3 = bvh::v2::Vec; + // these are the corners of the bounding box + const std::array, 8> corners{ + Vec3{minCorner[0], minCorner[1], minCorner[2]}, Vec3{minCorner[0], minCorner[1], maxCorner[2]}, + Vec3{minCorner[0], maxCorner[1], minCorner[2]}, Vec3{minCorner[0], maxCorner[1], maxCorner[2]}, + Vec3{maxCorner[0], minCorner[1], minCorner[2]}, Vec3{maxCorner[0], minCorner[1], maxCorner[2]}, + Vec3{maxCorner[0], maxCorner[1], minCorner[2]}, Vec3{maxCorner[0], maxCorner[1], maxCorner[2]}}; + + T Rmax_sq{0}; + for (const auto& corner : corners) { + float R_sq = 0.; + const auto dx = corner[0] - p[0]; + R_sq += dx * dx; + const auto dy = corner[1] - p[1]; + R_sq += dy * dy; + const auto dz = corner[2] - p[2]; + R_sq += dz * dz; + Rmax_sq = std::max(Rmax_sq, R_sq); + } + return Rmax_sq; +}; + +// determines the minimum squared distance of point to a bounding box ("safey square") +template +auto SafetySqToNode(bvh::v2::BBox const& box, bvh::v2::Vec const& p) +{ + T sqDist{0.0}; + for (int i = 0; i < 3; i++) { + T v = p[i]; + if (v < box.min[i]) { + sqDist += (box.min[i] - v) * (box.min[i] - v); + } else if (v > box.max[i]) { + sqDist += (v - box.max[i]) * (v - box.max[i]); + } + } + return sqDist; +}; + +} // namespace bvh::v2::extra + +#endif \ No newline at end of file diff --git a/Detectors/Base/src/bvh2_third_party.h b/Detectors/Base/src/bvh2_third_party.h new file mode 100644 index 0000000000000..0611997d39508 --- /dev/null +++ b/Detectors/Base/src/bvh2_third_party.h @@ -0,0 +1,36 @@ +#ifndef ROOT_GEOM_BVH2_THIRD_PARTY + +// A single entry header into third-party BVH2 +// Good place to manage compiler warnings etc. + +#if defined(__clang__) +#pragma clang diagnostic push +#pragma clang diagnostic ignored "-Wshadow" +#pragma clang diagnostic ignored "-Wpsabi" +#elif defined(__GNUC__) +#pragma GCC diagnostic push +#pragma GCC diagnostic ignored "-Wshadow" +#pragma GCC diagnostic ignored "-Wpsabi" +#pragma GCC diagnostic ignored "-Wall" +#pragma GCC diagnostic ignored "-Wshadow" +#pragma GCC diagnostic ignored "-Wunknown-pragmas" +#pragma GCC diagnostic ignored "-Wattributes" +#elif defined(_MSC_VER) +#pragma warning(push) +#pragma warning(disable : 5051) +#endif + +#include +#include +#include +#include +#include +#include + +#if defined(__clang__) +#pragma clang diagnostic pop +#elif defined(__GNUC__) +#pragma GCC diagnostic pop +#endif + +#endif \ No newline at end of file diff --git a/Detectors/Passive/CMakeLists.txt b/Detectors/Passive/CMakeLists.txt index 0976530bc6571..a24954ad10539 100644 --- a/Detectors/Passive/CMakeLists.txt +++ b/Detectors/Passive/CMakeLists.txt @@ -23,6 +23,7 @@ o2_add_library(DetectorsPassive src/Hall.cxx src/HallSimParam.cxx src/PassiveBase.cxx + src/ExternalModule.cxx PUBLIC_LINK_LIBRARIES O2::Field O2::DetectorsBase O2::SimConfig) o2_target_root_dictionary(DetectorsPassive @@ -39,6 +40,7 @@ o2_target_root_dictionary(DetectorsPassive include/DetectorsPassive/Hall.h include/DetectorsPassive/HallSimParam.h include/DetectorsPassive/PassiveBase.h + include/DetectorsPassive/ExternalModule.h LINKDEF src/PassiveLinkDef.h) # FIXME: if PutFrameInTop really depends on TRD, then the following can not work diff --git a/Detectors/Passive/include/DetectorsPassive/ExternalModule.h b/Detectors/Passive/include/DetectorsPassive/ExternalModule.h new file mode 100644 index 0000000000000..155870ae42a6d --- /dev/null +++ b/Detectors/Passive/include/DetectorsPassive/ExternalModule.h @@ -0,0 +1,64 @@ +// 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. + +#ifndef ALICEO2_PASSIVE_EXTERNALMODULE_H +#define ALICEO2_PASSIVE_EXTERNALMODULE_H + +#include "DetectorsPassive/PassiveBase.h" // base class of passive modules +#include "Rtypes.h" // for Pipe::Class, ClassDef, Pipe::Streamer + +class TGeoVolume; +class TGeoTransformation; + +namespace o2 +{ +namespace passive +{ + +// options used to configure a generic plug and play external module +struct ExternalModuleOptions { + std::string root_macro_file; // the file where to lookup the ROOT geometry building macro + std::string top_volume; // the volume to be added + std::string anchor_volume; // the volume into which the detector will be hooked + TGeoMatrix const* placement = nullptr; // how to place the module inside anchor_volume +}; + +// a module (passive material) defined externally (ROOT macro / GDML / TGeo geometry) +class ExternalModule : public PassiveBase +{ + public: + ExternalModule(const char* name, const char* long_title, ExternalModuleOptions options); + ExternalModule() = default; // default constructor + + ~ExternalModule() override = default; + void ConstructGeometry() override; + + /// 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 +} // namespace o2 +#endif diff --git a/Detectors/Passive/src/ExternalModule.cxx b/Detectors/Passive/src/ExternalModule.cxx new file mode 100644 index 0000000000000..414863883b714 --- /dev/null +++ b/Detectors/Passive/src/ExternalModule.cxx @@ -0,0 +1,162 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// ClassImp(o2::passive::ExternalModule) + +namespace o2::passive +{ + +ExternalModule::ExternalModule(const char* name, const char* long_title, ExternalModuleOptions options) : PassiveBase(name, long_title), mOptions(options) +{ +} + +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(); + if (!module_top) { + LOG(error) << "No module found\n"; + return; + } + + remapMedia(const_cast(module_top)); + + // 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(const_cast(module_top), 1, const_cast(mOptions.placement)); +} + +bool ExternalModule::initGeomBuilderHook() +{ + 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; + } + } + return false; +} + +} // namespace o2::passive \ No newline at end of file diff --git a/Steer/include/Steer/O2MCApplicationBase.h b/Steer/include/Steer/O2MCApplicationBase.h index 36966be9bde62..d61199baba0ae 100644 --- a/Steer/include/Steer/O2MCApplicationBase.h +++ b/Steer/include/Steer/O2MCApplicationBase.h @@ -58,6 +58,8 @@ class O2MCApplicationBase : public FairMCApplication typedef std::function TrackRefFcn; + void fixTGeoRuntimeShapes(); + protected: o2::conf::SimCutParams const& mCutParams; // reference to parameter system unsigned long long mStepCounter{0}; diff --git a/Steer/src/O2MCApplication.cxx b/Steer/src/O2MCApplication.cxx index 584598d350581..f832ab70ab121 100644 --- a/Steer/src/O2MCApplication.cxx +++ b/Steer/src/O2MCApplication.cxx @@ -37,6 +37,11 @@ #include "SimConfig/GlobalProcessCutSimParam.h" #include "DetectorsBase/GeometryManagerParam.h" #include +#include +#include +#include +#include +#include namespace o2 { @@ -209,10 +214,65 @@ bool O2MCApplicationBase::MisalignGeometry() gGeoManager->SetUseParallelWorldNav(true); } + // performs possible optimizations (shape replacements on the runtime geometry) + fixTGeoRuntimeShapes(); + // return original return value of misalignment procedure return true; } +void O2MCApplicationBase::fixTGeoRuntimeShapes() +{ + // Replace TGeo shapes by other ones for performance or other reasons. + // Should only affect runtime of simulation. + + // TODO: make this configurable via external JSON rules/macro + + // Also delete original shapes for memory reasons + + // We follow a visitor pattern on a geom hierarchy + // for now replace a TGeoTessellate by our own implementation + std::unordered_set volumes_visited; + std::unordered_set old_shape_pointers; + + std::function visit; + visit = [&](TGeoNode* node) -> void { + if (!node) { + return; + } + auto vol = node->GetVolume(); + if (volumes_visited.find(vol) != volumes_visited.end()) { + return; + } + volumes_visited.insert(vol); + + // transform the shape of this volume + auto shape = vol->GetShape(); + if (shape->IsA() == TGeoTessellated::Class()) { + auto tsl = static_cast(shape); + + // make a new O2Tessellated until ROOT has proper support for navigation in TGeoTessellated + std::cout << "Converting to O2Tessellated for vol " << vol->GetName() << "\n"; + auto replacement_shape = new o2::base::O2Tessellated(*tsl, false); + vol->SetShape(replacement_shape); + old_shape_pointers.insert(shape); + } + // other cases could come here + + for (int i = 0; i < vol->GetNdaughters(); ++i) { + auto child_node = vol->GetNode(i); + visit(child_node); + } + }; + + visit(gGeoManager->GetTopNode()); + + for (auto ptr : old_shape_pointers) { + delete ptr; + ptr = nullptr; + } +} + void O2MCApplicationBase::finishEventCommon() { LOG(info) << "This event/chunk did " << mStepCounter << " steps"; diff --git a/macro/build_geometry.C b/macro/build_geometry.C index 6b13f2eac2766..b1266cd58b5eb 100644 --- a/macro/build_geometry.C +++ b/macro/build_geometry.C @@ -63,6 +63,8 @@ #include #endif +#include + using Return = o2::base::Detector*; void finalize_geometry(FairRunSim* run); @@ -182,6 +184,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)); + } + // the absorber if (isActivated("ABSO")) { // the frame structure to support other detectors diff --git a/scripts/geometry/O2_CADtoTGeo.py b/scripts/geometry/O2_CADtoTGeo.py new file mode 100644 index 0000000000000..2686d7a24e3d9 --- /dev/null +++ b/scripts/geometry/O2_CADtoTGeo.py @@ -0,0 +1,602 @@ +#!/usr/bin/env python3 +""" +A Python script, doing a deep STEP/XCAF -> ROOT TGeo conversion. +For now, all CAD solids are simply meshed. The ROOT geometry is build as a C++ ROOT macro +and facet data is stored in binary form to keep disc space minimal. + +Generates (into --output-folder): + - geom.C (small ROOT macro) + - facets__.bin for each leaf logical volume (float32 triangles) + +Facet file format (little-endian): + uint32 nTriangles + then nTriangles * 9 * float32: + ax ay az bx by bz cx cy cz + +VOLNAME is a filename-safe version of the XCAF label name when available (e.g. "nut"), +and LID is the XCAF label entry (e.g. "0:1:1:7" -> "0_1_1_7") to keep filenames unique. + +Naming: + - C++ variable names stay based on XCAF label entry (e.g. 0:1:1:7) for uniqueness. + - ROOT object names (TGeoVolume / TGeoTessellated / TGeoVolumeAssembly) use the label's + human name when available (e.g. "nut", "rod-assembly"), falling back to the entry. + +Units: + - By default, the script tries to detect the STEP LENGTH unit by scanning the STEP file + header/contents (common patterns like .MILLI. / .CENTI. / .METRE. / INCH / FOOT). + - You can override with --step-unit {auto,mm,cm,m,in,ft}. TGeo expects cm. + +Author: + - Sandro Wenzel, CERN (02/2026) +""" + +import warnings +warnings.filterwarnings("ignore", message=".*all to deprecated function.*", category=DeprecationWarning) + +import argparse +import re +import struct +from pathlib import Path as _Path + +from OCC.Core.Bnd import Bnd_Box +from OCC.Core.BRepBndLib import brepbndlib +from OCC.Core.BRepMesh import BRepMesh_IncrementalMesh +from OCC.Core.BRep import BRep_Tool +from OCC.Core.TopLoc import TopLoc_Location +from OCC.Core.TopAbs import TopAbs_REVERSED +from OCC.Extend.TopologyUtils import TopologyExplorer + +from OCC.Core.STEPCAFControl import STEPCAFControl_Reader +from OCC.Core.TDocStd import TDocStd_Document +from OCC.Core.XCAFDoc import XCAFDoc_DocumentTool +from OCC.Core.IFSelect import IFSelect_RetDone + +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 + + +# ------------------------------- +# STEP/XCAF loading +# ------------------------------- + +def load_step_with_xcaf(path: str): + doc = TDocStd_Document("pythonocc-doc") + reader = STEPCAFControl_Reader() + reader.SetColorMode(True) + reader.SetNameMode(True) + reader.SetLayerMode(True) + + status = reader.ReadFile(path) + if status != IFSelect_RetDone: + raise RuntimeError(f"STEP read failed for: {path}") + + reader.Transfer(doc) + shape_tool = XCAFDoc_DocumentTool.ShapeTool(doc.Main()) + return doc, shape_tool + + +def label_id(label: TDF_Label) -> str: + s = TCollection_AsciiString() + TDF_Tool.Entry(label, s) + return s.ToCString() + + +def label_name(label: TDF_Label) -> str: + # Uses the XCAF/STEP name when present; can be empty. + try: + n = label.GetLabelName() + if n: + return str(n) + except Exception: + pass + return "" + + +# ------------------------------- +# Units +# ------------------------------- + +def step_unit_scale_to_cm(step_unit: str) -> float: + step_unit = (step_unit or "auto").lower() + if step_unit == "mm": + return 0.1 + if step_unit == "cm": + return 1.0 + if step_unit == "m": + return 100.0 + if step_unit == "in": + return 2.54 + if step_unit == "ft": + return 30.48 + raise ValueError(f"Unknown --step-unit {step_unit} (use auto, mm, cm, m, in, ft)") + + +def detect_step_length_unit(step_path: str) -> str: + """ + Heuristic unit detection by scanning STEP file text for common unit tokens. + This avoids relying on OCCT APIs that can vary across pythonOCC builds. + + Returns one of: mm, cm, m, in, ft. Defaults to mm if uncertain. + """ + p = _Path(step_path) + # STEP can be huge: read only the first few MB; units are near the header. + max_bytes = 4 * 1024 * 1024 + data = p.open("rb").read(max_bytes).decode("latin-1", errors="ignore").upper() + + if ".MILLI." in data: + return "mm" + if ".CENTI." in data: + return "cm" + if ".METRE." in data or ".METER." in data: + return "m" + if "INCH" in data: + return "in" + if "FOOT" in data or "FEET" in data: + return "ft" + + # Conservative default for mechanical CAD STEP is mm + return "mm" + + +# ------------------------------- +# Triangulation helpers +# ------------------------------- + +def _scale_triangles(triangles, s: float): + if s == 1.0: + return triangles + out = [] + for (a, b, c) in triangles: + out.append(( + (a[0] * s, a[1] * s, a[2] * s), + (b[0] * s, b[1] * s, b[2] * s), + (c[0] * s, c[1] * s, c[2] * s), + )) + return out + + +def triangulate_asbbox(shape, scale_to_cm: float = 1.0): + box = Bnd_Box() + brepbndlib.Add(shape, box) + xmin, ymin, zmin, xmax, ymax, zmax = box.Get() + + p000 = (xmin, ymin, zmin) + p001 = (xmin, ymin, zmax) + p010 = (xmin, ymax, zmin) + p011 = (xmin, ymax, zmax) + p100 = (xmax, ymin, zmin) + p101 = (xmax, ymin, zmax) + p110 = (xmax, ymax, zmin) + p111 = (xmax, ymax, zmax) + + triangles = [ + (p000, p100, p110), (p000, p110, p010), + (p001, p111, p101), (p001, p011, p111), + (p000, p101, p100), (p000, p001, p101), + (p010, p110, p111), (p010, p111, p011), + (p000, p010, p011), (p000, p011, p001), + (p100, p101, p111), (p100, p111, p110), + ] + return _scale_triangles(triangles, scale_to_cm) + + +def triangulate_CAD_solid(my_solid, meshparam, scale_to_cm: float = 1.0): + lin_defl = float(meshparam.get("lin_defl", 0.1)) + ang_defl = float(meshparam.get("ang_defl", 0.1)) + + parallel = True + try: + BRepMesh_IncrementalMesh(my_solid, lin_defl, False, ang_defl, bool(parallel)) + except TypeError: + BRepMesh_IncrementalMesh(my_solid, lin_defl, False, ang_defl) + + triangles = [] + for face in TopologyExplorer(my_solid).faces(): + loc = TopLoc_Location() + triangulation = BRep_Tool.Triangulation(face, loc) + if triangulation is None: + continue + + trsf = loc.Transformation() + reverse = (face.Orientation() == TopAbs_REVERSED) + + for i in range(1, triangulation.NbTriangles() + 1): + tri = triangulation.Triangle(i) + n1, n2, n3 = tri.Get() + + p1 = triangulation.Node(n1).Transformed(trsf) + p2 = triangulation.Node(n2).Transformed(trsf) + p3 = triangulation.Node(n3).Transformed(trsf) + + if reverse: + p2, p3 = p3, p2 + + triangles.append(( + (p1.X(), p1.Y(), p1.Z()), + (p2.X(), p2.Y(), p2.Z()), + (p3.X(), p3.Y(), p3.Z()), + )) + + return _scale_triangles(triangles, scale_to_cm) + + +# ------------------------------- +# Naming helpers +# ------------------------------- + +def sanitize_cpp_name(s: str) -> str: + safe = re.sub(r"[^0-9a-zA-Z]", "_", s) + if not safe: + safe = "x" + if not (safe[0].isalpha() or safe[0] == "_"): + safe = "_" + safe + return safe + + +def sanitize_filename(s: str) -> str: + safe = re.sub(r"[^0-9a-zA-Z]", "_", s) + return safe or "x" + + +# ------------------------------- +# Binary facet IO +# ------------------------------- + +def write_facets_bin(path: _Path, triangles): + path.parent.mkdir(parents=True, exist_ok=True) + with open(path, "wb") as f: + f.write(struct.pack(" str: + m = trsf.GetRotation().GetMatrix() + t = trsf.TranslationPart() + return f""" + Double_t {name}_m[9] = {{ + {m.Value(1,1)}, {m.Value(1,2)}, {m.Value(1,3)}, + {m.Value(2,1)}, {m.Value(2,2)}, {m.Value(2,3)}, + {m.Value(3,1)}, {m.Value(3,2)}, {m.Value(3,3)} + }}; + TGeoRotation *{name}_rot = new TGeoRotation(); + {name}_rot->SetMatrix({name}_m); + TGeoCombiTrans *{name} = new TGeoCombiTrans({t.X()*scale_to_cm}, {t.Y()*scale_to_cm}, {t.Z()*scale_to_cm}, {name}_rot); +""" + + +def emit_cpp_prelude() -> str: + return """#include +#include +#include +#include +#include +#include + +static void LoadFacets(const std::string& file, TGeoTessellated* solid, bool check=false) +{ + std::ifstream in(file, std::ios::binary); + if (!in) throw std::runtime_error("Cannot open facet file: " + file); + + uint32_t nTri = 0; + in.read(reinterpret_cast(&nTri), sizeof(nTri)); + if (!in) throw std::runtime_error("Bad facet header in: " + file); + + for (uint32_t i=0;i(v), sizeof(v)); + if (!in) throw std::runtime_error("Unexpected EOF in: " + file); + + solid->AddFacet(TGeoTessellated::Vertex_t(v[0],v[1],v[2]), + TGeoTessellated::Vertex_t(v[3],v[4],v[5]), + TGeoTessellated::Vertex_t(v[6],v[7],v[8])); + } + solid->CloseShape(check, true); +} +""" + + +def emit_materials_cpp() -> str: + return """ // Default material/medium (placeholder; can be replaced later) + TGeoMaterial *mat_Default = new TGeoMaterial("Default", 0., 0., 0.); + TGeoMedium *med_Default = new TGeoMedium("Default", 1, mat_Default); +""" + + +def emit_tessellated_cpp(lid: str, vol_display_name: str, facet_abspath: str, ntriangles: int) -> str: + safe = sanitize_cpp_name(lid) + shape_name = vol_display_name if vol_display_name else lid + + if ntriangles <= 0: + out = [] + out.append(f' TGeoBBox *solid_{safe} = new TGeoBBox("{shape_name}", 0.001, 0.001, 0.001);') + out.append(f' TGeoVolume *vol_{safe} = new TGeoVolume("{shape_name}", solid_{safe}, med_Default);') + return "\n".join(out) + + out = [] + out.append(f' TGeoTessellated *solid_{safe} = new TGeoTessellated("{shape_name}", {ntriangles});') + out.append(f' LoadFacets("{facet_abspath}", solid_{safe}, check);') + out.append(f' TGeoVolume *vol_{safe} = new TGeoVolume("{shape_name}", solid_{safe}, med_Default);') + return "\n".join(out) + + +def emit_assembly_cpp(lid: str, asm_display_name: str) -> str: + safe = sanitize_cpp_name(lid) + name = asm_display_name if asm_display_name else lid + return f' TGeoVolumeAssembly *asm_{safe} = new TGeoVolumeAssembly("{name}");' + + +# ------------------------------- +# Definition graph extraction +# ------------------------------- + +logical_volumes = {} # def_lid -> triangles +def_names = {} # def_lid -> human display name (may be "") +assemblies = set() # def_lid +placements = [] # (parent_def_lid, child_def_lid, gp_Trsf local) +top_defs = set() # top definition lids +visited_defs = set() # expanded defs + + +def cpp_var_for_def(lid: str) -> str: + safe = sanitize_cpp_name(lid) + 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) + + 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] = "" + + 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) + + for i in range(children.Length()): + child = children.Value(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) + 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 + + if shape_tool.IsSimpleShape(def_label): + if def_lid not in logical_volumes: + shape = shape_tool.GetShape(def_label) + 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 + + assemblies.add(def_lid) + + +def extract_graph(step_path: str, meshparam=None, scale_to_cm: float = 1.0): + global logical_volumes, def_names, assemblies, placements, top_defs, visited_defs + logical_volumes = {} + def_names = {} + assemblies = set() + placements = [] + top_defs = set() + visited_defs = set() + + doc, shape_tool = load_step_with_xcaf(step_path) + + roots = TDF_LabelSequence() + shape_tool.GetFreeShapes(roots) + + for i in range(roots.Length()): + root = roots.Value(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) + else: + top_defs.add(label_id(root)) + expand_definition(root, shape_tool, meshparam=meshparam, scale_to_cm=scale_to_cm) + + return doc, shape_tool + + +# ------------------------------- +# ROOT macro emission +# ------------------------------- + +def emit_placement_cpp(parent_def: str, child_def: str, trsf: gp_Trsf, copy_no: int, scale_to_cm: float) -> str: + parent_cpp = cpp_var_for_def(parent_def) + child_cpp = cpp_var_for_def(child_def) + tr_name = f"tr_{sanitize_cpp_name(parent_def)}_{sanitize_cpp_name(child_def)}_{copy_no}" + return trsf_to_tgeo(trsf, tr_name, scale_to_cm) + f" {parent_cpp}->AddNode({child_cpp}, {copy_no}, {tr_name});\n" + + +def emit_root_macro(step_path: str, out_folder: _Path, meshparam=None, step_unit: str = "auto"): + if (step_unit or "auto").lower() == "auto": + detected = detect_step_length_unit(step_path) + scale_to_cm = step_unit_scale_to_cm(detected) + print(f"Detected STEP length unit: {detected} (scale to cm = {scale_to_cm})") + else: + 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) + + out_folder = out_folder.expanduser().resolve() + out_folder.mkdir(parents=True, exist_ok=True) + + facet_files = {} # def_lid -> absolute path string + for lid, tris in logical_volumes.items(): + disp = def_names.get(lid, "") + volname = sanitize_filename(disp) if disp else "vol" + lidname = sanitize_filename(lid) + fname = f"facets_{volname}_{lidname}.bin" + fpath = (out_folder / fname).resolve() + write_facets_bin(fpath, tris) + facet_files[lid] = str(fpath).replace("\\", "\\\\") # C++ string literal safety + + cpp = [] + cpp.append(emit_cpp_prelude()) + + cpp.append("TGeoVolume* build(bool check=true) {") + cpp.append(' if (!gGeoManager) { throw std::runtime_error("gGeoManager is null. Call build_and_export() or create a TGeoManager first."); }') + cpp.append(emit_materials_cpp()) + + for lid in logical_volumes.keys(): + ntriangles = len(logical_volumes[lid]) + cpp.append(emit_tessellated_cpp(lid, def_names.get(lid, ""), facet_files[lid], ntriangles)) + + for lid in sorted(assemblies): + cpp.append(emit_assembly_cpp(lid, def_names.get(lid, ""))) + + for idx, (parent, child, trsf) in enumerate(placements, start=1): + cpp.append(emit_placement_cpp(parent, child, trsf, idx, scale_to_cm)) + + if len(top_defs) == 1: + top = next(iter(top_defs)) + cpp.append(f" return {cpp_var_for_def(top)};") + else: + cpp.append(' TGeoVolumeAssembly *asm_WORLD = new TGeoVolumeAssembly("WORLD");') + for i, node in enumerate(sorted(top_defs), start=1): + cpp.append(f" asm_WORLD->AddNode({cpp_var_for_def(node)}, {i});") + cpp.append(" return asm_WORLD;") + + cpp.append("}") + + # exports a function allowing to export the geometry to TGeo file + cpp.append('void build_and_export(const char* out_root = "geom.root", bool check=true) {') + cpp.append(' if (!gGeoManager) { new TGeoManager("geom","geom"); }') + cpp.append(' TGeoVolume* top = build(check);') + cpp.append(' gGeoManager->SetTopVolume(top);') + cpp.append(' gGeoManager->CloseGeometry();') + cpp.append(' gGeoManager->CheckOverlaps();') + cpp.append(' gGeoManager->Export(out_root);') + cpp.append('}') + + # exports a function to get get hold of the builder function in ALICE O2 + cpp.append('std::function get_builder_hook_checked() {') + cpp.append(' return []() { return build(true); };') + cpp.append('}') + # exports a function to get get hold of the builder function in ALICE O2 + cpp.append('std::function get_builder_hook_unchecked() {') + cpp.append(' return []() { return build(false); };') + cpp.append('}') + + return "\n".join(cpp) + + +# ------------------------------- +# Geometry Tree printing (debug) +# ------------------------------- + +def label_entry(label): + s = TCollection_AsciiString() + TDF_Tool.Entry(label, s) + return s.ToCString() + + +def traverse_print(label, shape_tool, depth=0): + indent = " " * depth + name = label.GetLabelName() + entry = label_entry(label) + print(f"{indent}- {name} =>[{entry}]") + + if shape_tool.IsReference(label): + ref_label = TDF_Label() + shape_tool.GetReferredShape(label, ref_label) + traverse_print(ref_label, shape_tool, depth + 1) + return + + children = TDF_LabelSequence() + shape_tool.GetComponents(label, children) + if children.Length() > 0 or shape_tool.IsAssembly(label): + for i in range(children.Length()): + traverse_print(children.Value(i + 1), shape_tool, depth + 1) + return + + if shape_tool.IsSimpleShape(label): + shape = shape_tool.GetShape(label) + print(f"{indent} [LogicalShape id={id(shape)}]") + + +def print_geom(step_file): + print(f"Printing GEOM hierarchy for {step_file}") + doc, shape_tool = load_step_with_xcaf(step_file) + roots = TDF_LabelSequence() + shape_tool.GetFreeShapes(roots) + for i in range(roots.Length()): + traverse_print(roots.Value(i + 1), shape_tool) + + +# ------------------------------- +# CLI +# ------------------------------- + +def main(): + ap = argparse.ArgumentParser(description="Convert STEP/XCAF to ROOT TGeo macro, facets in per-volume binary files.") + ap.add_argument("step", help="Input STEP file") + ap.add_argument("-o", "--out", default="geom.C", help="Output ROOT macro file name (default: geom.C)") + ap.add_argument("--output-folder", default="./", help="Output folder for macro + facet files") + ap.add_argument("--out-path", default=None, help="(deprecated) Alias for --output-folder") + ap.add_argument("--mesh", action="store_true", help="Use full BRepMesh triangulation instead of bounding boxes") + 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)") + + args = ap.parse_args() + + step_path = str(_Path(args.step).expanduser().resolve()) + if args.print_tree: + print_geom(step_path) + return + + out_folder = _Path(args.output_folder) + if args.out_path is not None: + out_folder = _Path(args.out_path) + + meshparam = {"do_meshing": args.mesh, "lin_defl": args.mesh_prec, "ang_defl": args.mesh_prec} + + out_folder = out_folder.expanduser().resolve() + out_folder.mkdir(parents=True, exist_ok=True) + + out_macro = (out_folder / _Path(args.out).name).resolve() + code = emit_root_macro(step_path, out_folder, meshparam=meshparam, step_unit=args.step_unit) + out_macro.write_text(code) + + print(f"Wrote ROOT macro: {out_macro}") + print(f"Wrote facet files into: {out_folder}") + print("In ROOT you can do:") + print(f" root -l {out_macro}") + print(' build_and_export("geom.root");') + + +if __name__ == "__main__": + main() diff --git a/scripts/geometry/README.md b/scripts/geometry/README.md new file mode 100644 index 0000000000000..1815536e1af73 --- /dev/null +++ b/scripts/geometry/README.md @@ -0,0 +1,27 @@ +This is the tool O2_CADtoTGeo.py which translates from geometries in STEP format (CAD export) to +TGeo. + +To use the tool, setup a conda environment with python-occ core installed. +The following should work on standard linux x86: + +``` +# -) 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 + +# -) Create an OCC environment (for OpenCacade) +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 diff --git a/scripts/geometry/simulating_CAD_modules.md b/scripts/geometry/simulating_CAD_modules.md new file mode 100644 index 0000000000000..67a9c5136bc91 --- /dev/null +++ b/scripts/geometry/simulating_CAD_modules.md @@ -0,0 +1,72 @@ +# ALICE-O2 GEANT Simulation of CAD Geometries + +These are a few notes related to the inclusion of external (CAD-described) detector modules into the O2 simulation framework. + +## Description of the Workflow + +In principle, such integration is now possible and requires the following steps: + +1. The CAD geometry needs to be exported to STEP format and must contain only the final geometry (no artificial eta-cut elements). Ideally, the geometry should be fully hierarchical with proper solid reuse. The solids should retain their proper surface representation for detailed analysis. + +2. A tool `O2-CADtoTGeo.py` is provided to convert the STEP geometry into TGeo format. The tool is part of AliceO2 and is based on Python bindings (OCC) for OpenCascade. The tool can be used as follows: + + ```bash + python O2-CADtoTGeo.py STEP_FILE --output-folder my_detector -o geom.C --mesh \ + --mesh-prec 0.2 + ``` + + This will create a ROOT macro file `geom.C` containing the geometry description in ROOT format, as well as several binary files describing the TGeo solids. The `geom.C` file can either be used directly in ROOT to inspect the geometry or be provided to ALICE-O2 for inclusion in the geometry. + +3. Introduction of materials/media in the file `geom.C`. Currently, the file `geom.C` needs to be patched or edited to properly include `TGeoMaterial`/`TGeoMedium` definitions and connect them to the relevant `TGeoVolume` objects. At present, every solid has the same dummy material attached, which is not realistic. It may be a good idea to create a new file `geom_withMaterials.C`, which differs from `geom.C` by the addition of these material definitions. + +4. Once the conversion is complete, the module can be inserted into the O2 geometry via the `ExternalModule` class. To do so, follow this pattern in `build_geometry.C`: + + ```cpp + if (isActivated("EXT")) { + o2::passive::ExternalModuleOptions options; + options.root_macro_file = "PATH_TO_MY_DETECTOR/my_detector/geom_withMaterials.C"; + options.anchor_volume = "barrel"; // hook this into barrel + auto rot = new TGeoCombiTrans(); + rot->RotateX(90); + rot->SetDy(30); // compensate for a shift of the barrel with respect to zero + options.placement = rot; + run->AddModule(new o2::passive::ExternalModule("A3VTX", "ALICE3 beam pipe", options)); + } + ``` + +5. Create a custom detector geometry list file `my_det.json` in JSON format that includes the external detector (and any other required components, such as the L3 magnet in this example): + + ```json + { + "MY_DET": [ + "EXT", + "MAG" + ] + } + ``` + +6. Run the Geant simulation with: + + ```bash + o2-sim --detectorList MY_DET:my_det.json -g pythia8pp .... + ``` + +## Known Limitations + +- The `O2-CADtoTGeo.py` tool currently converts geometries only into TGeoTessellated solids. This may be suboptimal for primitive shapes or only an approximation for shapes with exact second-order surfaces (e.g., tubes). The precision (and therefore the number of surface triangles) can be controlled with the `--mesh-prec` parameter. The smaller the value, the more precise the mesh. + +- Meshed solids created by the tool may have issues, such as topological errors or non-watertight surfaces. It is planned to include "healing" steps via additional processing with well-known geometry kernels (e.g., CGAL). + +- The tool does not currently export materials or TGeoMedia. These must be inserted or edited manually. It is planned to make this process more automatic and user-friendly. + +- The Python tool requires the OCC Python module, which is currently not part of our software distribution. We have found it most practical to run the tool in a separate conda environment (fully decoupled from the ALICE software stack). + +- The tool currently generates a `geom.C` macro file. In the future, it may be possible to directly create an in-memory TGeo representation for deeper integration. + +- Currently, only passive modules can be integrated. Treatment of sensitive volumes or parts will be addressed in a future step. + +## Software Installation + +- The simulation must be run in the standard O2 environment built with alibuild. + +- The CAD conversion tool must currently be run in a dedicated conda environment, as described in scripts/geometry/README.md in the AliceO2 source code. \ No newline at end of file From 7d0e2e1b08ad09a2157e73797db63789af65e3a2 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Thu, 12 Feb 2026 13:06:50 +0000 Subject: [PATCH 2/2] Please consider the following formatting changes --- macro/build_geometry.C | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/macro/build_geometry.C b/macro/build_geometry.C index b1266cd58b5eb..ccc3b13fe728d 100644 --- a/macro/build_geometry.C +++ b/macro/build_geometry.C @@ -188,7 +188,7 @@ void build_geometry(FairRunSim* run = nullptr) // 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 + 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