Introduce shape & volume sorters

GeoModelTools/GeoModelFuncSnippets/GeoModelFuncSnippets/GeoPhysVolSorter.h

+
+/*
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
+*/
+#ifndef GEOMODELFUNCSNIPPETS_GEOPHYSVOLSORTER_H
+#define GEOMODELFUNCSNIPPETS_GEOPHYSVOLSORTER_H
+
+#include "GeoModelKernel/GeoVPhysVol.h"
+
+/**
+ *  Helper struct to compare two physical volumes for the insertion into a set
+ *  Two volumes are considered to be equal if they have:
+ * 
+ *      1) Equivalent shapes
+ *      2) Same material pointer
+ *      3) Same number of children
+ *      4) Each child is placed at the same position inside the volume & 1-3 holds *          
+*/
+struct GeoPhysVolSorter {
+
+    bool operator()(const PVConstLink& a, const PVConstLink& b) const{ 
+       return compare(a, b) < 0;
+    }
+    
+    bool operator()(const GeoVPhysVol* a, const GeoVPhysVol* b) const {
+        return compare(a, b) < 0;
+    }
+
+    int compare(const GeoVPhysVol*a, const GeoVPhysVol* b) const;
+};
+#endif
\ No newline at end of file

GeoModelTools/GeoModelFuncSnippets/GeoModelFuncSnippets/GeoShapeSorter.h

+/*
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
+*/
+#ifndef GEOMODELFUNCSNIPPETS_GeoShapeSorter_H
+#define GEOMODELFUNCSNIPPETS_GeoShapeSorter_H
+
+#include "GeoModelKernel/GeoShape.h"
+#include "GeoModelKernel/GeoIntrusivePtr.h"
+#include <set>
+
+/// @brief  Helper struct to construct set of GeoShapes where the defining parameters of
+///         each element are unique within 1 micrometer
+struct GeoShapeSorter {
+   template<class ShapeType1, class ShapeType2>
+    bool operator()(const GeoIntrusivePtr<ShapeType1>& a, 
+                    const GeoIntrusivePtr<ShapeType2>& b) const {
+            return (*this)(a.get(), b.get());
+    }
+
+    bool operator()(const GeoShape* a, const GeoShape* b) const;
+    /// @brief Compares 2 GeoShapes
+    /// @param a : Pointer to shape A
+    /// @param b : Pointer to shape B
+    /// @return  Returns 0 if no defining shape property could be found that differs
+    ///          Returns -1 if the first defining & differing property of A is smaller
+    ///          Returns 1 otherwise
+    int compare(const GeoShape* a, const GeoShape* b) const;
+
+};
+
+struct GeoComposedShapeSorter {
+    template<class ShapeType1, class ShapeType2>
+    bool operator()(const GeoIntrusivePtr<ShapeType1>& a, 
+                    const GeoIntrusivePtr<ShapeType2>& b) const {
+            return (*this)(a.get(), b.get());
+    }
+    bool operator()(const GeoShape* a, const GeoShape* b) const;
+
+};
+/// @brief 
+/// @tparam ShapeType 
+template<class ShapeType>
+using GeoShapeSet = std::set<GeoIntrusivePtr<const ShapeType>, GeoShapeSorter>;
+
+#endif
\ No newline at end of file

GeoModelTools/GeoModelFuncSnippets/GeoModelFuncSnippets/GeoShapeUtils.h

+/*
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
+*/
+#ifndef GEOMDOELFUNCSNPPETS_PRINTSHAPE_H
+#define GEOMDOELFUNCSNPPETS_PRINTSHAPE_H
+#include "GeoModelKernel/GeoShape.h"
+#include "GeoModelKernel/GeoIntrusivePtr.h"
+#include <string>
+
+std::string printGeoShape(const GeoShape* shape);
+
+
+/// Returns the sub operators of a GeoShape (GeoShapeUnion, GeoShapeSubtraction, GeoShapeUnion, GeoShapeShift)
+std::pair<const GeoShape* , const GeoShape*> getOps(const GeoShape* composed);
+/// Counts out of how many shapes the shape object is made up
+unsigned int countComposedShapes(const GeoShape* shape);
+
+
+/// @brief  If a the operating shape of a ShapeShift is another ShapeShift a
+///         new shapeShift is created with the combined transformation 
+GeoIntrusivePtr<const GeoShape> compressShift(const GeoShape* shift);
+
+/// @brief Returns all acting of coomponents of a BooleanVolume. 
+///        I.e. in case Of GeoShapeInterSection & GeoShapeUnion all building shape volumes
+///        In the case of the GeoSubtraction all subtracting shapes
+std::vector<const GeoShape*> getBooleanComponents(const GeoShape* booleanShape);
+
+
+
+#endif
\ No newline at end of file

GeoModelTools/GeoModelFuncSnippets/GeoModelFuncSnippets/TransformSorter.h

+/*
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
+*/
+#ifndef GEOMODELFUNCSNIPPETS_TRANSFORMSORTER_H
+#define GEOMODELFUNCSNIPPETS_TRANSFORMSORTER_H
+
+#include "GeoModelKernel/GeoDefinitions.h"
+#include "GeoModelKernel/Units.h"
+
+#include <memory>
+
+
+namespace GeoTrf {
+      /// @brief: Helper struct to construct sets of Tranform3D or sets of
+      ///         smart pointers of Transform3D. Two Transforms are considered to
+      ///         be equivalent if their translation difference is less 
+      ///         than 0.1 microns in each vector component and if their angles
+      ///         coincide within 0.01 degrees.  
+      struct TransformSorter {
+      
+          bool operator()(const std::unique_ptr<Transform3D>& a, 
+                          const std::unique_ptr<Transform3D>& b) const {
+              return (*this) (*a, *b);
+          }
+          bool operator()(const std::shared_ptr<Transform3D>& a,
+                          const std::shared_ptr<Transform3D>& b) const {
+              return (*this)(*a, *b);          
+          }
+
+          bool operator()(const Transform3D& a, const Transform3D& b) const;
+          bool operator()(const RotationMatrix3D&a, const RotationMatrix3D& b) const;
+
+          int compare(const Transform3D&a, const Transform3D& b) const;
+          int compare(const RotationMatrix3D&a, const RotationMatrix3D& b) const;
+
+          /// @brief  Compares two N-dimension vectors component wise
+          /// @param vecA Vector a to compare
+          /// @param vecB Vector b to compare
+          /// @return 0 if all components are the same within the tolerance
+          ///        -1 if the first differing component of vecA is smaller than the one of vecB
+          ///         1 otherwise
+          template <int N> int compare(const VectorN<N>& vecA, const VectorN<N>& vecB) const {
+              constexpr double transTolerance = 0.1 * GeoModelKernelUnits::micrometer;
+              for (int d =0 ; d < N ; ++d) {
+                  const double diffValue = vecA[d] - vecB[d];
+                  if (std::abs(diffValue) > transTolerance) {
+                     return diffValue < 0. ? -1 : 1;
+                  }
+              }
+              return 0;
+          }
+          template <int N> bool operator()(const VectorN<N>& vecA, const VectorN<N>& vecB) const {
+              return compare(vecA, vecB) < 0;
+          }
+  };
+}
+
+#endif
\ No newline at end of file

