Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

feat: Adding G4Trap converter in Geant4Converters #3775

Open
wants to merge 9 commits into
base: main
Choose a base branch
from
Open

feat: Adding G4Trap converter in Geant4Converters #3775

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
0d3d440
Adding Trap converter in Geant4Converters
Oct 23, 2024
f04d7fa
Adding Trap converter in Geant4Converters
Oct 23, 2024
495dd67
Adding Trap converter in Geant4Converters
Oct 23, 2024
71d1646
Adding Trap converter in Geant4Converters
Oct 23, 2024
aff6611
Adding Trap converter in Geant4Converters
Oct 23, 2024
f16ac93
Adding G4Trap converter in Geant4Converters and its test-case
Oct 25, 2024
d1bc9fe
Adding G4Trap converter in Geant4Converters and its test-case
Nov 5, 2024
4bd4cbc
Adding G4Trap converter in Geant4Converters and its test-case
Dec 5, 2024
3eb26c5
Adding G4Trap converter in Geant4Converters and its test-case
Dec 5, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 8 additions & 0 deletions Plugins/Geant4/include/Acts/Plugins/Geant4/Geant4Converters.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -133,6 +133,14 @@ struct Geant4ShapeConverter {
std::tuple<std::shared_ptr<TrapezoidBounds>, std::array<int, 2u>, ActsScalar>
trapezoidBounds(const G4Trd& g4Trd);

/// @brief Convert to trapezoid bounds - from Trap
///
/// @param g4Trap a Geant4 trapezoid shape
///
/// @return an ACTS Trapezoid bounds object, axis orientation, and thickness
std::tuple<std::shared_ptr<TrapezoidBounds>, std::array<int, 2u>, ActsScalar>
trapezoidBounds(const G4Trap& g4Trap);

/// @brief Convert to general solid into a planar shape
///
/// @param g4Solid a Geant4 solid shape
Expand Down
86 changes: 86 additions & 0 deletions Plugins/Geant4/src/Geant4Converters.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,7 @@
#include "G4RotationMatrix.hh"
#include "G4ThreeVector.hh"
#include "G4Transform3D.hh"
#include "G4Trap.hh"
asalzburger marked this conversation as resolved.
Show resolved Hide resolved
#include "G4Trd.hh"
#include "G4Tubs.hh"
#include "G4VPhysicalVolume.hh"
Expand Down Expand Up @@ -234,6 +235,74 @@ Acts::Geant4ShapeConverter::trapezoidBounds(const G4Trd& g4Trd) {
return std::make_tuple(std::move(tBounds), rAxes, thickness);
}

std::tuple<std::shared_ptr<Acts::TrapezoidBounds>, std::array<int, 2u>,
Acts::ActsScalar>
Acts::Geant4ShapeConverter::trapezoidBounds(const G4Trap& g4Trap) {
// The conversion considers:
// - Theta = Tilt angle of the trapezoid along the Y-axis
// - Phi = Rotation angle around the Z-axis to determine the Y-axis
// - Half-lengths are adjusted based on these angles to maintain the shape
ActsScalar y1 = static_cast<ActsScalar>(g4Trap.GetYHalfLength1());
ActsScalar y2 = static_cast<ActsScalar>(g4Trap.GetYHalfLength2());
ActsScalar x1 = static_cast<ActsScalar>(g4Trap.GetXHalfLength1());
ActsScalar x2 = static_cast<ActsScalar>(g4Trap.GetXHalfLength2());
ActsScalar x3 = static_cast<ActsScalar>(g4Trap.GetXHalfLength3());
ActsScalar x4 = static_cast<ActsScalar>(g4Trap.GetXHalfLength4());
ActsScalar phi = static_cast<ActsScalar>(g4Trap.GetPhi());
ActsScalar theta = static_cast<ActsScalar>(g4Trap.GetTheta());
ActsScalar z = static_cast<ActsScalar>(g4Trap.GetZHalfLength());

Comment on lines +245 to +254
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

💡 Codebase verification

Validate the parameters, we must. Prevent the dark side of undefined behavior, this will.

For G4Trap parameters, essential validation is, hmm:

  • Negative half-lengths, forbidden they are
  • Theta angle, between 0 and π it must be
  • Phi angle, between -π and π it must lie
 ActsScalar y1 = static_cast<ActsScalar>(g4Trap.GetYHalfLength1());
+if (y1 < 0 || y2 < 0 || x1 < 0 || x2 < 0 || x3 < 0 || x4 < 0 || z < 0) {
+  throw std::invalid_argument("Negative half-lengths in G4Trap not allowed, they are");
+}
+if (theta < 0 || theta > M_PI) {
+  throw std::invalid_argument("Theta angle must between 0 and π be");
+}
+if (phi < -M_PI || phi > M_PI) {
+  throw std::invalid_argument("Phi angle must between -π and π be");
+}

From the test cases seen, validate these parameters we must. Prevent invalid geometric configurations that lead to dark undefined behavior, this will. Strong with the Force, this validation is.

🔗 Analysis chain

Consider adding validation for input parameters, we must.

Validate we must that input parameters from G4Trap are within expected ranges. Prevent undefined behavior in calculations, this will.

 ActsScalar y1 = static_cast<ActsScalar>(g4Trap.GetYHalfLength1());
+if (y1 < 0 || y2 < 0 || x1 < 0 || x2 < 0 || x3 < 0 || x4 < 0 || z < 0) {
+  throw std::invalid_argument("Negative half-lengths in G4Trap are not allowed");
+}
+if (!std::isfinite(theta) || !std::isfinite(phi)) {
+  throw std::invalid_argument("Invalid angles in G4Trap");
+}

Let's search for existing validation patterns and trap-related checks in the codebase, hmm? Help us understand the current practices, this will.

Hmm, deeper into the codebase we must look. Validation patterns for G4Trap and similar geometric conversions, seek we shall.

🏁 Scripts executed

The following scripts were executed for the analysis:

Script:

#!/bin/bash
# Search for existing parameter validation patterns in the codebase
rg -A 2 "GetTheta|GetPhi" --type cpp

Length of output: 1168

Script:

#!/bin/bash
# Search for validation patterns in geometric converters
rg -A 3 "invalid_argument|isfinite|epsilon" --type cpp

# Look for other trap-related validations
rg -A 3 "Trap.*Half" --type cpp

Length of output: 69663

Script:

#!/bin/bash
# Search for existing validation patterns in Geant4 converters
rg -A 5 "G4Trap" --type cpp

# Look for existing validation patterns in similar geometric converters
rg -A 5 "static_cast.*GetTheta|static_cast.*GetPhi" --type cpp

Length of output: 6629

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;
Comment on lines +255 to +259
Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

💡 Codebase verification

Safeguards for numerical stability, add we must.

Through the Force, I sense many trigonometric calculations in the codebase, but no explicit checks for numerical stability in this context, I do. A pattern of using std::tan without validation, widespread it is. Protect against infinite or undefined results, we should.

  • In Geant4Converters.cpp, add validation for theta angle before tangent calculation:
const ActsScalar tanTheta = std::tan(theta);
if (!std::isfinite(tanTheta)) {
    throw std::runtime_error("Invalid theta angle in G4Trap resulting in infinite tangent");
}
ActsScalar hlX1 = 2 * z * tanTheta * std::cos(phi) + (x3 + x4) * 0.5;
ActsScalar hlY1 = y2 + 2 * z * tanTheta * std::sin(phi);
🔗 Analysis chain

Careful with trigonometric calculations, we must be.

Potential numerical instability in calculations with extreme angles, there might be. Consider adding safeguards against numerical overflow.

-ActsScalar hlX1 = 2 * z * std::tan(theta) * std::cos(phi) + (x3 + x4) * 0.5;
+const ActsScalar tanTheta = std::tan(theta);
+if (!std::isfinite(tanTheta)) {
+  throw std::runtime_error("Invalid theta angle in G4Trap resulting in infinite tangent");
+}
+ActsScalar hlX1 = 2 * z * tanTheta * std::cos(phi) + (x3 + x4) * 0.5;
🏁 Scripts executed

The following scripts were executed for the analysis:

Script:

#!/bin/bash
# Search for similar trigonometric calculations in the codebase
ast-grep --pattern 'std::tan($x)'

Length of output: 35

Script:

#!/bin/bash
# Search for trigonometric calculations with pattern variations
rg "std::tan|tan\(" -A 2 -B 2

# Also search for angle-related variable names and calculations
rg "theta|angle|phi" -A 2 -B 2

# Look for error handling around trigonometric calculations
rg "isfinite|isinf|isnan" -A 2 -B 2

Length of output: 70292


std::vector<ActsScalar> 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<int, 2u> 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<TrapezoidBounds>(
halfLengthXminY, halfLengthXmaxY, halfLengthY);
return std::make_tuple(std::move(tBounds), rAxes, thickness);
}

std::tuple<std::shared_ptr<Acts::PlanarBounds>, std::array<int, 2u>,
Acts::ActsScalar>
Acts::Geant4ShapeConverter::planarBounds(const G4VSolid& g4Solid) {
Expand Down Expand Up @@ -332,6 +401,23 @@ std::shared_ptr<Acts::Surface> Acts::Geant4PhysicalVolumeConverter::surface(
}
}

// Into a Trapezoid (G4Trap)
auto g4Trap = dynamic_cast<const G4Trap*>(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<PlaneSurface>(orientedToGlobal,
std::move(bounds));
assignMaterial(*surface.get(), original, compressed);
return surface;
} else {
throw std::runtime_error("Cannot convert 'G4Trap' into forced shape.");
}
}

