Plugins/TGeo/TGeoBBoxConversionTests.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2020-2023 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/RectangleBounds.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 <memory>
0025 #include <utility>
0026
0027 #include "TGeoBBox.h"
0028 #include "TGeoManager.h"
0029 #include "TGeoMaterial.h"
0030 #include "TGeoMatrix.h"
0031 #include "TGeoMedium.h"
0032 #include "TGeoVolume.h"
0033 #include "TView.h"
0034
0035 namespace Acts::Test {
0036
0037 GeometryContext tgContext = GeometryContext();
0038
0039 ViewConfig red({200, 0, 0});
0040 ViewConfig green({0, 200, 0});
0041 ViewConfig blue({0, 0, 200});
0042
0043 /// @brief Unit test to convert a Bbox into a Plane
0044 ///
0045 /// This test also tests:
0046 /// * the "(x/X)(y/Y)(z/Z)" orientations
0047 /// * the scaling functionality
0048 BOOST_AUTO_TEST_CASE(TGeoBBox_to_PlaneSurface) {
0049   ObjVisualization3D objVis;
0050
0051   // BBox is defined [-dX,dX] x [-dY,dY] x [-dZ,dZ]
0052   double dX = 10.;
0053   double dY = 30.;
0054   double dZ = 1.;
0055
0056   new TGeoManager("box", "poza1");
0057   TGeoMaterial *mat = new TGeoMaterial("Al", 26.98, 13, 2.7);
0058   TGeoMedium *med = new TGeoMedium("MED", 1, mat);
0059   TGeoVolume *top = gGeoManager->MakeBox("TOP", med, 100, 100, 100);
0060   gGeoManager->SetTopVolume(top);
0061   TGeoVolume *vol = gGeoManager->MakeBox("BOX", med, dX, dY, dZ);
0062   vol->SetLineWidth(2);
0063   top->AddNode(vol, 1);
0064   gGeoManager->CloseGeometry();
0065
0066   // Upper case ---------------------------------
0067   auto [plane_XYZ, thickness_XYZ] = TGeoSurfaceConverter::toSurface(
0068       *vol->GetShape(), *gGeoIdentity, "XY*", 1);
0069   BOOST_REQUIRE_NE(plane_XYZ, nullptr);
0070   BOOST_CHECK_EQUAL(plane_XYZ->type(), Surface::Plane);
0071   CHECK_CLOSE_ABS(thickness_XYZ, 2 * dZ, s_epsilon);
0072
0073   auto bounds_XYZ =
0074       dynamic_cast<const RectangleBounds *>(&(plane_XYZ->bounds()));
0075   BOOST_REQUIRE_NE(bounds_XYZ, nullptr);
0076   double maxX = bounds_XYZ->get(RectangleBounds::eMaxX);
0077   double minX = bounds_XYZ->get(RectangleBounds::eMinX);
0078   double maxY = bounds_XYZ->get(RectangleBounds::eMaxY);
0079   double minY = bounds_XYZ->get(RectangleBounds::eMinY);
0080   CHECK_CLOSE_ABS(maxX - minX, 2 * dX, s_epsilon);
0081   CHECK_CLOSE_ABS(maxY - minY, 2 * dY, s_epsilon);
0082
0083   // Check if the surface is the (negative) identity
0084   auto transform_XYZ = plane_XYZ->transform(tgContext);
0085   auto rotation_XYZ = transform_XYZ.rotation();
0086   BOOST_CHECK(transform_XYZ.isApprox(Transform3::Identity()));
0087
0088   const Vector3 offset_XYZ{-5.5 * dX, 0., 0.};
0089   GeometryView3D::drawSurface(objVis, *plane_XYZ, tgContext,
0090                               Transform3(Translation3{offset_XYZ}));
0091   const Vector3 center_XYZ = plane_XYZ->center(tgContext) + offset_XYZ;
0092   GeometryView3D::drawArrowForward(
0093       objVis, center_XYZ,
0094       center_XYZ + 0.6 * (maxX - minX) * rotation_XYZ.col(0), 4., 2.5, red);
0095   GeometryView3D::drawArrowForward(
0096       objVis, center_XYZ,
0097       center_XYZ + 0.6 * (maxY - minY) * rotation_XYZ.col(1), 4., 2.5, green);
0098   GeometryView3D::drawArrowForward(
0099       objVis, center_XYZ, center_XYZ + 2 * rotation_XYZ.col(2), 4., 2.5, blue);
0100
0101   // Lower case ---------------------------------
0102   auto [plane_xyz, thickness_xyz] = TGeoSurfaceConverter::toSurface(
0103       *vol->GetShape(), *gGeoIdentity, "xy*", 1);
0104   BOOST_CHECK_NE(plane_xyz, nullptr);
0105   BOOST_CHECK_EQUAL(plane_xyz->type(), Surface::Plane);
0106   CHECK_CLOSE_ABS(thickness_xyz, 2 * dZ, s_epsilon);
0107
0108   auto bounds_xyz =
0109       dynamic_cast<const RectangleBounds *>(&(plane_XYZ->bounds()));
0110   BOOST_REQUIRE_NE(bounds_xyz, nullptr);
0111   BOOST_CHECK_EQUAL(bounds_xyz, bounds_XYZ);
0112   auto transform_xyz = plane_xyz->transform(tgContext);
0113   auto rotation_xyz = transform_xyz.rotation();
0114   BOOST_CHECK(rotation_xyz.col(0).isApprox(-1 * rotation_XYZ.col(0)));
0115   BOOST_CHECK(rotation_xyz.col(1).isApprox(-1 * rotation_XYZ.col(1)));
0116   BOOST_CHECK(rotation_xyz.col(2).isApprox(rotation_XYZ.col(2)));
0117
0118   const Vector3 offset_xyz{-2 * dX, 0., 0.};
0119   GeometryView3D::drawSurface(objVis, *plane_xyz, tgContext,
0120                               Transform3(Translation3{offset_xyz}));
0121   const Vector3 center_xyz = plane_xyz->center(tgContext) + offset_xyz;
0122   GeometryView3D::drawArrowForward(
0123       objVis, center_xyz,
0124       center_xyz + 0.6 * (maxX - minX) * rotation_xyz.col(0), 4., 2.5, red);
0125   GeometryView3D::drawArrowForward(
0126       objVis, center_xyz,
0127       center_xyz + 0.6 * (maxY - minY) * rotation_xyz.col(1), 4., 2.5, green);
0128   GeometryView3D::drawArrowForward(
0129       objVis, center_xyz, center_xyz + 2 * rotation_xyz.col(2), 4., 2.5, blue);
0130
0131   // Mixed case ---------------------------------
0132   auto [plane_xYz, thickness_xYz] = TGeoSurfaceConverter::toSurface(
0133       *vol->GetShape(), *gGeoIdentity, "xY*", 1);
0134   BOOST_REQUIRE_NE(plane_xYz, nullptr);
0135   BOOST_CHECK_EQUAL(plane_xYz->type(), Surface::Plane);
0136   CHECK_CLOSE_ABS(thickness_xYz, 2 * dZ, s_epsilon);
0137
0138   auto bounds_xYz =
0139       dynamic_cast<const RectangleBounds *>(&(plane_xYz->bounds()));
0140   BOOST_CHECK_NE(bounds_xYz, nullptr);
0141   BOOST_CHECK_EQUAL(bounds_xYz, bounds_xYz);
0142   auto transform_xYz = plane_xYz->transform(tgContext);
0143   auto rotation_xYz = transform_xYz.rotation();
0144   BOOST_CHECK(rotation_xYz.col(0).isApprox(-1 * rotation_XYZ.col(0)));
0145   BOOST_CHECK(rotation_xYz.col(1).isApprox(rotation_XYZ.col(1)));
0146   BOOST_CHECK(rotation_xYz.col(2).isApprox(-1. * rotation_XYZ.col(2)));
0147
0148   const Vector3 offset_xYz{2 * dX, 0., 0.};
0149   GeometryView3D::drawSurface(
0150       objVis, *plane_xYz, tgContext,
0151       Translation3{offset_xYz} * Transform3::Identity());
0152   const Vector3 center_xYz = plane_xYz->center(tgContext) + offset_xYz;
0153   GeometryView3D::drawArrowForward(
0154       objVis, center_xYz,
0155       center_xYz + 0.6 * (maxX - minX) * rotation_xYz.col(0), 4., 2.5, red);
0156   GeometryView3D::drawArrowForward(
0157       objVis, center_xYz,
0158       center_xYz + 0.6 * (maxY - minY) * rotation_xYz.col(1), 4., 2.5, green);
0159   GeometryView3D::drawArrowForward(
0160       objVis, center_xYz, center_xYz + 2 * rotation_xYz.col(2), 4., 2.5, blue);
0161
0162   // Swap case --------------------------------- (x/y) here
0163   auto [plane_YXz, thickness_YXz] = TGeoSurfaceConverter::toSurface(
0164       *vol->GetShape(), *gGeoIdentity, "YX*", 1);
0165   BOOST_REQUIRE_NE(plane_YXz, nullptr);
0166   BOOST_CHECK_EQUAL(plane_YXz->type(), Surface::Plane);
0167   CHECK_CLOSE_ABS(thickness_YXz, 2 * dZ, s_epsilon);
0168
0169   auto bounds_YXz =
0170       dynamic_cast<const RectangleBounds *>(&(plane_YXz->bounds()));
0171   maxX = bounds_YXz->get(RectangleBounds::eMaxX);
0172   minX = bounds_YXz->get(RectangleBounds::eMinX);
0173   maxY = bounds_YXz->get(RectangleBounds::eMaxY);
0174   minY = bounds_YXz->get(RectangleBounds::eMinY);
0175   CHECK_CLOSE_ABS(maxX - minX, 2 * dY, s_epsilon);
0176   CHECK_CLOSE_ABS(maxY - minY, 2 * dX, s_epsilon);
0177
0178   auto transform_YXz = plane_YXz->transform(tgContext);
0179   auto rotation_YXz = transform_YXz.rotation();
0180   BOOST_CHECK(rotation_YXz.col(0).isApprox(rotation_XYZ.col(1)));
0181   BOOST_CHECK(rotation_YXz.col(1).isApprox(rotation_XYZ.col(0)));
0182   BOOST_CHECK(rotation_YXz.col(2).isApprox(-1. * rotation_XYZ.col(2)));
0183
0184   const Vector3 offset_YXz{5.5 * dX, 0., 0.};
0185   GeometryView3D::drawSurface(objVis, *plane_YXz, tgContext,
0186                               Transform3(Translation3{offset_YXz}));
0187   const Vector3 center_YXz = plane_YXz->center(tgContext) + offset_YXz;
0188   GeometryView3D::drawArrowForward(
0189       objVis, center_YXz,
0190       center_YXz + 0.6 * (maxX - minX) * rotation_YXz.col(0), 4., 2.5, red);
0191   GeometryView3D::drawArrowForward(
0192       objVis, center_YXz,
0193       center_YXz + 0.6 * (maxY - minY) * rotation_YXz.col(1), 4., 2.5, green);
0194   GeometryView3D::drawArrowForward(
0195       objVis, center_YXz, center_YXz + 2 * rotation_YXz.col(2), 4., 2.5, blue);
0196
0197   // Scaling test ---------------------------------
0198   auto [plane_XYZ10, thickness_XYZ10] = TGeoSurfaceConverter::toSurface(
0199       *vol->GetShape(), *gGeoIdentity, "xY*", 10);
0200   BOOST_CHECK_NE(plane_XYZ10, nullptr);
0201   CHECK_CLOSE_ABS(thickness_XYZ10, 20 * dZ, s_epsilon);
0202
0203   auto bounds_XYZ10 =
0204       dynamic_cast<const RectangleBounds *>(&(plane_XYZ10->bounds()));
0205   double maxX10 = bounds_XYZ10->get(RectangleBounds::eMaxX);
0206   double minX10 = bounds_XYZ10->get(RectangleBounds::eMinX);
0207   double maxY10 = bounds_XYZ10->get(RectangleBounds::eMaxY);
0208   double minY10 = bounds_XYZ10->get(RectangleBounds::eMinY);
0209   CHECK_CLOSE_ABS(maxX10 - minX10, 20 * dX, s_epsilon);
0210   CHECK_CLOSE_ABS(maxY10 - minY10, 20 * dY, s_epsilon);
0211
0212   objVis.write("TGeoConversion_TGeoBBox_PlaneSurface");
0213 }
0214
0215 }  // namespace Acts::Test