GeoModelTools/GeoModelFuncSnippets/GeoModelFuncSnippets/TransformToStringConverter.h

+/*
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
+*/
+#ifndef GEOMODELFUNCSNIPPETS_TRANSFORMTOSTRINGCONVERTER_H
+#define GEOMODELFUNCSNIPPETS_TRANSFORMTOSTRINGCONVERTER_H
+#include "GeoModelKernel/GeoDefinitions.h"
+
+#include <string>
+#include <sstream>
+#include <iomanip>
+namespace GeoTrf{
+    
+    template <int N>
+    inline std::string toString(const VectorN<N>& vector, int precision = 4) {
+        std::stringstream sout;
+        sout << std::setiosflags(std::ios::fixed) << std::setprecision(precision);
+        sout << "(";
+        for (int i = 0; i < N; ++i) {
+            sout << vector[i];
+            if (i != N - 1) {
+                sout << ", ";
+            }
+        }
+        sout << ")";
+        return sout.str();
+    }
+
+    std::string toString(const Transform3D& trans, bool useCoordAngles = false, int precision = 4);
+    std::string toString(const EulerAngles& angles, int precision = 4);
+    std::string toString(const CoordEulerAngles& angles, int precision = 4);
+}
+
+#endif
\ No newline at end of file

GeoModelTools/GeoModelFuncSnippets/GeoModelFuncSnippets/getChildNodesWithTrf.h

+/*
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
+*/
+#ifndef GEOMDOELFUNCSNPPETS_GETCHILDNODESWITHTF_H
+#define GEOMDOELFUNCSNPPETS_GETCHILDNODESWITHTF_H
+
+#include "GeoModelKernel/GeoVPhysVol.h"
+#include "GeoModelKernel/GeoDefinitions.h"
+
+
+class GeoVolumeCursor;
+
+struct GeoChildNodeWithTrf {
+        
+        /// @brief  Transformation to the child node
+        GeoTrf::Transform3D transform{GeoTrf::Transform3D::Identity()};
+        /// @brief In cases of multiple copies. What's the transformation 
+        ///        to go from n --> n+1
+        GeoTrf::Transform3D inductionRule{GeoTrf::Transform3D::Identity()};
+        /// @brief Physical volume pointer to the child
+        PVConstLink volume{};
+        /// @brief Name of the physical volume empty if cursor finds ANON
+        std::string nodeName{};
+        /// @brief How many times is the volume placed in the child tree
+        unsigned int nCopies{1};
+        /// @brief tag whether the transformation is alignable
+        bool isAlignable{false};
+        /// @brief tag whether the physical volume is a full physVol
+        bool isSensitive{false};
+
+        GeoChildNodeWithTrf() = default;
+        GeoChildNodeWithTrf(GeoVolumeCursor& curs);
+};
+
+std::vector <GeoChildNodeWithTrf> getChildrenWithRef (PVConstLink& physVol,
+                                                      bool summarizeEqualVol = true);
+
+
+#endif
\ No newline at end of file

GeoModelTools/GeoModelFuncSnippets/src/GeoPhysVolSorter.cxx

+/*
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
+*/
+#include "GeoModelFuncSnippets/GeoPhysVolSorter.h"
+#include "GeoModelFuncSnippets/TransformSorter.h"
+#include "GeoModelFuncSnippets/GeoShapeSorter.h"
+#include "GeoModelKernel/GeoVolumeCursor.h"
+#include "GeoModelFuncSnippets/getChildNodesWithTrf.h"
+
+int GeoPhysVolSorter::compare(const GeoVPhysVol* a, const GeoVPhysVol* b) const {
+
+    const GeoLogVol* logVolA{a->getLogVol()}, *logVolB{b->getLogVol()};
+    if (logVolA->getMaterial() != logVolA->getMaterial()) {
+        return logVolA->getMaterial()->getName() <
+               logVolB->getMaterial()->getName() ? -1 : 1;
+    }
+    /// Check whether the shape differ
+    {
+        static const GeoShapeSorter sorter{};
+        const int shapeCmp = sorter.compare(logVolA->getShape(), logVolB->getShape());
+        if (shapeCmp) return shapeCmp;
+    }
+
+    static const GeoTrf::TransformSorter sorter{};
+
+    GeoVolumeCursor cursA{a}, cursB{b};
+    while (!cursA.atEnd() && !cursB.atEnd()) {
+        
+        GeoChildNodeWithTrf childA{cursA};
+        GeoChildNodeWithTrf childB{cursB};
+        cursA.next();
+        cursB.next();
+        /// Check whether there's an alignable transform somewhere
+        if (childA.isAlignable != childB.isAlignable) {
+            return childA.isAlignable;
+        }
+        /// Check whether the voumes are full physical volumes
+        if (childA.isSensitive != childB.isSensitive) {
+            return childA.isSensitive;
+        }
+        /// Check equivalance of the transformations
+        const int transCmp = sorter.compare(childA.transform,
+                                            childB.transform);
+        if (transCmp) return transCmp;
+        /// Every day greets the marmot
+        const int childCmp = compare(childA.volume, childB.volume);
+        if (childCmp) return childCmp;
+    }
+    if (cursA.atEnd() != cursB.atEnd()) {
+        if (cursA.atEnd()) return -1;
+        return 1;
+    }
+
+    return 0;
+}
\ No newline at end of file

GeoModelTools/GeoModelFuncSnippets/src/GeoShapeSorter.cxx

