Skip to content

Commit

Permalink
feat: Geant4 line/straw tube conversion capability (#1811)
Browse files Browse the repository at this point in the history 
This PR is a prerequisite for the the MS mockup station, it allows to do
straw tube creation from Geant4.
  • Loading branch information
@asalzburger
asalzburger committed Jan 31, 2023
1 parent c6fcdb8 commit b1e5cf2
Show file tree
Hide file tree
Showing 3 changed files with 229 additions and 25 deletions.
14 changes: 12 additions & 2 deletions Plugins/Geant4/include/Acts/Plugins/Geant4/Geant4Converters.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@

#include "Acts/Definitions/Algebra.hpp"
#include "Acts/Definitions/Common.hpp"
#include "Acts/Surfaces/Surface.hpp"

#include <memory>
#include <tuple>
Expand Down Expand Up @@ -62,6 +63,7 @@ class RadialBounds;
class RectangleBounds;
class TrapezoidBounds;
class PlanarBounds;
class LineBounds;

// The following set of converters convert a Geant4 volume shape
// to an ACTS surface bounds object, this is for converting the volume
Expand Down Expand Up @@ -93,6 +95,13 @@ struct Geant4ShapeConverter {
std::tuple<std::shared_ptr<RadialBounds>, ActsScalar> radialBounds(
const G4Tubs& g4Tubs);

/// @brief Convert to line/straw bounds
///
/// @param g4Tubs a Geant4 tube shape
///
/// @return an Acts line bounds object and thickness
std::shared_ptr<LineBounds> lineBounds(const G4Tubs& g4Tubs);

/// @brief Convert to rectangle bounds
///
/// @param g4Box a Geant4 box shape
Expand All @@ -119,9 +128,10 @@ struct Geant4ShapeConverter {
planarBounds(const G4VSolid& g4Solid);
};

class Surface;

struct Geant4PhysicalVolumeConverter {
/// Optionally allow to foce a type, throws exception if not possbile
Surface::SurfaceType forcedType = Surface::SurfaceType::Other;

/// @brief Convert a Geant4 phsyical volume to a surface
///
/// @param g4PhysVol the physical volume to be constructed
Expand Down
66 changes: 45 additions & 21 deletions Plugins/Geant4/src/Geant4Converters.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,10 +15,11 @@
#include "Acts/Surfaces/CylinderBounds.hpp"
#include "Acts/Surfaces/CylinderSurface.hpp"
#include "Acts/Surfaces/DiscSurface.hpp"
#include "Acts/Surfaces/LineBounds.hpp"
#include "Acts/Surfaces/PlaneSurface.hpp"
#include "Acts/Surfaces/RadialBounds.hpp"
#include "Acts/Surfaces/RectangleBounds.hpp"
#include "Acts/Surfaces/Surface.hpp"
#include "Acts/Surfaces/StrawSurface.hpp"
#include "Acts/Surfaces/TrapezoidBounds.hpp"
#include "Acts/Utilities/Helpers.hpp"

Expand Down Expand Up @@ -95,6 +96,13 @@ Acts::Geant4ShapeConverter::radialBounds(const G4Tubs& g4Tubs) {
return std::tie(rBounds, thickness);
}

std::shared_ptr<Acts::LineBounds> Acts::Geant4ShapeConverter::lineBounds(
const G4Tubs& g4Tubs) {
auto r = static_cast<ActsScalar>(g4Tubs.GetOuterRadius());
auto hlZ = static_cast<ActsScalar>(g4Tubs.GetZHalfLength());
return std::make_shared<LineBounds>(r, hlZ);
}

std::tuple<std::shared_ptr<Acts::RectangleBounds>, std::array<int, 2u>,
Acts::ActsScalar>
Acts::Geant4ShapeConverter::rectangleBounds(const G4Box& g4Box) {
Expand Down Expand Up @@ -250,47 +258,63 @@ std::shared_ptr<Acts::Surface> Acts::Geant4PhysicalVolumeConverter::surface(
// Into a rectangle
auto g4Box = dynamic_cast<const G4Box*>(g4Solid);
if (g4Box != nullptr) {
auto [bounds, axes, original] =
Geant4ShapeConverter{}.rectangleBounds(*g4Box);
auto orientedToGlobal = axesOriented(toGlobal, axes);
surface = Acts::Surface::makeShared<PlaneSurface>(orientedToGlobal,
std::move(bounds));
assignMaterial(*surface.get(), original, compressed);
return surface;
if (forcedType == Surface::SurfaceType::Other or
forcedType == Surface::SurfaceType::Plane) {
auto [bounds, axes, original] =
Geant4ShapeConverter{}.rectangleBounds(*g4Box);
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("Can not convert 'G4Box' into forced shape.");
}
}

// Into a Trapezoid
auto g4Trd = dynamic_cast<const G4Trd*>(g4Solid);
if (g4Trd != nullptr) {
auto [bounds, axes, original] =
Geant4ShapeConverter{}.trapezoidBounds(*g4Trd);
auto orientedToGlobal = axesOriented(toGlobal, axes);
surface = Acts::Surface::makeShared<PlaneSurface>(orientedToGlobal,
std::move(bounds));
assignMaterial(*surface.get(), original, compressed);
return surface;
if (forcedType == Surface::SurfaceType::Other or
forcedType == Surface::SurfaceType::Plane) {
auto [bounds, axes, original] =
Geant4ShapeConverter{}.trapezoidBounds(*g4Trd);
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("Can not convert 'G4Trd' into forced shape.");
}
}

// Into a Cylinder or disc
// Into a Cylinder, disc or line
auto g4Tubs = dynamic_cast<const G4Tubs*>(g4Solid);
if (g4Tubs != nullptr) {
ActsScalar diffR = g4Tubs->GetOuterRadius() - g4Tubs->GetInnerRadius();
ActsScalar diffZ = 2 * g4Tubs->GetZHalfLength();
// Detect if cylinder or disc case
ActsScalar original = 0.;
if (diffR < diffZ) {
if (forcedType == Surface::SurfaceType::Cylinder or
(diffR < diffZ and forcedType == Surface::SurfaceType::Other)) {
auto [bounds, originalT] = Geant4ShapeConverter{}.cylinderBounds(*g4Tubs);

std::cout << "Creating cylinder with " << *bounds << std::endl;

original = originalT;
surface = Acts::Surface::makeShared<CylinderSurface>(toGlobal,
std::move(bounds));
} else {
} else if (forcedType == Surface::SurfaceType::Disc or
forcedType == Surface::SurfaceType::Other) {
auto [bounds, originalT] = Geant4ShapeConverter{}.radialBounds(*g4Tubs);
original = originalT;
surface =
Acts::Surface::makeShared<DiscSurface>(toGlobal, std::move(bounds));
} else if (forcedType == Surface::SurfaceType::Straw) {
auto bounds = Geant4ShapeConverter{}.lineBounds(*g4Tubs);
surface =
Acts::Surface::makeShared<StrawSurface>(toGlobal, std::move(bounds));

} else {
throw std::runtime_error("Can not convert 'G4Tubs' into forced shape.");
}
assignMaterial(*surface.get(), original, compressed);
return surface;
Expand Down
174 changes: 172 additions & 2 deletions Tests/UnitTests/Plugins/Geant4/Geant4ConvertersTests.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,20 +10,31 @@

#include "Acts/Definitions/Algebra.hpp"
#include "Acts/Geometry/GeometryContext.hpp"
#include "Acts/Material/HomogeneousSurfaceMaterial.hpp"
#include "Acts/Material/MaterialSlab.hpp"
#include "Acts/Plugins/Geant4/Geant4Converters.hpp"
#include "Acts/Surfaces/CylinderBounds.hpp"
#include "Acts/Surfaces/LineBounds.hpp"
#include "Acts/Surfaces/RadialBounds.hpp"
#include "Acts/Surfaces/RectangleBounds.hpp"
#include "Acts/Surfaces/TrapezoidBounds.hpp"
#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"

#include "G4Box.hh"
#include "G4LogicalVolume.hh"
#include "G4Material.hh"
#include "G4PVPlacement.hh"
#include "G4RotationMatrix.hh"
#include "G4ThreeVector.hh"
#include "G4Trd.hh"
#include "G4Tubs.hh"
#include "G4VPhysicalVolume.hh"

Acts::GeometryContext tContext;

Acts::ActsScalar rho = 1.2345;
G4Material* g4Material = new G4Material("Material", 6., 12., rho);

BOOST_AUTO_TEST_SUITE(Geant4Plugin)

BOOST_AUTO_TEST_CASE(Geant4AlgebraConversion) {
Expand Down Expand Up @@ -60,8 +71,8 @@ BOOST_AUTO_TEST_CASE(Geant4CylinderConversion) {
}

BOOST_AUTO_TEST_CASE(Geant4RadialConversion) {
G4Tubs disk("disk", 40., 400., 2., 0., 2 * M_PI);
auto [bounds, thickness] = Acts::Geant4ShapeConverter{}.radialBounds(disk);
G4Tubs disc("disc", 40., 400., 2., 0., 2 * M_PI);
auto [bounds, thickness] = Acts::Geant4ShapeConverter{}.radialBounds(disc);
CHECK_CLOSE_ABS(bounds->get(Acts::RadialBounds::BoundValues::eMinR), 40.,
10e-10);
CHECK_CLOSE_ABS(bounds->get(Acts::RadialBounds::BoundValues::eMaxR), 400.,
Expand All @@ -73,6 +84,14 @@ BOOST_AUTO_TEST_CASE(Geant4RadialConversion) {
CHECK_CLOSE_ABS(thickness, 4., 10e-10);
}

BOOST_AUTO_TEST_CASE(Geant4LineConversion) {
G4Tubs line("line", 0., 20., 400., 0., 2 * M_PI);
auto bounds = Acts::Geant4ShapeConverter{}.lineBounds(line);
CHECK_CLOSE_ABS(bounds->get(Acts::LineBounds::BoundValues::eR), 20., 10e-10);
CHECK_CLOSE_ABS(bounds->get(Acts::LineBounds::BoundValues::eHalfLengthZ),
400., 10e-10);
}

BOOST_AUTO_TEST_CASE(Geant4BoxConversion) {
// Test the standard orientations
G4Box sensorXY("SensorXY", 23., 34., 1.);
Expand Down Expand Up @@ -197,4 +216,155 @@ BOOST_AUTO_TEST_CASE(Geant4PlanarConversion) {
BOOST_CHECK(tBounds != nullptr);
}

BOOST_AUTO_TEST_CASE(Geant4BoxVPhysConversion) {
Acts::ActsScalar thickness = 2.;

G4Box* g4Box = new G4Box("Box", 23., 34., 0.5 * thickness);
G4RotationMatrix* g4Rot = new G4RotationMatrix({0., 0., 1.}, 1.2);
G4LogicalVolume* g4BoxLog = new G4LogicalVolume(g4Box, g4Material, "BoxLog");

G4ThreeVector g4Trans(0., 0., 100.);
G4PVPlacement g4BoxPhys(g4Rot, g4Trans, g4BoxLog, "BoxPhys", nullptr, false,
1);

auto planeSurface = Acts::Geant4PhysicalVolumeConverter{}.surface(
g4BoxPhys, Acts::Transform3::Identity(), true, thickness);
BOOST_CHECK(planeSurface != nullptr);
BOOST_CHECK(planeSurface->type() == Acts::Surface::SurfaceType::Plane);

auto material = planeSurface->surfaceMaterial();
BOOST_CHECK(material != nullptr);

auto materialSlab = material->materialSlab(Acts::Vector3{0., 0., 0.});
// Here it should be uncompressed material
CHECK_CLOSE_ABS(materialSlab.material().massDensity(), rho, 0.001);
CHECK_CLOSE_REL(thickness / g4Material->GetRadlen(),
materialSlab.thicknessInX0(), 0.1);

// Convert with compression
Acts::ActsScalar compression = 4.;
planeSurface = Acts::Geant4PhysicalVolumeConverter{}.surface(
g4BoxPhys, Acts::Transform3::Identity(), true, thickness / compression);
BOOST_CHECK(planeSurface != nullptr);
BOOST_CHECK(planeSurface->type() == Acts::Surface::SurfaceType::Plane);

material = planeSurface->surfaceMaterial();
BOOST_CHECK(material != nullptr);
materialSlab = material->materialSlab(Acts::Vector3{0., 0., 0.});

// Here it should be uncompressed material
CHECK_CLOSE_ABS(materialSlab.material().massDensity(), compression * rho,
0.001);
CHECK_CLOSE_REL(thickness / g4Material->GetRadlen(),
materialSlab.thicknessInX0(), 0.01);

CHECK_CLOSE_ABS(materialSlab.thickness(), thickness / compression, 0.01);
CHECK_CLOSE_REL(materialSlab.material().X0() * compression,
g4Material->GetRadlen(), 0.01);

delete g4Box;
delete g4Rot;
delete g4BoxLog;
}

BOOST_AUTO_TEST_CASE(Geant4CylVPhysConversion) {
Acts::ActsScalar radius = 45.;
Acts::ActsScalar thickness = 1.;
Acts::ActsScalar halfLengthZ = 200;

G4Tubs* g4Tube =
new G4Tubs("Tube", radius, radius + thickness, halfLengthZ, 0., 2 * M_PI);

G4RotationMatrix* g4Rot = new G4RotationMatrix({0., 0., 1.}, 0.);
G4LogicalVolume* g4TubeLog =
new G4LogicalVolume(g4Tube, g4Material, "TubeLog");
G4ThreeVector g4Trans(0., 0., 100.);
G4PVPlacement g4CylinderPhys(g4Rot, g4Trans, g4TubeLog, "TubePhys", nullptr,
false, 1);

auto cylinderSurface = Acts::Geant4PhysicalVolumeConverter{}.surface(
g4CylinderPhys, Acts::Transform3::Identity(), true, thickness);
BOOST_CHECK(cylinderSurface != nullptr);
BOOST_CHECK(cylinderSurface->type() == Acts::Surface::SurfaceType::Cylinder);

auto material = cylinderSurface->surfaceMaterial();
BOOST_CHECK(material != nullptr);

auto materialSlab = material->materialSlab(Acts::Vector3{0., 0., 0.});
CHECK_CLOSE_REL(thickness / g4Material->GetRadlen(),
materialSlab.thicknessInX0(), 0.1);

// Here it should be uncompressed material
CHECK_CLOSE_ABS(materialSlab.material().massDensity(), rho, 0.001);

/// CHECK exception throwing
BOOST_CHECK_THROW(
Acts::Geant4PhysicalVolumeConverter{Acts::Surface::SurfaceType::Plane}
.surface(g4CylinderPhys, Acts::Transform3::Identity(), true,
thickness),
std::runtime_error);

delete g4Tube;
delete g4Rot;
delete g4TubeLog;
}

BOOST_AUTO_TEST_CASE(Geant4VDiscVPhysConversion) {
Acts::ActsScalar innerRadius = 45.;
Acts::ActsScalar outerRadius = 75.;
Acts::ActsScalar thickness = 2.;

G4Tubs* g4Tube = new G4Tubs("Disc", innerRadius, outerRadius, 0.5 * thickness,
0., 2 * M_PI);

G4RotationMatrix* g4Rot = new G4RotationMatrix({0., 0., 1.}, 0.);
G4LogicalVolume* g4TubeLog =
new G4LogicalVolume(g4Tube, g4Material, "TubeLog");
G4ThreeVector g4Trans(0., 0., 100.);
G4PVPlacement g4discPhys(g4Rot, g4Trans, g4TubeLog, "TubePhys", nullptr,
false, 1);

auto discSurface = Acts::Geant4PhysicalVolumeConverter{}.surface(
g4discPhys, Acts::Transform3::Identity(), true, thickness);
BOOST_CHECK(discSurface != nullptr);
BOOST_CHECK(discSurface->type() == Acts::Surface::SurfaceType::Disc);

auto material = discSurface->surfaceMaterial();
BOOST_CHECK(material != nullptr);

auto materialSlab = material->materialSlab(Acts::Vector3{0., 0., 0.});
// Here it should be uncompressed material
CHECK_CLOSE_ABS(materialSlab.material().massDensity(), rho, 0.001);

delete g4Tube;
delete g4Rot;
delete g4TubeLog;
}

BOOST_AUTO_TEST_CASE(Geant4LineVPhysConversion) {
Acts::ActsScalar innerRadius = 0.;
Acts::ActsScalar outerRadius = 20.;
Acts::ActsScalar thickness = 400.;

G4Tubs* g4Tube = new G4Tubs("Line", innerRadius, outerRadius, 0.5 * thickness,
0., 2 * M_PI);

G4RotationMatrix* g4Rot = new G4RotationMatrix({0., 0., 1.}, 0.);
G4LogicalVolume* g4TubeLog =
new G4LogicalVolume(g4Tube, g4Material, "LineLog");
G4ThreeVector g4Trans(0., 0., 100.);
G4PVPlacement g4linePhys(g4Rot, g4Trans, g4TubeLog, "LinePhys", nullptr,
false, 1);

auto lineSurface =
Acts::Geant4PhysicalVolumeConverter{Acts::Surface::SurfaceType::Straw}
.surface(g4linePhys, Acts::Transform3::Identity(), true, thickness);
BOOST_CHECK(lineSurface != nullptr);
BOOST_CHECK(lineSurface->type() == Acts::Surface::SurfaceType::Straw);

delete g4Tube;
delete g4Rot;
delete g4TubeLog;
}

BOOST_AUTO_TEST_SUITE_END()

0 comments on commit b1e5cf2

Please sign in to comment.