Skip to content

Commit

Permalink
Add tests
Browse files Browse the repository at this point in the history
  • Loading branch information
@sethrj
sethrj committed Apr 3, 2024
1 parent d72356f commit 9bf8ab0
Show file tree
Hide file tree
Showing 7 changed files with 414 additions and 44 deletions.
1 change: 1 addition & 0 deletions src/geocel/g4vg/Scaler.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
#pragma once

#include <utility>
#include <G4ThreeVector.hh>
#include <G4TwoVector.hh>

#include "geocel/detail/LengthUnits.hh"
Expand Down
22 changes: 10 additions & 12 deletions src/orange/g4org/Scaler.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,8 +8,10 @@
#pragma once

#include <utility>
#include <G4ThreeVector.hh>
#include <G4TwoVector.hh>

#include "corecel/Assert.hh"
#include "corecel/cont/Array.hh"
#include "geocel/detail/LengthUnits.hh"

Expand All @@ -20,26 +22,23 @@ namespace g4org
//---------------------------------------------------------------------------//
/*!
* Convert a unit from Geant4 scale to another.
*
* Currently the scale is hardcoded as mm (i.e., CLHEP units) but could easily
* be a class attribute.
*/
class Scaler
{
public:
//!@{
//! \name Type aliases
using Real3 = Array<double, 3>;
//!@}
//! Default scale to CLHEP units (mm)
Scaler() : scale_{celeritas::lengthunits::millimeter} {}

//! Scale with an explicit factor, probably for testing
explicit Scaler(double sc) : scale_{sc} { CELER_EXPECT(scale_ > 0); }

public:
//! Convert a positional scalar
double operator()(double val) const { return val * scale_; }

//! Convert a two-vector
std::pair<double, double> operator()(G4TwoVector const& vec) const
Array<double, 2> operator()(G4TwoVector const& vec) const
{
return this->to<std::pair<double, double>>(vec.x(), vec.y());
return this->to<Array<double, 2>>(vec.x(), vec.y());
}

//! Convert a three-vector
Expand All @@ -56,8 +55,7 @@ class Scaler
}

private:
inline static constexpr double scale_
= ::celeritas::lengthunits::millimeter;
double scale_;
};

//---------------------------------------------------------------------------//
Expand Down
168 changes: 139 additions & 29 deletions src/orange/g4org/SolidConverter.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,8 +24,6 @@
#include <G4GenericTrap.hh>
#include <G4Hype.hh>
#include <G4IntersectionSolid.hh>
#include <G4LogicalVolume.hh>
#include <G4LogicalVolumeStore.hh>
#include <G4Navigator.hh>
#include <G4Orb.hh>
#include <G4PVDivision.hh>
Expand All @@ -35,17 +33,14 @@
#include <G4Paraboloid.hh>
#include <G4Polycone.hh>
#include <G4Polyhedra.hh>
#include <G4PropagatorInField.hh>
#include <G4ReflectedSolid.hh>
#include <G4ReflectionFactory.hh>
#include <G4RotationMatrix.hh>
#include <G4Sphere.hh>
#include <G4SubtractionSolid.hh>
#include <G4TessellatedSolid.hh>
#include <G4Tet.hh>
#include <G4ThreeVector.hh>
#include <G4Torus.hh>
#include <G4Transform3D.hh>
#include <G4Trap.hh>
#include <G4Trd.hh>
#include <G4Tubs.hh>
Expand All @@ -63,6 +58,7 @@
#include "orange/orangeinp/CsgObject.hh"
#include "orange/orangeinp/Shape.hh"
#include "orange/orangeinp/Solid.hh"
#include "orange/orangeinp/Transformed.hh"

#include "Scaler.hh"
#include "Transformer.hh"
Expand Down Expand Up @@ -132,15 +128,15 @@ SolidEnclosedAngle get_polar_wedge(S const& solid)