+/*
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
+*/
+#include "GeoModelFuncSnippets/GeoShapeSorter.h"
+#include "GeoModelFuncSnippets/TransformSorter.h"
+#include "GeoModelFuncSnippets/TransformToStringConverter.h"
+#include "GeoModelFuncSnippets/GeoShapeUtils.h"
+#include "GeoModelFuncSnippets/throwExcept.h"
+
+#include "GeoModelKernel/Units.h"
+
+#include "GeoModelKernel/GeoShapeUnion.h"
+#include "GeoModelKernel/GeoShapeSubtraction.h"
+#include "GeoModelKernel/GeoShapeIntersection.h"
+#include "GeoModelKernel/GeoShapeShift.h"
+
+
+#include "GeoModelKernel/GeoBox.h"
+#include "GeoModelKernel/GeoTrd.h"
+#include "GeoModelKernel/GeoTube.h"
+#include "GeoModelKernel/GeoTubs.h"
+#include "GeoModelKernel/GeoCons.h"
+#include "GeoModelKernel/GeoEllipticalTube.h"
+#include "GeoModelKernel/GeoPara.h"
+#include "GeoModelKernel/GeoPcon.h"
+#include "GeoModelKernel/GeoGenericTrap.h"
+#include "GeoModelKernel/GeoPgon.h"
+#include "GeoModelKernel/GeoSimplePolygonBrep.h"
+#include "GeoModelKernel/GeoTorus.h"
+#include "GeoModelKernel/GeoTrap.h"
+#include "GeoModelKernel/GeoTwistedTrap.h"
+
+#include <exception>
+#include <sstream>
+#include <iostream>
+
+namespace {
+   constexpr double tolerance = 1.* GeoModelKernelUnits::micrometer;
+   static const GeoTrf::TransformSorter transCmp{};
+}
+
+#define CHECK_PROPERTY(shapeA, shapeB, PROP_NAME)         \
+   {                                                      \
+      const double diffValue = 1.*shapeA->PROP_NAME() -   \
+                               1.*shapeB->PROP_NAME();    \
+      if (std::abs(diffValue) > tolerance) {              \
+        return diffValue < 0 ? - 1 : 1;                   \
+      }                                                   \
+   }                                                      \
+
+#define CHECK_VEC_PROPERTY(shapeA, shapeB, PROP_NAME, nChecks)  \
+    {                                                           \
+      for(unsigned int n = 0 ; n < nChecks; ++n) {              \
+        const double diffValue = 1.*shapeA->PROP_NAME(n) -      \
+                                 1.*shapeB->PROP_NAME(n);       \
+                                                                \
+        if (std::abs(diffValue) > tolerance) {                  \
+          return diffValue < 0 ? - 1 : 1;                       \
+        }                                                       \
+      }                                                         \
+    }
+
+#define CALL_SORTER(shapeA, shapeB)             \
+    {                                           \
+      const int cmpAB{compare(shapeA, shapeB)}; \
+      if (cmpAB) return cmpAB;                  \
+    }
+
+bool GeoShapeSorter::operator()(const GeoShape* objA, const GeoShape* objB) const {  
+  if (!objA || !objB) {
+     THROW_EXCEPTION("Nullptr given to the comparator");
+  }
+  return compare(objA, objB) < 0;
+}
+
+int GeoShapeSorter::compare(const GeoShape* objA, const GeoShape* objB) const {
+  CHECK_PROPERTY(objA, objB, typeID);
+  /// Check the defining parameters of each shape one by one
+  const int typeID = objA->typeID(); 
+  if (typeID == GeoShapeUnion::getClassTypeID() ||
+      typeID == GeoShapeIntersection::getClassTypeID() ||
+      typeID == GeoShapeSubtraction::getClassTypeID()) {
+    
+    std::pair<const GeoShape*, const GeoShape*> shapeOpsA{getOps(objA)};
+    std::pair<const GeoShape*, const GeoShape*> shapeOpsB{getOps(objB)};
+    CALL_SORTER(shapeOpsA.first, shapeOpsB.first);
+    CALL_SORTER(shapeOpsA.second, shapeOpsB.second);
+  } 
+  ///
+  else if (typeID == GeoShapeShift::getClassTypeID()) {
+    const GeoShapeShift* shapeA = dynamic_cast<const GeoShapeShift*>(objA);
+    const GeoShapeShift* shapeB = dynamic_cast<const GeoShapeShift*>(objB);
+    const int  xCmp = transCmp.compare(shapeA->getX(), shapeB->getX());
+    if (xCmp) {
+        return xCmp;
+    }
+    CALL_SORTER(shapeA->getOp(), shapeB->getOp());
+  } 
+  else if (typeID == GeoBox::getClassTypeID()) {
+    const GeoBox* boxA = dynamic_cast<const GeoBox*>(objA);
+    const GeoBox* boxB = dynamic_cast<const GeoBox*>(objB);
+    CHECK_PROPERTY(boxA, boxB, getXHalfLength);
+    CHECK_PROPERTY(boxA, boxB, getYHalfLength);
+    CHECK_PROPERTY(boxA, boxB, getZHalfLength);
+  } else if (typeID == GeoTrd::getClassTypeID()) {
+    const GeoTrd* trdA = dynamic_cast<const GeoTrd*>(objA);
+    const GeoTrd* trdB = dynamic_cast<const GeoTrd*>(objB);
+    CHECK_PROPERTY(trdA, trdB, getXHalfLength1);
+    CHECK_PROPERTY(trdA, trdB, getXHalfLength2);
+    CHECK_PROPERTY(trdA, trdB, getYHalfLength1);
+    CHECK_PROPERTY(trdA, trdB, getYHalfLength2);
+    CHECK_PROPERTY(trdA, trdB, getZHalfLength);
+  } else if (typeID == GeoTube::getClassTypeID()) {
+    const GeoTube* tubeA = dynamic_cast<const GeoTube*>(objA);
+    const GeoTube* tubeB = dynamic_cast<const GeoTube*>(objB);
+    CHECK_PROPERTY(tubeA, tubeB, getRMin);
+    CHECK_PROPERTY(tubeA, tubeB, getRMax);
+    CHECK_PROPERTY(tubeA, tubeB, getZHalfLength);
+  } else if (typeID == GeoTubs::getClassTypeID()) {
+    const GeoTubs* tubeA = dynamic_cast<const GeoTubs*>(objA);
+    const GeoTubs* tubeB = dynamic_cast<const GeoTubs*>(objB);
+    CHECK_PROPERTY(tubeA, tubeB, getRMin);
+    CHECK_PROPERTY(tubeA, tubeB, getRMax);
+    CHECK_PROPERTY(tubeA, tubeB, getZHalfLength);
+    CHECK_PROPERTY(tubeA, tubeB, getSPhi);
+    CHECK_PROPERTY(tubeA, tubeB, getDPhi);
+  } else if (typeID == GeoCons::getClassTypeID()) {
+    const GeoCons* consA = dynamic_cast<const GeoCons*>(objA);
+    const GeoCons* consB = dynamic_cast<const GeoCons*>(objB);
+    CHECK_PROPERTY(consA, consB, getRMin1);
+    CHECK_PROPERTY(consA, consB, getRMin2);
+    CHECK_PROPERTY(consA, consB, getRMin1);
+    CHECK_PROPERTY(consA, consB, getRMax1);
+    CHECK_PROPERTY(consA, consB, getRMax2);
+    CHECK_PROPERTY(consA, consB, getDZ);
+    CHECK_PROPERTY(consA, consB, getSPhi);
+    CHECK_PROPERTY(consA, consB, getDPhi);
+  } else if (typeID == GeoEllipticalTube::getClassTypeID()) {
+    const GeoEllipticalTube* tubeA = dynamic_cast<const GeoEllipticalTube*>(objA);
+    const GeoEllipticalTube* tubeB = dynamic_cast<const GeoEllipticalTube*>(objB);
+    CHECK_PROPERTY(tubeA, tubeB, getXHalfLength);
+    CHECK_PROPERTY(tubeA, tubeB, getYHalfLength);
+    CHECK_PROPERTY(tubeA, tubeB, getZHalfLength);
+   } else if (typeID == GeoPara::getClassTypeID()) {
+     const GeoPara* paraA = dynamic_cast<const GeoPara*>(objA);
+     const GeoPara* paraB = dynamic_cast<const GeoPara*>(objB);
+     CHECK_PROPERTY(paraA, paraB, getXHalfLength);
+     CHECK_PROPERTY(paraA, paraB, getYHalfLength);
+     CHECK_PROPERTY(paraA, paraB, getZHalfLength);   
+     CHECK_PROPERTY(paraA, paraB, getTheta);
+     CHECK_PROPERTY(paraA, paraB, getAlpha);
+     CHECK_PROPERTY(paraA, paraB, getPhi); 
+   } else if (typeID == GeoGenericTrap::getClassTypeID()) {
+      const GeoGenericTrap* trapA = dynamic_cast<const GeoGenericTrap*>(objA);
+      const GeoGenericTrap* trapB = dynamic_cast<const GeoGenericTrap*>(objB);
+      CHECK_PROPERTY(trapA, trapB, getZHalfLength);
+      CHECK_PROPERTY(trapA, trapB, getVertices().size);
+      for (size_t v = 0; v < trapA->getVertices().size(); ++v) {
+          if (transCmp(trapA->getVertices()[v], trapB->getVertices()[v])) return true;
+      }
+   } else if (typeID == GeoPcon::getClassTypeID()) {
+     const GeoPcon* pconA = dynamic_cast<const GeoPcon*>(objA);
+     const GeoPcon* pconB = dynamic_cast<const GeoPcon*>(objB);
+     CHECK_PROPERTY(pconA, pconB, getNPlanes);
+     CHECK_PROPERTY(pconA, pconB, getSPhi);
+     CHECK_PROPERTY(pconA, pconB, getDPhi);
+     CHECK_VEC_PROPERTY(pconA, pconB, getZPlane, pconA->getNPlanes());
+     CHECK_VEC_PROPERTY(pconA, pconB, getRMinPlane, pconA->getNPlanes());
+     CHECK_VEC_PROPERTY(pconA, pconB, getRMaxPlane, pconA->getNPlanes());
+   } else if (typeID == GeoPgon::getClassTypeID()) {
+     const GeoPgon* pgonA = dynamic_cast<const GeoPgon*>(objA);
+     const GeoPgon* pgonB = dynamic_cast<const GeoPgon*>(objB);
+     CHECK_PROPERTY(pgonA, pgonB, getNPlanes);
+     CHECK_PROPERTY(pgonA, pgonB, getNSides);
+     CHECK_PROPERTY(pgonA, pgonB, getSPhi);
+     CHECK_PROPERTY(pgonA, pgonB, getDPhi);
+     CHECK_VEC_PROPERTY(pgonA, pgonB, getZPlane, pgonA->getNPlanes());
+     CHECK_VEC_PROPERTY(pgonA, pgonB, getRMinPlane, pgonA->getNPlanes());
+     CHECK_VEC_PROPERTY(pgonA, pgonB, getRMaxPlane, pgonA->getNPlanes());
+   } else if (typeID == GeoSimplePolygonBrep::getClassTypeID()) {
+      const GeoSimplePolygonBrep* brepA = dynamic_cast<const GeoSimplePolygonBrep*>(objA);
+      const GeoSimplePolygonBrep* brepB = dynamic_cast<const GeoSimplePolygonBrep*>(objB);
+      CHECK_PROPERTY(brepA, brepB, getNVertices);
+      CHECK_PROPERTY(brepA, brepB, getDZ);
+      CHECK_VEC_PROPERTY(brepA, brepB, getXVertex, brepA->getNVertices());
+      CHECK_VEC_PROPERTY(brepA, brepB, getYVertex, brepA->getNVertices());
+   } else if (typeID == GeoTorus::getClassTypeID()) {
+      const GeoTorus* torA = dynamic_cast<const GeoTorus*>(objA);
+      const GeoTorus* torB = dynamic_cast<const GeoTorus*>(objB);
+      CHECK_PROPERTY(torA, torB, getRMin);
+      CHECK_PROPERTY(torA, torB, getRMax);
+      CHECK_PROPERTY(torA, torB, getRTor);
+      CHECK_PROPERTY(torA, torB, getSPhi);
+      CHECK_PROPERTY(torA, torB, getDPhi);
+   } else if (typeID == GeoTrap::getClassTypeID()) {
+      const GeoTrap* trapA = dynamic_cast<const GeoTrap*>(objA);
+      const GeoTrap* trapB = dynamic_cast<const GeoTrap*>(objB);
+      CHECK_PROPERTY(trapA, trapB, getZHalfLength);
+      CHECK_PROPERTY(trapA, trapB, getTheta);
+      CHECK_PROPERTY(trapA, trapB, getPhi);
+      CHECK_PROPERTY(trapA, trapB, getDydzn);
+      CHECK_PROPERTY(trapA, trapB, getDxdyndzn);      
+      CHECK_PROPERTY(trapA, trapB, getDxdypdzn);
+      CHECK_PROPERTY(trapA, trapB, getAngleydzn);
+      CHECK_PROPERTY(trapA, trapB, getDydzp);
+      CHECK_PROPERTY(trapA, trapB, getDxdyndzp);
+      CHECK_PROPERTY(trapA, trapB, getDxdypdzp);
+      CHECK_PROPERTY(trapA, trapB, getAngleydzp);
+   } else if (typeID == GeoTwistedTrap::getClassTypeID()) {
+      const GeoTwistedTrap* trapA = dynamic_cast<const GeoTwistedTrap*>(objA);
+      const GeoTwistedTrap* trapB = dynamic_cast<const GeoTwistedTrap*>(objB);
+      CHECK_PROPERTY(trapA, trapB, getY1HalfLength);
+      CHECK_PROPERTY(trapA, trapB, getX1HalfLength);
+      CHECK_PROPERTY(trapA, trapB, getX2HalfLength);
+      CHECK_PROPERTY(trapA, trapB, getY2HalfLength);
+      CHECK_PROPERTY(trapA, trapB, getX3HalfLength);
+      CHECK_PROPERTY(trapA, trapB, getX4HalfLength);
+      CHECK_PROPERTY(trapA, trapB, getZHalfLength);
+      CHECK_PROPERTY(trapA, trapB, getPhiTwist);
+      CHECK_PROPERTY(trapA, trapB, getTheta);
+      CHECK_PROPERTY(trapA, trapB, getPhi);
+      CHECK_PROPERTY(trapA, trapB, getTiltAngleAlpha);
+   } else {
+    THROW_EXCEPTION("The shape "<<objA->type()<<" has not yet been implemented into the sorter");
+  }
+  return 0;
+}
+
+bool GeoComposedShapeSorter::operator()(const GeoShape* a, const GeoShape* b) const {
+    const unsigned int  aComposed{countComposedShapes(a)};
+    const unsigned int  bComposed{countComposedShapes(b)};
+    if (aComposed != bComposed) return aComposed < bComposed;
+    return a->typeID() < b->typeID();
+}
+#undef CHECK_PROPERTY
+#undef CHECK_VEC_PROPERTY
+#undef COMPARE_GEOVEC
+#undef CALL_SORTER
\ No newline at end of file

