Root/src/TGeoSurfaceConverter.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "ActsPlugins/Root/TGeoSurfaceConverter.hpp"
0010
0011 #include "Acts/Definitions/Tolerance.hpp"
0012 #include "Acts/Surfaces/AnnulusBounds.hpp"
0013 #include "Acts/Surfaces/ConvexPolygonBounds.hpp"
0014 #include "Acts/Surfaces/CylinderBounds.hpp"
0015 #include "Acts/Surfaces/CylinderSurface.hpp"
0016 #include "Acts/Surfaces/DiscSurface.hpp"
0017 #include "Acts/Surfaces/PlaneSurface.hpp"
0018 #include "Acts/Surfaces/RadialBounds.hpp"
0019 #include "Acts/Surfaces/RectangleBounds.hpp"
0020 #include "Acts/Surfaces/Surface.hpp"
0021 #include "Acts/Surfaces/TrapezoidBounds.hpp"
0022 #include "Acts/Utilities/Helpers.hpp"
0023 #include "ActsPlugins/Root/TGeoPrimitivesHelper.hpp"
0024
0025 #include <algorithm>
0026 #include <array>
0027 #include <cctype>
0028 #include <cstddef>
0029 #include <memory>
0030 #include <numbers>
0031 #include <stdexcept>
0032 #include <tuple>
0033 #include <utility>
0034 #include <vector>
0035
0036 #include <boost/algorithm/string.hpp>
0037 #include <boost/algorithm/string/predicate.hpp>
0038
0039 #include "RtypesCore.h"
0040 #include "TGeoArb8.h"
0041 #include "TGeoBBox.h"
0042 #include "TGeoBoolNode.h"
0043 #include "TGeoCompositeShape.h"
0044 #include "TGeoMatrix.h"
0045 #include "TGeoShape.h"
0046 #include "TGeoTrd1.h"
0047 #include "TGeoTrd2.h"
0048 #include "TGeoTube.h"
0049
0050 using namespace Acts;
0051
0052 std::tuple<std::shared_ptr<const CylinderBounds>, const Transform3, double>
0053 ActsPlugins::TGeoSurfaceConverter::cylinderComponents(
0054     const TGeoShape& tgShape, const Double_t* rotation,
0055     const Double_t* translation, ActsPlugins::TGeoAxes axes,
0056     double scalor) noexcept(false) {
0057   auto a = axes.value();
0058   std::shared_ptr<const CylinderBounds> bounds = nullptr;
0059   Transform3 transform = Transform3::Identity();
0060   double thickness = 0.;
0061
0062   // Check if it's a tube (segment)
0063   auto tube = dynamic_cast<const TGeoTube*>(&tgShape);
0064   if (tube != nullptr) {
0065     if (!boost::istarts_with(a, "XY") && !boost::istarts_with(a, "YX")) {
0066       throw std::invalid_argument(
0067           "TGeoShape -> CylinderSurface (full): can only be converted with "
0068           "'(x/X)(y/Y)(*)' or '(y/Y)(x/X)(*) axes.");
0069     }
0070
0071     // The sign of the axes
0072     int xs = std::islower(a[0]) != 0 ? -1 : 1;
0073     int ys = std::islower(a[1]) != 0 ? -1 : 1;
0074
0075     // Create translation and rotation
0076     Vector3 t(scalor * translation[0], scalor * translation[1],
0077               scalor * translation[2]);
0078     bool flipxy = !boost::istarts_with(a, "X");
0079     Vector3 ax = flipxy ? xs * Vector3(rotation[1], rotation[4], rotation[7])
0080                         : xs * Vector3(rotation[0], rotation[3], rotation[6]);
0081     Vector3 ay = flipxy ? ys * Vector3(rotation[0], rotation[3], rotation[6])
0082                         : ys * Vector3(rotation[1], rotation[4], rotation[7]);
0083     Vector3 az = ax.cross(ay);
0084
0085     double minR = tube->GetRmin() * scalor;
0086     double maxR = tube->GetRmax() * scalor;
0087     double deltaR = maxR - minR;
0088     double medR = 0.5 * (minR + maxR);
0089     double halfZ = tube->GetDz() * scalor;
0090     if (halfZ > deltaR) {
0091       transform = TGeoPrimitivesHelper::makeTransform(ax, ay, az, t);
0092       double halfPhi = std::numbers::pi;
0093       double avgPhi = 0.;
0094       // Check if it's a segment
0095       auto tubeSeg = dynamic_cast<const TGeoTubeSeg*>(tube);
0096       if (tubeSeg != nullptr) {
0097         double phi1 = toRadian(tubeSeg->GetPhi1());
0098         double phi2 = toRadian(tubeSeg->GetPhi2());
0099         if (std::abs(phi2 - phi1) < std::numbers::pi * (1. - s_epsilon)) {
0100           if (!boost::starts_with(a, "X")) {
0101             throw std::invalid_argument(
0102                 "TGeoShape -> CylinderSurface (sectorial): can only be "
0103                 "converted "
0104                 "with "
0105                 "'(X)(y/Y)(*)' axes.");
0106           }
0107           halfPhi = 0.5 * (std::max(phi1, phi2) - std::min(phi1, phi2));
0108           avgPhi = 0.5 * (phi1 + phi2);
0109         }
0110       }
0111       bounds = std::make_shared<CylinderBounds>(medR, halfZ, halfPhi, avgPhi);
0112       thickness = deltaR;
0113     }
0114   }
0115   return {bounds, transform, thickness};
0116 }
0117
0118 std::tuple<std::shared_ptr<const DiscBounds>, const Transform3, double>
0119 ActsPlugins::TGeoSurfaceConverter::discComponents(
0120     const TGeoShape& tgShape, const Double_t* rotation,
0121     const Double_t* translation, ActsPlugins::TGeoAxes axes,
0122     double scalor) noexcept(false) {
0123   auto a = axes.value();
0124   using Line2D = Eigen::Hyperplane<double, 2>;
0125   std::shared_ptr<const DiscBounds> bounds = nullptr;
0126   Transform3 transform = Transform3::Identity();
0127
0128   double thickness = 0.;
0129   // Special test for composite shape of silicon
0130   auto compShape = dynamic_cast<const TGeoCompositeShape*>(&tgShape);
0131   if (compShape != nullptr) {
0132     if (!boost::istarts_with(a, "XY")) {
0133       throw std::invalid_argument(
0134           "TGeoShape -> DiscSurface (Annulus): can only be converted with "
0135           "'(x/X)(y/Y)(*)' "
0136           "axes");
0137     }
0138
0139     // Create translation and rotation
0140     Vector3 t(scalor * translation[0], scalor * translation[1],
0141               scalor * translation[2]);
0142     Vector3 ax(rotation[0], rotation[3], rotation[6]);
0143     Vector3 ay(rotation[1], rotation[4], rotation[7]);
0144     Vector3 az(rotation[2], rotation[5], rotation[8]);
0145
0146     transform = TGeoPrimitivesHelper::makeTransform(ax, ay, az, t);
0147
0148     auto interNode = dynamic_cast<TGeoIntersection*>(compShape->GetBoolNode());
0149     if (interNode != nullptr) {
0150       auto baseTube = dynamic_cast<TGeoTubeSeg*>(interNode->GetLeftShape());
0151       if (baseTube != nullptr) {
0152         double rMin = baseTube->GetRmin() * scalor;
0153         double rMax = baseTube->GetRmax() * scalor;
0154         auto maskShape = dynamic_cast<TGeoArb8*>(interNode->GetRightShape());
0155         if (maskShape != nullptr) {
0156           auto maskTransform = interNode->GetRightMatrix();
0157           // Get the only vertices
0158           const Double_t* polyVrt = maskShape->GetVertices();
0159           // Apply the whole transformation stored for the
0160           // polyhedron, since there is a translation and
0161           // also a side flip that needs to be applied.
0162           // @TODO check that 3rd coordinate is not altered by
0163           // the transformation ?
0164           std::vector<Vector2> vertices;
0165           for (unsigned int v = 0; v < 8; v += 2) {
0166             std::array<double, 3> local{polyVrt[v + 0], polyVrt[v + 1], 0.};
0167             std::array<double, 3> global{};
0168             maskTransform->LocalToMaster(local.data(), global.data());
0169             Vector2 vtx = Vector2(global[0] * scalor, global[1] * scalor);
0170             vertices.push_back(vtx);
0171           }
0172
0173           std::vector<std::pair<Vector2, Vector2>> boundLines;
0174           for (std::size_t i = 0; i < vertices.size(); ++i) {
0175             Vector2 va = vertices.at(i);
0176             Vector2 vb = vertices.at((i + 1) % vertices.size());
0177             Vector2 ab = vb - va;
0178             double phi = VectorHelpers::phi(ab);
0179
0180             if (std::abs(phi) > 3 * std::numbers::pi / 4. ||
0181                 std::abs(phi) < std::numbers::pi / 4.) {
0182               if (va.norm() < vb.norm()) {
0183                 boundLines.push_back(std::make_pair(va, vb));
0184               } else {
0185                 boundLines.push_back(std::make_pair(vb, va));
0186               }
0187             }
0188           }
0189
0190           if (boundLines.size() != 2) {
0191             throw std::logic_error(
0192                 "Input DiscPoly bounds type does not have sensible edges.");
0193           }
0194
0195           Line2D lA =
0196               Line2D::Through(boundLines[0].first, boundLines[0].second);
0197           Line2D lB =
0198               Line2D::Through(boundLines[1].first, boundLines[1].second);
0199           Vector2 ix = lA.intersection(lB);
0200
0201           const Eigen::Translation3d originTranslation(ix.x(), ix.y(), 0.);
0202           const Vector2 originShift = -ix;
0203
0204           // Update transform by prepending the origin shift translation
0205           transform = transform * originTranslation;
0206           // Transform phi line point to new origin and get phi
0207           double phi1 =
0208               VectorHelpers::phi(boundLines[0].second - boundLines[0].first);
0209           double phi2 =
0210               VectorHelpers::phi(boundLines[1].second - boundLines[1].first);
0211           double phiMax = std::max(phi1, phi2);
0212           double phiMin = std::min(phi1, phi2);
0213           double phiShift = 0.;
0214
0215           // Create the bounds
0216           auto annulusBounds = std::make_shared<const AnnulusBounds>(
0217               rMin, rMax, phiMin, phiMax, originShift, phiShift);
0218
0219           thickness = maskShape->GetDZ() * scalor;
0220           bounds = annulusBounds;
0221         }
0222       }
0223     }
0224   } else {
0225     // Check if it's a tube
0226     auto tube = dynamic_cast<const TGeoTube*>(&tgShape);
0227     if (tube != nullptr) {
0228       if (!boost::istarts_with(a, "XY") && !boost::istarts_with(a, "YX")) {
0229         throw std::invalid_argument(
0230             "TGeoShape -> DiscSurface: can only be converted with "
0231             "'(x/X)(y/Y)(*)' or '(y/Y)(x/X)(*) axes.");
0232       }
0233
0234       // The sign of the axes
0235       int xs = std::islower(a[0]) != 0 ? -1 : 1;
0236       int ys = std::islower(a[1]) != 0 ? -1 : 1;
0237
0238       // Create translation and rotation
0239       Vector3 t(scalor * translation[0], scalor * translation[1],
0240                 scalor * translation[2]);
0241       Vector3 ax = xs * Vector3(rotation[0], rotation[3], rotation[6]);
0242       Vector3 ay = ys * Vector3(rotation[1], rotation[4], rotation[7]);
0243       Vector3 az = ax.cross(ay);
0244       transform = TGeoPrimitivesHelper::makeTransform(ax, ay, az, t);
0245
0246       double minR = tube->GetRmin() * scalor;
0247       double maxR = tube->GetRmax() * scalor;
0248       double halfZ = tube->GetDz() * scalor;
0249       double halfPhi = std::numbers::pi;
0250       double avgPhi = 0.;
0251       // Check if it's a segment
0252       auto tubeSeg = dynamic_cast<const TGeoTubeSeg*>(tube);
0253       if (tubeSeg != nullptr) {
0254         double phi1 = toRadian(tubeSeg->GetPhi1());
0255         double phi2 = toRadian(tubeSeg->GetPhi2());
0256         if (std::abs(phi2 - phi1) < 2 * std::numbers::pi * (1. - s_epsilon)) {
0257           if (!boost::starts_with(a, "X")) {
0258             throw std::invalid_argument(
0259                 "TGeoShape -> CylinderSurface (sectorial): can only be "
0260                 "converted "
0261                 "with "
0262                 "'(X)(y/Y)(*)' axes.");
0263           }
0264           halfPhi = 0.5 * (std::max(phi1, phi2) - std::min(phi1, phi2));
0265           avgPhi = 0.5 * (phi1 + phi2);
0266         }
0267       }
0268       bounds = std::make_shared<RadialBounds>(minR, maxR, halfPhi, avgPhi);
0269       thickness = 2 * halfZ;
0270     }
0271   }
0272   return {bounds, transform, thickness};
0273 }
0274
0275 std::tuple<std::shared_ptr<const PlanarBounds>, const Transform3, double>
0276 ActsPlugins::TGeoSurfaceConverter::planeComponents(
0277     const TGeoShape& tgShape, const Double_t* rotation,
0278     const Double_t* translation, TGeoAxes axes, double scalor) noexcept(false) {
0279   auto a = axes.value();
0280   // Create translation and rotation
0281   Vector3 t(scalor * translation[0], scalor * translation[1],
0282             scalor * translation[2]);
0283   Vector3 ax(rotation[0], rotation[3], rotation[6]);
0284   Vector3 ay(rotation[1], rotation[4], rotation[7]);
0285   Vector3 az(rotation[2], rotation[5], rotation[8]);
0286
0287   std::shared_ptr<const PlanarBounds> bounds = nullptr;
0288
0289   // Check if it's a box - always true, hence last ressort
0290   auto box = dynamic_cast<const TGeoBBox*>(&tgShape);
0291
0292   // Check if it's a trapezoid2
0293   auto trapezoid1 = dynamic_cast<const TGeoTrd1*>(&tgShape);
0294   if ((trapezoid1 != nullptr) && !boost::istarts_with(a, "XZ")) {
0295     throw std::invalid_argument(
0296         "TGeoTrd1 -> PlaneSurface: can only be converted with '(x/X)(z/Z)(*)' "
0297         "axes");
0298   }
0299
0300   // Check if it's a trapezoid2
0301   auto trapezoid2 = dynamic_cast<const TGeoTrd2*>(&tgShape);
0302   if (trapezoid2 != nullptr) {
0303     if (!boost::istarts_with(a, "X") &&
0304         std::abs(trapezoid2->GetDx1() - trapezoid2->GetDx2()) > s_epsilon) {
0305       throw std::invalid_argument(
0306           "TGeoTrd2 -> PlaneSurface: dx1 must be be equal to dx2 if not taken "
0307           "as trapezoidal side.");
0308     } else if (!boost::istarts_with(a, "Y") &&
0309                std::abs(trapezoid2->GetDy1() - trapezoid2->GetDy2()) >
0310                    s_epsilon) {
0311       throw std::invalid_argument(
0312           "TGeoTrd2 -> PlaneSurface: dy1 must be be equal to dy2 if not taken "
0313           "as trapezoidal side.");
0314     }
0315     // Not allowed
0316     if (boost::istarts_with(a, "XY") || boost::istarts_with(a, "YX")) {
0317       throw std::invalid_argument(
0318           "TGeoTrd2 -> PlaneSurface: only works with (x/X)(z/Z) and "
0319           "(y/Y)(z/Z).");
0320     }
0321   }
0322
0323   // Check if it's a Arb8
0324   auto polygon8c = dynamic_cast<const TGeoArb8*>(&tgShape);
0325   TGeoArb8* polygon8 = nullptr;
0326   if (polygon8c != nullptr) {
0327     // Needed otherwise you can access GetVertices
0328     polygon8 = const_cast<TGeoArb8*>(polygon8c);
0329   }
0330
0331   if ((polygon8c != nullptr) &&
0332       !(boost::istarts_with(a, "XY") || boost::istarts_with(a, "YX"))) {
0333     throw std::invalid_argument(
0334         "TGeoArb8 -> PlaneSurface: dz must be normal component of Surface.");
0335   }
0336
0337   // The thickness will be filled
0338   double thickness = 0.;
0339
0340   // The sign of the axes
0341   int xs = std::islower(a[0]) != 0 ? -1 : 1;
0342   int ys = std::islower(a[1]) != 0 ? -1 : 1;
0343
0344   // Set up the columns : only cyclic iterations are allowed
0345   Vector3 cx = xs * ax;
0346   Vector3 cy = ys * ay;
0347   if (boost::istarts_with(a, "XY")) {
0348     if (trapezoid2 != nullptr) {
0349       double dx1 = (ys < 0) ? trapezoid1->GetDx2() : trapezoid1->GetDx1();
0350       double dx2 = (ys < 0) ? trapezoid1->GetDx1() : trapezoid1->GetDx2();
0351       bounds = std::make_shared<const TrapezoidBounds>(
0352           scalor * dx1, scalor * dx2, scalor * trapezoid2->GetDy1());
0353       thickness = 2 * scalor * trapezoid2->GetDz();
0354     } else if (polygon8 != nullptr) {
0355       Double_t* tgverts = polygon8->GetVertices();
0356       std::vector<Vector2> pVertices;
0357       for (unsigned int ivtx = 0; ivtx < 4; ++ivtx) {
0358         pVertices.push_back(Vector2(scalor * xs * tgverts[ivtx * 2],
0359                                     scalor * ys * tgverts[ivtx * 2 + 1]));
0360       }
0361       bounds = std::make_shared<ConvexPolygonBounds<4>>(pVertices);
0362       thickness = 2 * scalor * polygon8->GetDz();
0363     } else if (box != nullptr) {
0364       bounds = std::make_shared<const RectangleBounds>(scalor * box->GetDX(),
0365                                                        scalor * box->GetDY());
0366       thickness = 2 * scalor * box->GetDZ();
0367     }
0368   } else if (boost::istarts_with(a, "YZ")) {
0369     cx = xs * ay;
0370     cy = ys * az;
0371     if (trapezoid1 != nullptr) {
0372       throw std::invalid_argument(
0373           "TGeoTrd1 can only be converted with '(x/X)(z/Z)(y/Y)' axes");
0374     } else if (trapezoid2 != nullptr) {
0375       double dx1 = (ys < 0) ? trapezoid2->GetDy2() : trapezoid2->GetDy1();
0376       double dx2 = (ys < 0) ? trapezoid2->GetDy1() : trapezoid2->GetDy2();
0377       bounds = std::make_shared<const TrapezoidBounds>(
0378           scalor * dx1, scalor * dx2, scalor * trapezoid2->GetDz());
0379       thickness = 2 * scalor * trapezoid2->GetDx1();
0380     } else if (box != nullptr) {
0381       bounds = std::make_shared<const RectangleBounds>(scalor * box->GetDY(),
0382                                                        scalor * box->GetDZ());
0383       thickness = 2 * scalor * box->GetDX();
0384     }
0385   } else if (boost::istarts_with(a, "ZX")) {
0386     cx = xs * az;
0387     cy = ys * ax;
0388     if (box != nullptr) {
0389       bounds = std::make_shared<const RectangleBounds>(scalor * box->GetDZ(),
0390                                                        scalor * box->GetDX());
0391       thickness = 2 * scalor * box->GetDY();
0392     }
0393   } else if (boost::istarts_with(a, "XZ")) {
0394     cx = xs * ax;
0395     cy = ys * az;
0396     if (trapezoid1 != nullptr) {
0397       double dx1 = (ys < 0) ? trapezoid1->GetDx2() : trapezoid1->GetDx1();
0398       double dx2 = (ys < 0) ? trapezoid1->GetDx1() : trapezoid1->GetDx2();
0399       bounds = std::make_shared<const TrapezoidBounds>(
0400           scalor * dx1, scalor * dx2, scalor * trapezoid1->GetDz());
0401       thickness = 2 * scalor * trapezoid1->GetDy();
0402     } else if (trapezoid2 != nullptr) {
0403       double dx1 = (ys < 0) ? trapezoid2->GetDx2() : trapezoid2->GetDx1();
0404       double dx2 = (ys < 0) ? trapezoid2->GetDx1() : trapezoid2->GetDx2();
0405       bounds = std::make_shared<const TrapezoidBounds>(
0406           scalor * dx1, scalor * dx2, scalor * trapezoid2->GetDz());
0407       thickness = 2 * scalor * trapezoid2->GetDy1();
0408     } else if (box != nullptr) {
0409       bounds = std::make_shared<const RectangleBounds>(scalor * box->GetDX(),
0410                                                        scalor * box->GetDZ());
0411       thickness = 2 * scalor * box->GetDY();
0412     }
0413   } else if (boost::istarts_with(a, "YX")) {
0414     cx = xs * ay;
0415     cy = ys * ax;
0416     if (trapezoid2 != nullptr) {
0417       double dx1 = (ys < 0) ? trapezoid2->GetDy2() : trapezoid2->GetDy1();
0418       double dx2 = (ys < 0) ? trapezoid2->GetDy1() : trapezoid2->GetDy2();
0419       bounds = std::make_shared<const TrapezoidBounds>(
0420           scalor * dx1, scalor * dx2, scalor * trapezoid2->GetDx1());
0421       thickness = 2 * scalor * trapezoid2->GetDz();
0422     } else if (polygon8 != nullptr) {
0423       const Double_t* tgverts = polygon8->GetVertices();
0424       std::vector<Vector2> pVertices;
0425       for (unsigned int ivtx = 0; ivtx < 4; ++ivtx) {
0426         pVertices.push_back(Vector2(scalor * xs * tgverts[ivtx * 2 + 1],
0427                                     scalor * ys * tgverts[ivtx * 2]));
0428       }
0429       bounds = std::make_shared<ConvexPolygonBounds<4>>(pVertices);
0430       thickness = 2 * scalor * polygon8->GetDz();
0431     } else if (box != nullptr) {
0432       bounds = std::make_shared<const RectangleBounds>(scalor * box->GetDY(),
0433                                                        scalor * box->GetDX());
0434       thickness = 2 * scalor * box->GetDZ();
0435     }
0436   } else if (boost::istarts_with(a, "ZY")) {
0437     cx = xs * az;
0438     cy = ys * ay;
0439     if (box != nullptr) {
0440       bounds = std::make_shared<const RectangleBounds>(scalor * box->GetDZ(),
0441                                                        scalor * box->GetDY());
0442       thickness = 2 * scalor * box->GetDX();
0443     }
0444   } else {
0445     throw std::invalid_argument(
0446         "TGeoConverter: axes definition must be permutation of "
0447         "'(x/X)(y/Y)(z/Z)'");
0448   }
0449
0450   // Create the normal vector & the transform
0451   auto cz = cx.cross(cy);
0452   auto transform = TGeoPrimitivesHelper::makeTransform(cx, cy, cz, t);
0453
0454   return {bounds, transform, thickness};
0455 }
0456
0457 std::tuple<std::shared_ptr<Surface>, double>
0458 ActsPlugins::TGeoSurfaceConverter::toSurface(const TGeoShape& tgShape,
0459                                              const TGeoMatrix& tgMatrix,
0460                                              TGeoAxes axes,
0461                                              double scalor) noexcept(false) {
0462   // Get the placement and orientation in respect to its mother
0463   const Double_t* rotation = tgMatrix.GetRotationMatrix();
0464   const Double_t* translation = tgMatrix.GetTranslation();
0465
0466   auto [cBounds, cTransform, cThickness] =
0467       cylinderComponents(tgShape, rotation, translation, axes, scalor);
0468   if (cBounds != nullptr) {
0469     return {Surface::makeShared<CylinderSurface>(cTransform, cBounds),
0470             cThickness};
0471   }
0472
0473   auto [dBounds, dTransform, dThickness] =
0474       discComponents(tgShape, rotation, translation, axes, scalor);
0475   if (dBounds != nullptr) {
0476     return {Surface::makeShared<DiscSurface>(dTransform, dBounds), dThickness};
0477   }
0478
0479   auto [pBounds, pTransform, pThickness] =
0480       planeComponents(tgShape, rotation, translation, axes, scalor);
0481   if (pBounds != nullptr) {
0482     return {Surface::makeShared<PlaneSurface>(pTransform, pBounds), pThickness};
0483   }
0484
0485   return {nullptr, 0.};
0486 }