Skip to content

Commit

Permalink
fix: Geant4 transformation conversion (acts-project#3051)
Browse files Browse the repository at this point in the history
This PR addresses the issue raised in acts-project#3042 and acts-project#3038.

The underlying problem in both cases stemmed from the fact that the `Geant4DetectorSurfaceFactory` did not respect the Geant4 transformation relations between the mother and the daughter volumes. 

For example, in the picture below the `Blue` volume is the daughter of the `Red` volume and the `Red` volume is the daughter of the `Magenta` volume. The `Green` volume is the daughter of the `World` and its transformation is the composition of the `Magenta`, `Red`, and `Blue` volume transformations implemented as it is currently done in the `Geant4DetectorSurfaceFactory`. If the composition is realized correctly, the `Green` volume should completely overlap the `Blue` one, which is not the case here.  
![bad_conversion](https://github.com/acts-project/acts/assets/113530373/10f661c7-6592-4a38-8915-fad610db74b5)

The changes introduced here: 

a) Fix the transform composition in the `Geant4DetectorSurfaceFactory`;
b) Fix the transformation conversion in the `Geant4Converters` so that the Geant4 rotation/translation order (first the daughter is translated with respect to the mother-volume frame and then rotated) is respected;
c) Remove the unnecessary transposition of the rotation matrix in the `Geant4Converters` that resulted in the wrong normal direction of the converted surface.
d) Add a Unit Test that will insure the correct conversion. 