GeoModelTools/GeoModelFuncSnippets/src/GeoShapeUtils.cxx

+/*
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
+*/
+#include "GeoModelFuncSnippets/GeoShapeUtils.h"
+#include "GeoModelFuncSnippets/TransformToStringConverter.h"
+
+
+/// Boolean volume shapes
+#include "GeoModelKernel/GeoShapeUnion.h"
+#include "GeoModelKernel/GeoShapeSubtraction.h"
+#include "GeoModelKernel/GeoShapeIntersection.h"
+#include "GeoModelKernel/GeoShapeShift.h"
+/// Ordinary shapes
+#include "GeoModelKernel/GeoBox.h"
+#include "GeoModelKernel/GeoTrd.h"
+#include "GeoModelKernel/GeoTube.h"
+#include "GeoModelKernel/GeoTubs.h"
+#include "GeoModelKernel/GeoEllipticalTube.h"
+#include "GeoModelKernel/GeoCons.h"
+#include "GeoModelKernel/GeoPara.h"
+#include "GeoModelKernel/GeoTorus.h"
+
+#include "GeoModelKernel/Units.h"
+
+std::pair<const GeoShape* , const GeoShape*> getOps(const GeoShape* composed) {
+    if (composed->typeID() == GeoShapeUnion::getClassTypeID()) {
+        const GeoShapeUnion* unionShape = dynamic_cast<const GeoShapeUnion*>(composed);
+        return std::make_pair(unionShape->getOpA(), unionShape->getOpB());
+    } else if (composed->typeID() == GeoShapeSubtraction::getClassTypeID()) {
+        const GeoShapeSubtraction* shapeSubtract = dynamic_cast<const GeoShapeSubtraction*>(composed);
+        return std::make_pair(shapeSubtract->getOpA(), shapeSubtract->getOpB());
+    } else if (composed->typeID() == GeoShapeUnion::getClassTypeID()) {
+        const GeoShapeIntersection* shapeIntersect = dynamic_cast<const GeoShapeIntersection*>(composed);
+        return std::make_pair(shapeIntersect->getOpA(), shapeIntersect->getOpB());
+    } else if (composed->typeID() == GeoShapeShift::getClassTypeID()) {
+        const GeoShapeShift* shapeShift = dynamic_cast<const GeoShapeShift*>(composed);
+        return std::make_pair(shapeShift->getOp(), nullptr);
+    }
+    return std::make_pair(nullptr, nullptr);
+}
+
+unsigned int countComposedShapes(const GeoShape* shape) {
+    std::pair<const GeoShape*, const GeoShape*> ops = getOps(shape);
+    return 1 + (ops.first ? countComposedShapes(ops.first) : 0) +
+                (ops.second ? countComposedShapes(ops.second) : 0);
+}
+
+GeoIntrusivePtr<const GeoShape> compressShift(const GeoShape* shift) {    
+    if (shift->typeID() != GeoShapeShift::getClassTypeID()) return GeoIntrusivePtr<const GeoShape>{shift};
+    const GeoShapeShift* shapeShift = dynamic_cast<const GeoShapeShift*>(shift);
+    if (shapeShift->getOp()->typeID() != GeoShapeShift::getClassTypeID()) return GeoIntrusivePtr<const GeoShape>{shift};
+    GeoIntrusivePtr<const GeoShape> subShape{compressShift(shapeShift->getOp())};
+    const GeoShapeShift* subShift = dynamic_cast<const GeoShapeShift*>(subShape.get());
+    return GeoIntrusivePtr<const GeoShape>{new GeoShapeShift(subShift->getOp(), subShift->getX() * shapeShift->getX())};
+}
+std::vector<const GeoShape*> getBooleanComponents(const GeoShape* booleanShape) {
+     std::pair<const GeoShape*, const GeoShape*> operands = getOps(booleanShape);
+     if (!operands.first || !operands.second){
+        return {compressShift(booleanShape)};
+     }
+     std::vector<const GeoShape*> components{};
+     if (booleanShape->typeID() != GeoShapeSubtraction::getClassTypeID()) {
+        components = getBooleanComponents(operands.first);
+        std::vector<const GeoShape*> secCmp = getBooleanComponents(operands.second);
+        components.insert(components.end(), std::make_move_iterator(secCmp.begin()), 
+                                            std::make_move_iterator(secCmp.end()));
+     } else {
+        if (operands.first->typeID() == GeoShapeSubtraction::getClassTypeID()) {
+            components = getBooleanComponents(operands.first);
+        }
+        components.push_back(compressShift(operands.second));             
+     }
+     return components;
+}
+
+
+
+std::string printGeoShape(const GeoShape* shape) {
+    std::stringstream ostr{};
+    constexpr double toDeg{1./GeoModelKernelUnits::deg};
+    ostr<<shape->type()<<" ("<<shape<<") ";
+    const int typeID = shape->typeID();
+    if (typeID == GeoShapeUnion::getClassTypeID()) {
+        const GeoShapeUnion* shapeUnion = dynamic_cast<const GeoShapeUnion*>(shape);
+        ostr<<"of {"<<printGeoShape(shapeUnion->getOpA())<<"} & {"<<printGeoShape(shapeUnion->getOpB())<<"}";        
+    } else if (typeID == GeoShapeSubtraction::getClassTypeID()) {
+        const GeoShapeSubtraction* subtractShape = dynamic_cast<const GeoShapeSubtraction*>(shape);
+        ostr<<"of {"<<printGeoShape(subtractShape->getOpB())<<"} from {"<<printGeoShape(subtractShape->getOpA())<<"}";
+    } else if (typeID == GeoShapeIntersection::getClassTypeID()) {
+        const GeoShapeIntersection* intersectShape = dynamic_cast<const GeoShapeIntersection*>(shape);
+        ostr<<"between {"<<printGeoShape(intersectShape->getOpA())<<"} & {"<<printGeoShape(intersectShape->getOpB())<<"}";        
+    } else if (typeID == GeoShapeShift::getClassTypeID()) {
+        const GeoShapeShift* shiftShape = dynamic_cast<const GeoShapeShift*>(shape);
+        ostr<<"of "<<printGeoShape(shiftShape->getOp())<<" by "<<GeoTrf::toString(shiftShape->getX()); 
+    }
+    /// Elementary shape types
+    else if (typeID == GeoBox::getClassTypeID()) {
+        const GeoBox* boxShape = dynamic_cast<const GeoBox*>(shape);
+        ostr<<" halfX="<<boxShape->getXHalfLength()<<", halfY="<<boxShape->getYHalfLength()<<", halfZ="<<boxShape->getZHalfLength();
+    } else if (typeID == GeoTrd::getClassTypeID()) {
+        const GeoTrd* trd = dynamic_cast<const GeoTrd*>(shape);
+        ostr<<"halfX="<<trd->getXHalfLength1()<<"/"<<trd->getXHalfLength2()<<", ";
+        ostr<<"halfY="<<trd->getYHalfLength1()<<"/"<<trd->getYHalfLength2()<<", ";
+        ostr<<"halfZ="<<trd->getZHalfLength();
+   } else if (typeID == GeoTube::getClassTypeID()) {
+        const GeoTube* tube = dynamic_cast<const GeoTube*>(shape);
+        ostr<<"rMin="<<tube->getRMin()<<", ";
+        ostr<<"rMax="<<tube->getRMax()<<", ";
+        ostr<<"halfZLength="<<tube->getZHalfLength();
+    } else if (typeID == GeoTubs::getClassTypeID()) {
+        const GeoTubs* tube = dynamic_cast<const GeoTubs*>(shape);
+        ostr<<"rMin="<<tube->getRMin()<<", ";
+        ostr<<"rMax="<<tube->getRMax()<<", ";
+        ostr<<"start phi="<<tube->getSPhi()*toDeg<<", ";
+        ostr<<"dPhi="<<tube->getDPhi()*toDeg<<", ";
+        ostr<<"halfZ="<<tube->getZHalfLength();
+    } else if (typeID == GeoCons::getClassTypeID()) {
+        const GeoCons* cons = dynamic_cast<const GeoCons*>(shape);
+        ostr<<"r1="<<cons->getRMin1()<<"/"<<cons->getRMax1()<<", ";
+        ostr<<"r2="<<cons->getRMin2()<<"/"<<cons->getRMax2()<<", ";
+        ostr<<"dZ="<<cons->getDZ()<<", ";
+        ostr<<"start phi="<<cons->getSPhi()*toDeg<<", ";
+        ostr<<"dPhi="<<cons->getDPhi()*toDeg;
+   } else if (typeID == GeoEllipticalTube::getClassTypeID()) {
+        const GeoEllipticalTube* tube = dynamic_cast<const GeoEllipticalTube*>(shape);
+        ostr<<"halfX="<<tube->getXHalfLength()<<", halfY="<<tube->getYHalfLength()<<", halfZ="<<tube->getZHalfLength();
+   } else if (typeID == GeoPara::getClassTypeID()) {
+        const GeoPara* para = dynamic_cast<const GeoPara*>(shape);
+        ostr<<"halfX="<<para->getXHalfLength()<<", halfY="<<para->getYHalfLength()<<", halfZ="<<para->getZHalfLength();
+        ostr<<"theta="<<para->getTheta()*toDeg<<", alpha="<<para->getAlpha()*toDeg<<", phi="<<para->getPhi()*toDeg; 
+   } else if (typeID == GeoTorus::getClassTypeID()) {
+      const GeoTorus* torus = dynamic_cast<const GeoTorus*>(shape);
+      ostr<<"r="<<torus->getRMin()<<"/"<<torus->getRMax()<<", ";
+      ostr<<"Torus R="<<torus->getRTor()<<", ";
+      ostr<<"sPhi="<<torus->getSPhi()*toDeg<<", ";
+      ostr<<"dPhi="<<torus->getDPhi()*toDeg;
+   }
+    return ostr.str();
+}
\ No newline at end of file

