diff --git a/Plugins/Geant4/include/Acts/Plugins/Geant4/Geant4Converters.hpp b/Plugins/Geant4/include/Acts/Plugins/Geant4/Geant4Converters.hpp index a65b9d5c711..d2f43044f7b 100644 --- a/Plugins/Geant4/include/Acts/Plugins/Geant4/Geant4Converters.hpp +++ b/Plugins/Geant4/include/Acts/Plugins/Geant4/Geant4Converters.hpp @@ -133,6 +133,14 @@ struct Geant4ShapeConverter { std::tuple, std::array, ActsScalar> trapezoidBounds(const G4Trd& g4Trd); + /// @brief Convert to trapezoid bounds - from Trap + /// + /// @param g4Trd a Geant4 trapezoid shape + /// + /// @return an ACTS Trapezoid bounds object, axis orientation, and thickness + std::tuple, std::array, ActsScalar> + trapezoidBounds(const G4Trap& g4Trap); + /// @brief Convert to general solid into a planar shape /// /// @param g4Solid a Geant4 solid shape diff --git a/Plugins/Geant4/src/Geant4Converters.cpp b/Plugins/Geant4/src/Geant4Converters.cpp index db7d9f93f1e..4e0ba107893 100644 --- a/Plugins/Geant4/src/Geant4Converters.cpp +++ b/Plugins/Geant4/src/Geant4Converters.cpp @@ -37,6 +37,7 @@ #include "G4RotationMatrix.hh" #include "G4ThreeVector.hh" #include "G4Transform3D.hh" +#include "G4Trap.hh" #include "G4Trd.hh" #include "G4Tubs.hh" #include "G4VPhysicalVolume.hh" @@ -234,6 +235,71 @@ Acts::Geant4ShapeConverter::trapezoidBounds(const G4Trd& g4Trd) { return std::make_tuple(std::move(tBounds), rAxes, thickness); } +std::tuple, std::array, + Acts::ActsScalar> +Acts::Geant4ShapeConverter::trapezoidBounds(const G4Trap& g4Trap) { + // primary parameters + ActsScalar y1 = static_cast(g4Trap.GetYHalfLength1()); + ActsScalar y2 = static_cast(g4Trap.GetYHalfLength2()); + ActsScalar x1 = static_cast(g4Trap.GetXHalfLength1()); + ActsScalar x2 = static_cast(g4Trap.GetXHalfLength2()); + ActsScalar x3 = static_cast(g4Trap.GetXHalfLength3()); + ActsScalar x4 = static_cast(g4Trap.GetXHalfLength4()); + ActsScalar phi = static_cast(g4Trap.GetPhi()); + ActsScalar theta = static_cast(g4Trap.GetTheta()); + ActsScalar z = static_cast(g4Trap.GetZHalfLength()); + + ActsScalar hlX0 = (x1 + x2) * 0.5; + ActsScalar hlX1 = 2 * z * std::tan(theta) * std::cos(phi) + (x3 + x4) * 0.5; + ActsScalar hlY0 = y1; + ActsScalar hlY1 = y2 + 2 * z * std::tan(theta) * std::sin(phi); + ActsScalar hlZ = z; + + std::vector dXYZ = {(hlX0 + hlX1) * 0.5, (hlY0 + hlY1) * 0.5, + hlZ}; + + auto minAt = std::min_element(dXYZ.begin(), dXYZ.end()); + std::size_t minPos = std::distance(dXYZ.begin(), minAt); + ActsScalar thickness = 2. * dXYZ[minPos]; + + ActsScalar halfLengthXminY = 0.; + ActsScalar halfLengthXmaxY = 0.; + ActsScalar halfLengthY = 0.; + + std::array rAxes = {}; + switch (minPos) { + case 0: { + halfLengthXminY = std::min(hlY0, hlY1); + halfLengthXmaxY = std::max(hlY0, hlY1); + halfLengthY = hlZ; + rAxes = {1, 2}; + } break; + case 1: { + halfLengthXminY = std::min(hlX0, hlX1); + halfLengthXmaxY = std::max(hlX0, hlX1); + halfLengthY = hlZ; + rAxes = {0, -2}; + } break; + case 2: { + if (std::abs(hlY0 - hlY1) < std::abs(hlX0 - hlX1)) { + halfLengthXminY = std::min(hlX0, hlX1); + halfLengthXmaxY = std::max(hlX0, hlX1); + halfLengthY = (hlY0 + hlY1) * 0.5; + rAxes = {0, 1}; + } else { + halfLengthXminY = std::min(hlY0, hlY1); + halfLengthXmaxY = std::max(hlY0, hlY1); + halfLengthY = (hlX0 + hlX1) * 0.5; + rAxes = {-1, 0}; + } + } break; + } + + auto tBounds = std::make_shared( + halfLengthXminY, halfLengthXmaxY, halfLengthY); + return std::make_tuple(std::move(tBounds), rAxes, thickness); +} + std::tuple, std::array, Acts::ActsScalar> Acts::Geant4ShapeConverter::planarBounds(const G4VSolid& g4Solid) { @@ -332,6 +398,23 @@ std::shared_ptr Acts::Geant4PhysicalVolumeConverter::surface( } } + // Into a Trapezoid (G4Trap) + auto g4Trap = dynamic_cast(g4Solid); + if (g4Trap != nullptr) { + if (forcedType == Surface::SurfaceType::Other || + forcedType == Surface::SurfaceType::Plane) { + auto [bounds, axes, original] = + Geant4ShapeConverter{}.trapezoidBounds(*g4Trap); + auto orientedToGlobal = axesOriented(toGlobal, axes); + surface = Acts::Surface::makeShared(orientedToGlobal, + std::move(bounds)); + assignMaterial(*surface.get(), original, compressed); + return surface; + } else { + throw std::runtime_error("Can not convert 'G4Trap' into forced shape."); + } + } + // Into a Cylinder, disc or line auto g4Tubs = dynamic_cast(g4Solid); if (g4Tubs != nullptr) { diff --git a/Tests/UnitTests/Plugins/Geant4/Geant4ConvertersTests.cpp b/Tests/UnitTests/Plugins/Geant4/Geant4ConvertersTests.cpp index c748921b4d0..80083144e80 100644 --- a/Tests/UnitTests/Plugins/Geant4/Geant4ConvertersTests.cpp +++ b/Tests/UnitTests/Plugins/Geant4/Geant4ConvertersTests.cpp @@ -210,6 +210,80 @@ BOOST_AUTO_TEST_CASE(Geant4TrapzoidConversion) { CHECK_CLOSE_ABS(thicknessY, 2., 10e-10); } +BOOST_AUTO_TEST_CASE(Geant4TrapzoidConversion) { + // Standard Trap: XY are already well defined + G4Trap trapXY("trapXY", 2, 0.523599, 0.785398, 125, 200, 125, 0.174533, 50, + 125, 50, 0.174533); + auto [boundsXY, axesXY, thicknessZ] = + Acts::Geant4ShapeConverter{}.TrapezoidBounds(trapXY); + CHECK_CLOSE_ABS( + boundsXY->get(Acts::TrapezoidBounds::BoundValues::eHalfLengthXnegY), 51, + 10e-10); + CHECK_CLOSE_ABS( + boundsXY->get(Acts::TrapezoidBounds::BoundValues::eHalfLengthXposY), 125, + 10e-10); + CHECK_CLOSE_ABS( + boundsXY->get(Acts::TrapezoidBounds::BoundValues::eHalfLengthY), 125, + 10e-10); + auto refXY = std::array{0, 1}; + BOOST_CHECK(axesXY == refXY); + CHECK_CLOSE_ABS(thicknessZ, 4., 10e-10); + + // Flipped, yX are the coordinates + G4Trap trapyX("trapyX", 2, 0.523599, 0.785398, 50, 125, 50, 0.174533, 125, + 200, 125, 0.174533); + auto [boundsyX, axesyX, thicknessZ2] = + Acts::Geant4ShapeConverter{}.trapezoidBounds(trapyX); + CHECK_CLOSE_ABS( + boundsyX->get(Acts::TrapezoidBounds::BoundValues::eHalfLengthXnegY), 87, + 10e-10); + CHECK_CLOSE_ABS( + boundsyX->get(Acts::TrapezoidBounds::BoundValues::eHalfLengthXposY), 164, + 10e-10); + CHECK_CLOSE_ABS( + boundsyX->get(Acts::TrapezoidBounds::BoundValues::eHalfLengthY), 88, + 10e-10); + auto refyX = std::array{-1, 0}; + BOOST_CHECK(axesyX == refyX); + CHECK_CLOSE_ABS(thicknessZ2, 4., 10e-10); + + // YZ span the trapezoid + G4Trap trapYZ("trapYZ", 200, 0.523599, 0.785398, 140, 2, 2, 0.174533, 120, 2, + 2, 0.174533); + auto [boundsYZ, axesYZ, thicknessX] = + Acts::Geant4ShapeConverter{}.trapezoidBounds(trapYZ); + CHECK_CLOSE_ABS( + boundsYZ->get(Acts::TrapezoidBounds::BoundValues::eHalfLengthXnegY), 140., + 10e-10); + CHECK_CLOSE_ABS( + boundsYZ->get(Acts::TrapezoidBounds::BoundValues::eHalfLengthXposY), 283., + 10e-10); + CHECK_CLOSE_ABS( + boundsYZ->get(Acts::TrapezoidBounds::BoundValues::eHalfLengthY), 200., + 10e-10); + auto refYZ = std::array{1, 2}; + BOOST_CHECK(axesYZ == refYZ); + CHECK_CLOSE_ABS(thicknessX, 166, 10e-10); + + // Xz span the trapezoid + G4Trap trapXz("trapXz", 200, 0.523599, 0.785398, 2, 150, 100, 0.174533, 2, + 150, 100, 0.174533); + auto [boundsXz, axesXz, thicknessY] = + Acts::Geant4ShapeConverter{}.trapezoidBounds(trapXz); + CHECK_CLOSE_ABS( + boundsXz->get(Acts::TrapezoidBounds::BoundValues::eHalfLengthXnegY), 125., + 10e-10); + CHECK_CLOSE_ABS( + boundsXz->get(Acts::TrapezoidBounds::BoundValues::eHalfLengthXposY), 288., + 10e-10); + CHECK_CLOSE_ABS( + boundsXz->get(Acts::TrapezoidBounds::BoundValues::eHalfLengthY), 200., + 10e-10); + auto refXz = std::array{0, -2}; + BOOST_CHECK(axesXz == refXz); + CHECK_CLOSE_ABS(thicknessY, 166, 10e-10); +} + BOOST_AUTO_TEST_CASE(Geant4PlanarConversion) { G4Box boxXY("boxXY", 23., 34., 1.); auto pBoundsBox =