Plugins/Root/TGeoTubeConversionTests.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 <boost/test/unit_test.hpp>
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Tolerance.hpp"
0013 #include "Acts/Geometry/GeometryContext.hpp"
0014 #include "Acts/Surfaces/CylinderBounds.hpp"
0015 #include "Acts/Surfaces/RadialBounds.hpp"
0016 #include "Acts/Surfaces/Surface.hpp"
0017 #include "Acts/Surfaces/SurfaceBounds.hpp"
0018 #include "Acts/Visualization/GeometryView3D.hpp"
0019 #include "Acts/Visualization/ObjVisualization3D.hpp"
0020 #include "Acts/Visualization/ViewConfig.hpp"
0021 #include "ActsPlugins/Root/TGeoSurfaceConverter.hpp"
0022 #include "ActsTests/CommonHelpers/FloatComparisons.hpp"
0023
0024 #include <cmath>
0025 #include <cstddef>
0026 #include <memory>
0027 #include <numbers>
0028 #include <stdexcept>
0029 #include <string>
0030 #include <vector>
0031
0032 #include "TGeoManager.h"
0033 #include "TGeoMaterial.h"
0034 #include "TGeoMatrix.h"
0035 #include "TGeoMedium.h"
0036 #include "TGeoTube.h"
0037 #include "TGeoVolume.h"
0038 #include "TView.h"
0039
0040 using namespace Acts;
0041 using namespace ActsPlugins;
0042
0043 namespace ActsTests {
0044
0045 GeometryContext tgContext = GeometryContext();
0046
0047 ViewConfig red{.color = {200, 0, 0}};
0048 ViewConfig green{.color = {0, 200, 0}};
0049 ViewConfig blue{.color = {0, 0, 200}};
0050
0051 std::vector<std::string> allowedAxes = {"XY*", "Xy*", "xy*", "xY*",
0052                                         "YX*", "yx*", "yX*", "Yx*"};
0053
0054 std::vector<std::string> notAllowedAxes = {"YZ*", "ZX*", "ZY*"};
0055
0056 BOOST_AUTO_TEST_SUITE(RootSuite)
0057
0058 /// @brief Unit test to convert a TGeoTube into a CylinderSurface
0059 ///
0060 /// * The TGeoTrd1 can only have (x/X)(y/Y) orientation
0061 BOOST_AUTO_TEST_CASE(TGeoTube_to_CylinderSurface) {
0062   ObjVisualization3D objVis;
0063
0064   double rmin = 10.;
0065   double rmax = 11;
0066   double hz = 40.;
0067   double phimin = 45.;
0068   double phimax = -45.;
0069
0070   new TGeoManager("trd1", "poza9");
0071   TGeoMaterial *mat = new TGeoMaterial("Al", 26.98, 13, 2.7);
0072   TGeoMedium *med = new TGeoMedium("MED", 1, mat);
0073   TGeoVolume *top = gGeoManager->MakeBox("TOP", med, 100, 100, 100);
0074   gGeoManager->SetTopVolume(top);
0075   TGeoVolume *vol = gGeoManager->MakeTube("Tube", med, rmin, rmax, hz);
0076   TGeoVolume *vols =
0077       gGeoManager->MakeTubs("Tube", med, rmin, rmax, hz, phimin, phimax);
0078   gGeoManager->CloseGeometry();
0079
0080   std::size_t icyl = 0;
0081   for (const auto &axes : allowedAxes) {
0082     auto [cylinder, thickness] = TGeoSurfaceConverter::toSurface(
0083         *vol->GetShape(), *gGeoIdentity, axes, 1);
0084     BOOST_REQUIRE_NE(cylinder, nullptr);
0085     BOOST_CHECK_EQUAL(cylinder->type(), Surface::Cylinder);
0086     CHECK_CLOSE_ABS(thickness, rmax - rmin, s_epsilon);
0087
0088     auto bounds = dynamic_cast<const CylinderBounds *>(&(cylinder->bounds()));
0089     BOOST_REQUIRE_NE(bounds, nullptr);
0090     double bR = bounds->get(CylinderBounds::eR);
0091     double bhZ = bounds->get(CylinderBounds::eHalfLengthZ);
0092
0093     CHECK_CLOSE_ABS(bR, 10.5, s_epsilon);
0094     CHECK_CLOSE_ABS(bhZ, hz, s_epsilon);
0095
0096     auto transform = cylinder->transform(tgContext);
0097     auto rotation = transform.rotation();
0098
0099     // Check if the surface is the (negative) identity
0100     GeometryView3D::drawSurface(objVis, *cylinder, tgContext);
0101     const Vector3 center = cylinder->center(tgContext);
0102     GeometryView3D::drawArrowForward(
0103         objVis, center, center + 1.2 * bR * rotation.col(0), 4., 2.5, red);
0104     GeometryView3D::drawArrowForward(
0105         objVis, center, center + 1.2 * bR * rotation.col(1), 4., 2.5, green);
0106     GeometryView3D::drawArrowForward(
0107         objVis, center, center + 1.2 * bhZ * rotation.col(2), 4., 2.5, blue);
0108     objVis.write("TGeoConversion_TGeoTube_CylinderSurface_" +
0109                  std::to_string(icyl));
0110     objVis.clear();
0111
0112     if (icyl < 2) {
0113       auto [cylinderSegment, cThickness] = TGeoSurfaceConverter::toSurface(
0114           *vols->GetShape(), *gGeoIdentity, axes, 1);
0115       BOOST_REQUIRE_NE(cylinderSegment, nullptr);
0116       BOOST_CHECK_EQUAL(cylinderSegment->type(), Surface::Cylinder);
0117       CHECK_CLOSE_ABS(cThickness, rmax - rmin, s_epsilon);
0118
0119       auto boundsSegment =
0120           dynamic_cast<const CylinderBounds *>(&(cylinderSegment->bounds()));
0121       BOOST_REQUIRE_NE(boundsSegment, nullptr);
0122       bR = boundsSegment->get(CylinderBounds::eR);
0123       bhZ = boundsSegment->get(CylinderBounds::eHalfLengthZ);
0124       double hphi = boundsSegment->get(CylinderBounds::eHalfPhiSector);
0125       double mphi = boundsSegment->get(CylinderBounds::eAveragePhi);
0126       CHECK_CLOSE_ABS(bR, 10.5, s_epsilon);
0127       CHECK_CLOSE_ABS(bhZ, hz, s_epsilon);
0128       CHECK_CLOSE_ABS(hphi, std::numbers::pi / 4., s_epsilon);
0129       CHECK_CLOSE_ABS(mphi, 0., s_epsilon);
0130       GeometryView3D::drawSurface(objVis, *cylinderSegment, tgContext);
0131       GeometryView3D::drawArrowForward(
0132           objVis, center, center + 1.2 * bR * rotation.col(0), 4., 2.5, red);
0133       GeometryView3D::drawArrowForward(
0134           objVis, center, center + 1.2 * bR * rotation.col(1), 4., 2.5, green);
0135       GeometryView3D::drawArrowForward(
0136           objVis, center, center + 1.2 * bhZ * rotation.col(2), 4., 2.5, blue);
0137       objVis.write("TGeoConversion_TGeoTube_CylinderSegmentSurface_" +
0138                    std::to_string(icyl));
0139       objVis.clear();
0140     } else {
0141       BOOST_CHECK_THROW(TGeoSurfaceConverter::toSurface(*vols->GetShape(),
0142                                                         *gGeoIdentity, axes, 1),
0143                         std::invalid_argument);
0144     }
0145     ++icyl;
0146   }
0147
0148   // Check exceptions for not allowed axis definition
0149   for (const auto &naxes : notAllowedAxes) {
0150     BOOST_CHECK_THROW(TGeoSurfaceConverter::toSurface(*vol->GetShape(),
0151                                                       *gGeoIdentity, naxes, 1),
0152                       std::invalid_argument);
0153   }
0154 }
0155
0156 /// @brief Unit test to convert a TGeoTube into a DiscSurface
0157 ///
0158 /// * The TGeoTrd1 can only have (x/X)(y/Y) orientation
0159 BOOST_AUTO_TEST_CASE(TGeoTube_to_DiscSurface) {
0160   ObjVisualization3D objVis;
0161
0162   double rmin = 5.;
0163   double rmax = 25;
0164   double hz = 2.;
0165   double phimin = 45.;
0166   double phimax = -45.;
0167
0168   new TGeoManager("trd1", "poza9");
0169   TGeoMaterial *mat = new TGeoMaterial("Al", 26.98, 13, 2.7);
0170   TGeoMedium *med = new TGeoMedium("MED", 1, mat);
0171   TGeoVolume *top = gGeoManager->MakeBox("TOP", med, 100, 100, 100);
0172   gGeoManager->SetTopVolume(top);
0173   TGeoVolume *vol = gGeoManager->MakeTube("Tube", med, rmin, rmax, hz);
0174   vol->SetLineWidth(2);
0175   TGeoVolume *vols =
0176       gGeoManager->MakeTubs("Tube", med, rmin, rmax, hz, phimin, phimax);
0177   gGeoManager->CloseGeometry();
0178
0179   std::size_t idisc = 0;
0180   for (const auto &axes : allowedAxes) {
0181     auto [disc, thickness] = TGeoSurfaceConverter::toSurface(
0182         *vol->GetShape(), *gGeoIdentity, axes, 1);
0183     BOOST_REQUIRE_NE(disc, nullptr);
0184     BOOST_CHECK_EQUAL(disc->type(), Surface::Disc);
0185     CHECK_CLOSE_ABS(thickness, 2 * hz, s_epsilon);
0186
0187     auto bounds = dynamic_cast<const RadialBounds *>(&(disc->bounds()));
0188     BOOST_REQUIRE_NE(bounds, nullptr);
0189     double bminr = bounds->get(RadialBounds::eMinR);
0190     double bmaxr = bounds->get(RadialBounds::eMaxR);
0191
0192     CHECK_CLOSE_ABS(bminr, rmin, s_epsilon);
0193     CHECK_CLOSE_ABS(bmaxr, rmax, s_epsilon);
0194
0195     // Check if the surface is the (negative) identity
0196     GeometryView3D::drawSurface(objVis, *disc, tgContext);
0197     const Vector3 center = disc->center(tgContext);
0198     GeometryView3D::drawArrowForward(
0199         objVis, center, center + 1.2 * rmax * Vector3::UnitX(), 4., 2.5, red);
0200     GeometryView3D::drawArrowForward(
0201         objVis, center, center + 1.2 * rmax * Vector3::UnitY(), 4., 2.5, green);
0202     GeometryView3D::drawArrowForward(
0203         objVis, center, center + 1.2 * hz * Vector3::UnitZ(), 4., 2.5, blue);
0204     objVis.write("TGeoConversion_TGeoTube_DiscSurface_" +
0205                  std::to_string(idisc));
0206     objVis.clear();
0207
0208     if (idisc < 2) {
0209       auto [discSegment, dThickness] = TGeoSurfaceConverter::toSurface(
0210           *vols->GetShape(), *gGeoIdentity, axes, 1);
0211       BOOST_REQUIRE_NE(discSegment, nullptr);
0212       BOOST_CHECK_EQUAL(discSegment->type(), Surface::Disc);
0213       CHECK_CLOSE_ABS(dThickness, 2 * hz, s_epsilon);
0214
0215       auto boundsSegment =
0216           dynamic_cast<const RadialBounds *>(&(discSegment->bounds()));
0217       BOOST_REQUIRE_NE(boundsSegment, nullptr);
0218       bminr = boundsSegment->get(RadialBounds::eMinR);
0219       bmaxr = boundsSegment->get(RadialBounds::eMaxR);
0220       double hphi = boundsSegment->get(RadialBounds::eHalfPhiSector);
0221       double mphi = boundsSegment->get(RadialBounds::eAveragePhi);
0222       CHECK_CLOSE_ABS(bminr, rmin, s_epsilon);
0223       CHECK_CLOSE_ABS(bmaxr, rmax, s_epsilon);
0224       CHECK_CLOSE_ABS(hphi, std::numbers::pi / 4., s_epsilon);
0225       CHECK_CLOSE_ABS(mphi, 0., s_epsilon);
0226       GeometryView3D::drawSurface(objVis, *discSegment, tgContext);
0227       GeometryView3D::drawArrowForward(objVis, center,
0228                                        center + 1.2 * bmaxr * Vector3::UnitX(),
0229                                        4., 2.5, red);
0230       GeometryView3D::drawArrowForward(objVis, center,
0231                                        center + 1.2 * bmaxr * Vector3::UnitY(),
0232                                        4., 2.5, green);
0233       GeometryView3D::drawArrowForward(
0234           objVis, center, center + 1.2 * hz * Vector3::UnitZ(), 4., 2.5, blue);
0235       objVis.write("TGeoConversion_TGeoTube_DiscSegmentSurface_" +
0236                    std::to_string(idisc));
0237       objVis.clear();
0238
0239     } else {
0240       BOOST_CHECK_THROW(TGeoSurfaceConverter::toSurface(*vols->GetShape(),
0241                                                         *gGeoIdentity, axes, 1),
0242                         std::invalid_argument);
0243     }
0244     ++idisc;
0245   }
0246
0247   // Check exceptions for not allowed axis definition
0248   for (const auto &naxes : notAllowedAxes) {
0249     BOOST_CHECK_THROW(TGeoSurfaceConverter::toSurface(*vol->GetShape(),
0250                                                       *gGeoIdentity, naxes, 1),
0251                       std::invalid_argument);
0252   }
0253 }
0254
0255 BOOST_AUTO_TEST_SUITE_END()
0256
0257 }  // namespace ActsTests