GeoModelTools/GeoModelFuncSnippets/src/TransformSorter.cxx

+/*
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
+*/
+#include "GeoModelFuncSnippets/TransformSorter.h"
+
+
+namespace GeoTrf {
+    bool TransformSorter::operator()(const Transform3D& a, const Transform3D& b) const{
+        return compare(a, b) < 0;
+    }
+    bool TransformSorter::operator()(const RotationMatrix3D&a, const RotationMatrix3D& b) const{
+        return compare(a, b) < 0;
+    }
+
+    int TransformSorter::compare(const Transform3D&a, const Transform3D& b) const{ 
+        const Vector3D transA{a.translation()}, transB{b.translation()};
+        const int transCmp = compare(transA, transB);
+        if (transCmp) return transCmp;
+        return compare(a.linear(), b.linear());
+    }
+    int TransformSorter::compare(const RotationMatrix3D&a, const RotationMatrix3D& b) const{ 
+        return getGeoRotationAngles(a).compare(getGeoRotationAngles(b));
+    }
+}
\ No newline at end of file

GeoModelTools/GeoModelFuncSnippets/src/TransformToStringConverter.cxx

+/*
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
+*/
+#include "GeoModelFuncSnippets/TransformToStringConverter.h"
+#include "GeoModelKernel/Units.h"
+
+namespace{
+    constexpr double toDeg = 1./ GeoModelKernelUnits::deg;
+}
+
+namespace GeoTrf{
+    std::string toString(const Transform3D& transform, bool useCoordAngles, int precision) {
+        std::stringstream ostr{};
+        ostr << std::setiosflags(std::ios::fixed) << std::setprecision(precision);
+        const Vector3D trans = transform.translation();
+        bool printed{false};
+        if (std::sqrt(trans.dot(trans)) > 0.1 *GeoModelKernelUnits::micrometer) {
+            ostr<<"translation -- ";
+            ostr<<toString(trans, precision);
+            printed = true;
+        }
+        if (useCoordAngles) {
+            CoordEulerAngles angles = getCoordRotationAngles(transform);
+            if (angles) {
+                if (printed) ostr<<", ";
+                ostr<<"Euler angles -- ";
+                ostr<<toString(angles, precision); 
+                printed = true;               
+            }
+        } else {
+            EulerAngles angles = getGeoRotationAngles(transform);
+            if (angles){
+                 if(printed) ostr<<", ";
+                 ostr<<"Euler angles -- ";
+                 ostr<<toString(angles, precision);
+                 printed = true;                 
+            }
+        }
+
+        if (!printed) ostr<<"Identity transformation";
+        return ostr.str();
+    }  
+    std::string toString(const EulerAngles& angles, int precision) {
+        std::stringstream ostr{};
+        ostr << std::setiosflags(std::ios::fixed) << std::setprecision(precision);
+        ostr<<"phi: "<<angles.phi *toDeg<<", ";
+        ostr<<"theta: "<<angles.theta *toDeg<<", ";
+        ostr<<"psi: "<<angles.psi *toDeg;
+        return ostr.str();
+    }
+    std::string toString(const CoordEulerAngles& angles, int precision) {
+        std::stringstream ostr{};
+        ostr << std::setiosflags(std::ios::fixed) << std::setprecision(precision);
+        ostr <<"alpha: "<<angles.alpha * toDeg<<", ";
+        ostr <<"beta: "<<angles.beta * toDeg<<", ";
+        ostr <<"gamma: "<<angles.gamma * toDeg;        
+        return ostr.str();
+    }
+}

