Plugins/Geant4/Geant4ConvertersTests.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/Material/ISurfaceMaterial.hpp"
0013 #include "Acts/Material/Material.hpp"
0014 #include "Acts/Material/MaterialSlab.hpp"
0015 #include "Acts/Surfaces/CylinderBounds.hpp"
0016 #include "Acts/Surfaces/LineBounds.hpp"
0017 #include "Acts/Surfaces/RadialBounds.hpp"
0018 #include "Acts/Surfaces/RectangleBounds.hpp"
0019 #include "Acts/Surfaces/Surface.hpp"
0020 #include "Acts/Surfaces/TrapezoidBounds.hpp"
0021 #include "ActsPlugins/Geant4/Geant4Converters.hpp"
0022 #include "ActsTests/CommonHelpers/FloatComparisons.hpp"
0023
0024 #include <array>
0025 #include <cmath>
0026 #include <memory>
0027 #include <numbers>
0028 #include <stdexcept>
0029 #include <tuple>
0030
0031 #include "G4Box.hh"
0032 #include "G4LogicalVolume.hh"
0033 #include "G4Material.hh"
0034 #include "G4PVPlacement.hh"
0035 #include "G4RotationMatrix.hh"
0036 #include "G4SystemOfUnits.hh"
0037 #include "G4ThreeVector.hh"
0038 #include "G4Trap.hh"
0039 #include "G4Trd.hh"
0040 #include "G4Tubs.hh"
0041 #include "G4VPhysicalVolume.hh"
0042
0043 using namespace Acts;
0044 using namespace ActsPlugins;
0045
0046 double rho = 1.2345;
0047 G4Material* g4Material = new G4Material("Material", 6., 12., rho);
0048
0049 namespace ActsTests {
0050
0051 BOOST_AUTO_TEST_SUITE(Geant4Suite)
0052
0053 BOOST_AUTO_TEST_CASE(Geant4AlgebraConversion) {
0054   G4ThreeVector g4Translation(10., 20., 30.);
0055
0056   auto translated = Geant4AlgebraConverter{}.transform(g4Translation);
0057   auto actsTranslation = translated.translation();
0058   BOOST_CHECK_EQUAL(actsTranslation[0], 10.);
0059   BOOST_CHECK_EQUAL(actsTranslation[1], 20.);
0060   BOOST_CHECK_EQUAL(actsTranslation[2], 30.);
0061
0062   auto translatedScaled = Geant4AlgebraConverter{10.}.transform(g4Translation);
0063   auto actsTranslationScaled = translatedScaled.translation();
0064   BOOST_CHECK_EQUAL(actsTranslationScaled[0], 100.);
0065   BOOST_CHECK_EQUAL(actsTranslationScaled[1], 200.);
0066   BOOST_CHECK_EQUAL(actsTranslationScaled[2], 300.);
0067 }
0068
0069 BOOST_AUTO_TEST_CASE(Geant4CylinderConversion) {
0070   G4Tubs cylinder("Cylinder", 399., 401., 800.,
0071                   -std::numbers::pi * CLHEP::radian,
0072                   2 * std::numbers::pi * CLHEP::radian);
0073   auto [bounds, thickness] = Geant4ShapeConverter{}.cylinderBounds(cylinder);
0074   CHECK_CLOSE_ABS(bounds->get(CylinderBounds::BoundValues::eR), 400., 10e-10);
0075   CHECK_CLOSE_ABS(bounds->get(CylinderBounds::BoundValues::eHalfLengthZ), 800.,
0076                   10e-10);
0077   CHECK_CLOSE_ABS(bounds->get(CylinderBounds::BoundValues::eHalfPhiSector),
0078                   std::numbers::pi, 10e-10);
0079   CHECK_CLOSE_ABS(bounds->get(CylinderBounds::BoundValues::eAveragePhi), 0.,
0080                   10e-10);
0081   CHECK_CLOSE_ABS(thickness, 2., 10e-10);
0082 }
0083
0084 BOOST_AUTO_TEST_CASE(Geant4RadialConversion) {
0085   G4Tubs disc("disc", 40., 400., 2., -std::numbers::pi * CLHEP::radian,
0086               2 * std::numbers::pi * CLHEP::radian);
0087   auto [bounds, thickness] = Geant4ShapeConverter{}.radialBounds(disc);
0088   CHECK_CLOSE_ABS(bounds->get(RadialBounds::BoundValues::eMinR), 40., 10e-10);
0089   CHECK_CLOSE_ABS(bounds->get(RadialBounds::BoundValues::eMaxR), 400., 10e-10);
0090   CHECK_CLOSE_ABS(bounds->get(RadialBounds::BoundValues::eHalfPhiSector),
0091                   std::numbers::pi, 10e-10);
0092   CHECK_CLOSE_ABS(bounds->get(RadialBounds::BoundValues::eAveragePhi), 0.,
0093                   10e-10);
0094   CHECK_CLOSE_ABS(thickness, 4., 10e-10);
0095 }
0096
0097 BOOST_AUTO_TEST_CASE(Geant4LineConversion) {
0098   G4Tubs line("line", 0., 20., 400., 0., 2 * std::numbers::pi);
0099   auto bounds = Geant4ShapeConverter{}.lineBounds(line);
0100   CHECK_CLOSE_ABS(bounds->get(LineBounds::BoundValues::eR), 20., 10e-10);
0101   CHECK_CLOSE_ABS(bounds->get(LineBounds::BoundValues::eHalfLengthZ), 400.,
0102                   10e-10);
0103 }
0104
0105 BOOST_AUTO_TEST_CASE(Geant4BoxConversion) {
0106   // Test the standard orientations
0107   G4Box sensorXY("SensorXY", 23., 34., 1.);
0108   auto [boundsXY, axesXY, thicknessZ] =
0109       Geant4ShapeConverter{}.rectangleBounds(sensorXY);
0110   CHECK_CLOSE_ABS(boundsXY->halfLengthX(), 23., 10e-10);
0111   CHECK_CLOSE_ABS(boundsXY->halfLengthY(), 34., 10e-10);
0112   auto refXY = std::array<int, 2u>{0, 1};
0113   BOOST_CHECK(axesXY == refXY);
0114   CHECK_CLOSE_ABS(thicknessZ, 2., 10e-10);
0115
0116   G4Box sensorYZ("SensorYZ", 2., 45., 56.);
0117   auto [boundsYZ, axesYZ, thicknessX] =
0118       Geant4ShapeConverter{}.rectangleBounds(sensorYZ);
0119   CHECK_CLOSE_ABS(boundsYZ->halfLengthX(), 45., 10e-10);
0120   CHECK_CLOSE_ABS(boundsYZ->halfLengthY(), 56., 10e-10);
0121   auto refYZ = std::array<int, 2u>{1, 2};
0122   BOOST_CHECK(axesYZ == refYZ);
0123   CHECK_CLOSE_ABS(thicknessX, 4., 10e-10);
0124
0125   G4Box sensorZX("SensorZX", 78., 2., 67.);
0126   auto [boundsZX, axesZX, thicknessY] =
0127       Geant4ShapeConverter{}.rectangleBounds(sensorZX);
0128   CHECK_CLOSE_ABS(boundsZX->halfLengthX(), 67., 10e-10);
0129   CHECK_CLOSE_ABS(boundsZX->halfLengthY(), 78., 10e-10);
0130   auto refZX = std::array<int, 2u>{2, 0};
0131   BOOST_CHECK(axesZX == refZX);
0132   CHECK_CLOSE_ABS(thicknessY, 4., 10e-10);
0133
0134   // Test the flipped axis
0135   G4Box sensorXz("SensorXz", 78., 2., 67.);
0136   auto [boundsXz, axesXz, thicknessY2] =
0137       Geant4ShapeConverter{1, true}.rectangleBounds(sensorXz);
0138   CHECK_CLOSE_ABS(boundsXz->halfLengthX(), 78., 10e-10);
0139   CHECK_CLOSE_ABS(boundsXz->halfLengthY(), 67., 10e-10);
0140   auto refXz = std::array<int, 2u>{0, -2};
0141   BOOST_CHECK(axesXz == refXz);
0142   CHECK_CLOSE_ABS(thicknessY2, 4., 10e-10);
0143 }
0144
0145 BOOST_AUTO_TEST_CASE(Geant4TrapzoidConversionTrd) {
0146   // Standard TRD: XY are already well defined
0147   G4Trd trdXY("trdXY", 100, 150, 200, 200, 2);
0148   auto [boundsXY, axesXY, thicknessZ] =
0149       Geant4ShapeConverter{}.trapezoidBounds(trdXY);
0150   CHECK_CLOSE_ABS(boundsXY->get(TrapezoidBounds::BoundValues::eHalfLengthXnegY),
0151                   100, 10e-10);
0152   CHECK_CLOSE_ABS(boundsXY->get(TrapezoidBounds::BoundValues::eHalfLengthXposY),
0153                   150, 10e-10);
0154   CHECK_CLOSE_ABS(boundsXY->get(TrapezoidBounds::BoundValues::eHalfLengthY),
0155                   200, 10e-10);
0156   auto refXY = std::array<int, 2u>{0, 1};
0157   BOOST_CHECK(axesXY == refXY);
0158   CHECK_CLOSE_ABS(thicknessZ, 4., 10e-10);
0159
0160   // Flipped, yX are the coordinates
0161   G4Trd trdyX("trdyX", 200, 200, 100, 150, 2);
0162   auto [boundsyX, axesyX, thicknessZ2] =
0163       Geant4ShapeConverter{}.trapezoidBounds(trdyX);
0164   CHECK_CLOSE_ABS(boundsyX->get(TrapezoidBounds::BoundValues::eHalfLengthXnegY),
0165                   100, 10e-10);
0166   CHECK_CLOSE_ABS(boundsyX->get(TrapezoidBounds::BoundValues::eHalfLengthXposY),
0167                   150, 10e-10);
0168   CHECK_CLOSE_ABS(boundsyX->get(TrapezoidBounds::BoundValues::eHalfLengthY),
0169                   200, 10e-10);
0170   auto refyX = std::array<int, 2u>{-1, 0};
0171   BOOST_CHECK(axesyX == refyX);
0172   CHECK_CLOSE_ABS(thicknessZ2, 4., 10e-10);
0173
0174   // YZ span the trapezoid
0175   G4Trd trdYZ("trdYZ", 2, 2, 120, 140, 200);
0176   auto [boundsYZ, axesYZ, thicknessX] =
0177       Geant4ShapeConverter{}.trapezoidBounds(trdYZ);
0178   CHECK_CLOSE_ABS(boundsYZ->get(TrapezoidBounds::BoundValues::eHalfLengthXnegY),
0179                   120., 10e-10);
0180   CHECK_CLOSE_ABS(boundsYZ->get(TrapezoidBounds::BoundValues::eHalfLengthXposY),
0181                   140., 10e-10);
0182   CHECK_CLOSE_ABS(boundsYZ->get(TrapezoidBounds::BoundValues::eHalfLengthY),
0183                   200., 10e-10);
0184   auto refYZ = std::array<int, 2u>{1, 2};
0185   BOOST_CHECK(axesYZ == refYZ);
0186   CHECK_CLOSE_ABS(thicknessX, 4., 10e-10);
0187
0188   // Xz span the trapezoid
0189   G4Trd trdXz("trdXz", 50, 75, 1, 1, 200);
0190   auto [boundsXz, axesXz, thicknessY] =
0191       Geant4ShapeConverter{}.trapezoidBounds(trdXz);
0192   CHECK_CLOSE_ABS(boundsXz->get(TrapezoidBounds::BoundValues::eHalfLengthXnegY),
0193                   50., 10e-10);
0194   CHECK_CLOSE_ABS(boundsXz->get(TrapezoidBounds::BoundValues::eHalfLengthXposY),
0195                   75., 10e-10);
0196   CHECK_CLOSE_ABS(boundsXz->get(TrapezoidBounds::BoundValues::eHalfLengthY),
0197                   200., 10e-10);
0198   auto refXz = std::array<int, 2u>{0, -2};
0199   BOOST_CHECK(axesXz == refXz);
0200   CHECK_CLOSE_ABS(thicknessY, 2., 10e-10);
0201 }
0202
0203 BOOST_AUTO_TEST_CASE(Geant4TrapzoidConversionTrap) {
0204   // x value changes, y is the symmetric side, z is the thickness
0205   G4Trap trapXY("trapXY", 10., 15., 20., 20., 0.125);
0206   auto [boundsXY, axesXY, thicknessZ] =
0207       Geant4ShapeConverter{}.trapezoidBounds(trapXY);
0208   CHECK_CLOSE_ABS(boundsXY->get(TrapezoidBounds::BoundValues::eHalfLengthXnegY),
0209                   10, 10e-10);
0210   CHECK_CLOSE_ABS(boundsXY->get(TrapezoidBounds::BoundValues::eHalfLengthXposY),
0211                   15, 10e-10);
0212   CHECK_CLOSE_ABS(boundsXY->get(TrapezoidBounds::BoundValues::eHalfLengthY), 20,
0213                   10e-10);
0214   auto refXY = std::array<int, 2u>{0, 1};
0215   BOOST_CHECK(axesXY == refXY);
0216   CHECK_CLOSE_ABS(thicknessZ, 0.25, 10e-10);
0217
0218   // y value changes, x is the symmetric side, z is the thickness
0219   G4Trap trapYx("trapYx", 22., 22., 11., 16., 0.250);
0220   auto [boundsYx, axesYx, thicknessYx] =
0221       Geant4ShapeConverter{}.trapezoidBounds(trapYx);
0222   CHECK_CLOSE_ABS(boundsYx->get(TrapezoidBounds::BoundValues::eHalfLengthXnegY),
0223                   11, 10e-10);
0224   CHECK_CLOSE_ABS(boundsYx->get(TrapezoidBounds::BoundValues::eHalfLengthXposY),
0225                   16, 10e-10);
0226   CHECK_CLOSE_ABS(boundsYx->get(TrapezoidBounds::BoundValues::eHalfLengthY), 22,
0227                   10e-10);
0228   auto refYx = std::array<int, 2u>{-1, 0};
0229   BOOST_CHECK(axesYx == refYx);
0230   CHECK_CLOSE_ABS(thicknessYx, 0.5, 10e-10);
0231
0232   // x is the thickness, y changes, z is symmetric
0233   G4Trap trapXz("trapXz", 0.5, 0.5, 8., 16., 10.);
0234   auto [boundsYZ, axesYZ, thicknessXZ] =
0235       Geant4ShapeConverter{}.trapezoidBounds(trapXz);
0236   CHECK_CLOSE_ABS(boundsYZ->get(TrapezoidBounds::BoundValues::eHalfLengthXnegY),
0237                   8, 10e-10);
0238   CHECK_CLOSE_ABS(boundsYZ->get(TrapezoidBounds::BoundValues::eHalfLengthXposY),
0239                   16, 10e-10);
0240   CHECK_CLOSE_ABS(boundsYZ->get(TrapezoidBounds::BoundValues::eHalfLengthY),
0241                   10., 10e-10);
0242   auto refYZ = std::array<int, 2u>{1, 2};
0243   BOOST_CHECK(axesYZ == refYZ);
0244   CHECK_CLOSE_ABS(thicknessXZ, 1., 10e-10);
0245 }
0246
0247 BOOST_AUTO_TEST_CASE(Geant4PlanarConversion) {
0248   G4Box boxXY("boxXY", 23., 34., 1.);
0249   auto pBoundsBox = std::get<0u>(Geant4ShapeConverter{}.planarBounds(boxXY));
0250   auto rBounds = dynamic_cast<const RectangleBounds*>(pBoundsBox.get());
0251   BOOST_CHECK_NE(rBounds, nullptr);
0252
0253   G4Trd trdXY("trdXY", 100, 150, 200, 200, 2);
0254   auto pBoundsTrd = std::get<0u>(Geant4ShapeConverter{}.planarBounds(trdXY));
0255   auto tBounds = dynamic_cast<const TrapezoidBounds*>(pBoundsTrd.get());
0256   BOOST_CHECK_NE(tBounds, nullptr);
0257 }
0258
0259 BOOST_AUTO_TEST_CASE(Geant4BoxVPhysConversion) {
0260   double thickness = 2.;
0261
0262   G4Box* g4Box = new G4Box("Box", 23., 34., 0.5 * thickness);
0263   G4RotationMatrix* g4Rot = new G4RotationMatrix({0., 0., 1.}, 1.2);
0264   G4LogicalVolume* g4BoxLog = new G4LogicalVolume(g4Box, g4Material, "BoxLog");
0265
0266   G4ThreeVector g4Trans(0., 0., 100.);
0267   G4PVPlacement g4BoxPhys(g4Rot, g4Trans, g4BoxLog, "BoxPhys", nullptr, false,
0268                           1);
0269
0270   auto planeSurface = Geant4PhysicalVolumeConverter{}.surface(
0271       g4BoxPhys, Transform3::Identity(), true, thickness);
0272   BOOST_REQUIRE_NE(planeSurface, nullptr);
0273   BOOST_CHECK_EQUAL(planeSurface->type(), Surface::SurfaceType::Plane);
0274
0275   auto material = planeSurface->surfaceMaterial();
0276   BOOST_REQUIRE_NE(material, nullptr);
0277
0278   auto materialSlab = material->materialSlab(Vector3{0., 0., 0.});
0279   // Here it should be uncompressed material
0280   CHECK_CLOSE_ABS(materialSlab.material().massDensity(), rho, 0.001);
0281   CHECK_CLOSE_REL(thickness / g4Material->GetRadlen(),
0282                   materialSlab.thicknessInX0(), 0.1);
0283
0284   // Convert with compression
0285   double compression = 4.;
0286   planeSurface = Geant4PhysicalVolumeConverter{}.surface(
0287       g4BoxPhys, Transform3::Identity(), true, thickness / compression);
0288   BOOST_REQUIRE_NE(planeSurface, nullptr);
0289   BOOST_CHECK_EQUAL(planeSurface->type(), Surface::SurfaceType::Plane);
0290
0291   material = planeSurface->surfaceMaterial();
0292   BOOST_REQUIRE_NE(material, nullptr);
0293   materialSlab = material->materialSlab(Vector3{0., 0., 0.});
0294
0295   // Here it should be uncompressed material
0296   CHECK_CLOSE_ABS(materialSlab.material().massDensity(), compression * rho,
0297                   0.001);
0298   CHECK_CLOSE_REL(thickness / g4Material->GetRadlen(),
0299                   materialSlab.thicknessInX0(), 0.01);
0300
0301   CHECK_CLOSE_ABS(materialSlab.thickness(), thickness / compression, 0.01);
0302   CHECK_CLOSE_REL(materialSlab.material().X0() * compression,
0303                   g4Material->GetRadlen(), 0.01);
0304
0305   delete g4Box;
0306   delete g4Rot;
0307   delete g4BoxLog;
0308 }
0309
0310 BOOST_AUTO_TEST_CASE(Geant4CylVPhysConversion) {
0311   double radius = 45.;
0312   double thickness = 1.;
0313   double halfLengthZ = 200;
0314
0315   G4Tubs* g4Tube = new G4Tubs("Tube", radius, radius + thickness, halfLengthZ,
0316                               -std::numbers::pi * CLHEP::radian,
0317                               2 * std::numbers::pi * CLHEP::radian);
0318
0319   G4RotationMatrix* g4Rot = new G4RotationMatrix({0., 0., 1.}, 0.);
0320   G4LogicalVolume* g4TubeLog =
0321       new G4LogicalVolume(g4Tube, g4Material, "TubeLog");
0322   G4ThreeVector g4Trans(0., 0., 100.);
0323   G4PVPlacement g4CylinderPhys(g4Rot, g4Trans, g4TubeLog, "TubePhys", nullptr,
0324                                false, 1);
0325
0326   auto cylinderSurface = Geant4PhysicalVolumeConverter{}.surface(
0327       g4CylinderPhys, Transform3::Identity(), true, thickness);
0328   BOOST_REQUIRE_NE(cylinderSurface, nullptr);
0329   BOOST_CHECK_EQUAL(cylinderSurface->type(), Surface::SurfaceType::Cylinder);
0330
0331   auto material = cylinderSurface->surfaceMaterial();
0332   BOOST_REQUIRE_NE(material, nullptr);
0333
0334   auto materialSlab = material->materialSlab(Vector3{0., 0., 0.});
0335   CHECK_CLOSE_REL(thickness / g4Material->GetRadlen(),
0336                   materialSlab.thicknessInX0(), 0.1);
0337
0338   // Here it should be uncompressed material
0339   CHECK_CLOSE_ABS(materialSlab.material().massDensity(), rho, 0.001);
0340
0341   /// CHECK exception throwing
0342   BOOST_CHECK_THROW(
0343       Geant4PhysicalVolumeConverter{Surface::SurfaceType::Plane}.surface(
0344           g4CylinderPhys, Transform3::Identity(), true, thickness),
0345       std::runtime_error);
0346
0347   delete g4Tube;
0348   delete g4Rot;
0349   delete g4TubeLog;
0350 }
0351
0352 BOOST_AUTO_TEST_CASE(Geant4VDiscVPhysConversion) {
0353   double innerRadius = 45.;
0354   double outerRadius = 75.;
0355   double thickness = 2.;
0356
0357   G4Tubs* g4Tube = new G4Tubs("Disc", innerRadius, outerRadius, 0.5 * thickness,
0358                               -std::numbers::pi * CLHEP::radian,
0359                               2 * std::numbers::pi * CLHEP::radian);
0360
0361   G4RotationMatrix* g4Rot = new G4RotationMatrix({0., 0., 1.}, 0.);
0362   G4LogicalVolume* g4TubeLog =
0363       new G4LogicalVolume(g4Tube, g4Material, "TubeLog");
0364   G4ThreeVector g4Trans(0., 0., 100.);
0365   G4PVPlacement g4discPhys(g4Rot, g4Trans, g4TubeLog, "TubePhys", nullptr,
0366                            false, 1);
0367
0368   auto discSurface = Geant4PhysicalVolumeConverter{}.surface(
0369       g4discPhys, Transform3::Identity(), true, thickness);
0370   BOOST_REQUIRE_NE(discSurface, nullptr);
0371   BOOST_CHECK_EQUAL(discSurface->type(), Surface::SurfaceType::Disc);
0372
0373   auto material = discSurface->surfaceMaterial();
0374   BOOST_REQUIRE_NE(material, nullptr);
0375
0376   auto materialSlab = material->materialSlab(Vector3{0., 0., 0.});
0377   // Here it should be uncompressed material
0378   CHECK_CLOSE_ABS(materialSlab.material().massDensity(), rho, 0.001);
0379
0380   delete g4Tube;
0381   delete g4Rot;
0382   delete g4TubeLog;
0383 }
0384
0385 BOOST_AUTO_TEST_CASE(Geant4LineVPhysConversion) {
0386   double innerRadius = 0.;
0387   double outerRadius = 20.;
0388   double thickness = 400.;
0389
0390   G4Tubs* g4Tube = new G4Tubs("Line", innerRadius, outerRadius, 0.5 * thickness,
0391                               -std::numbers::pi * CLHEP::radian,
0392                               2 * std::numbers::pi * CLHEP::radian);
0393
0394   G4RotationMatrix* g4Rot = new G4RotationMatrix({0., 0., 1.}, 0.);
0395   G4LogicalVolume* g4TubeLog =
0396       new G4LogicalVolume(g4Tube, g4Material, "LineLog");
0397   G4ThreeVector g4Trans(0., 0., 100.);
0398   G4PVPlacement g4linePhys(g4Rot, g4Trans, g4TubeLog, "LinePhys", nullptr,
0399                            false, 1);
0400
0401   auto lineSurface =
0402       Geant4PhysicalVolumeConverter{Surface::SurfaceType::Straw}.surface(
0403           g4linePhys, Transform3::Identity(), true, thickness);
0404   BOOST_REQUIRE_NE(lineSurface, nullptr);
0405   BOOST_CHECK_EQUAL(lineSurface->type(), Surface::SurfaceType::Straw);
0406
0407   delete g4Tube;
0408   delete g4Rot;
0409   delete g4TubeLog;
0410 }
0411
0412 BOOST_AUTO_TEST_SUITE_END()
0413
0414 }  // namespace ActsTests