Plugins/TGeo/TGeoTubeConversionTests.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2020 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 http://mozilla.org/MPL/2.0/.
0008
0009 #include <boost/test/data/test_case.hpp>
0010 #include <boost/test/unit_test.hpp>
0011
0012 #include "Acts/Definitions/Algebra.hpp"
0013 #include "Acts/Definitions/Tolerance.hpp"
0014 #include "Acts/Geometry/GeometryContext.hpp"
0015 #include "Acts/Plugins/TGeo/TGeoSurfaceConverter.hpp"
0016 #include "Acts/Surfaces/CylinderBounds.hpp"
0017 #include "Acts/Surfaces/RadialBounds.hpp"
0018 #include "Acts/Surfaces/Surface.hpp"
0019 #include "Acts/Surfaces/SurfaceBounds.hpp"
0020 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0021 #include "Acts/Visualization/GeometryView3D.hpp"
0022 #include "Acts/Visualization/ObjVisualization3D.hpp"
0023 #include "Acts/Visualization/ViewConfig.hpp"
0024
0025 #include <cmath>
0026 #include <cstddef>
0027 #include <memory>
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({200, 0, 0});
0045 ViewConfig green({0, 200, 0});
0046 ViewConfig blue({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
0104     objVis.write("TGeoConversion_TGeoTube_CylinderSurface_" +
0105                  std::to_string(icyl));
0106     objVis.clear();
0107
0108     if (icyl < 2) {
0109       auto [cylinderSegment, cThickness] = TGeoSurfaceConverter::toSurface(
0110           *vols->GetShape(), *gGeoIdentity, axes, 1);
0111       BOOST_REQUIRE_NE(cylinderSegment, nullptr);
0112       BOOST_CHECK_EQUAL(cylinderSegment->type(), Surface::Cylinder);
0113       CHECK_CLOSE_ABS(cThickness, rmax - rmin, s_epsilon);
0114
0115       auto boundsSegment =
0116           dynamic_cast<const CylinderBounds *>(&(cylinderSegment->bounds()));
0117       BOOST_REQUIRE_NE(boundsSegment, nullptr);
0118       bR = boundsSegment->get(CylinderBounds::eR);
0119       bhZ = boundsSegment->get(CylinderBounds::eHalfLengthZ);
0120       double hphi = boundsSegment->get(CylinderBounds::eHalfPhiSector);
0121       double mphi = boundsSegment->get(CylinderBounds::eAveragePhi);
0122       CHECK_CLOSE_ABS(bR, 10.5, s_epsilon);
0123       CHECK_CLOSE_ABS(bhZ, hz, s_epsilon);
0124       CHECK_CLOSE_ABS(hphi, 0.25 * M_PI, s_epsilon);
0125       CHECK_CLOSE_ABS(mphi, 0., s_epsilon);
0126       GeometryView3D::drawSurface(objVis, *cylinderSegment, tgContext);
0127       GeometryView3D::drawArrowForward(
0128           objVis, center, center + 1.2 * bR * rotation.col(0), 4., 2.5, red);
0129       GeometryView3D::drawArrowForward(
0130           objVis, center, center + 1.2 * bR * rotation.col(1), 4., 2.5, green);
0131       GeometryView3D::drawArrowForward(
0132           objVis, center, center + 1.2 * bhZ * rotation.col(2), 4., 2.5, blue);
0133       objVis.write("TGeoConversion_TGeoTube_CylinderSegmentSurface_" +
0134                    std::to_string(icyl));
0135       objVis.clear();
0136     } else {
0137       BOOST_CHECK_THROW(TGeoSurfaceConverter::toSurface(*vols->GetShape(),
0138                                                         *gGeoIdentity, axes, 1),
0139                         std::invalid_argument);
0140     }
0141     ++icyl;
0142   }
0143
0144   // Check exceptions for not allowed axis definition
0145   for (const auto &naxes : notAllowedAxes) {
0146     BOOST_CHECK_THROW(TGeoSurfaceConverter::toSurface(*vol->GetShape(),
0147                                                       *gGeoIdentity, naxes, 1),
0148                       std::invalid_argument);
0149   }
0150 }
0151
0152 /// @brief Unit test to convert a TGeoTube into a DiscSurface
0153 ///
0154 /// * The TGeoTrd1 can only have (x/X)(y/Y) orientation
0155 BOOST_AUTO_TEST_CASE(TGeoTube_to_DiscSurface) {
0156   ObjVisualization3D objVis;
0157
0158   double rmin = 5.;
0159   double rmax = 25;
0160   double hz = 2.;
0161   double phimin = 45.;
0162   double phimax = -45.;
0163
0164   new TGeoManager("trd1", "poza9");
0165   TGeoMaterial *mat = new TGeoMaterial("Al", 26.98, 13, 2.7);
0166   TGeoMedium *med = new TGeoMedium("MED", 1, mat);
0167   TGeoVolume *top = gGeoManager->MakeBox("TOP", med, 100, 100, 100);
0168   gGeoManager->SetTopVolume(top);
0169   TGeoVolume *vol = gGeoManager->MakeTube("Tube", med, rmin, rmax, hz);
0170   vol->SetLineWidth(2);
0171   TGeoVolume *vols =
0172       gGeoManager->MakeTubs("Tube", med, rmin, rmax, hz, phimin, phimax);
0173   gGeoManager->CloseGeometry();
0174
0175   std::size_t idisc = 0;
0176   for (const auto &axes : allowedAxes) {
0177     auto [disc, thickness] = TGeoSurfaceConverter::toSurface(
0178         *vol->GetShape(), *gGeoIdentity, axes, 1);
0179     BOOST_REQUIRE_NE(disc, nullptr);
0180     BOOST_CHECK_EQUAL(disc->type(), Surface::Disc);
0181     CHECK_CLOSE_ABS(thickness, 2 * hz, s_epsilon);
0182
0183     auto bounds = dynamic_cast<const RadialBounds *>(&(disc->bounds()));
0184     BOOST_REQUIRE_NE(bounds, nullptr);
0185     double bminr = bounds->get(RadialBounds::eMinR);
0186     double bmaxr = bounds->get(RadialBounds::eMaxR);
0187
0188     CHECK_CLOSE_ABS(bminr, rmin, s_epsilon);
0189     CHECK_CLOSE_ABS(bmaxr, rmax, s_epsilon);
0190
0191     // Check if the surface is the (negative) identity
0192     GeometryView3D::drawSurface(objVis, *disc, tgContext);
0193     const Vector3 center = disc->center(tgContext);
0194     GeometryView3D::drawArrowForward(
0195         objVis, center, center + 1.2 * rmax * Vector3::UnitX(), 4., 2.5, red);
0196     GeometryView3D::drawArrowForward(
0197         objVis, center, center + 1.2 * rmax * Vector3::UnitY(), 4., 2.5, green);
0198     GeometryView3D::drawArrowForward(
0199         objVis, center, center + 1.2 * hz * Vector3::UnitZ(), 4., 2.5, blue);
0200     objVis.write("TGeoConversion_TGeoTube_DiscSurface_" +
0201                  std::to_string(idisc));
0202     objVis.clear();
0203
0204     if (idisc < 2) {
0205       auto [discSegment, dThickness] = TGeoSurfaceConverter::toSurface(
0206           *vols->GetShape(), *gGeoIdentity, axes, 1);
0207       BOOST_REQUIRE_NE(discSegment, nullptr);
0208       BOOST_CHECK_EQUAL(discSegment->type(), Surface::Disc);
0209       CHECK_CLOSE_ABS(dThickness, 2 * hz, s_epsilon);
0210
0211       auto boundsSegment =
0212           dynamic_cast<const RadialBounds *>(&(discSegment->bounds()));
0213       BOOST_REQUIRE_NE(boundsSegment, nullptr);
0214       bminr = boundsSegment->get(RadialBounds::eMinR);
0215       bmaxr = boundsSegment->get(RadialBounds::eMaxR);
0216       double hphi = boundsSegment->get(RadialBounds::eHalfPhiSector);
0217       double mphi = boundsSegment->get(RadialBounds::eAveragePhi);
0218       CHECK_CLOSE_ABS(bminr, rmin, s_epsilon);
0219       CHECK_CLOSE_ABS(bmaxr, rmax, s_epsilon);
0220       CHECK_CLOSE_ABS(hphi, 0.25 * M_PI, s_epsilon);
0221       CHECK_CLOSE_ABS(mphi, 0., s_epsilon);
0222       GeometryView3D::drawSurface(objVis, *discSegment, tgContext);
0223       GeometryView3D::drawArrowForward(objVis, center,
0224                                        center + 1.2 * bmaxr * Vector3::UnitX(),
0225                                        4., 2.5, red);
0226       GeometryView3D::drawArrowForward(objVis, center,
0227                                        center + 1.2 * bmaxr * Vector3::UnitY(),
0228                                        4., 2.5, green);
0229       GeometryView3D::drawArrowForward(
0230           objVis, center, center + 1.2 * hz * Vector3::UnitZ(), 4., 2.5, blue);
0231       objVis.write("TGeoConversion_TGeoTube_DiscSegmentSurface_" +
0232                    std::to_string(idisc));
0233       objVis.clear();
0234
0235     } else {
0236       BOOST_CHECK_THROW(TGeoSurfaceConverter::toSurface(*vols->GetShape(),
0237                                                         *gGeoIdentity, axes, 1),
0238                         std::invalid_argument);
0239     }
0240     ++idisc;
0241   }
0242
0243   // Check exceptions for not allowed axis definition
0244   for (const auto &naxes : notAllowedAxes) {
0245     BOOST_CHECK_THROW(TGeoSurfaceConverter::toSurface(*vol->GetShape(),
0246                                                       *gGeoIdentity, naxes, 1),
0247                       std::invalid_argument);
0248   }
0249 }
0250
0251 }  // namespace Acts::Test