GeoModelTools/GeoModelFuncSnippets/src/getChildNodesWithTrf.cxx

+/*
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
+*/
+#include "GeoModelFuncSnippets/getChildNodesWithTrf.h"
+
+#include "GeoModelFuncSnippets/GeoPhysVolSorter.h"
+#include "GeoModelFuncSnippets/TransformSorter.h"
+#include "GeoModelKernel/GeoVolumeCursor.h"
+#include "GeoModelKernel/GeoFullPhysVol.h"
+
+namespace {
+    constexpr std::string_view dummyNodeName{"ANON"};
+}
+
+
+ GeoChildNodeWithTrf::GeoChildNodeWithTrf(GeoVolumeCursor& curs):
+    transform{curs.getTransform()}, 
+    volume{curs.getVolume()},
+    nodeName{curs.getName()},
+    isAlignable{curs.hasAlignableTransform()},
+    isSensitive{typeid(*volume) == typeid(GeoFullPhysVol)} {
+    //// Do not specify a node name if it's a dummy one
+    if (nodeName == dummyNodeName) {
+        nodeName = volume->getLogVol()->getName();
+
+    }
+}
+
+std::vector <GeoChildNodeWithTrf> getChildrenWithRef(PVConstLink& physVol,
+                                                     bool summarizeEqualVol) {
+    std::vector<GeoChildNodeWithTrf> children{};
+ 
+    static const GeoPhysVolSorter physVolSort{};
+    static const GeoTrf::TransformSorter transSort{};
+
+    
+    GeoVolumeCursor cursor{physVol};
+    GeoTrf::Transform3D lastNodeTrf{GeoTrf::Transform3D::Identity()};
+
+    while (!cursor.atEnd()) {
+        if (children.empty() || !summarizeEqualVol) {
+            children.emplace_back(cursor);
+            cursor.next();
+            continue;
+         }
+         GeoChildNodeWithTrf& prevChild{children.back()};
+         GeoChildNodeWithTrf currentChild{cursor};
+         if (prevChild.isAlignable != currentChild.isAlignable ||
+             prevChild.isSensitive != currentChild.isSensitive ||
+             physVolSort.compare(prevChild.volume, currentChild.volume)) {
+            children.emplace_back(std::move(currentChild));
+         } else if (prevChild.nCopies == 1) {
+            ++prevChild.nCopies;
+            prevChild.inductionRule = prevChild.transform.inverse() * 
+                                      currentChild.transform;            
+         } else if (!transSort.compare(prevChild.inductionRule, 
+                                      lastNodeTrf.inverse() * currentChild.transform)) {
+            ++prevChild.nCopies;
+        } else {
+            children.emplace_back(std::move(currentChild));
+        }
+        lastNodeTrf = cursor.getTransform();   
+        cursor.next();
+    }
+    return children;
+}