With the new composition applied this is how the picture from above changes (i.e. we see the complete overlap):    
![good_conversion](https://github.com/acts-project/acts/assets/113530373/91baddec-23cf-47fb-afde-dba193635ce4)

And this is how it looks inside Acts:
![good_conversion_acts](https://github.com/acts-project/acts/assets/113530373/c5d095ba-42d3-441f-bf81-f721b6df504c)

Using the opportunity, I'd like to raise a question of the convention for the rotation/translation order inside the Acts codebase, as it looks like right now there are no strict guidelines, which may be pretty confusing.
  • Loading branch information
ssdetlab authored and EleniXoch committed May 6, 2024
1 parent 3febc9f commit 566d07f
Show file tree
Hide file tree
Showing 4 changed files with 108 additions and 9 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -127,8 +127,10 @@ class Geant4SurfaceProvider : public Acts::Experimental::ISurfacesProvider {
Acts::Geant4DetectorSurfaceFactory::Cache g4SurfaceCache;

/// Find and store surfaces in the cache object
auto g4ToWorldConsistent = G4Transform3D(
m_g4ToWorld.getRotation().inverse(), m_g4ToWorld.getTranslation());
Acts::Geant4DetectorSurfaceFactory{}.construct(
g4SurfaceCache, m_g4ToWorld, *m_g4World, g4SurfaceOptions);
g4SurfaceCache, g4ToWorldConsistent, *m_g4World, g4SurfaceOptions);

auto surfaces = g4SurfaceCache.passiveSurfaces;

Expand Down
6 changes: 3 additions & 3 deletions Plugins/Geant4/src/Geant4Converters.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -56,11 +56,11 @@ Acts::Transform3 Acts::Geant4AlgebraConverter::transform(
scale * g4Trans[2]);
// And the rotation to it
RotationMatrix3 rotation;
rotation << g4Rot.xx(), g4Rot.yx(), g4Rot.zx(), g4Rot.xy(), g4Rot.yy(),
g4Rot.zy(), g4Rot.xz(), g4Rot.yz(), g4Rot.zz();
rotation << g4Rot.xx(), g4Rot.xy(), g4Rot.xz(), g4Rot.yx(), g4Rot.yy(),
g4Rot.yz(), g4Rot.zx(), g4Rot.zy(), g4Rot.zz();
Transform3 transform = Transform3::Identity();
transform.rotate(rotation);
transform.pretranslate(translation);
transform.prerotate(rotation);
return transform;
}

Expand Down
11 changes: 6 additions & 5 deletions Plugins/Geant4/src/Geant4DetectorSurfaceFactory.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -26,12 +26,13 @@ void Acts::Geant4DetectorSurfaceFactory::construct(
auto g4Translation = g4PhysVol.GetTranslation();
auto g4Rotation = g4PhysVol.GetRotation();

G4Transform3D g4Transform =
(g4Rotation == nullptr)
? G4Transform3D(CLHEP::HepRotation(), g4Translation)
: G4Transform3D(*g4Rotation, g4Translation);
auto newTranslation =
g4ToGlobal.getTranslation() + g4ToGlobal.getRotation() * g4Translation;
auto newRotation = (g4Rotation == nullptr)
? g4ToGlobal.getRotation() * CLHEP::HepRotation()
: g4ToGlobal.getRotation() * g4Rotation->inverse();

G4Transform3D newToGlobal = g4ToGlobal * g4Transform;
G4Transform3D newToGlobal(newRotation, newTranslation);

// Get the logical volume
auto g4LogicalVolume = g4PhysVol.GetLogicalVolume();
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,9 @@

#include "Acts/Plugins/Geant4/Geant4DetectorSurfaceFactory.hpp"
#include "Acts/Plugins/Geant4/Geant4PhysicalVolumeSelectors.hpp"
#include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
#include "Acts/Visualization/GeometryView3D.hpp"
#include "Acts/Visualization/ObjVisualization3D.hpp"

#include <memory>
#include <string>
Expand Down Expand Up @@ -94,4 +97,97 @@ BOOST_AUTO_TEST_CASE(Geant4DetecturSurfaceFactory_Cylinder) {
BOOST_CHECK_EQUAL(surface->type(), Acts::Surface::SurfaceType::Cylinder);
}

BOOST_AUTO_TEST_CASE(Geant4DetecturSurfaceFactory_Transforms) {
Acts::GeometryContext gctx;

G4Box* worldS = new G4Box("world", 1000, 1000, 1000);
G4LogicalVolume* worldLV = new G4LogicalVolume(worldS, nullptr, "World");
G4VPhysicalVolume* worldPV = new G4PVPlacement(
nullptr, G4ThreeVector(), worldLV, "World", nullptr, false, 0, false);

auto vol1S = new G4Box("volume1", 25, 10, 50);
auto vol1L = new G4LogicalVolume(vol1S, nullptr, "Volume1");

G4Transform3D transformVol1(CLHEP::HepRotationX(M_PI / 4),
G4ThreeVector(20, 0, 0));

[[maybe_unused]] auto vol1PV = new G4PVPlacement(
transformVol1, vol1L, "Volume1", worldLV, false, 0, false);

auto vol2S = new G4Box("volume2", 25, 10, 50);
auto vol2L = new G4LogicalVolume(vol2S, nullptr, "Volume2");

G4Transform3D transformVol2(CLHEP::HepRotationY(M_PI / 6),
G4ThreeVector(0, 100, 20));

[[maybe_unused]] auto vol2PV = new G4PVPlacement(
transformVol2, vol2L, "Volume2", vol1L, false, 0, false);

auto vol3S = new G4Box("volume3", 25, 10, 50);
auto vol3L = new G4LogicalVolume(vol3S, nullptr, "Volume3");

G4Transform3D transformVol3(CLHEP::HepRotationZ(M_PI / 12),
G4ThreeVector(30, 100, 0));

[[maybe_unused]] auto vol3PV = new G4PVPlacement(
transformVol3, vol3L, "Volume3", vol2L, false, 0, false);

// Get the lowest volume
auto nameSelector =
std::make_shared<Acts::Geant4PhysicalVolumeSelectors::NameSelector>(
std::vector<std::string>{"olume"}, false);

Acts::Geant4DetectorSurfaceFactory::Cache cache;
Acts::Geant4DetectorSurfaceFactory::Options options;
options.sensitiveSurfaceSelector = nameSelector;

G4Transform3D nominal;

Acts::Geant4DetectorSurfaceFactory factory;
factory.construct(cache, nominal, *worldPV, options);

auto [element, surface] = cache.sensitiveSurfaces.front();
BOOST_CHECK_EQUAL(surface->type(), Acts::Surface::SurfaceType::Plane);

auto center = surface->center(gctx);
auto normal = surface->normal(gctx, center, Acts::Vector3(1, 0, 0));
CHECK_CLOSE_ABS(center.x(), 45.981, 1e-3);
CHECK_CLOSE_ABS(center.y(), 137.886, 1e-3);
CHECK_CLOSE_ABS(center.z(), 144.957, 1e-3);
CHECK_CLOSE_ABS(normal.x(), -0.224, 1e-3);
CHECK_CLOSE_ABS(normal.y(), 0.592, 1e-3);
CHECK_CLOSE_ABS(normal.z(), 0.775, 1e-3);

Acts::ObjVisualization3D obj;
Acts::Vector3 origin(0, 0, 0);
Acts::GeometryView3D::drawArrowForward(obj, origin, Acts::Vector3(100, 0, 0),
1000, 10,
Acts::ViewConfig({255, 0, 0}));
Acts::GeometryView3D::drawArrowForward(obj, origin, Acts::Vector3(0, 100, 0),
1000, 10,
Acts::ViewConfig({0, 255, 0}));
Acts::GeometryView3D::drawArrowForward(obj, origin, Acts::Vector3(0, 0, 100),
1000, 10,
Acts::ViewConfig({0, 0, 255}));
Acts::GeometryView3D::drawArrowForward(
obj, surface->center(gctx), surface->center(gctx) + 100 * normal, 1000,
10, Acts::ViewConfig({0, 255, 0}));
auto surfaces = cache.sensitiveSurfaces;
for (const auto& [k, val] : Acts::enumerate(cache.sensitiveSurfaces)) {
const auto& [el, surf] = val;
Acts::ViewConfig vCfg;
if (k == 0) {
vCfg.color = Acts::ColorRGB({0, 255, 0});
} else if (k == 1) {
vCfg.color = Acts::ColorRGB({255, 0, 0});
} else if (k == 2) {
vCfg.color = Acts::ColorRGB({0, 255, 255});
}
Acts::GeometryView3D::drawSurface(obj, *surf, gctx,
Acts::Transform3::Identity(), vCfg);
}

obj.write("RotatedSurface.obj");
}

BOOST_AUTO_TEST_SUITE_END()

0 comments on commit 566d07f

Please sign in to comment.