Plugins/TGeo/TGeoTrd2ConversionTests.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/Surface.hpp"
0017 #include "Acts/Surfaces/SurfaceBounds.hpp"
0018 #include "Acts/Surfaces/TrapezoidBounds.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 <algorithm>
0025 #include <cstddef>
0026 #include <memory>
0027 #include <stdexcept>
0028 #include <string>
0029 #include <utility>
0030 #include <vector>
0031
0032 #include "TGeoManager.h"
0033 #include "TGeoMaterial.h"
0034 #include "TGeoMatrix.h"
0035 #include "TGeoMedium.h"
0036 #include "TGeoTrd1.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 /// @brief Unit test to convert a TGeoTrd2 into a Plane
0049 ///
0050 /// * The TGeoTrd2 has x/z orientation
0051 BOOST_AUTO_TEST_CASE(TGeoTrd2_xz_to_PlaneSurface) {
0052   ObjVisualization3D objVis;
0053
0054   double hxmin = 10.;
0055   double hxmax = 30.;
0056   double ht = 1.;  // this is the half thickness
0057   double hy = 40.;
0058
0059   new TGeoManager("trd1", "poza9");
0060   TGeoMaterial *mat = new TGeoMaterial("Al", 26.98, 13, 2.7);
0061   TGeoMedium *med = new TGeoMedium("MED", 1, mat);
0062   TGeoVolume *top = gGeoManager->MakeBox("TOP", med, 100, 100, 100);
0063   gGeoManager->SetTopVolume(top);
0064   TGeoVolume *vol =
0065       gGeoManager->MakeTrd2("Trd2", med, hxmin, hxmax, ht, ht, hy);
0066   gGeoManager->CloseGeometry();
0067
0068   // Check the 4 possible ways
0069   std::vector<std::string> axesTypes = {"XZ*", "xZ*", "xz*", "Xz*"};
0070
0071   std::size_t itrd = 0;
0072   for (const auto &axes : axesTypes) {
0073     auto [plane, thickness] = TGeoSurfaceConverter::toSurface(
0074         *vol->GetShape(), *gGeoIdentity, axes, 1);
0075     BOOST_REQUIRE_NE(plane, nullptr);
0076     BOOST_CHECK_EQUAL(plane->type(), Surface::Plane);
0077     CHECK_CLOSE_ABS(thickness, 2 * ht, s_epsilon);
0078
0079     auto bounds = dynamic_cast<const TrapezoidBounds *>(&(plane->bounds()));
0080     BOOST_REQUIRE_NE(bounds, nullptr);
0081     double hXminY = bounds->get(TrapezoidBounds::eHalfLengthXnegY);
0082     double hXmaxY = bounds->get(TrapezoidBounds::eHalfLengthXposY);
0083     double hY = bounds->get(TrapezoidBounds::eHalfLengthY);
0084
0085     CHECK_CLOSE_ABS(hxmin, std::min(hXminY, hXmaxY), s_epsilon);
0086     CHECK_CLOSE_ABS(hxmax, std::max(hXminY, hXmaxY), s_epsilon);
0087     CHECK_CLOSE_ABS(hy, hY, s_epsilon);
0088
0089     // Check if the surface is the (negative) identity
0090     auto transform = plane->transform(tgContext);
0091     auto rotation = transform.rotation();
0092     const Vector3 offset{(-5.5 + (itrd++) * 2.5) * hxmax, 0., 0.};
0093     GeometryView3D::drawSurface(objVis, *plane, tgContext,
0094                                 Translation3{offset} * Transform3::Identity());
0095     const Vector3 center = plane->center(tgContext) + offset;
0096     GeometryView3D::drawArrowForward(
0097         objVis, center, center + 1.2 * (hXminY + hXmaxY) * rotation.col(0), 4.,
0098         2.5, red);
0099     GeometryView3D::drawArrowForward(
0100         objVis, center, center + 1.2 * hY * rotation.col(1), 4., 2.5, green);
0101     GeometryView3D::drawArrowForward(
0102         objVis, center, center + 2 * rotation.col(2), 4., 2.5, blue);
0103   }
0104   objVis.write("TGeoConversion_TGeoTrd2_xz_PlaneSurface");
0105
0106   // Check exceptions for not allowed axis definition
0107   std::vector<std::string> notAllowed = {"XY*", "xy*", "Xy*", "xY*"};
0108   for (const auto &naxis : notAllowed) {
0109     BOOST_CHECK_THROW(TGeoSurfaceConverter::toSurface(*vol->GetShape(),
0110                                                       *gGeoIdentity, naxis, 1),
0111                       std::invalid_argument);
0112   }
0113 }
0114
0115 /// @brief Unit test to convert a TGeoTrd2 into a Plane
0116 ///
0117 /// * The TGeoTrd2 has y/z orientation
0118 BOOST_AUTO_TEST_CASE(TGeoTrd2_yz_to_PlaneSurface) {
0119   ObjVisualization3D objVis;
0120
0121   double hxmin = 10.;
0122   double hxmax = 30.;
0123   double ht = 1.;  // this is the half thickness
0124   double hy = 40.;
0125
0126   new TGeoManager("trd1", "poza9");
0127   TGeoMaterial *mat = new TGeoMaterial("Al", 26.98, 13, 2.7);
0128   TGeoMedium *med = new TGeoMedium("MED", 1, mat);
0129   TGeoVolume *top = gGeoManager->MakeBox("TOP", med, 100, 100, 100);
0130   gGeoManager->SetTopVolume(top);
0131   TGeoVolume *vol =
0132       gGeoManager->MakeTrd2("Trd2", med, ht, ht, hxmin, hxmax, hy);
0133   gGeoManager->CloseGeometry();
0134
0135   // Check the 4 possible ways
0136   std::vector<std::string> axesTypes = {"YZ*", "yZ*", "yz*", "Yz*"};
0137
0138   std::size_t itrd = 0;
0139   for (const auto &axes : axesTypes) {
0140     auto [plane, thickness] = TGeoSurfaceConverter::toSurface(
0141         *vol->GetShape(), *gGeoIdentity, axes, 1);
0142     BOOST_REQUIRE_NE(plane, nullptr);
0143     BOOST_CHECK_EQUAL(plane->type(), Surface::Plane);
0144     CHECK_CLOSE_ABS(thickness, 2 * ht, s_epsilon);
0145
0146     auto bounds = dynamic_cast<const TrapezoidBounds *>(&(plane->bounds()));
0147     BOOST_REQUIRE_NE(bounds, nullptr);
0148     double hXminY = bounds->get(TrapezoidBounds::eHalfLengthXnegY);
0149     double hXmaxY = bounds->get(TrapezoidBounds::eHalfLengthXposY);
0150     double hY = bounds->get(TrapezoidBounds::eHalfLengthY);
0151
0152     CHECK_CLOSE_ABS(hxmin, std::min(hXminY, hXmaxY), s_epsilon);
0153     CHECK_CLOSE_ABS(hxmax, std::max(hXminY, hXmaxY), s_epsilon);
0154     CHECK_CLOSE_ABS(hy, hY, s_epsilon);
0155
0156     // Check if the surface is the (negative) identity
0157     auto transform = plane->transform(tgContext);
0158     auto rotation = transform.rotation();
0159     const Vector3 offset{(-5.5 + (itrd++) * 2.5) * hxmax, 0., 0.};
0160     GeometryView3D::drawSurface(objVis, *plane, tgContext,
0161                                 Translation3{offset} * Transform3::Identity());
0162     const Vector3 center = plane->center(tgContext) + offset;
0163     GeometryView3D::drawArrowForward(
0164         objVis, center, center + 1.2 * (hXminY + hXmaxY) * rotation.col(0), 4.,
0165         2.5, red);
0166     GeometryView3D::drawArrowForward(
0167         objVis, center, center + 1.2 * hY * rotation.col(1), 4., 2.5, green);
0168     GeometryView3D::drawArrowForward(
0169         objVis, center, center + 2 * rotation.col(2), 4., 2.5, blue);
0170   }
0171   objVis.write("TGeoConversion_TGeoTrd2_yz_PlaneSurface");
0172
0173   // Check exceptions for not allowed axis definition
0174   std::vector<std::string> notAllowed = {"YX*", "yx*", "yX*", "Yx*"};
0175   for (const auto &naxis : notAllowed) {
0176     BOOST_CHECK_THROW(TGeoSurfaceConverter::toSurface(*vol->GetShape(),
0177                                                       *gGeoIdentity, naxis, 1),
0178                       std::invalid_argument);
0179   }
0180 }
0181
0182 }  // namespace Acts::Test