GeoModelTools/GeoModelFuncSnippets/tests/testGeoShapeSet.cxx

+/*
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
+*/
+#include "GeoModelKernel/Units.h"
+#include "GeoModelKernel/GeoDefinitions.h"
+#include "GeoModelFuncSnippets/GeoShapeSorter.h"
+#include "GeoModelFuncSnippets/GeoShapeUtils.h"
+
+#include "GeoModelKernel/GeoBox.h"
+#include "GeoModelKernel/GeoTube.h"
+#include "GeoModelKernel/GeoShapeShift.h"
+
+
+#include <stdlib.h>
+#include <iostream>
+
+template<class T, typename... args> bool insertVolume(const bool printFailure,
+                                                     GeoShapeSet<GeoShape>& set,
+                                                     args... a) {
+    GeoIntrusivePtr<T> newPtr{new T(a...)};
+    
+    const auto insert_itr = set.insert(newPtr);
+    if (insert_itr.second && newPtr->refCount() == 2) return true;
+    if (printFailure) {
+        std::cerr<<"Shape insertion failed "<<printGeoShape(newPtr)<<std::endl;
+        std::cerr<<"Shape "<<printGeoShape(*insert_itr.first)<<" is taking its place "<<std::endl;
+    }
+    return false;
+}
+
+int main() {
+    GeoShapeSet<GeoShape> shapeColl{};
+    GeoShapeSorter shapeSorter{};
+
+    if (!insertVolume<GeoBox>(true, shapeColl, 10.,10.,10.)) {
+        std::cerr<<"testGeoShapeSorter() "<<__LINE__<<" test failed. "<<std::endl;
+        return EXIT_FAILURE;
+    }
+    if (insertVolume<GeoBox>(false, shapeColl, 10.,10.,10.)) {
+        std::cerr<<"testGeoShapeSorter() "<<__LINE__<<" GeoBox A has been inserted twice. "<<std::endl;
+        return EXIT_FAILURE;
+    }
+    if (!insertVolume<GeoBox>(true, shapeColl, 11.,10., 10.)) {
+        std::cerr<<"testGeoShapeSorter() "<<__LINE__<<" test failed. "<<std::endl;
+        return EXIT_FAILURE;
+    }
+    if (!insertVolume<GeoBox>(true, shapeColl, 10., 11., 10.)) {
+        std::cerr<<"testGeoShapeSorter() "<<__LINE__<<" test failed. "<<std::endl;
+        return EXIT_FAILURE;
+    }
+
+    if (!insertVolume<GeoBox>(true, shapeColl, 10., 10., 11.)) {
+        std::cerr<<"testGeoShapeSorter() "<<__LINE__<<" test failed. "<<std::endl;
+        return EXIT_FAILURE;
+    }
+
+    if (!insertVolume<GeoBox>(true, shapeColl, 10., 11., 11.)) {
+        std::cerr<<"testGeoShapeSorter() "<<__LINE__<<" test failed. "<<std::endl;
+        return EXIT_FAILURE;
+    }
+    ///##############################################################
+    ///                 GeoTube
+    ///##############################################################
+    if (!insertVolume<GeoTube>(true, shapeColl, 10.,10.,10.)) {
+        std::cerr<<"testGeoShapeSorter() "<<__LINE__<<" test failed. "<<std::endl;
+        return EXIT_FAILURE;
+    }
+    if (insertVolume<GeoTube>(false, shapeColl, 10.,10.,10.)) {
+        std::cerr<<"testGeoShapeSorter() "<<__LINE__<<" GeoTube A has been inserted twice. "<<std::endl;
+        return EXIT_FAILURE;
+    }
+    if (!insertVolume<GeoTube>(true, shapeColl, 11.,10., 10.)) {
+        std::cerr<<"testGeoShapeSorter() "<<__LINE__<<" test failed. "<<std::endl;
+        return EXIT_FAILURE;
+    }
+    if (!insertVolume<GeoTube>(true, shapeColl, 10., 11., 10.)) {
+        std::cerr<<"testGeoShapeSorter() "<<__LINE__<<" test failed. "<<std::endl;
+        return EXIT_FAILURE;
+    }
+
+    if (!insertVolume<GeoTube>(true, shapeColl, 10., 10., 11.)) {
+        std::cerr<<"testGeoShapeSorter() "<<__LINE__<<" test failed. "<<std::endl;
+        return EXIT_FAILURE;
+    }
+
+    if (!insertVolume<GeoTube>(true, shapeColl, 10., 11., 11.)) {
+        std::cerr<<"testGeoShapeSorter() "<<__LINE__<<" test failed. "<<std::endl;
+        return EXIT_FAILURE;
+    }
+
+    {
+        GeoIntrusivePtr<GeoTube> tubeShift{new GeoTube(2,3,4)};
+        GeoTrf::Transform3D trans1{GeoTrf::GeoTransformRT{GeoTrf::GeoRotation(0., 2., 2.), GeoTrf::Vector3D{1.,0.,3}}};
+        GeoTrf::Transform3D trans2{GeoTrf::GeoTransformRT{GeoTrf::GeoRotation(2.,0., 2.), GeoTrf::Vector3D{2.,0.,3}}};
+        if (!insertVolume<GeoShapeShift>(true, shapeColl, tubeShift, trans1)){
+            std::cerr<<"testGeoShapeSorter() "<<__LINE__<<" test failed. "<<std::endl;
+            return EXIT_FAILURE;
+        }
+        if (!insertVolume<GeoShapeShift>(true, shapeColl, tubeShift, trans2)){
+            std::cerr<<"testGeoShapeSorter() "<<__LINE__<<" test failed. "<<std::endl;
+            return EXIT_FAILURE;
+        }
+    }
+    return EXIT_SUCCESS;
+}
\ No newline at end of file

