[geom] Keep reflections through a GDML round trip

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,14 +2512,42 @@ 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);
          }
 
          // 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())

geom/gdml/src/TGDMLWrite.cxx

@@ -699,32 +699,40 @@ void TGDMLWrite::ExtractVolumes(TGeoNode *node)
          fGdmlE->AddChild(fDefineNode, childN); // adding node to <define> 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 <scale> and <rotation> 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 <define> 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 <scale> and <rotation> 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);

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

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 <gtest/gtest.h>
+
+#include <TGeoManager.h>
+#include <TGeoMaterial.h>
+#include <TGeoMatrix.h>
+#include <TGeoMedium.h>
+#include <TGeoNode.h>
+#include <TGeoVolume.h>
+#include <TString.h>
+
+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");
+}