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"); +}