GeoModelTools/GeoModelFuncSnippets/tests/testTransformSorter.cxx

+/*
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "GeoModelFuncSnippets/TransformSorter.h"
+#include "GeoModelFuncSnippets/TransformToStringConverter.h"
+#include "GeoModelKernel/Units.h"
+#include <stdlib.h>
+#include <iostream>
+
+std::ostream& operator<<(std::ostream& ostr, const GeoTrf::Transform3D& trans) {
+    ostr<<GeoTrf::toString(trans, true, 2);
+    return ostr;
+}
+
+#define TEST_TRANSFORM(TRANS) \
+    std::cout<<"testTransformSorter() "<<__LINE__<<"  Test transformation "<<TRANS<<std::endl; \
+    if (sorter(TRANS, TRANS)) { \
+        std::cerr<<"testTransformSorter() "<<__LINE__<<" transformation "<<TRANS<<" is less than itself "<<std::endl; \
+        return EXIT_FAILURE; \
+    }
+#define COMPARE_SORTER(TRANSA, TRANSB)         \
+    if (sorter(TRANSA, TRANSB) == sorter(TRANSB, TRANSA)) {                                   \
+        std::cerr<<"testTransformSorter() "<<__LINE__<<" transformations "                    \
+                  <<TRANSA<<" & "<<TRANSB<<" are identical. "<<std::endl;                     \
+        const GeoTrf::Vector3D vecA{TRANSA.translation()}, vecB{TRANSB.translation()};        \
+        std::cerr<<"testTransformSorter() comparison result (A,B): "<<sorter(TRANSA, TRANSB)  \
+                 << " compare result: "<<sorter.compare(TRANSA, TRANSB)                       \
+                 << " -- linear part: "<<sorter(TRANSA.linear(), TRANSB.linear()) << " / "    \
+                 << sorter.compare(TRANSA.linear(), TRANSB.linear())                          \
+                 << " -- translation part: "<<sorter(vecA, vecB)                              \
+                 << " / "<<sorter.compare(vecA, vecB)<<std::endl;                             \
+        std::cerr<<"testTransformSorter() comparison result (B,A): "<<sorter(TRANSB, TRANSA)  \
+                 << " compare result: "<<sorter.compare(TRANSB, TRANSA)                       \
+                 << " -- linear part: "<<sorter(TRANSB.linear(), TRANSA.linear()) << " / "    \
+                 << sorter.compare(TRANSB.linear(), TRANSA.linear())                          \
+                 << " -- translation part: "<<sorter(vecB, vecA)                              \
+                 << " / "<<sorter.compare(vecB, vecA)<<std::endl;                             \
+        return EXIT_FAILURE;                                                                  \
+    }
+int main() {
+    constexpr double deg = GeoModelKernelUnits::deg;
+    GeoTrf::TransformSorter sorter{};
+    GeoTrf::CoordEulerAngles anglesA{45.*deg, 0., 0.};
+    GeoTrf::Transform3D transA{GeoTrf::GeoTransformRT{GeoTrf::GeoRotation{anglesA}, 10. *GeoTrf::Vector3D::UnitX()}};
+    TEST_TRANSFORM(transA);
+
+    GeoTrf::CoordEulerAngles anglesB{0.,0.,0.};
+    GeoTrf::Transform3D transB{GeoTrf::GeoTransformRT{GeoTrf::GeoRotation{anglesB}, 10. *GeoTrf::Vector3D::UnitX()}};
+    TEST_TRANSFORM(transB);
+
+    COMPARE_SORTER(transA, transB);
+   
+
+    GeoTrf::Transform3D transC{GeoTrf::Translate3D{1, -1, 0.}};
+    TEST_TRANSFORM(transC)
+    COMPARE_SORTER(transA, transC);
+    COMPARE_SORTER(transB, transC);
+   
+    return EXIT_SUCCESS;
+}
\ No newline at end of file