// Into a Cylinder, disc or line
auto g4Tubs = dynamic_cast<const G4Tubs*>(g4Solid);
if (g4Tubs != nullptr) {
Expand Down
102 changes: 88 additions & 14 deletions Tests/UnitTests/Plugins/Geant4/Geant4ConvertersTests.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,7 @@
#include "Acts/Surfaces/RadialBounds.hpp"
#include "Acts/Surfaces/RectangleBounds.hpp"
#include "Acts/Surfaces/Surface.hpp"
#include "Acts/Surfaces/TrapezoidBounds.hpp"
#include "Acts/Surfaces/trapezoidBounds.hpp"
#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"

#include <array>
Expand Down Expand Up @@ -146,13 +146,13 @@ BOOST_AUTO_TEST_CASE(Geant4TrapzoidConversion) {
auto [boundsXY, axesXY, thicknessZ] =
Acts::Geant4ShapeConverter{}.trapezoidBounds(trdXY);
CHECK_CLOSE_ABS(
boundsXY->get(Acts::TrapezoidBounds::BoundValues::eHalfLengthXnegY), 100,
boundsXY->get(Acts::trapezoidBounds::BoundValues::eHalfLengthXnegY), 100,
10e-10);
CHECK_CLOSE_ABS(
boundsXY->get(Acts::TrapezoidBounds::BoundValues::eHalfLengthXposY), 150,
boundsXY->get(Acts::trapezoidBounds::BoundValues::eHalfLengthXposY), 150,
10e-10);
CHECK_CLOSE_ABS(
boundsXY->get(Acts::TrapezoidBounds::BoundValues::eHalfLengthY), 200,
boundsXY->get(Acts::trapezoidBounds::BoundValues::eHalfLengthY), 200,
10e-10);
auto refXY = std::array<int, 2u>{0, 1};
BOOST_CHECK(axesXY == refXY);
Expand All @@ -163,13 +163,13 @@ BOOST_AUTO_TEST_CASE(Geant4TrapzoidConversion) {
auto [boundsyX, axesyX, thicknessZ2] =
Acts::Geant4ShapeConverter{}.trapezoidBounds(trdyX);
CHECK_CLOSE_ABS(
boundsyX->get(Acts::TrapezoidBounds::BoundValues::eHalfLengthXnegY), 100,
boundsyX->get(Acts::trapezoidBounds::BoundValues::eHalfLengthXnegY), 100,
10e-10);
CHECK_CLOSE_ABS(
boundsyX->get(Acts::TrapezoidBounds::BoundValues::eHalfLengthXposY), 150,
boundsyX->get(Acts::trapezoidBounds::BoundValues::eHalfLengthXposY), 150,
10e-10);
CHECK_CLOSE_ABS(
boundsyX->get(Acts::TrapezoidBounds::BoundValues::eHalfLengthY), 200,
boundsyX->get(Acts::trapezoidBounds::BoundValues::eHalfLengthY), 200,
10e-10);
auto refyX = std::array<int, 2u>{-1, 0};
BOOST_CHECK(axesyX == refyX);
Expand All @@ -180,13 +180,13 @@ BOOST_AUTO_TEST_CASE(Geant4TrapzoidConversion) {
auto [boundsYZ, axesYZ, thicknessX] =
Acts::Geant4ShapeConverter{}.trapezoidBounds(trdYZ);
CHECK_CLOSE_ABS(
boundsYZ->get(Acts::TrapezoidBounds::BoundValues::eHalfLengthXnegY), 120.,
boundsYZ->get(Acts::trapezoidBounds::BoundValues::eHalfLengthXnegY), 120.,
10e-10);
CHECK_CLOSE_ABS(
boundsYZ->get(Acts::TrapezoidBounds::BoundValues::eHalfLengthXposY), 140.,
boundsYZ->get(Acts::trapezoidBounds::BoundValues::eHalfLengthXposY), 140.,
10e-10);
CHECK_CLOSE_ABS(
boundsYZ->get(Acts::TrapezoidBounds::BoundValues::eHalfLengthY), 200.,
boundsYZ->get(Acts::trapezoidBounds::BoundValues::eHalfLengthY), 200.,
10e-10);
auto refYZ = std::array<int, 2u>{1, 2};
BOOST_CHECK(axesYZ == refYZ);
Expand All @@ -197,19 +197,93 @@ BOOST_AUTO_TEST_CASE(Geant4TrapzoidConversion) {
auto [boundsXz, axesXz, thicknessY] =
Acts::Geant4ShapeConverter{}.trapezoidBounds(trdXz);
CHECK_CLOSE_ABS(
boundsXz->get(Acts::TrapezoidBounds::BoundValues::eHalfLengthXnegY), 50.,
boundsXz->get(Acts::trapezoidBounds::BoundValues::eHalfLengthXnegY), 50.,
10e-10);
CHECK_CLOSE_ABS(
boundsXz->get(Acts::TrapezoidBounds::BoundValues::eHalfLengthXposY), 75.,
boundsXz->get(Acts::trapezoidBounds::BoundValues::eHalfLengthXposY), 75.,
10e-10);
CHECK_CLOSE_ABS(
boundsXz->get(Acts::TrapezoidBounds::BoundValues::eHalfLengthY), 200.,
boundsXz->get(Acts::trapezoidBounds::BoundValues::eHalfLengthY), 200.,
10e-10);
auto refXz = std::array<int, 2u>{0, -2};
BOOST_CHECK(axesXz == refXz);
CHECK_CLOSE_ABS(thicknessY, 2., 10e-10);
}

BOOST_AUTO_TEST_CASE(Geant4TrapConversion_G4Trap) {
// 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<int, 2u>{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<int, 2u>{-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<int, 2u>{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<int, 2u>{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 =
Expand All @@ -220,7 +294,7 @@ BOOST_AUTO_TEST_CASE(Geant4PlanarConversion) {
G4Trd trdXY("trdXY", 100, 150, 200, 200, 2);
auto pBoundsTrd =
std::get<0u>(Acts::Geant4ShapeConverter{}.planarBounds(trdXY));
auto tBounds = dynamic_cast<const Acts::TrapezoidBounds*>(pBoundsTrd.get());
auto tBounds = dynamic_cast<const Acts::trapezoidBounds*>(pBoundsTrd.get());
BOOST_CHECK_NE(tBounds, nullptr);
}

Expand Down
Loading