TGeo/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 "Acts/Plugins/TGeo/TGeoSurfaceConverter.hpp"
0010
0011 #include "Acts/Definitions/Tolerance.hpp"
0012 #include "Acts/Plugins/TGeo/TGeoPrimitivesHelper.hpp"
0013 #include "Acts/Surfaces/AnnulusBounds.hpp"
0014 #include "Acts/Surfaces/ConvexPolygonBounds.hpp"
0015 #include "Acts/Surfaces/CylinderBounds.hpp"
0016 #include "Acts/Surfaces/CylinderSurface.hpp"
0017 #include "Acts/Surfaces/DiscSurface.hpp"
0018 #include "Acts/Surfaces/PlaneSurface.hpp"
0019 #include "Acts/Surfaces/RadialBounds.hpp"
0020 #include "Acts/Surfaces/RectangleBounds.hpp"
0021 #include "Acts/Surfaces/Surface.hpp"
0022 #include "Acts/Surfaces/TrapezoidBounds.hpp"
0023 #include "Acts/Utilities/Helpers.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 std::tuple<std::shared_ptr<const Acts::CylinderBounds>, const Acts::Transform3,
0051            double>
0052 Acts::TGeoSurfaceConverter::cylinderComponents(const TGeoShape& tgShape,
0053                                                const Double_t* rotation,
0054                                                const Double_t* translation,
0055                                                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 Acts::DiscBounds>, const Acts::Transform3,
0118            double>
0119 Acts::TGeoSurfaceConverter::discComponents(const TGeoShape& tgShape,
0120                                            const Double_t* rotation,
0121                                            const Double_t* translation,
0122                                            const std::string& axes,
0123                                            double scalor) noexcept(false) {
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(axes, "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 a = vertices.at(i);
0176             Vector2 b = vertices.at((i + 1) % vertices.size());
0177             Vector2 ab = b - a;
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 (a.norm() < b.norm()) {
0183                 boundLines.push_back(std::make_pair(a, b));
0184               } else {
0185                 boundLines.push_back(std::make_pair(b, a));
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(axes, "XY") &&
0229           !boost::istarts_with(axes, "YX")) {
0230         throw std::invalid_argument(
0231             "TGeoShape -> DiscSurface: can only be converted with "
0232             "'(x/X)(y/Y)(*)' or '(y/Y)(x/X)(*) axes.");
0233       }
0234
0235       // The sign of the axes
0236       int xs = std::islower(axes.at(0)) != 0 ? -1 : 1;
0237       int ys = std::islower(axes.at(1)) != 0 ? -1 : 1;
0238
0239       // Create translation and rotation
0240       Vector3 t(scalor * translation[0], scalor * translation[1],
0241                 scalor * translation[2]);
0242       Vector3 ax = xs * Vector3(rotation[0], rotation[3], rotation[6]);
0243       Vector3 ay = ys * Vector3(rotation[1], rotation[4], rotation[7]);
0244       Vector3 az = ax.cross(ay);
0245       transform = TGeoPrimitivesHelper::makeTransform(ax, ay, az, t);
0246
0247       double minR = tube->GetRmin() * scalor;
0248       double maxR = tube->GetRmax() * scalor;
0249       double halfZ = tube->GetDz() * scalor;
0250       double halfPhi = std::numbers::pi;
0251       double avgPhi = 0.;
0252       // Check if it's a segment
0253       auto tubeSeg = dynamic_cast<const TGeoTubeSeg*>(tube);
0254       if (tubeSeg != nullptr) {
0255         double phi1 = toRadian(tubeSeg->GetPhi1());
0256         double phi2 = toRadian(tubeSeg->GetPhi2());
0257         if (std::abs(phi2 - phi1) < 2 * std::numbers::pi * (1. - s_epsilon)) {
0258           if (!boost::starts_with(axes, "X")) {
0259             throw std::invalid_argument(
0260                 "TGeoShape -> CylinderSurface (sectorial): can only be "
0261                 "converted "
0262                 "with "
0263                 "'(X)(y/Y)(*)' axes.");
0264           }
0265           halfPhi = 0.5 * (std::max(phi1, phi2) - std::min(phi1, phi2));
0266           avgPhi = 0.5 * (phi1 + phi2);
0267         }
0268       }
0269       bounds = std::make_shared<RadialBounds>(minR, maxR, halfPhi, avgPhi);
0270       thickness = 2 * halfZ;
0271     }
0272   }
0273   return {bounds, transform, thickness};
0274 }
0275
0276 std::tuple<std::shared_ptr<const Acts::PlanarBounds>, const Acts::Transform3,
0277            double>
0278 Acts::TGeoSurfaceConverter::planeComponents(const TGeoShape& tgShape,
0279                                             const Double_t* rotation,
0280                                             const Double_t* translation,
0281                                             const std::string& axes,
0282                                             double scalor) noexcept(false) {
0283   // Create translation and rotation
0284   Vector3 t(scalor * translation[0], scalor * translation[1],
0285             scalor * translation[2]);
0286   Vector3 ax(rotation[0], rotation[3], rotation[6]);
0287   Vector3 ay(rotation[1], rotation[4], rotation[7]);
0288   Vector3 az(rotation[2], rotation[5], rotation[8]);
0289
0290   std::shared_ptr<const PlanarBounds> bounds = nullptr;
0291
0292   // Check if it's a box - always true, hence last ressort
0293   const TGeoBBox* box = dynamic_cast<const TGeoBBox*>(&tgShape);
0294
0295   // Check if it's a trapezoid2
0296   const TGeoTrd1* trapezoid1 = dynamic_cast<const TGeoTrd1*>(&tgShape);
0297   if ((trapezoid1 != nullptr) && !boost::istarts_with(axes, "XZ")) {
0298     throw std::invalid_argument(
0299         "TGeoTrd1 -> PlaneSurface: can only be converted with '(x/X)(z/Z)(*)' "
0300         "axes");
0301   }
0302
0303   // Check if it's a trapezoid2
0304   const TGeoTrd2* trapezoid2 = dynamic_cast<const TGeoTrd2*>(&tgShape);
0305   if (trapezoid2 != nullptr) {
0306     if (!boost::istarts_with(axes, "X") &&
0307         std::abs(trapezoid2->GetDx1() - trapezoid2->GetDx2()) > s_epsilon) {
0308       throw std::invalid_argument(
0309           "TGeoTrd2 -> PlaneSurface: dx1 must be be equal to dx2 if not taken "
0310           "as trapezoidal side.");
0311     } else if (!boost::istarts_with(axes, "Y") &&
0312                std::abs(trapezoid2->GetDy1() - trapezoid2->GetDy2()) >
0313                    s_epsilon) {
0314       throw std::invalid_argument(
0315           "TGeoTrd2 -> PlaneSurface: dy1 must be be equal to dy2 if not taken "
0316           "as trapezoidal side.");
0317     }
0318     // Not allowed
0319     if (boost::istarts_with(axes, "XY") || boost::istarts_with(axes, "YX")) {
0320       throw std::invalid_argument(
0321           "TGeoTrd2 -> PlaneSurface: only works with (x/X)(z/Z) and "
0322           "(y/Y)(z/Z).");
0323     }
0324   }
0325
0326   // Check if it's a Arb8
0327   const TGeoArb8* polygon8c = dynamic_cast<const TGeoArb8*>(&tgShape);
0328   TGeoArb8* polygon8 = nullptr;
0329   if (polygon8c != nullptr) {
0330     // Needed otherwise you can access GetVertices
0331     polygon8 = const_cast<TGeoArb8*>(polygon8c);
0332   }
0333
0334   if ((polygon8c != nullptr) &&
0335       !(boost::istarts_with(axes, "XY") || boost::istarts_with(axes, "YX"))) {
0336     throw std::invalid_argument(
0337         "TGeoArb8 -> PlaneSurface: dz must be normal component of Surface.");
0338   }
0339
0340   // The thickness will be filled
0341   double thickness = 0.;
0342
0343   // The sign of the axes
0344   int xs = std::islower(axes.at(0)) != 0 ? -1 : 1;
0345   int ys = std::islower(axes.at(1)) != 0 ? -1 : 1;
0346
0347   // Set up the columns : only cyclic iterations are allowed
0348   Vector3 cx = xs * ax;
0349   Vector3 cy = ys * ay;
0350   if (boost::istarts_with(axes, "XY")) {
0351     if (trapezoid2 != nullptr) {
0352       double dx1 = (ys < 0) ? trapezoid1->GetDx2() : trapezoid1->GetDx1();
0353       double dx2 = (ys < 0) ? trapezoid1->GetDx1() : trapezoid1->GetDx2();
0354       bounds = std::make_shared<const TrapezoidBounds>(
0355           scalor * dx1, scalor * dx2, scalor * trapezoid2->GetDy1());
0356       thickness = 2 * scalor * trapezoid2->GetDz();
0357     } else if (polygon8 != nullptr) {
0358       Double_t* tgverts = polygon8->GetVertices();
0359       std::vector<Vector2> pVertices;
0360       for (unsigned int ivtx = 0; ivtx < 4; ++ivtx) {
0361         pVertices.push_back(Vector2(scalor * xs * tgverts[ivtx * 2],
0362                                     scalor * ys * tgverts[ivtx * 2 + 1]));
0363       }
0364       bounds = std::make_shared<ConvexPolygonBounds<4>>(pVertices);
0365       thickness = 2 * scalor * polygon8->GetDz();
0366     } else if (box != nullptr) {
0367       bounds = std::make_shared<const RectangleBounds>(scalor * box->GetDX(),
0368                                                        scalor * box->GetDY());
0369       thickness = 2 * scalor * box->GetDZ();
0370     }
0371   } else if (boost::istarts_with(axes, "YZ")) {
0372     cx = xs * ay;
0373     cy = ys * az;
0374     if (trapezoid1 != nullptr) {
0375       throw std::invalid_argument(
0376           "TGeoTrd1 can only be converted with '(x/X)(z/Z)(y/Y)' axes");
0377     } else if (trapezoid2 != nullptr) {
0378       double dx1 = (ys < 0) ? trapezoid2->GetDy2() : trapezoid2->GetDy1();
0379       double dx2 = (ys < 0) ? trapezoid2->GetDy1() : trapezoid2->GetDy2();
0380       bounds = std::make_shared<const TrapezoidBounds>(
0381           scalor * dx1, scalor * dx2, scalor * trapezoid2->GetDz());
0382       thickness = 2 * scalor * trapezoid2->GetDx1();
0383     } else if (box != nullptr) {
0384       bounds = std::make_shared<const RectangleBounds>(scalor * box->GetDY(),
0385                                                        scalor * box->GetDZ());
0386       thickness = 2 * scalor * box->GetDX();
0387     }
0388   } else if (boost::istarts_with(axes, "ZX")) {
0389     cx = xs * az;
0390     cy = ys * ax;
0391     if (box != nullptr) {
0392       bounds = std::make_shared<const RectangleBounds>(scalor * box->GetDZ(),
0393                                                        scalor * box->GetDX());
0394       thickness = 2 * scalor * box->GetDY();
0395     }
0396   } else if (boost::istarts_with(axes, "XZ")) {
0397     cx = xs * ax;
0398     cy = ys * az;
0399     if (trapezoid1 != nullptr) {
0400       double dx1 = (ys < 0) ? trapezoid1->GetDx2() : trapezoid1->GetDx1();
0401       double dx2 = (ys < 0) ? trapezoid1->GetDx1() : trapezoid1->GetDx2();
0402       bounds = std::make_shared<const TrapezoidBounds>(
0403           scalor * dx1, scalor * dx2, scalor * trapezoid1->GetDz());
0404       thickness = 2 * scalor * trapezoid1->GetDy();
0405     } else if (trapezoid2 != nullptr) {
0406       double dx1 = (ys < 0) ? trapezoid2->GetDx2() : trapezoid2->GetDx1();
0407       double dx2 = (ys < 0) ? trapezoid2->GetDx1() : trapezoid2->GetDx2();
0408       bounds = std::make_shared<const TrapezoidBounds>(
0409           scalor * dx1, scalor * dx2, scalor * trapezoid2->GetDz());
0410       thickness = 2 * scalor * trapezoid2->GetDy1();
0411     } else if (box != nullptr) {
0412       bounds = std::make_shared<const RectangleBounds>(scalor * box->GetDX(),
0413                                                        scalor * box->GetDZ());
0414       thickness = 2 * scalor * box->GetDY();
0415     }
0416   } else if (boost::istarts_with(axes, "YX")) {
0417     cx = xs * ay;
0418     cy = ys * ax;
0419     if (trapezoid2 != nullptr) {
0420       double dx1 = (ys < 0) ? trapezoid2->GetDy2() : trapezoid2->GetDy1();
0421       double dx2 = (ys < 0) ? trapezoid2->GetDy1() : trapezoid2->GetDy2();
0422       bounds = std::make_shared<const TrapezoidBounds>(
0423           scalor * dx1, scalor * dx2, scalor * trapezoid2->GetDx1());
0424       thickness = 2 * scalor * trapezoid2->GetDz();
0425     } else if (polygon8 != nullptr) {
0426       const Double_t* tgverts = polygon8->GetVertices();
0427       std::vector<Vector2> pVertices;
0428       for (unsigned int ivtx = 0; ivtx < 4; ++ivtx) {
0429         pVertices.push_back(Vector2(scalor * xs * tgverts[ivtx * 2 + 1],
0430                                     scalor * ys * tgverts[ivtx * 2]));
0431       }
0432       bounds = std::make_shared<ConvexPolygonBounds<4>>(pVertices);
0433       thickness = 2 * scalor * polygon8->GetDz();
0434     } else if (box != nullptr) {
0435       bounds = std::make_shared<const RectangleBounds>(scalor * box->GetDY(),
0436                                                        scalor * box->GetDX());
0437       thickness = 2 * scalor * box->GetDZ();
0438     }
0439   } else if (boost::istarts_with(axes, "ZY")) {
0440     cx = xs * az;
0441     cy = ys * ay;
0442     if (box != nullptr) {
0443       bounds = std::make_shared<const RectangleBounds>(scalor * box->GetDZ(),
0444                                                        scalor * box->GetDY());
0445       thickness = 2 * scalor * box->GetDX();
0446     }
0447   } else {
0448     throw std::invalid_argument(
0449         "TGeoConverter: axes definition must be permutation of "
0450         "'(x/X)(y/Y)(z/Z)'");
0451   }
0452
0453   // Create the normal vector & the transform
0454   auto cz = cx.cross(cy);
0455   auto transform = TGeoPrimitivesHelper::makeTransform(cx, cy, cz, t);
0456
0457   return {bounds, transform, thickness};
0458 }
0459
0460 std::tuple<std::shared_ptr<Acts::Surface>, double>
0461 Acts::TGeoSurfaceConverter::toSurface(const TGeoShape& tgShape,
0462                                       const TGeoMatrix& tgMatrix,
0463                                       const std::string& axes,
0464                                       double scalor) noexcept(false) {
0465   // Get the placement and orientation in respect to its mother
0466   const Double_t* rotation = tgMatrix.GetRotationMatrix();
0467   const Double_t* translation = tgMatrix.GetTranslation();
0468
0469   auto [cBounds, cTransform, cThickness] =
0470       cylinderComponents(tgShape, rotation, translation, axes, scalor);
0471   if (cBounds != nullptr) {
0472     return {Surface::makeShared<CylinderSurface>(cTransform, cBounds),
0473             cThickness};
0474   }
0475
0476   auto [dBounds, dTransform, dThickness] =
0477       discComponents(tgShape, rotation, translation, axes, scalor);
0478   if (dBounds != nullptr) {
0479     return {Surface::makeShared<DiscSurface>(dTransform, dBounds), dThickness};
0480   }
0481
0482   auto [pBounds, pTransform, pThickness] =
0483       planeComponents(tgShape, rotation, translation, axes, scalor);
0484   if (pBounds != nullptr) {
0485     return {Surface::makeShared<PlaneSurface>(pTransform, pBounds), pThickness};
0486   }
0487
0488   return {nullptr, 0.};
0489 }