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