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