Skip to content

Commit

Permalink
Add ORANGE converter for G4Trd (#1218)
Browse files Browse the repository at this point in the history 
* Add ORANGE converter for G4Trd

* Add a TRD-like GenTrap constructor for reusability

* Add helper functions to calculate corners for GenTrap construction

The helper functions have arguments that correspond to the main G4Trap
and G4Trd constructors, but arranges those parameters into z-face
rectangles or trapezes.

Testing cases were added for the new helper functions.

* Apply reviewer suggestions

* Fix comparison tolerance for planarity

* Rearrange constructors for consistency

* Remove failing test to be deleted soon
  • Loading branch information
@mrguilima
mrguilima committed May 7, 2024
1 parent 6684723 commit a8c69d2
Show file tree
Hide file tree
Showing 7 changed files with 277 additions and 86 deletions.
47 changes: 24 additions & 23 deletions src/orange/g4org/SolidConverter.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -609,55 +609,56 @@ auto SolidConverter::trap(arg_type solid_base) -> result_type
{
auto const& solid = dynamic_cast<G4Trap const&>(solid_base);

real_type sin_phi{}, cos_phi{}, tan_theta{};
real_type tan_theta{};
Turn phi{};
#if G4VERSION_NUMBER < 1100
// Geant4 10.7 and earlier - axis = (sinth.cosphi, sinth.sinphi, costh)
// Note: both sinth,costh >= 0 since theta is in [0, pi/2] for a G4Trap
auto axis = solid.GetSymAxis();
auto sin_theta = std::sqrt(real_type(1.0) - axis.z() * axis.z());
auto sin_theta = std::max(0.0, std::sqrt(1.0 - ipow<2>(axis.z())));
tan_theta = SoftZero<double>{1.e-8}(axis.z())
? sin_theta / axis.z()
: numeric_limits<real_type>::infinity();
cos_phi = sin_theta > 0 ? axis.x() / sin_theta : 1.0;
sin_phi = sin_theta > 0 ? axis.y() / sin_theta : 1.0;
real_type cos_phi = sin_theta > 0 ? axis.x() / sin_theta : 1.0;
real_type sin_phi = sin_theta > 0 ? axis.y() / sin_theta : 1.0;
phi = native_value_to<Turn>(std::atan2(sin_phi, cos_phi));
#else
// Geant4 11 and later
sincos(solid.GetPhi(), &sin_phi, &cos_phi);
tan_theta = std::tan(solid.GetTheta());
phi = native_value_to<Turn>(solid.GetPhi());
#endif

auto hz = scale_(solid.GetZHalfLength());
auto hy1 = scale_(solid.GetYHalfLength1());
auto hy2 = scale_(solid.GetYHalfLength2());
auto hx1 = scale_(solid.GetXHalfLength1());
auto hx2 = scale_(solid.GetXHalfLength2());
auto hy2 = scale_(solid.GetYHalfLength2());
auto hx3 = scale_(solid.GetXHalfLength3());
auto hx4 = scale_(solid.GetXHalfLength4());
auto dxdz_hz = tan_theta * cos_phi * hz;
auto dydz_hz = tan_theta * sin_phi * hz;
auto dxdy_hy1 = solid.GetTanAlpha1() * hy1;
auto dxdy_hy2 = solid.GetTanAlpha2() * hy2;

std::vector<GenTrap::Real2> lower(4), upper(4);
lower[0] = GenTrap::Real2{-dxdz_hz - dxdy_hy1 - hx1, -dydz_hz - hy1};
lower[1] = GenTrap::Real2{-dxdz_hz - dxdy_hy1 + hx1, -dydz_hz - hy1};
lower[2] = GenTrap::Real2{-dxdz_hz + dxdy_hy1 + hx2, -dydz_hz + hy1};
lower[3] = GenTrap::Real2{-dxdz_hz + dxdy_hy1 - hx2, -dydz_hz + hy1};
upper[0] = GenTrap::Real2{+dxdz_hz - dxdy_hy2 - hx3, +dydz_hz - hy2};
upper[1] = GenTrap::Real2{+dxdz_hz - dxdy_hy2 + hx3, +dydz_hz - hy2};
upper[2] = GenTrap::Real2{+dxdz_hz + dxdy_hy2 + hx4, +dydz_hz + hy2};
upper[3] = GenTrap::Real2{+dxdz_hz + dxdy_hy2 - hx4, +dydz_hz + hy2};

return make_shape<GenTrap>(solid, hz, lower, upper);
return make_shape<GenTrap>(
solid,
GenTrap::from_trap(hz,
tan_theta,
phi,
{hy1, hx1, hx2, solid.GetTanAlpha1()},
{hy2, hx3, hx4, solid.GetTanAlpha2()}));
}

//---------------------------------------------------------------------------//
//! Convert a simple trapezoid
auto SolidConverter::trd(arg_type solid_base) -> result_type
{
auto const& solid = dynamic_cast<G4Trd const&>(solid_base);
CELER_DISCARD(solid);
CELER_NOT_IMPLEMENTED("trd");

auto hz = scale_(solid.GetZHalfLength());
auto hy1 = scale_(solid.GetYHalfLength1());
auto hy2 = scale_(solid.GetYHalfLength2());
auto hx1 = scale_(solid.GetXHalfLength1());
auto hx2 = scale_(solid.GetXHalfLength2());

return make_shape<GenTrap>(solid,
GenTrap::from_trd(hz, {hx1, hy1}, {hx2, hy2}));
}

//---------------------------------------------------------------------------//
Expand Down
76 changes: 76 additions & 0 deletions src/orange/orangeinp/ConvexRegion.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -314,6 +314,82 @@ void Ellipsoid::output(JsonPimpl* j) const

//---------------------------------------------------------------------------//
// GENTRAP
//---------------------------------------------------------------------------//
/*!
* Construct a TRD from 5 half-lengths: one half-height, {hx1,hy1} for a in -z,
* and {hx2,hy2} in +z
*/
GenTrap GenTrap::from_trd(real_type halfz, Real2 const& lo, Real2 const& hi)
{
CELER_VALIDATE(lo[0] > 0, << "nonpositive lower x half-edge: " << lo[0]);
CELER_VALIDATE(hi[0] > 0, << "nonpositive upper x half-edge: " << hi[0]);
CELER_VALIDATE(lo[1] > 0, << "nonpositive lower y half-edge: " << lo[1]);
CELER_VALIDATE(hi[1] > 0, << "nonpositive upper y half-edge: " << hi[1]);
CELER_VALIDATE(halfz > 0, << "nonpositive half-height: " << halfz);

// Construct points counterclockwise from lower left
VecReal2 lower
= {{-lo[0], -lo[1]}, {lo[0], -lo[1]}, {lo[0], lo[1]}, {-lo[0], lo[1]}};
VecReal2 upper
= {{-hi[0], -hi[1]}, {hi[0], -hi[1]}, {hi[0], hi[1]}, {-hi[0], hi[1]}};

return GenTrap{halfz, std::move(lower), std::move(upper)};
}

//---------------------------------------------------------------------------//
/*!
* Construct from half Z height and 1-4 vertices for top and bottom planes.
*/
GenTrap GenTrap::from_trap(real_type hz,
real_type tan_theta,
Turn const& phi,
TrapFace const& lo,
TrapFace const& hi)
{
// TrapFace is validated in its constructor
CELER_VALIDATE(hz > 0, << "nonpositive half-height: " << hz);

CELER_VALIDATE(tan_theta >= 0, << "negative tan(theta): " << tan_theta);

CELER_VALIDATE(phi >= zero_quantity() && phi < Turn{1.},
<< "invalid angle " << phi.value()
<< " [turns]: must be in the range [0, 1)");

for (TrapFace const* tf : {&lo, &hi})
{
CELER_VALIDATE(tf->hx_lo_ > 0,
<< "nonpositive lower x half-edge: " << tf->hx_lo_);
CELER_VALIDATE(tf->hx_hi_ > 0,
<< "nonpositive upper x half-edge: " << tf->hx_hi_);
CELER_VALIDATE(tf->hy_ > 0,
<< "nonpositive y half-distance: " << tf->hy_);
}
real_type cos_phi{}, sin_phi{};
sincos(phi, &sin_phi, &cos_phi);
auto dxdz_hz = tan_theta * cos_phi * hz;
auto dydz_hz = tan_theta * sin_phi * hz;

auto hy1 = lo.hy_;
auto hx1 = lo.hx_lo_;
auto hx2 = lo.hx_hi_;
auto hy2 = hi.hy_;
auto hx3 = hi.hx_lo_;
auto hx4 = hi.hx_hi_;
auto dxdy_hy1 = lo.tan_alpha_ * lo.hy_;
auto dxdy_hy2 = hi.tan_alpha_ * hi.hy_;

VecReal2 lower = {{-dxdz_hz - dxdy_hy1 - hx1, -dydz_hz - hy1},
{-dxdz_hz - dxdy_hy1 + hx1, -dydz_hz - hy1},
{-dxdz_hz + dxdy_hy1 + hx2, -dydz_hz + hy1},
{-dxdz_hz + dxdy_hy1 - hx2, -dydz_hz + hy1}};
VecReal2 upper = {{+dxdz_hz - dxdy_hy2 - hx3, +dydz_hz - hy2},
{+dxdz_hz - dxdy_hy2 + hx3, +dydz_hz - hy2},
{+dxdz_hz + dxdy_hy2 + hx4, +dydz_hz + hy2},
{+dxdz_hz + dxdy_hy2 - hx4, +dydz_hz + hy2}};

return GenTrap{hz, std::move(lower), std::move(upper)};
}

//---------------------------------------------------------------------------//
/*!
* Construct from half Z height and 1-4 vertices for top and bottom planes.
Expand Down
28 changes: 27 additions & 1 deletion src/orange/orangeinp/ConvexRegion.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -216,8 +216,34 @@ class GenTrap final : public ConvexRegionInterface
using VecReal2 = std::vector<Real2>;
//!@}

//! Argument for building from regular trapezoidal top/bottom polygons
struct TrapFace
{
real_type hy_{}; //!< half the vertical distance between horizontal
//!< edges
real_type hx_lo_{}; //!< top horizontal edge half-length
real_type hx_hi_{}; //!< botton horizontal edge half-length

// tan(alpha), where alpha is the clockwise angle between the
// _centers_ of horizontal edges, with respect to the vertical
// (alpha=0)
real_type tan_alpha_{};
};

public:
// Construct from half Z height and 1-4 vertices for top and bottom planes
// Helper function to construct a Trd shape from hz and two rectangles,
// one for each z-face
static GenTrap from_trd(real_type halfz, Real2 const& lo, Real2 const& hi);

// Helper function to construct a general trap from its half-height and
// the two trapezoids defining its lower and upper faces
static GenTrap from_trap(real_type hz,
real_type tan_theta,
Turn const& phi,
TrapFace const& lo,
TrapFace const& hi);

// Construct from half Z height and 4 vertices for top and bottom planes
GenTrap(real_type halfz, VecReal2 const& lo, VecReal2 const& hi);

// Build surfaces
Expand Down
4 changes: 2 additions & 2 deletions src/orange/orangeinp/detail/PolygonUtils.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,8 +122,8 @@ is_planar(Real3 const& a, Real3 const& b, Real3 const& c, Real3 const& d)
auto norm = make_unit_vector(cross_product(b - a, c - a));
auto val = dot_product(norm, d - a);

// FIXME: SoftEqual and SoftZero should have rel = abs
return SoftZero{SoftEqual<>{}.rel()}(val);
return SoftZero{
::celeritas::detail::SoftEqualTraits<real_type>::sqrt_prec()}(val);
}

//---------------------------------------------------------------------------//
Expand Down
69 changes: 48 additions & 21 deletions test/orange/g4org/SolidConverter.test.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -398,49 +398,76 @@ TEST_F(SolidConverterTest, subtractionsolid)
{{0, 0, 0}, {0, 0, 10}, {1, 0, 0}, {0, 1, 0}, {0, 0, 1}});
}

TEST_F(SolidConverterTest, trap)
TEST_F(SolidConverterTest, trap_box)
{
{
this->build_and_test(
G4Trap("boxTrap", 30, 0, 0, 20, 10, 10, 0, 20, 10, 10, 0),
G4Trap("trap_box", 30, 0, 0, 20, 10, 10, 0, 20, 10, 10, 0),
R"json({"_type":"shape","interior":{"_type":"gentrap",
"halfheight":3.0,
"lower":[[-1.0,-2.0],[1.0,-2.0],[1.0,2.0],[-1.0,2.0]],
"upper":[[-1.0,-2.0],[1.0,-2.0],[1.0,2.0],[-1.0,2.0]]},
"label":"boxTrap"})json",
"label":"trap_box"})json",
{{-1, -2, -3}, {1, 2, 3}, {1, 2, 4}});
}
}

TEST_F(SolidConverterTest, trap_trd)
{
{
this->build_and_test(
G4Trap("trdGenTrap", 50, 100, 100, 200, 300),
G4Trap("trap_trd", 50, 100, 100, 200, 300),
R"json({"_type":"shape","interior":{"_type":"gentrap",
"halfheight":30.0,
"lower":[[-5.0,-10.0],[5.0,-10.0],[5.0,10.0],[-5.0,10.0]],
"upper":[[-10.0,-20.0],[10.0,-20.0],[10.0,20.0],[-10.0,20.0]]},
"label":"trdGenTrap"})json",
"label":"trap_trd"})json",
{{-10, -20, -40},
{-10, -20, -30 + 1.e-6},
{5, 10, 30},
{10, 10, 30}});
}
}

{
this->build_and_test(
G4Trap("trap", 40, 0.02, 0.05, 20, 10, 10, 0.01, 30, 15, 15, 0.01),
R"json({"_type":"shape","interior":{"_type":"gentrap",
"halfheight":4.0,
"lower":[[-1.0999113425658524,-2.0039988667381046],
[0.9000886574341476,-2.0039988667381046],
[0.9400899908208165,1.9960011332618952],
[-1.0599100091791835,1.9960011332618952]],
"upper":[[-1.4500903241674836,-2.9960011332618954],
[1.5499096758325164,-2.9960011332618954],
[1.6099116759125196,3.0039988667381046],
[-1.3900883240874804,3.0039988667381046]]},
"label":"trap"})json",
{{-1, -2, -4 - 1.e-6}, {-1, -2, -3}, {0.5, 1, 3}, {1, 1, 3}});
}
TEST_F(SolidConverterTest, trap)
{
this->build_and_test(
G4Trap("trap", 40, 0.02, 0.05, 20, 10, 10, 0.01, 30, 15, 15, 0.01),
R"json({"_type":"shape","interior":{"_type":"gentrap",
"halfheight":4.0,
"lower":[[-1.0999113425658524,-2.0039988667381046],
[0.9000886574341476,-2.0039988667381046],
[0.9400899908208165,1.9960011332618952],
[-1.0599100091791835,1.9960011332618952]],
"upper":[[-1.4500903241674836,-2.9960011332618954],
[1.5499096758325164,-2.9960011332618954],
[1.6099116759125196,3.0039988667381046],
[-1.3900883240874804,3.0039988667381046]]},
"label":"trap"})json",
{{-1, -2, -4 - 1.e-6}, {-1, -2, -3}, {0.5, 1, 3}, {1, 1, 3}});
}

TEST_F(SolidConverterTest, trd_box)
{
this->build_and_test(G4Trd("trd_box", 10, 10, 20, 20, 30),
R"json({"_type":"shape","interior":{"_type":"gentrap",
"halfheight":3.0,
"lower":[[-1.0,-2.0],[1.0,-2.0],[1.0,2.0],[-1.0,2.0]],
"upper":[[-1.0,-2.0],[1.0,-2.0],[1.0,2.0],[-1.0,2.0]]},
"label":"trd_box"})json",
{{-1, -2, -3}, {1, 2, 3}, {1, 2, 4}});
}

TEST_F(SolidConverterTest, trd)
{
this->build_and_test(
G4Trd("trd", 50, 100, 100, 200, 300),
R"json({"_type":"shape","interior":{"_type":"gentrap",
"halfheight":30.0,
"lower":[[-5.0,-10.0],[5.0,-10.0],[5.0,10.0],[-5.0,10.0]],
"upper":[[-10.0,-20.0],[10.0,-20.0],[10.0,20.0],[-10.0,20.0]]},
"label":"trd"})json",
{{-10, -20, -40}, {-10, -20, -30 + 1.e-6}, {5, 10, 30}, {10, 10, 30}});
}

TEST_F(SolidConverterTest, tubs)
Expand Down

0 comments on commit a8c69d2

Please sign in to comment.