//---------------------------------------------------------------------------//
template<class CR>
SPConstObject make_solid(G4VSolid const& solid,
CR&& interior,
std::optional<CR>&& excluded,
SolidEnclosedAngle&& enclosed)
auto make_solid(G4VSolid const& solid,
CR&& interior,
std::optional<CR>&& excluded,
SolidEnclosedAngle&& enclosed)
{
return Solid<CR>::from_optional(std::string{solid.GetName()},
std::move(interior),
std::move(excluded),
std::move(enclosed));
return Solid<CR>::or_shape(std::string{solid.GetName()},
std::move(interior),
std::move(excluded),
std::move(enclosed));
}

//---------------------------------------------------------------------------//
Expand Down Expand Up @@ -191,6 +187,7 @@ auto SolidConverter::convert_impl(arg_type solid_base) -> result_type
SC_TYPE_FUNC(Box , box),
SC_TYPE_FUNC(Cons , cons),
SC_TYPE_FUNC(CutTubs , cuttubs),
SC_TYPE_FUNC(DisplacedSolid , displaced),
SC_TYPE_FUNC(Ellipsoid , ellipsoid),
SC_TYPE_FUNC(EllipticalCone , ellipticalcone),
SC_TYPE_FUNC(EllipticalTube , ellipticaltube),
Expand Down Expand Up @@ -259,6 +256,20 @@ auto SolidConverter::cons(arg_type solid_base) -> result_type
solid.GetInnerRadiusPlusZ());
auto hh = scale_(solid.GetZHalfLength());

if (outer_r[0] == outer_r[1])
{
std::optional<Cylinder> inner;
if (inner_r[0] || inner_r[1])
{
inner = Cylinder{inner_r[0], hh};
}

return make_solid(solid,
Cylinder{outer_r[0], hh},
std::move(inner),
get_azimuthal_wedge(solid));
}

std::optional<Cone> inner;
if (inner_r[0] || inner_r[1])
{
Expand All @@ -278,6 +289,18 @@ auto SolidConverter::cuttubs(arg_type solid_base) -> result_type
CELER_NOT_IMPLEMENTED("cuttubs");
}

//---------------------------------------------------------------------------//
//! Convert a displaced solid
auto SolidConverter::displaced(arg_type solid_base) -> result_type
{
auto const& solid = dynamic_cast<G4DisplacedSolid const&>(solid_base);
G4VSolid* g4daughter = solid.GetConstituentMovedSolid();
CELER_ASSERT(g4daughter);
auto daughter = (*this)(*g4daughter);
return std::make_shared<Transformed>(
daughter, transform_(solid.GetDirectTransform()));
}

//---------------------------------------------------------------------------//
//! Convert an ellipsoid
auto SolidConverter::ellipsoid(arg_type solid_base) -> result_type
Expand Down Expand Up @@ -345,26 +368,40 @@ auto SolidConverter::hype(arg_type solid_base) -> result_type
//! Convert an intersection solid
auto SolidConverter::intersectionsolid(arg_type solid_base) -> result_type
{
CELER_DISCARD(solid_base);
CELER_NOT_IMPLEMENTED("intersectionsolid");
auto vols = this->make_bool_solids(
dynamic_cast<G4BooleanSolid const&>(solid_base));
return std::make_shared<AllObjects>(
std::string{solid_base.GetName()},
std::vector{std::move(vols[0]), std::move(vols[1])});
}

//---------------------------------------------------------------------------//
//! Convert an orb
auto SolidConverter::orb(arg_type solid_base) -> result_type
{
auto const& solid = dynamic_cast<G4Orb const&>(solid_base);
CELER_DISCARD(solid);
CELER_NOT_IMPLEMENTED("orb");
return make_shape<Sphere>(solid, scale_(solid.GetRadius()));
}

