From 81043cfc46b10a1d5d3848a48992300c8f0882f4 Mon Sep 17 00:00:00 2001 From: Sandro Wenzel Date: Wed, 2 Sep 2026 17:11:39 +0200 Subject: [PATCH] [geom] Keep reflections through a GDML round trip This fixes https://github.com/root-project/root/issues/23237, a bug in the GDML writer and reader that could lose reflections when a geometry was exported to GDML and read back. The writer now preserves reflections combined with rotations by factoring the reflection into a z-mirror and writing the remaining proper rotation. Diagonal reflection matrices keep their existing representation. The reader now also handles the scale of placements inside assemblies, so reflections are preserved there as well. A new test exports and re-imports a reflected box with an off-centre daughter, both in the world and in an assembly. It checks the global position, frame handedness and TGeoHMatrix::IsReflection() before and after the round trip. The ALICE geometry is now correctly exported to GDML, with reflections preserved through the round trip. --- geom/gdml/src/TGDMLParse.cxx | 33 +++++++- geom/gdml/src/TGDMLWrite.cxx | 84 +++++++++++--------- geom/test/CMakeLists.txt | 3 + geom/test/test_gdml_reflection.cxx | 118 +++++++++++++++++++++++++++++ 4 files changed, 202 insertions(+), 36 deletions(-) create mode 100644 geom/test/test_gdml_reflection.cxx diff --git a/geom/gdml/src/TGDMLParse.cxx b/geom/gdml/src/TGDMLParse.cxx index 5c3f9e5792671..b17af2441fa2b 100644 --- a/geom/gdml/src/TGDMLParse.cxx +++ b/geom/gdml/src/TGDMLParse.cxx @@ -2466,7 +2466,7 @@ XMLNodePointer_t TGDMLParse::AssProcess(TXMLEngine *gdml, XMLNodePointer_t node) TGeoVolume *lv = nullptr; TGeoTranslation *pos = nullptr; TGeoRotation *rot = nullptr; - TGeoCombiTrans *matr; + TGeoScale *scl = nullptr; TGeoVolumeAssembly *assem = new TGeoVolumeAssembly(NameShort(name)); @@ -2483,6 +2483,7 @@ XMLNodePointer_t TGDMLParse::AssProcess(TXMLEngine *gdml, XMLNodePointer_t node) subchild = gdml->GetChild(child); pos = new TGeoTranslation(0, 0, 0); rot = new TGeoRotation(); + scl = nullptr; while (subchild != nullptr) { tempattr = gdml->GetNodeName(subchild); @@ -2511,6 +2512,16 @@ XMLNodePointer_t TGDMLParse::AssProcess(TXMLEngine *gdml, XMLNodePointer_t node) RotProcess(gdml, subchild, attr); reftemp = gdml->GetAttr(subchild, "name"); rot = GetRotation(reftemp.Data()); + } else if (tempattr == "scale") { + attr = gdml->GetFirstAttr(subchild); + SclProcess(gdml, subchild, attr); + reftemp = gdml->GetAttr(subchild, "name"); + scl = GetScaleObj(reftemp.Data()); + } else if (tempattr == "scaleref") { + reftemp = gdml->GetAttr(subchild, "ref"); + scl = GetScaleObj(reftemp.Data()); + if (!scl) + Fatal("AssProcess", "Physvol's scale %s not found", reftemp.Data()); } subchild = gdml->GetNext(subchild); @@ -2518,7 +2529,25 @@ XMLNodePointer_t TGDMLParse::AssProcess(TXMLEngine *gdml, XMLNodePointer_t node) // ADD PHYSVOL TO GEOMETRY fVolID = fVolID + 1; - matr = new TGeoCombiTrans(*pos, *rot); + TGeoHMatrix *matr = new TGeoHMatrix(); + matr->SetTranslation(pos->GetTranslation()); + matr->SetRotation(rot->GetRotationMatrix()); + + if (scl != nullptr) { // Scaling must be added to the rotation matrix! + + Double_t scale3x3[9]; + memset(scale3x3, 0, 9 * sizeof(Double_t)); + const Double_t *diagonal = scl->GetScale(); + + scale3x3[0] = diagonal[0]; + scale3x3[4] = diagonal[1]; + scale3x3[8] = diagonal[2]; + + TGeoRotation scaleMatrix; + scaleMatrix.SetMatrix(scale3x3); + matr->Multiply(&scaleMatrix); + } + assem->AddNode(lv, copynum, matr); TGeoNode *lastnode = (TGeoNode *)assem->GetNodes()->Last(); if (!pnodename.IsNull()) diff --git a/geom/gdml/src/TGDMLWrite.cxx b/geom/gdml/src/TGDMLWrite.cxx index d2599020f612e..af753dbd3443d 100644 --- a/geom/gdml/src/TGDMLWrite.cxx +++ b/geom/gdml/src/TGDMLWrite.cxx @@ -699,32 +699,40 @@ void TGDMLWrite::ExtractVolumes(TGeoNode *node) fGdmlE->AddChild(fDefineNode, childN); // adding node to node // Deal with reflection XMLNodePointer_t scaleN = nullptr; - Double_t lx, ly, lz; - Double_t xangle = 0; - Double_t zangle = 0; - lx = geoNode->GetMatrix()->GetRotationMatrix()[0]; - ly = geoNode->GetMatrix()->GetRotationMatrix()[4]; - lz = geoNode->GetMatrix()->GetRotationMatrix()[8]; - if (geoNode->GetMatrix()->IsReflection() && TMath::Abs(lx) == 1 && TMath::Abs(ly) == 1 && - TMath::Abs(lz) == 1) { + Double_t rot3x3[9]; + for (Int_t i = 0; i < 9; ++i) { + rot3x3[i] = geoNode->GetMatrix()->GetRotationMatrix()[i]; + } + Double_t lx = rot3x3[0]; + Double_t ly = rot3x3[4]; + Double_t lz = rot3x3[8]; + if (geoNode->GetMatrix()->IsReflection()) { + // A physvol is read back as translation * rotation * scale, so the reflection + // is carried by and gets the matrix with the scale taken + // out again, which is a proper rotation. A matrix that is already diagonal + // keeps its own signs and leaves an identity rotation; every other reflection + // is written as a mirror on z plus a rotation. + if (TMath::Abs(lx) != 1 || TMath::Abs(ly) != 1 || TMath::Abs(lz) != 1) { + lx = 1.; + ly = 1.; + lz = -1.; + } scaleN = fGdmlE->NewChild(nullptr, nullptr, "scale", nullptr); fGdmlE->NewAttr(scaleN, nullptr, "name", (nodename + "scl").Data()); fGdmlE->NewAttr(scaleN, nullptr, "x", TString::Format(fltPrecision.Data(), lx)); fGdmlE->NewAttr(scaleN, nullptr, "y", TString::Format(fltPrecision.Data(), ly)); fGdmlE->NewAttr(scaleN, nullptr, "z", TString::Format(fltPrecision.Data(), lz)); - // experimentally found out, that rotation should be updated like this - if (lx == -1) { - zangle = 180; - } - if (lz == -1) { - xangle = 180; + // take the scale out column by column; the entries are +-1, so multiplying + // by them again is the same as dividing + for (Int_t i = 0; i < 3; ++i) { + rot3x3[3 * i] *= lx; + rot3x3[3 * i + 1] *= ly; + rot3x3[3 * i + 2] *= lz; } } // rotation - TGDMLWrite::Xyz lxyz = GetXYZangles(geoNode->GetMatrix()->GetRotationMatrix()); - lxyz.x -= xangle; - lxyz.z -= zangle; + TGDMLWrite::Xyz lxyz = GetXYZangles(rot3x3); if ((lxyz.x != 0.0) || (lxyz.y != 0.0) || (lxyz.z != 0.0)) { rotname = nodename + "rot"; childN = CreateRotationN(rotname.Data(), lxyz); @@ -2797,32 +2805,40 @@ void TGDMLWrite::ExtractVolumes(TGeoVolume *volume) fGdmlE->AddChild(fDefineNode, childN); // adding node to node // Deal with reflection XMLNodePointer_t scaleN = nullptr; - Double_t lx, ly, lz; - Double_t xangle = 0; - Double_t zangle = 0; - lx = geoNode->GetMatrix()->GetRotationMatrix()[0]; - ly = geoNode->GetMatrix()->GetRotationMatrix()[4]; - lz = geoNode->GetMatrix()->GetRotationMatrix()[8]; - if (geoNode->GetMatrix()->IsReflection() && TMath::Abs(lx) == 1 && TMath::Abs(ly) == 1 && - TMath::Abs(lz) == 1) { + Double_t rot3x3[9]; + for (Int_t i = 0; i < 9; ++i) { + rot3x3[i] = geoNode->GetMatrix()->GetRotationMatrix()[i]; + } + Double_t lx = rot3x3[0]; + Double_t ly = rot3x3[4]; + Double_t lz = rot3x3[8]; + if (geoNode->GetMatrix()->IsReflection()) { + // A physvol is read back as translation * rotation * scale, so the reflection + // is carried by and gets the matrix with the scale taken + // out again, which is a proper rotation. A matrix that is already diagonal + // keeps its own signs and leaves an identity rotation; every other reflection + // is written as a mirror on z plus a rotation. + if (TMath::Abs(lx) != 1 || TMath::Abs(ly) != 1 || TMath::Abs(lz) != 1) { + lx = 1.; + ly = 1.; + lz = -1.; + } scaleN = fGdmlE->NewChild(nullptr, nullptr, "scale", nullptr); fGdmlE->NewAttr(scaleN, nullptr, "name", (nodename + "scl").Data()); fGdmlE->NewAttr(scaleN, nullptr, "x", TString::Format(fltPrecision.Data(), lx)); fGdmlE->NewAttr(scaleN, nullptr, "y", TString::Format(fltPrecision.Data(), ly)); fGdmlE->NewAttr(scaleN, nullptr, "z", TString::Format(fltPrecision.Data(), lz)); - // experimentally found out, that rotation should be updated like this - if (lx == -1) { - zangle = 180; - } - if (lz == -1) { - xangle = 180; + // take the scale out column by column; the entries are +-1, so multiplying + // by them again is the same as dividing + for (Int_t i = 0; i < 3; ++i) { + rot3x3[3 * i] *= lx; + rot3x3[3 * i + 1] *= ly; + rot3x3[3 * i + 2] *= lz; } } // rotation - TGDMLWrite::Xyz lxyz = GetXYZangles(geoNode->GetMatrix()->GetRotationMatrix()); - lxyz.x -= xangle; - lxyz.z -= zangle; + TGDMLWrite::Xyz lxyz = GetXYZangles(rot3x3); if ((lxyz.x != 0.0) || (lxyz.y != 0.0) || (lxyz.z != 0.0)) { rotname = nodename + "rot"; childN = CreateRotationN(rotname.Data(), lxyz); diff --git a/geom/test/CMakeLists.txt b/geom/test/CMakeLists.txt index c52ed9610d743..55fb804ec755a 100644 --- a/geom/test/CMakeLists.txt +++ b/geom/test/CMakeLists.txt @@ -16,6 +16,9 @@ if(gdml) ROOT_ADD_GTEST(boolean_extrusion test_boolean_extrusion.cxx LIBRARIES Geom GeomChecker) + ROOT_ADD_GTEST(gdml_reflection + test_gdml_reflection.cxx + LIBRARIES Geom) endif() ROOT_ADD_GTEST(tessellated test_tessellated.cxx diff --git a/geom/test/test_gdml_reflection.cxx b/geom/test/test_gdml_reflection.cxx new file mode 100644 index 0000000000000..d0ffe1a2bd23b --- /dev/null +++ b/geom/test/test_gdml_reflection.cxx @@ -0,0 +1,118 @@ +// A reflected placement must survive a GDML round trip - also when the placement +// carries a rotation as well, and also when it sits inside an assembly. + +#include + +#include +#include +#include +#include +#include +#include +#include + +namespace { + +// The global placement of the volume named "inner" +struct Placement { + double translation[3] = {0., 0., 0.}; + double xaxis[3] = {0., 0., 0.}; + bool reflection = false; + bool found = false; +}; + +void FindInner(TGeoNode *node, const TGeoHMatrix &upto, Placement &p) +{ + TGeoHMatrix here = upto; + here.Multiply(node->GetMatrix()); + if (TString(node->GetVolume()->GetName()) == "inner") { + const double *tr = here.GetTranslation(); + for (int k = 0; k < 3; ++k) + p.translation[k] = tr[k]; + const double localx[3] = {1., 0., 0.}; + here.LocalToMasterVect(localx, p.xaxis); + p.reflection = here.IsReflection(); + p.found = true; + return; + } + for (int i = 0; i < node->GetNdaughters() && !p.found; ++i) + FindInner(node->GetDaughter(i), here, p); +} + +Placement InnerPlacement() +{ + Placement p; + TGeoHMatrix identity; + FindInner(gGeoManager->GetTopNode(), identity, p); + return p; +} + +// A world holding one reflected box with an off-centre daughter in it. The +// rotation passed to ReflectZ is what distinguishes the two cases, and +// inAssembly puts the reflected placement inside an assembly instead of +// directly in the world. +void BuildGeometry(double phi, double theta, double psi, bool inAssembly) +{ + delete gGeoManager; + new TGeoManager("refl", "reflected placement through GDML"); + TGeoMedium *medium = new TGeoMedium("m", 1, new TGeoMaterial("m", 26.98, 13, 2.7)); + TGeoVolume *top = gGeoManager->MakeBox("world", medium, 200., 200., 200.); + gGeoManager->SetTopVolume(top); + + TGeoVolume *outer = gGeoManager->MakeBox("outer", medium, 80., 80., 80.); + outer->AddNode(gGeoManager->MakeBox("inner", medium, 5., 5., 5.), 1, new TGeoTranslation(10., 20., 30.)); + + TGeoRotation *rotation = new TGeoRotation("r", phi, theta, psi); + rotation->ReflectZ(kTRUE); // determinant -1 + TGeoCombiTrans *placement = new TGeoCombiTrans(0., 0., 0., rotation); + if (inAssembly) { + TGeoVolume *assembly = new TGeoVolumeAssembly("asm"); + assembly->AddNode(outer, 1, placement); + top->AddNode(assembly, 1, new TGeoTranslation(0., 0., 0.)); + } else { + top->AddNode(outer, 1, placement); + } + gGeoManager->CloseGeometry(); +} + +void CheckRoundTrip(double phi, double theta, double psi, bool inAssembly, const char *filename) +{ + BuildGeometry(phi, theta, psi, inAssembly); + const Placement before = InnerPlacement(); + ASSERT_TRUE(before.found); + EXPECT_TRUE(before.reflection); + + gGeoManager->Export(filename); + delete gGeoManager; + ASSERT_NE(TGeoManager::Import(filename), nullptr); + + const Placement after = InnerPlacement(); + ASSERT_TRUE(after.found); + EXPECT_EQ(before.reflection, after.reflection); + for (int k = 0; k < 3; ++k) { + EXPECT_NEAR(before.translation[k], after.translation[k], 1e-9); + EXPECT_NEAR(before.xaxis[k], after.xaxis[k], 1e-9); + } +} + +} // namespace + +TEST(GDMLReflection, PureReflection) +{ + CheckRoundTrip(0., 0., 0., false, "gdml_reflection_pure.gdml"); +} + +TEST(GDMLReflection, ReflectionWithRotation) +{ + CheckRoundTrip(30., 40., 50., false, "gdml_reflection_rotated.gdml"); +} + +TEST(GDMLReflection, PureReflectionInAssembly) +{ + CheckRoundTrip(0., 0., 0., true, "gdml_reflection_pure_assembly.gdml"); +} + +TEST(GDMLReflection, ReflectionWithRotationInAssembly) +{ + CheckRoundTrip(30., 40., 50., true, "gdml_reflection_rotated_assembly.gdml"); +}