//---------------------------------------------------------------------------//
//! Convert a parallelepiped
auto SolidConverter::para(arg_type solid_base) -> result_type
{
auto const& solid = dynamic_cast<G4Para const&>(solid_base);
CELER_DISCARD(solid);
CELER_NOT_IMPLEMENTED("para");
#if G4VERSION_NUMBER >= 1100
double const theta = solid.GetTheta();
double const phi = solid.GetPhi();
#else
CELER_NOT_IMPLEMENTED("older geant4, don't duplicate code");
#endif
return make_shape<Parallelepiped>(
solid,
scale_.to<Real3>(solid.GetXHalfLength(),
solid.GetYHalfLength(),
solid.GetZHalfLength()),
native_value_to<Turn>(std::atan(solid.GetTanAlpha())),
native_value_to<Turn>(theta),
native_value_to<Turn>(phi));
}

//---------------------------------------------------------------------------//
Expand All @@ -382,7 +419,36 @@ auto SolidConverter::polycone(arg_type solid_base) -> result_type
{
auto const& solid = dynamic_cast<G4Polycone const&>(solid_base);
auto const& params = *solid.GetOriginalParameters();
CELER_DISCARD(params);

std::vector<double> zs(params.Num_z_planes);
std::vector<double> rmin(zs.size());
std::vector<double> rmax(zs.size());
for (auto i : range(zs.size()))
{
zs[i] = scale_(params.Z_values[i]);
rmin[i] = scale_(params.Rmin[i]);
rmax[i] = scale_(params.Rmax[i]);
}

if (zs.size() == 2)
{
// Special case: displaced cone/cylinder
double const hh = (zs[1] - zs[0]) / 2;
Translation trans{{0, 0, zs[0] - hh}};
if (rmin[0] == rmax[1])
{
// Cylinder
return std::make_shared<Transformed>(
make_shape<Cylinder>(solid, rmin[0], hh), std::move(trans));
}
else
{
return std::make_shared<Transformed>(
make_shape<Cone>(solid, Cone::Real2{rmin[0], rmin[1]}, hh),
std::move(trans));
}
}

CELER_NOT_IMPLEMENTED("polycone");
}

Expand All @@ -393,7 +459,38 @@ auto SolidConverter::polyhedra(arg_type solid_base) -> result_type
auto const& solid = dynamic_cast<G4Polyhedra const&>(solid_base);
auto const& params = *solid.GetOriginalParameters();

CELER_DISCARD(params);
// Opening angle: end - start phi
double const radius_factor
= std::cos(0.5 * params.Opening_angle / params.numSide);

std::vector<double> zs(params.Num_z_planes);
std::vector<double> rmin(zs.size());
std::vector<double> rmax(zs.size());
for (auto i : range(zs.size()))
{
zs[i] = scale_(params.Z_values[i]);
rmin[i] = scale_(params.Rmin[i]) * radius_factor;
rmax[i] = scale_(params.Rmax[i]) * radius_factor;
}

auto startphi = native_value_to<Turn>(solid.GetStartPhi());
SolidEnclosedAngle angle(
startphi, native_value_to<Turn>(solid.GetEndPhi()) - startphi);

if (zs.size() == 2)
{
// Just a prism
double const hh = (zs[1] - zs[0]) / 2;
Translation trans{{0, 0, zs[0] - hh}};
return std::make_shared<Transformed>(
make_shape<Prism>(solid,
params.numSide,
rmin.front(),
hh,
startphi.value() * params.numSide),
trans);
}

CELER_NOT_IMPLEMENTED("polyhedra");
}

Expand All @@ -419,8 +516,11 @@ auto SolidConverter::sphere(arg_type solid_base) -> result_type
//! Convert a subtraction solid
auto SolidConverter::subtractionsolid(arg_type solid_base) -> result_type
{
CELER_DISCARD(solid_base);
CELER_NOT_IMPLEMENTED("subtractionsolid");
auto vols = this->make_bool_solids(
dynamic_cast<G4BooleanSolid const&>(solid_base));
return make_subtraction(std::string{solid_base.GetName()},
std::move(vols[0]),
std::move(vols[1]));
}

//---------------------------------------------------------------------------//
Expand Down Expand Up @@ -465,7 +565,7 @@ auto SolidConverter::trd(arg_type solid_base) -> result_type
{
auto const& solid = dynamic_cast<G4Trd const&>(solid_base);
CELER_DISCARD(solid);
CELER_NOT_IMPLEMENTED("tubs ");
CELER_NOT_IMPLEMENTED("trd");
}

//---------------------------------------------------------------------------//
Expand All @@ -491,18 +591,28 @@ auto SolidConverter::tubs(arg_type solid_base) -> result_type
//! Convert a union solid
auto SolidConverter::unionsolid(arg_type solid_base) -> result_type
{
CELER_DISCARD(solid_base);
CELER_NOT_IMPLEMENTED("unionsolid");
auto vols = this->make_bool_solids(
dynamic_cast<G4BooleanSolid const&>(solid_base));
return std::make_shared<AnyObjects>(
std::string{solid_base.GetName()},
std::vector{std::move(vols[0]), std::move(vols[1])});
}

//---------------------------------------------------------------------------//
// HELPERS
//---------------------------------------------------------------------------//
//! Create daughter volumes for a boolean solid
auto SolidConverter::convert_bool_impl(G4BooleanSolid const&) -> result_type
auto SolidConverter::make_bool_solids(G4BooleanSolid const& bs)
-> Array<result_type, 2>
{
CELER_DISCARD(transform_);
CELER_NOT_IMPLEMENTED("bool");
Array<result_type, 2> result;
for (auto i : range(result.size()))
{
G4VSolid const* solid = bs.GetConstituentSolid(i);
CELER_ASSERT(solid);
result[i] = (*this)(*solid);
}
return result;
}

//---------------------------------------------------------------------------//
Expand Down
3 changes: 2 additions & 1 deletion src/orange/g4org/SolidConverter.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -67,6 +67,7 @@ class SolidConverter
result_type box(arg_type);
result_type cons(arg_type);
result_type cuttubs(arg_type);
result_type displaced(arg_type);
result_type ellipsoid(arg_type);
result_type ellipticalcone(arg_type);
result_type ellipticaltube(arg_type);
Expand All @@ -92,7 +93,7 @@ class SolidConverter
result_type unionsolid(arg_type);

// Construct bool daughters
result_type convert_bool_impl(G4BooleanSolid const&);
Array<result_type, 2> make_bool_solids(G4BooleanSolid const&);
// Calculate solid capacity in native celeritas units
double calc_capacity(G4VSolid const&) const;
};
Expand Down
4 changes: 2 additions & 2 deletions src/orange/orangeinp/ConvexRegion.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -290,8 +290,8 @@ class InfWedge final : public ConvexRegionInterface
* The shape has the center in the origin and it is defined by:
*
* - `halfedges:` a 3-vector (dY, dY, dZ) with half-lengths of the
* projections of the edges on X, Y, Z. The lower Z face is positioned at
* `-dZ`, and the upper one at `+dZ`.
* projections of the edges on X, Y, Z. The lower Z face is positioned at
* `-dZ`, and the upper one at `+dZ`.
* - `alpha:` angle between the segment defined by the centers of the
* X-parallel edges and Y axis. Validity range is `(-1/4, 1/4)`;
* - `theta:` polar angle of the shape's main axis, e.g. the segment defined
Expand Down
1 change: 1 addition & 0 deletions test/orange/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -132,6 +132,7 @@ if(CELERITAS_USE_Geant4 AND CELERITAS_REAL_TYPE STREQUAL "double")
list(APPEND CELERITASTEST_LINK_LIBRARIES ${_g4_geo_libs})

celeritas_add_test(g4org/Transformer.test.cc)
celeritas_add_test(g4org/SolidConverter.test.cc)
endif()

#-----------------------------------------------------------------------------#

0 comments on commit 9bf8ab0

Please sign in to comment.