diff --git a/GeoModelCore/GeoModelHelpers/GeoModelHelpers/GeoShapeSorter.h b/GeoModelCore/GeoModelHelpers/GeoModelHelpers/GeoShapeSorter.h
index a18006c36a16d4088b6e08566a30fa19dfe0ebdb..a7214d78ff32c8e532a63fcaf6b280438aa1d4dd 100644
--- a/GeoModelCore/GeoModelHelpers/GeoModelHelpers/GeoShapeSorter.h
+++ b/GeoModelCore/GeoModelHelpers/GeoModelHelpers/GeoShapeSorter.h
@@ -28,15 +28,6 @@ struct GeoShapeSorter {
 
 };
 
-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>
diff --git a/GeoModelCore/GeoModelHelpers/src/GeoShapeSorter.cxx b/GeoModelCore/GeoModelHelpers/src/GeoShapeSorter.cxx
index 5d1aee7624151e8b5fc1b276836bd425052542e4..b103d3ef97055e549b87bb5f40bd706f71d2166e 100644
--- a/GeoModelCore/GeoModelHelpers/src/GeoShapeSorter.cxx
+++ b/GeoModelCore/GeoModelHelpers/src/GeoShapeSorter.cxx
@@ -83,6 +83,7 @@ int GeoShapeSorter::compare(const GeoShape* objA, const GeoShape* objB) const {
     
     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);
   } 
@@ -156,7 +157,8 @@ int GeoShapeSorter::compare(const GeoShape* objA, const GeoShape* objB) const {
       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;
+          const int compare = transCmp.compare(trapA->getVertices()[v], trapB->getVertices()[v]);
+          if (compare) return compare;
       }
    } else if (typeID == GeoPcon::getClassTypeID()) {
      const GeoPcon* pconA = dynamic_cast<const GeoPcon*>(objA);
@@ -226,12 +228,6 @@ int GeoShapeSorter::compare(const GeoShape* objA, const GeoShape* objB) const {
   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
diff --git a/GeoModelCore/GeoModelHelpers/src/GeoShapeUtils.cxx b/GeoModelCore/GeoModelHelpers/src/GeoShapeUtils.cxx
index 2117315dddf4126c109e9e629702de12f4ce8fad..b0a6f8d6f369220ec737edfd8b23dab95ada8cac 100644
--- a/GeoModelCore/GeoModelHelpers/src/GeoShapeUtils.cxx
+++ b/GeoModelCore/GeoModelHelpers/src/GeoShapeUtils.cxx
@@ -21,6 +21,7 @@
 #include "GeoModelKernel/GeoPara.h"
 #include "GeoModelKernel/GeoTorus.h"
 #include "GeoModelKernel/GeoSimplePolygonBrep.h"
+#include "GeoModelKernel/GeoGenericTrap.h"
 
 #include "GeoModelKernel/Units.h"
 
@@ -33,7 +34,7 @@ std::pair<const GeoShape* , const GeoShape*> getOps(const GeoShape* composed) {
     } 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()) {
+    } else if (composed->typeID() == GeoShapeIntersection::getClassTypeID()) {
         const GeoShapeIntersection* shapeIntersect = dynamic_cast<const GeoShapeIntersection*>(composed);
         return std::make_pair(shapeIntersect->getOpA(), shapeIntersect->getOpB());
     } else if (composed->typeID() == GeoShapeShift::getClassTypeID()) {
@@ -82,7 +83,7 @@ std::vector<const GeoShape*> getBooleanComponents(const GeoShape* booleanShape)
 std::string printGeoShape(const GeoShape* shape) {
     std::stringstream ostr{};
     constexpr double toDeg{1./GeoModelKernelUnits::deg};
-    ostr<<shape->type()<<" ("<<shape<<") ";
+    ostr<<shape->type()<<" ("<<shape<<"/"<<shape->refCount()<<") ";
     const int typeID = shape->typeID();
     if (typeID == GeoShapeUnion::getClassTypeID()) {
         const GeoShapeUnion* shapeUnion = dynamic_cast<const GeoShapeUnion*>(shape);
@@ -138,6 +139,12 @@ std::string printGeoShape(const GeoShape* shape) {
       ostr<<"Torus R="<<torus->getRTor()<<", ";
       ostr<<"sPhi="<<torus->getSPhi()*toDeg<<", ";
       ostr<<"dPhi="<<torus->getDPhi()*toDeg;
+   } else if (typeID == GeoGenericTrap::getClassTypeID()) {
+        const GeoGenericTrap* trap = dynamic_cast<const GeoGenericTrap*>(shape);
+        ostr<<"half Z: "<<trap->getZHalfLength()<<" n Vertices: "<<trap->getVertices().size()<<std::endl;
+        for ( const GeoTrf::Vector2D& vec : trap->getVertices()) {
+            ostr<<"  **** "<<GeoTrf::toString(vec)<<std::endl;
+        }
    }
     return ostr.str();
 }
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/AddPlane.h b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/AddPlane.h
index 5bdf2cf62e495f35dcfc2c6e42bf268d5168353a..b604566e81cde3a1aa19919f89a500ab82075245 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/AddPlane.h
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/AddPlane.h
@@ -14,8 +14,8 @@ class GmxUtil;
 
 class AddPlane {
 public:
-    AddPlane() {};
-    ~AddPlane() {};
+    AddPlane() = default;
+    ~AddPlane() = default;
     void process(const xercesc::DOMElement *element, double &zPlane, double &rMinPlane, double &rMaxPlane);
     GmxUtil* gmxUtil=nullptr;
 };
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeBox.h b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeBox.h
index 968bc48068bf509c0cf010480d124cee416c04ee..220b875c622379d5f6fdda0fd237153f23252850 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeBox.h
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeBox.h
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
 */
 
 //
@@ -12,10 +12,12 @@
 
 #include "GeoModelXml/Element2GeoItem.h"
 
+#include "GeoModelKernel/GeoBox.h"
 class MakeBox: public Element2GeoItem {
 public:
-    MakeBox();
+    MakeBox() = default;
     virtual RCBase * make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const override;
+
 };
 
 #endif // GEO_MODEL_XML_MAKE_BOX_H
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeCons.h b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeCons.h
index 7de87d2d0a02f67cd7e836c2ee273bacc37661a0..7fb7d0ebf5651973966906f93bc07d374a2680a2 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeCons.h
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeCons.h
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
 */
 
 //
@@ -13,9 +13,10 @@
 #include "GeoModelXml/Element2GeoItem.h"
 
 class MakeCons: public Element2GeoItem {
-public:
-    MakeCons();
-    virtual RCBase * make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const override;
+  public:
+      MakeCons() = default;
+      virtual RCBase * make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const override;
+
 };
 
 #endif // GEO_MODEL_XML_MAKE_CONS_H
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeEllipticalTube.h b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeEllipticalTube.h
index 11d894ca7274c2539d37bcc7805ff0e44b3e3f00..9609bc2e871dd47bbe0ef31f677d245fc1cac5c5 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeEllipticalTube.h
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeEllipticalTube.h
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
 */
 
 #ifndef GEO_MODEL_XML_MAKE_ELLIPTICAL_TUBE_H
@@ -10,8 +10,9 @@
 
 class MakeEllipticalTube: public Element2GeoItem {
 public:
-    MakeEllipticalTube();
+    MakeEllipticalTube() = default;;
     virtual RCBase * make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const override;
+
 };
 
 #endif // GEO_MODEL_XML_MAKE_ELLIPTICAL_TUBE_H
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeGenericTrap.h b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeGenericTrap.h
index 83b011534ea981d6f8bdde4f3c0cfb3bb8034536..4e33e9e081995d49641d71fa9ec7565735851a8c 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeGenericTrap.h
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeGenericTrap.h
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
 */
 
 //
@@ -12,9 +12,10 @@
 
 #include "GeoModelXml/Element2GeoItem.h"
 
+
 class MakeGenericTrap: public Element2GeoItem {
 public:
-    MakeGenericTrap();
+    MakeGenericTrap() = default;
     virtual RCBase * make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const override;
 };
 
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeIntersection.h b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeIntersection.h
index c8d7b71d9ac5b107a93d297afaa8bae1144dcc74..f3ebfdf3aa0494fee0cac3963737bfb88d8ebe5c 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeIntersection.h
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeIntersection.h
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
 */
 
 //
@@ -12,10 +12,12 @@
 
 #include "GeoModelXml/Element2GeoItem.h"
 
+
 class MakeIntersection: public Element2GeoItem {
 public:
-    MakeIntersection();
+    MakeIntersection() = default;
     virtual RCBase * make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const override;
+
 };
 
 #endif // GEO_MODEL_XML_MAKE_INTERSECTION_H
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakePara.h b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakePara.h
index c340ca70593702cf2fd9cf47b1c56793d29d8338..3cf210eb0dfec4df5cf3b8fe148a25c93635586e 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakePara.h
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakePara.h
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
 */
 
 //
@@ -13,9 +13,10 @@
 #include "GeoModelXml/Element2GeoItem.h"
 
 class MakePara: public Element2GeoItem {
-public:
-    MakePara();
+  public:
+    MakePara() = default;
     virtual RCBase * make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const override;
+   
 };
 
 #endif // GEO_MODEL_XML_MAKE_PARA_H
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakePcon.h b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakePcon.h
index 190badf08e30646822f366915ab6bc10c31744b1..62338329662c55ce19abb1cf033dc98abe791fa2 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakePcon.h
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakePcon.h
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
 */
 
 //
@@ -12,10 +12,13 @@
 
 #include "GeoModelXml/Element2GeoItem.h"
 
+
+
 class MakePcon: public Element2GeoItem {
 public:
-    MakePcon();
+    MakePcon() = default;
     virtual RCBase * make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const override;
+
 };
 
 #endif // GEO_MODEL_XML_MAKE_PCON_H
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakePgon.h b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakePgon.h
index cc48599aa884a7785c9ece17ba2bf2057bc72a8c..2b36bb38a68c458d9d72940609cc59e358b9dadc 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakePgon.h
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakePgon.h
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
 */
 
 //
@@ -12,10 +12,12 @@
 
 #include "GeoModelXml/Element2GeoItem.h"
 
+
 class MakePgon: public Element2GeoItem {
-public:
-    MakePgon();
-    virtual RCBase * make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const override;
+  public:
+      MakePgon() = default;
+      virtual RCBase * make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const override;
+
 };
 
 #endif // GEO_MODEL_XML_MAKE_PGON_H
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeSimplePolygonBrep.h b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeSimplePolygonBrep.h
index 92382bef94823fc60e070a6c810e369339c22bb5..9a4f2897bc9235bbb69bcb17927f5417c6970ac7 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeSimplePolygonBrep.h
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeSimplePolygonBrep.h
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
 */
 #include <xercesc/util/XercesDefs.hpp>
 //
@@ -13,7 +13,7 @@
 
 class MakeSimplePolygonBrep: public Element2GeoItem {
 public:
-    MakeSimplePolygonBrep();
+    MakeSimplePolygonBrep() = default;
     virtual RCBase * make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const override;
 };
 
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeSubtraction.h b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeSubtraction.h
index 359b43037de6c593f2f13581144f616c369a7e97..cb9ac7391d773e90245dbbc75b4ee27a918f4332 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeSubtraction.h
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeSubtraction.h
@@ -12,10 +12,12 @@
 
 #include "GeoModelXml/Element2GeoItem.h"
 
+
 class MakeSubtraction: public Element2GeoItem {
 public:
-    MakeSubtraction();
+    MakeSubtraction() = default;
     virtual RCBase * make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const override;
+
 };
 
 #endif // GEO_MODEL_XML_MAKE_SUBTRACTION_H
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeTorus.h b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeTorus.h
index 5d269cb7ed9ec9839d4afb33ca7803671c95e076..2425818b08623e29d9b909796ca5915796a6ca89 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeTorus.h
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeTorus.h
@@ -9,10 +9,12 @@
 
 #include "GeoModelXml/Element2GeoItem.h"
 
+
 class MakeTorus: public Element2GeoItem {
 public:
-    MakeTorus();
+    MakeTorus() = default;
     virtual RCBase * make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const override;
+
 };
 
 #endif // GEO_MODEL_XML_MAKE_TORUS_H
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeTrap.h b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeTrap.h
index 2832814597c75fb22f96aaa75c91d9ee3207dc6a..3c9d4435da320c2a4e1dbd15692c4462f283d5fd 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeTrap.h
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeTrap.h
@@ -14,7 +14,7 @@
 
 class MakeTrap: public Element2GeoItem {
 public:
-    MakeTrap();
+    MakeTrap() = default;
     virtual RCBase * make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const override;
 };
 
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeTrd.h b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeTrd.h
index ae05777cdd8e06dd21ad940118e07f62fe054710..cf758804b5df68e9c58e07ed8513f9c94c407f05 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeTrd.h
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeTrd.h
@@ -14,7 +14,7 @@
 
 class MakeTrd: public Element2GeoItem {
 public:
-    MakeTrd();
+    MakeTrd() = default;
     virtual RCBase * make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const override;
 };
 
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeTube.h b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeTube.h
index bdf81d3c207186ef38ba65dfba5ce840adff492b..eebad477cbadde207e2ef9c9dc905db0567b8fb3 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeTube.h
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeTube.h
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
 */
 
 //
@@ -12,10 +12,12 @@
 
 #include "GeoModelXml/Element2GeoItem.h"
 
+
 class MakeTube: public Element2GeoItem {
 public:
-    MakeTube();
+    MakeTube() = default;
     virtual RCBase * make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const override;
+
 };
 
 #endif // GEO_MODEL_XML_MAKE_TUBE_H
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeTubs.h b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeTubs.h
index a75b172635a4a5bd7a49570a91adf584e898de66..7d3520c18f1a8f9a9b8632d680a56cc1beb1ea6a 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeTubs.h
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeTubs.h
@@ -12,9 +12,10 @@
 
 #include "GeoModelXml/Element2GeoItem.h"
 
+
 class MakeTubs: public Element2GeoItem {
 public:
-    MakeTubs();
+    MakeTubs() = default;
     virtual RCBase * make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const override;
 };
 
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeTwistedTrap.h b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeTwistedTrap.h
index 984ae165a72ffd822c43b1811a4a07dc2e258da4..54a88ac91089d607461020afafdf4f902fc7e4c0 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeTwistedTrap.h
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeTwistedTrap.h
@@ -8,10 +8,12 @@
 
 #include "GeoModelXml/Element2GeoItem.h"
 
+
 class MakeTwistedTrap: public Element2GeoItem {
 public:
-    MakeTwistedTrap();
+    MakeTwistedTrap() = default;
     virtual RCBase * make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const override;
+
 };
 
 #endif // GEO_MODEL_XML_MAKE_TWISTED_TRAP_H
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeUnion.h b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeUnion.h
index cb5f32fc1ea58e9b140f3c212769ee41ae375acd..2606df77317dc4ac61f74bb9c656a7bf10bd96f7 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeUnion.h
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/GeoModelXml/shape/MakeUnion.h
@@ -14,8 +14,9 @@
 
 class MakeUnion: public Element2GeoItem {
 public:
-    MakeUnion();
+    MakeUnion() = default;
     virtual RCBase * make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const override;
+
 };
 
 #endif // GEO_MODEL_XML_MAKE_UNION_H
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeBox.cxx b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeBox.cxx
index bd8cf3e8bf4365e7a718fc6cb1ef08f37eff0fb0..854c5cd445060ff5a1d4f06e066d0b638806f8c1 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeBox.cxx
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeBox.cxx
@@ -10,21 +10,20 @@
 #include "xercesc/util/XMLString.hpp"
 #include "GeoModelXml/GmxUtil.h"
 
+#include <array>
 using namespace xercesc;
 
-MakeBox::MakeBox() {}
-
 RCBase * MakeBox::make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const {
-const int nParams = 3; 
-char const *parName[nParams] = {"xhalflength", "yhalflength", "zhalflength"};
-double p[nParams];
-char *toRelease;
+  constexpr int nParams = 3; 
+  static const std::array<std::string, nParams> parName{"xhalflength", "yhalflength", "zhalflength"};
+  std::array<double, nParams> p{};
+  char *toRelease;
 
-    for (int i = 0; i < nParams; ++i) {
-        toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[i])));
-        p[i] = gmxUtil.evaluate(toRelease);
-        XMLString::release(&toRelease);
-    }
+  for (int i = 0; i < nParams; ++i) {
+      toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[i].data())));
+      p[i] = gmxUtil.evaluate(toRelease);
+      XMLString::release(&toRelease);
+  }
 
-    return new GeoBox(p[0], p[1], p[2]);
+  return const_cast<GeoShape*>(cacheShape(new GeoBox(p[0], p[1], p[2])).get());
 }
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeCons.cxx b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeCons.cxx
index 092261f6db5081f8ef84940ff4b41c87e79e18bd..799b735d43fe8b05874e92f6b073a0c49bd4c823 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeCons.cxx
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeCons.cxx
@@ -10,21 +10,22 @@
 #include "xercesc/util/XMLString.hpp"
 #include "GeoModelXml/GmxUtil.h"
 
+#include <array>
+
 using namespace xercesc;
 
-MakeCons::MakeCons() {}
 
 RCBase * MakeCons::make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const {
-const int nParams = 7; 
-char const *parName[nParams] = {"rmin1", "rmin2", "rmax1", "rmax2", "dz", "sphi", "dphi"};
-double p[nParams];
-char *toRelease;
+    constexpr int nParams = 7; 
+    static const std::array<std::string, nParams> parName{"rmin1", "rmin2", "rmax1", "rmax2", "dz", "sphi", "dphi"};
+    std::array<double, nParams> p{};
+    char *toRelease;
 
     for (int i = 0; i < nParams; ++i) {
-        toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[i])));
+        toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[i].data())));
         p[i] = gmxUtil.evaluate(toRelease);
         XMLString::release(&toRelease);
     }
 
-    return new GeoCons(p[0], p[1], p[2], p[3], p[4], p[5], p[6]);
+    return const_cast<GeoShape*>(cacheShape(new GeoCons(p[0], p[1], p[2], p[3], p[4], p[5], p[6])).get());
 }
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeEllipticalTube.cxx b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeEllipticalTube.cxx
index 6b874b8792b9ebb4d7a02892dc645b4d1ed6a450..178f4b022df8671bf53722a150be1545e63954fe 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeEllipticalTube.cxx
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeEllipticalTube.cxx
@@ -8,22 +8,22 @@
 #include "GeoModelKernel/GeoEllipticalTube.h"
 #include "xercesc/util/XMLString.hpp"
 #include "GeoModelXml/GmxUtil.h"
+#include <array>
 
 using namespace xercesc;
 
-MakeEllipticalTube::MakeEllipticalTube() {}
 
 RCBase * MakeEllipticalTube::make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const {
-const int nParams = 3; 
-char const *parName[nParams] = {"xhalflength", "yhalflength", "zhalflength"};
-double p[nParams];
-char *toRelease;
+    constexpr int nParams = 3; 
+    static const std::array<std::string, nParams> parName{"xhalflength", "yhalflength", "zhalflength"};
+    std::array<double, nParams> p{};
+    char *toRelease;
 
     for (int i = 0; i < nParams; ++i) {
-        toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[i])));
+        toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[i].data())));
         p[i] = gmxUtil.evaluate(toRelease);
         XMLString::release(&toRelease);
     }
 
-    return new GeoEllipticalTube(p[0], p[1], p[2]);
+    return  const_cast<GeoShape*>(cacheShape(new GeoEllipticalTube(p[0], p[1], p[2])).get());
 }
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeGenericTrap.cxx b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeGenericTrap.cxx
index 4bf1df752462fb55565c92f4c2dfaf704d7b609e..5f9e1723c5692ebe05ac568093b6bb1f9866a5c0 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeGenericTrap.cxx
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeGenericTrap.cxx
@@ -16,33 +16,30 @@
 
 using namespace xercesc;
 
-MakeGenericTrap::MakeGenericTrap() {}
 
 RCBase * MakeGenericTrap::make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const {
-const int nParams = 17; 
-char const *parName[nParams] = {"x0", "y0", "x1", "y1", "x2", "y2", "x3", "y3", 
-                                "x4", "y4", "x5", "y5", "x6", "y6", "x7", "y7", "zhalflength"};
-double p;
-char *toRelease;
-std::vector<GeoTwoVector> vertices;
-double x;
-double y;
+    constexpr int nParams = 17; 
+    std::array<std::string, 17> parName{"x0", "y0", "x1", "y1", "x2", "y2", "x3", "y3", 
+                                        "x4", "y4", "x5", "y5", "x6", "y6", "x7", "y7", "zhalflength"};
+    char *toRelease;
+    std::vector<GeoTwoVector> vertices;
+    double x{0.}, y{0.}, p{0.};
 
     for (int i = 0; i < (nParams - 1) / 2; ++i) {
-        toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[2 * i])));
+        toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[2 * i].data())));
         x = gmxUtil.evaluate(toRelease);
         XMLString::release(&toRelease);
-        toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[2 * i + 1])));
+        toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[2 * i + 1].data())));
         y = gmxUtil.evaluate(toRelease);
         XMLString::release(&toRelease);
-        vertices.push_back(GeoTwoVector(x, y));
+        vertices.emplace_back(x, y);
     }
-//
-//    z-half-length
-//
-    toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[16])));
+    //
+    //    z-half-length
+    //
+    toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[16].data())));
     p = gmxUtil.evaluate(toRelease);
     XMLString::release(&toRelease);
 
-    return new GeoGenericTrap(p, vertices);
+    return const_cast<GeoShape*>(cacheShape(new GeoGenericTrap(p, vertices)).get());
 }
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeIntersection.cxx b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeIntersection.cxx
index e787255e4c525d0ca1db2f2608fdcdcf5f95c3e0..8348b666f78498f9a608efb90bfe5e67a106a3c4 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeIntersection.cxx
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeIntersection.cxx
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
 */
 
 // Automatically generated code from /home/hessey/prog/gmx2geo/makeshape
@@ -9,60 +9,64 @@
 #include "OutputDirector.h"
 #include <xercesc/dom/DOM.hpp>
 #include "GeoModelKernel/GeoDefinitions.h"
-#include "GeoModelKernel/RCBase.h"
-#include "GeoModelKernel/GeoShape.h"
-#include "GeoModelKernel/GeoTransform.h"
+#include "GeoModelKernel/GeoShapeIntersection.h"
+#include "GeoModelKernel/GeoShapeShift.h"
 
+#include "GeoModelHelpers/TransformSorter.h"
+#include "GeoModelHelpers/GeoShapeUtils.h"
+#include "GeoModelHelpers/throwExcept.h"
 #include "xercesc/util/XMLString.hpp"
 #include "GeoModelXml/GmxUtil.h"
 
 using namespace xercesc;
 using namespace std;
 
-MakeIntersection::MakeIntersection() {}
-
 RCBase * MakeIntersection::make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const {
 // 
 //    Process child elements; first is first shaperef; then transform; then second shaperef.
 //
-    const GeoShape *first = 0;
-    const GeoShape *second = 0;
+    GeoIntrusivePtr<const GeoShape> first{}, second{};
     GeoTrf::Transform3D hepXf=GeoTrf::Transform3D::Identity(); 
     int elementIndex = 0;
     for (DOMNode *child = element->getFirstChild(); child != 0; child = child->getNextSibling()) {
         if (child->getNodeType() == DOMNode::ELEMENT_NODE) { // Skips text nodes
-	  switch (elementIndex) {
-	  case 0: { // First element is first shaperef
-	    first = static_cast<const GeoShape *>(gmxUtil.tagHandler.shaperef.process(dynamic_cast<DOMElement *> (child), gmxUtil));
-	    break;
-	  }
-	  case 1:  { // Second element is transformation or transformationref
-	    char *toRelease = XMLString::transcode(child->getNodeName());
-	    string nodeName(toRelease);
-	    XMLString::release(&toRelease);
-	    const GeoTransform *geoXf = (nodeName == "transformation")
-	      ? static_cast<const GeoTransform *>( gmxUtil.tagHandler.transformation.process(dynamic_cast<DOMElement *>(child), gmxUtil))
-	      : static_cast<const GeoTransform *>( gmxUtil.tagHandler.transformationref.process(dynamic_cast<DOMElement *>(child), gmxUtil));
-	    hepXf = geoXf->getTransform();
-	    break;
-	  }
-	  case 2: { // Third element is second shaperef
-	    second = static_cast<const GeoShape *>( gmxUtil.tagHandler.shaperef.process(dynamic_cast<DOMElement *> (child), gmxUtil));
-	    break;
-	  }
-	  default: // More than 3 elements?
-	    msglog << MSG::FATAL << "MakeIntersection: Incompatible DTD? got more than 3 child elements\n";
-	  }
-	  elementIndex++;
+      switch (elementIndex) {
+      case 0: { // First element is first shaperef
+        first.reset(static_cast<const GeoShape *>(gmxUtil.tagHandler.shaperef.process(dynamic_cast<DOMElement *> (child), gmxUtil)));
+        break;
+      }
+      case 1:  { // Second element is transformation or transformationref
+        char *toRelease = XMLString::transcode(child->getNodeName());
+        string nodeName(toRelease);
+        XMLString::release(&toRelease);
+        const GeoTransform *geoXf = (nodeName == "transformation")
+          ? static_cast<const GeoTransform *>( gmxUtil.tagHandler.transformation.process(dynamic_cast<DOMElement *>(child), gmxUtil))
+          : static_cast<const GeoTransform *>( gmxUtil.tagHandler.transformationref.process(dynamic_cast<DOMElement *>(child), gmxUtil));
+        hepXf = geoXf->getTransform();
+        break;
+      }
+      case 2: { // Third element is second shaperef
+        second.reset(static_cast<const GeoShape *>(gmxUtil.tagHandler.shaperef.process(dynamic_cast<DOMElement *> (child), gmxUtil)));
+        break;
+      }
+      default: // More than 3 elements?
+        THROW_EXCEPTION("MakeIntersection: Incompatible DTD? got more than 3 child elements\n");
+      }
+      elementIndex++;
         }
     }
 
-    if (!first || !second) std::abort();
-
-    // FIXME: intersect() returns a new'd object --- should really be
-    // returning a `unique_ptr<GeoShapeIntersection>' not a
-    // `const GeoShapeIntersection'
-    GeoShapeIntersection *temp = const_cast<GeoShapeIntersection*>(&(first->intersect(*(GeoShape *) &(*(second) << hepXf))));
+    if (!first || !second) {
+        THROW_EXCEPTION("Not both shapes were defined");
+    } 
+    GeoIntrusivePtr<GeoShapeIntersection> isect{};
+    static const GeoTrf::TransformSorter trfSorter{};
+    if (trfSorter.compare(hepXf, GeoTrf::Transform3D::Identity())) {
+        isect = new GeoShapeIntersection(first, new GeoShapeShift(second, hepXf));
+    } else {
+        isect = new GeoShapeIntersection(first, second);
 
-    return (RCBase *) temp;
+    }
+    std::cout<<isect->refCount()<<std::endl;
+    return const_cast<GeoShape*>(cacheShape(isect).get());
 }
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakePara.cxx b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakePara.cxx
index c68e187611e95746c93a369f46a3149b15703a1d..edaf9c774ff0f2b562465f35112afeb8b0576cbd 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakePara.cxx
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakePara.cxx
@@ -10,21 +10,21 @@
 #include "xercesc/util/XMLString.hpp"
 #include "GeoModelXml/GmxUtil.h"
 
+#include <array>
 using namespace xercesc;
 
-MakePara::MakePara() {}
 
 RCBase * MakePara::make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const {
-const int nParams = 6; 
-char const *parName[nParams] = {"xhalflength", "yhalflength", "zhalflength", "alpha", "theta", "phi"};
-double p[nParams];
-char *toRelease;
+    constexpr int nParams = 6; 
+    static const std::array<std::string, nParams> parName{"xhalflength", "yhalflength", "zhalflength", "alpha", "theta", "phi"};
+    std::array<double, nParams> p{};
+    char *toRelease;
 
     for (int i = 0; i < nParams; ++i) {
-        toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[i])));
+        toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[i].data())));
         p[i] = gmxUtil.evaluate(toRelease);
         XMLString::release(&toRelease);
     }
 
-    return new GeoPara(p[0], p[1], p[2], p[3], p[4], p[5]);
+    return const_cast<GeoShape*>(cacheShape(new GeoPara(p[0], p[1], p[2], p[3], p[4], p[5])).get());
 }
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakePcon.cxx b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakePcon.cxx
index 4b98472a34cc87373c5c39110dacd0b6fdb78e25..ec658425a18fd66d84c3525c5add4258143022af 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakePcon.cxx
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakePcon.cxx
@@ -1,39 +1,37 @@
 /*
-  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
 */
 
 // Automatically generated code from /home/hessey/prog/gmx2geo/makeshape
 // But then edited for AddPlane stuff
 #include "GeoModelXml/shape/MakePcon.h"
 #include <xercesc/dom/DOM.hpp>
-#include "GeoModelKernel/RCBase.h"
 #include "GeoModelKernel/GeoPcon.h"
+#include "GeoModelHelpers/GeoShapeUtils.h"
+#include "GeoModelHelpers/throwExcept.h"
 #include "xercesc/util/XMLString.hpp"
 #include "GeoModelXml/GmxUtil.h"
 
+#include <array>
 using namespace xercesc;
 
-MakePcon::MakePcon() {}
 
 RCBase * MakePcon::make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const {
-const int nParams = 2; 
-char const *parName[nParams] = {"sphi", "dphi"};
-double p[nParams];
-char *toRelease;
-
+    constexpr int nParams = 2; 
+    static const std::array<std::string, nParams> parName{"sphi", "dphi"};
+    std::array<double, nParams> p{};
+    char *toRelease;
     for (int i = 0; i < nParams; ++i) {
-        toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[i])));
+        toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[i].data())));
         p[i] = gmxUtil.evaluate(toRelease);
         XMLString::release(&toRelease);
     }
 
-    GeoPcon *pcon = new GeoPcon(p[0], p[1]);
-//
-//    Add planes
-//
-    double zPlane = 0.;
-    double rMinPlane = 0.;
-    double rMaxPlane = 0.;
+    GeoIntrusivePtr<GeoPcon> pcon{new GeoPcon(p[0], p[1])};
+    //
+    //    Add planes
+    //
+    double zPlane{0.}, rMinPlane{0.}, rMaxPlane{0.};
     for (DOMNode *child = element->getFirstChild(); child != 0; child = child->getNextSibling()) {
         if (child->getNodeType() == DOMNode::ELEMENT_NODE) {
             toRelease = XMLString::transcode(child->getNodeName());
@@ -45,6 +43,8 @@ char *toRelease;
             }
         }
     }
-
-    return pcon;
+    if (!pcon->isValid()) {
+        THROW_EXCEPTION("Invalid Pcon defined "<<printGeoShape(pcon));
+    }
+    return const_cast<GeoShape*>(cacheShape(pcon).get());
 }
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakePgon.cxx b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakePgon.cxx
index a73b002c76279d6b03cca1bc818d26c9aa96dd6b..36d277c8726681b9dd70c0e1186cca17747a90e6 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakePgon.cxx
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakePgon.cxx
@@ -5,35 +5,32 @@
 // Automatically generated code from /home/hessey/prog/gmx2geo/makeshape
 // But then edited for AddPlane stuff
 #include "GeoModelXml/shape/MakePgon.h"
-#include <xercesc/dom/DOM.hpp>
-#include "GeoModelKernel/RCBase.h"
 #include "GeoModelKernel/GeoPgon.h"
+#include <xercesc/dom/DOM.hpp>
 #include "xercesc/util/XMLString.hpp"
 #include "GeoModelXml/GmxUtil.h"
 
+#include <array>
 using namespace xercesc;
 
-MakePgon::MakePgon() {}
 
 RCBase * MakePgon::make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const {
-const int nParams = 3; 
-char const *parName[nParams] = {"sphi", "dphi", "nsides"};
-double p[nParams];
-char *toRelease;
+    constexpr int nParams = 3; 
+    static const std::array<std::string, nParams> parName {"sphi", "dphi", "nsides"};
+    std::array<double, nParams> p{};
+    char *toRelease;
 
     for (int i = 0; i < nParams; ++i) {
-        toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[i])));
+        toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[i].data())));
         p[i] = gmxUtil.evaluate(toRelease);
         XMLString::release(&toRelease);
     }
 
-    GeoPgon *pgon = new GeoPgon(p[0], p[1], p[2]);
+    GeoIntrusivePtr<GeoPgon> pgon{new GeoPgon(p[0], p[1], p[2])};
 //
 //    Add planes
 //
-    double zPlane = 0.;
-    double rMinPlane = 0.;
-    double rMaxPlane = 0.;
+    double zPlane{0.}, rMinPlane{0.}, rMaxPlane{0.};
     for (DOMNode *child = element->getFirstChild(); child != 0; child = child->getNextSibling()) {
         if (child->getNodeType() == DOMNode::ELEMENT_NODE) {
             toRelease = XMLString::transcode(child->getNodeName());
@@ -45,6 +42,5 @@ char *toRelease;
             }
         }
     }
-
-    return pgon;
+    return const_cast<GeoShape*>(cacheShape(pgon).get());
 }
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeRotation.cxx b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeRotation.cxx
index db97f3f0159be6d74f89d148e3e23be5ec59c265..5473cbfafb7a70dc23f0ba2c4c907f6f3139ade8 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeRotation.cxx
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeRotation.cxx
@@ -18,12 +18,12 @@ using namespace xercesc;
 GeoTrf::Rotation3D MakeRotation::getTransform(const DOMElement *rotation, GmxUtil &gmxUtil) {
 
 const int nParams = 4; 
-char const *parName[nParams] = {"angle", "xcos", "ycos", "zcos"};
-double p[nParams];
+static const std::array<std::string, nParams> parName {"angle", "xcos", "ycos", "zcos"};
+std::array<double, nParams> p{};
 char *toRelease;
 
     for (int i = 0; i < nParams; ++i) {
-        toRelease = XMLString::transcode(rotation->getAttribute(XMLString::transcode(parName[i])));
+        toRelease = XMLString::transcode(rotation->getAttribute(XMLString::transcode(parName[i].data())));
         p[i] = gmxUtil.evaluate(toRelease);
         XMLString::release(&toRelease);
     }
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeSimplePolygonBrep.cxx b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeSimplePolygonBrep.cxx
index 80943dc738b323a85d6e468f84ba189a2be1cffd..4c6f22bc7acef03353e01ac7c0caf1aaaaf1a697 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeSimplePolygonBrep.cxx
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeSimplePolygonBrep.cxx
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
 */
 #include <string>
 #include <vector>
@@ -7,37 +7,35 @@
 #include <stdexcept>
 
 #include "GeoModelXml/shape/MakeSimplePolygonBrep.h"
-#include <xercesc/dom/DOM.hpp>
-#include "GeoModelKernel/RCBase.h"
 #include "GeoModelKernel/GeoSimplePolygonBrep.h"
+#include <xercesc/dom/DOM.hpp>
+#include "GeoModelHelpers/throwExcept.h"
 #include "xercesc/util/XMLString.hpp"
 #include "GeoModelXml/GmxUtil.h"
 using namespace xercesc;
 
-MakeSimplePolygonBrep::MakeSimplePolygonBrep() {}
 
 RCBase * MakeSimplePolygonBrep::make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const {
-const int nParams = 3; 
-char const *parName[nParams] = {"xpoints", "ypoints", "zhalflength"};
-double z;
-std::vector <double> x;
-std::vector <double> y;
+    constexpr int nParams = 3; 
+    static const std::array<std::string, nParams> parName {"xpoints", "ypoints", "zhalflength"};
+    double z{0.};
+    std::vector <double> x{}, y{};
 
-char *toRelease;
+    char *toRelease;
 
-    toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[0])));
+    toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[0].data())));
     std::string xPoints(toRelease);
     XMLString::release(&toRelease);
 
-    toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[1])));
+    toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[1].data())));
     std::string yPoints(toRelease);
     XMLString::release(&toRelease);
 
-    toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[2])));
+    toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[2].data())));
     z = gmxUtil.evaluate(toRelease);
     XMLString::release(&toRelease);
 
-    GeoSimplePolygonBrep * poly = new GeoSimplePolygonBrep(z);
+    GeoIntrusivePtr<GeoSimplePolygonBrep> poly{new GeoSimplePolygonBrep(z)};
 
     std::istringstream xSS(xPoints);
     while (!xSS.eof()) {
@@ -57,7 +55,7 @@ char *toRelease;
     int ny = y.size();
   
     if (nx < 3 || ny < 3 || nx != ny) {
-        throw std::runtime_error(std::string("MakeSimplePolygonBrep: Unequal number of x and y points, or less than 3\n\n") +
+        THROW_EXCEPTION("Unequal number of x and y points, or less than 3\n\n"<<
          "xpoints was:\n" + xPoints + "\nypoints was:\n" + yPoints + "\n\n");
     }
 
@@ -65,6 +63,6 @@ char *toRelease;
         poly->addVertex(x[i], y[i]);
     }
 
-    return poly;
+    return const_cast<GeoShape*>(cacheShape(poly).get());
 
 }
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeSubtraction.cxx b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeSubtraction.cxx
index 1e16c58973be49b4eea91dc742333d5a5c313608..4ee9be38ad74ae7f796a6cbb1e48715feae403d7 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeSubtraction.cxx
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeSubtraction.cxx
@@ -8,62 +8,64 @@
 #include "OutputDirector.h"
 #include <string>
 #include <xercesc/dom/DOM.hpp>
-//   #include <CLHEP/Geometry/Transform3D.h>
-#include "GeoModelKernel/GeoDefinitions.h"
-#include "GeoModelKernel/RCBase.h"
-#include "GeoModelKernel/GeoShape.h"
+#include "GeoModelKernel/GeoShapeSubtraction.h"
+#include "GeoModelKernel/GeoShapeShift.h"
 #include "GeoModelKernel/GeoTransform.h"
+#include "GeoModelHelpers/throwExcept.h"
 
 #include "xercesc/util/XMLString.hpp"
 #include "GeoModelXml/GmxUtil.h"
-
+#include <array>
 using namespace xercesc;
 using namespace std;
 
-MakeSubtraction::MakeSubtraction() {}
 
 RCBase * MakeSubtraction::make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const {
 // 
 //    Process child elements; first is first shaperef; then transformation; then second shaperef.
 //
-    const GeoShape *first = 0;
-    const GeoShape *second = 0;
+    GeoIntrusivePtr<const GeoShape> first{}, second{};
     GeoTrf::Transform3D hepXf=GeoTrf::Transform3D::Identity();
     int elementIndex = 0;
     for (DOMNode *child = element->getFirstChild(); child != 0; child = child->getNextSibling()) {
         if (child->getNodeType() == DOMNode::ELEMENT_NODE) { // Skips text nodes
-	  switch (elementIndex) {
-	  case 0: { // First element is first shaperef
-	    first = static_cast<const GeoShape *>(gmxUtil.tagHandler.shaperef.process(dynamic_cast<DOMElement *> (child), gmxUtil));
-	    break;
-	  }
-	  case 1: { // Second element is transformation or transformationref
-	    char *toRelease = XMLString::transcode(child->getNodeName());
-	    string nodeName(toRelease);
-	    XMLString::release(&toRelease);
-	    const GeoTransform *geoXf = (nodeName == "transformation")
-	      ? static_cast<const GeoTransform *>( gmxUtil.tagHandler.transformation.process(dynamic_cast<DOMElement *>(child), gmxUtil))
-	      : static_cast<const GeoTransform *>( gmxUtil.tagHandler.transformationref.process(dynamic_cast<DOMElement *>(child), gmxUtil));
-	    hepXf = geoXf->getTransform();
-	    break;
-	  }
-	  case 2: { // Third element is second shaperef
-	    second = static_cast<const GeoShape *>( gmxUtil.tagHandler.shaperef.process(dynamic_cast<DOMElement *> (child), gmxUtil));
-	    break;
-	  }
-	  default: // More than 3 elements?
-	    msglog << MSG::FATAL  << "MakeSubtraction: Incompatible DTD? got more than 3 child elements" << endmsg;
-	  }
-	  elementIndex++;
+      switch (elementIndex) {
+      case 0: { // First element is first shaperef
+        first = static_cast<const GeoShape *>(gmxUtil.tagHandler.shaperef.process(dynamic_cast<DOMElement *> (child), gmxUtil));
+        break;
+      }
+      case 1: { // Second element is transformation or transformationref
+        char *toRelease = XMLString::transcode(child->getNodeName());
+        string nodeName(toRelease);
+        XMLString::release(&toRelease);
+        const GeoTransform *geoXf = (nodeName == "transformation")
+          ? static_cast<const GeoTransform *>( gmxUtil.tagHandler.transformation.process(dynamic_cast<DOMElement *>(child), gmxUtil))
+          : static_cast<const GeoTransform *>( gmxUtil.tagHandler.transformationref.process(dynamic_cast<DOMElement *>(child), gmxUtil));
+        hepXf = geoXf->getTransform();
+        break;
+      }
+      case 2: { // Third element is second shaperef
+        second = static_cast<const GeoShape *>( gmxUtil.tagHandler.shaperef.process(dynamic_cast<DOMElement *> (child), gmxUtil));
+        break;
+      }
+      default: // More than 3 elements?
+        THROW_EXCEPTION("MakeSubtraction: Incompatible DTD? got more than 3 child elements");
+      }
+      elementIndex++;
         }
     }
 
-    if (!first || !second) std::abort();
-
-    // FIXME: subtract() returns a new'd object --- should really be
-    // returning a `unique_ptr<GeoShapeSubtraction>' not a
-    // `const GeoShapeSubtraction'
-    GeoShapeSubtraction *temp = const_cast<GeoShapeSubtraction*>(&(first->subtract(*(GeoShape *) &(*(second) << hepXf))));
+    if (!first || !second) {
+        THROW_EXCEPTION("Only one of the two subtraction operands was defined.");
+    }
 
-    return (RCBase *) temp;
+  
+    GeoIntrusivePtr<GeoShapeSubtraction> subtract{};
+    static const GeoTrf::TransformSorter sorter{};
+    if (sorter.compare(hepXf, GeoTrf::Transform3D::Identity())) {
+        subtract = new GeoShapeSubtraction(first, new GeoShapeShift(second, hepXf));
+    } else {
+        subtract = new GeoShapeSubtraction(first, second);
+    }
+    return const_cast<GeoShape*>(cacheShape(subtract).get());
 }
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTorus.cxx b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTorus.cxx
index 205a260e588f518dff02b0c142c7be2535405727..99bfa5dd63fbc6355f0b56f983277da8ee9d91a3 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTorus.cxx
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTorus.cxx
@@ -13,19 +13,18 @@
 
 using namespace xercesc;
 
-MakeTorus::MakeTorus() {}
 
 RCBase * MakeTorus::make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const {
-const int nParams = 5; 
-char const *parName[nParams] = {"rmin", "rmax", "rtor", "sphi", "dphi"};
-double p[nParams];
-char *toRelease;
+  constexpr int nParams = 5; 
+  static const std::array<std::string, nParams> parName {"rmin", "rmax", "rtor", "sphi", "dphi"};
+  std::array<double, nParams> p{};
+  char *toRelease;
 
-    for (int i = 0; i < nParams; ++i) {
-        toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[i])));
-        p[i] = gmxUtil.evaluate(toRelease);
-        XMLString::release(&toRelease);
-    }
+  for (int i = 0; i < nParams; ++i) {
+      toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[i].data())));
+      p[i] = gmxUtil.evaluate(toRelease);
+      XMLString::release(&toRelease);
+  }
 
-    return new GeoTorus(p[0], p[1], p[2], p[3], p[4]);
+  return  const_cast<GeoShape*>(cacheShape(new GeoTorus(p[0], p[1], p[2], p[3], p[4])).get());
 }
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTransformation.cxx b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTransformation.cxx
index 5c1e1cf30bbdf0f04e830de94579272389344f1c..200c1d7be771e73593ee6d680605b203dff558bb 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTransformation.cxx
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTransformation.cxx
@@ -55,7 +55,7 @@ char *name2release;
     if (alignable.compare(string("true")) == 0) {
         return (RCBase *) new GeoAlignableTransform(hepTransform);
     }
-    else {
-        return (RCBase *) new GeoTransform(hepTransform);
-    }
+    
+    return makeTransform(hepTransform);
+    
 }
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTranslation.cxx b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTranslation.cxx
index 7068dd222ffb8f1cfa57fe1c185eb8b26ba16235..3d73786d4b44f4f2cc499ea6666a9d072ca089ee 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTranslation.cxx
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTranslation.cxx
@@ -19,12 +19,12 @@ using namespace xercesc;
 GeoTrf::Translate3D MakeTranslation::getTransform(const DOMElement *translation, GmxUtil &gmxUtil) {
 
 const int nParams = 3; 
-char const *parName[nParams] = {"x", "y", "z"};
-double p[nParams];
+static const std::array<std::string, nParams> parName {"x", "y", "z"};
+std::array<double, nParams> p{};
 char *toRelease;
 
     for (int i = 0; i < nParams; ++i) {
-        toRelease = XMLString::transcode(translation->getAttribute(XMLString::transcode(parName[i])));
+        toRelease = XMLString::transcode(translation->getAttribute(XMLString::transcode(parName[i].data())));
         p[i] = gmxUtil.evaluate(toRelease);
         XMLString::release(&toRelease);
     }
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTrap.cxx b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTrap.cxx
index 6e881629f2def603fbbbf18218211ba2de957d9d..577f729d5cd5c86795dc117d164e0d2c55a5e4c3 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTrap.cxx
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTrap.cxx
@@ -1,30 +1,30 @@
 /*
-  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
 */
 
 // Automatically generated code from /home/hessey/prog/gmx2geo/makeshape
 #include "GeoModelXml/shape/MakeTrap.h"
 #include <xercesc/dom/DOM.hpp>
 #include "GeoModelKernel/RCBase.h"
-#include "GeoModelKernel/GeoTrap.h"
 #include "xercesc/util/XMLString.hpp"
 #include "GeoModelXml/GmxUtil.h"
-
+#include "GeoModelKernel/GeoTrap.h"
 using namespace xercesc;
 
-MakeTrap::MakeTrap() {}
 
 RCBase * MakeTrap::make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const {
-const int nParams = 11; 
-char const *parName[nParams] = {"zhalflength", "theta", "phi", "dydzn", "dxdyndzn", "dxdypdzn", "angleydzn", "dydzp", "dxdyndzp", "dxdypdzp", "angleydzp"};
-double p[nParams];
-char *toRelease;
+  constexpr int nParams = 11; 
+  static const std::array<std::string, nParams> parName {"zhalflength", "theta", "phi", "dydzn", "dxdyndzn", "dxdypdzn", "angleydzn", "dydzp", "dxdyndzp", "dxdypdzp", "angleydzp"};
+  std::array<double, nParams> p{};
+  char *toRelease;
 
     for (int i = 0; i < nParams; ++i) {
-        toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[i])));
+        toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[i].data())));
         p[i] = gmxUtil.evaluate(toRelease);
         XMLString::release(&toRelease);
     }
 
-    return new GeoTrap(p[0], p[1], p[2], p[3], p[4], p[5], p[6], p[7], p[8], p[9], p[10]);
+    return  const_cast<GeoShape*>(cacheShape(new GeoTrap(p[0], p[1], p[2], p[3], 
+                                                         p[4], p[5], p[6], p[7], 
+                                                         p[8], p[9], p[10])).get());
 }
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTrd.cxx b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTrd.cxx
index 7b326216dfe14327cdec8b4e54034ec65247be5a..5d8d64ff155a41e1f713bf54420acdb6c4cc32cb 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTrd.cxx
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTrd.cxx
@@ -12,19 +12,18 @@
 
 using namespace xercesc;
 
-MakeTrd::MakeTrd() {}
 
 RCBase * MakeTrd::make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const {
-const int nParams = 5; 
-char const *parName[nParams] = {"xhalflength1", "xhalflength2", "yhalflength1", "yhalflength2", "zhalflength"};
-double p[nParams];
-char *toRelease;
+    constexpr int nParams = 5; 
+    static const std::array<std::string, nParams> parName {"xhalflength1", "xhalflength2", "yhalflength1", "yhalflength2", "zhalflength"};
+    std::array<double, nParams> p{};
+    char *toRelease;
 
     for (int i = 0; i < nParams; ++i) {
-        toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[i])));
+        toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[i].data())));
         p[i] = gmxUtil.evaluate(toRelease);
         XMLString::release(&toRelease);
     }
 
-    return new GeoTrd(p[0], p[1], p[2], p[3], p[4]);
+    return  const_cast<GeoShape*>(cacheShape(new GeoTrd(p[0], p[1], p[2], p[3], p[4])).get());
 }
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTube.cxx b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTube.cxx
index e916cd98af4bccaf3425320fc48e1942a19e6c3c..00806d89af43b0b401cb220fedf377ec04481ad0 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTube.cxx
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTube.cxx
@@ -12,19 +12,18 @@
 
 using namespace xercesc;
 
-MakeTube::MakeTube() {}
 
 RCBase * MakeTube::make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const {
-const int nParams = 3; 
-char const *parName[nParams] = {"rmin", "rmax", "zhalflength"};
-double p[nParams];
-char *toRelease;
+    constexpr int nParams = 3; 
+    static const std::array<std::string, nParams> parName {"rmin", "rmax", "zhalflength"};
+    std::array<double, nParams> p{};
+    char *toRelease;
 
     for (int i = 0; i < nParams; ++i) {
-        toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[i])));
+        toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[i].data())));
         p[i] = gmxUtil.evaluate(toRelease);
         XMLString::release(&toRelease);
     }
 
-    return new GeoTube(p[0], p[1], p[2]);
+    return const_cast<GeoShape*>(cacheShape(new GeoTube(p[0], p[1], p[2])).get());;
 }
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTubs.cxx b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTubs.cxx
index f6c76a6a988ea2b329e84c3cc15db1dc1047a067..c8fe9cad93ae7a211fa9075e3e569d081927e9bc 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTubs.cxx
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTubs.cxx
@@ -12,19 +12,18 @@
 
 using namespace xercesc;
 
-MakeTubs::MakeTubs() {}
 
 RCBase * MakeTubs::make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const {
-const int nParams = 5; 
-char const *parName[nParams] = {"rmin", "rmax", "zhalflength", "sphi", "dphi"};
-double p[nParams];
-char *toRelease;
+    constexpr int nParams = 5; 
+    static const std::array<std::string, nParams> parName {"rmin", "rmax", "zhalflength", "sphi", "dphi"};
+    std::array<double, nParams> p{};
+    char *toRelease;
 
     for (int i = 0; i < nParams; ++i) {
-        toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[i])));
+        toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[i].data())));
         p[i] = gmxUtil.evaluate(toRelease);
         XMLString::release(&toRelease);
     }
 
-    return new GeoTubs(p[0], p[1], p[2], p[3], p[4]);
+    return  const_cast<GeoShape*>(cacheShape(new GeoTubs(p[0], p[1], p[2], p[3], p[4])).get());
 }
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTwistedTrap.cxx b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTwistedTrap.cxx
index 786e506a8f865383c39f434860b367a422dbbdee..ef1b479762bef17c40a48f2d2acf9978e5e50f00 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTwistedTrap.cxx
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeTwistedTrap.cxx
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
 */
 
 /*
@@ -31,19 +31,20 @@
 
 using namespace xercesc;
 
-MakeTwistedTrap::MakeTwistedTrap() {}
 
 RCBase * MakeTwistedTrap::make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const {
-const int nParams = 11; 
-char const *parName[nParams] = {"twist", "dz", "theta","phi","dy1","dx1","dx2","dy2","dx3","dx4","alpha"};
-double p[nParams];
-char *toRelease;
+  constexpr int nParams = 11; 
+  static const std::array<std::string, nParams> parName {"twist", "dz", "theta","phi","dy1","dx1","dx2","dy2","dx3","dx4","alpha"};
+  std::array<double, nParams> p{};
+  char *toRelease;
 
-    for (int i = 0; i < nParams; ++i) {
-        toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[i])));
+  for (int i = 0; i < nParams; ++i) {
+        toRelease = XMLString::transcode(element->getAttribute(XMLString::transcode(parName[i].data())));
         p[i] = gmxUtil.evaluate(toRelease);
         XMLString::release(&toRelease);
     }
 
-    return new GeoTwistedTrap(p[0], p[1], p[2],p[3], p[4], p[5],p[6], p[7], p[8],p[9], p[10]);
+    return const_cast<GeoShape*>(cacheShape(new GeoTwistedTrap(p[0], p[1], p[2], p[3], 
+                                                               p[4], p[5], p[6], p[7], 
+                                                               p[8], p[9], p[10])).get());
 }
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeUnion.cxx b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeUnion.cxx
index 2bcb0d530e2198f4fd379bcb78152c9e12a80e7f..fc949f0b2880c2a5a7cc5970a961f9cf8cdc37d6 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeUnion.cxx
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/src/MakeUnion.cxx
@@ -1,67 +1,67 @@
 /*
-  Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
 */
 
 #include "GeoModelXml/shape/MakeUnion.h"
 #include <string>
 #include "OutputDirector.h"
 #include <xercesc/dom/DOM.hpp>
-//   #include <CLHEP/Geometry/Transform3D.h>
-#include "GeoModelKernel/RCBase.h"
-#include "GeoModelKernel/GeoShape.h"
-#include "GeoModelKernel/GeoTransform.h"
-#include "GeoModelKernel/GeoDefinitions.h"
-
+#include "GeoModelKernel/GeoShapeUnion.h"
+#include "GeoModelKernel/GeoShapeShift.h"
+#include "GeoModelHelpers/throwExcept.h"
 #include "xercesc/util/XMLString.hpp"
 #include "GeoModelXml/GmxUtil.h"
 
 using namespace xercesc;
 using namespace std;
 
-MakeUnion::MakeUnion() {}
 
 RCBase * MakeUnion::make(const xercesc::DOMElement *element, GmxUtil &gmxUtil) const {
 // 
 //    Process child elements; first is first shaperef; then transform; then second shaperef.
 //
-    const GeoShape *first = 0;
-    const GeoShape *second = 0;
+    GeoIntrusivePtr<const GeoShape> first{}, second{};
     GeoTrf::Transform3D hepXf=GeoTrf::Transform3D::Identity();
     int elementIndex = 0;
     for (DOMNode *child = element->getFirstChild(); child != 0; child = child->getNextSibling()) {
         if (child->getNodeType() == DOMNode::ELEMENT_NODE) { // Skips text nodes
-	  switch (elementIndex) {
-	  case 0: { // First element is first shaperef
-	    first = static_cast<GeoShape *>( gmxUtil.tagHandler.shaperef.process(dynamic_cast<DOMElement *> (child), gmxUtil));
-	    break;
-	  }
-	  case 1: { // Second element is transformation or transformationref
-	    char *toRelease = XMLString::transcode(child->getNodeName());
-	    string nodeName(toRelease);
-	    XMLString::release(&toRelease);
-	    const GeoTransform *geoXf = (nodeName == "transformation")
-	      ? static_cast<const GeoTransform *>( gmxUtil.tagHandler.transformation.process(dynamic_cast<DOMElement *>(child), gmxUtil))
-	      : static_cast<const GeoTransform *>( gmxUtil.tagHandler.transformationref.process(dynamic_cast<DOMElement *>(child), gmxUtil));
-	    hepXf = geoXf->getTransform();
-	    break;
-	  }
-	  case 2: { // Third element is second shaperef
-	    second = static_cast<const GeoShape *>( gmxUtil.tagHandler.shaperef.process(dynamic_cast<DOMElement *> (child), gmxUtil));
-	    break;
-	  }
-	  default: // More than 3 elements?
-	    msglog << MSG::FATAL << "MakeUnion: Incompatible DTD? got more than 3 child elements" << endmsg;
-	  }
-	  elementIndex++;
+      switch (elementIndex) {
+      case 0: { // First element is first shaperef
+        first = static_cast<GeoShape *>( gmxUtil.tagHandler.shaperef.process(dynamic_cast<DOMElement *> (child), gmxUtil));
+        break;
+      }
+      case 1: { // Second element is transformation or transformationref
+        char *toRelease = XMLString::transcode(child->getNodeName());
+        string nodeName(toRelease);
+        XMLString::release(&toRelease);
+        const GeoTransform *geoXf = (nodeName == "transformation")
+          ? static_cast<const GeoTransform *>( gmxUtil.tagHandler.transformation.process(dynamic_cast<DOMElement *>(child), gmxUtil))
+          : static_cast<const GeoTransform *>( gmxUtil.tagHandler.transformationref.process(dynamic_cast<DOMElement *>(child), gmxUtil));
+        hepXf = geoXf->getTransform();
+        break;
+      }
+      case 2: { // Third element is second shaperef
+        second = static_cast<const GeoShape *>( gmxUtil.tagHandler.shaperef.process(dynamic_cast<DOMElement *> (child), gmxUtil));
+        break;
+      }
+      default: // More than 3 elements?
+        THROW_EXCEPTION("MakeUnion: Incompatible DTD? got more than 3 child elements");
+      }
+      elementIndex++;
         }
     }
 
-    if (!first || !second) std::abort();
-
-    // FIXME: add() returns a new'd object --- should really be
-    // returning a `unique_ptr<GeoShapeUnion>' not a
-    // `const GeoShapeUnion'
-    GeoShapeUnion *temp = const_cast<GeoShapeUnion*>(&(first->add(*(GeoShape *) &(*(second) << hepXf))));
+    if (!first || !second) {
+		THROW_EXCEPTION("One of the two operands is not defined");
+    }
+    GeoIntrusivePtr<const GeoShapeUnion> unionShape{};
 
-    return (RCBase *) temp;
+   static const GeoTrf::TransformSorter sorter{};
+    if (sorter.compare(hepXf, GeoTrf::Transform3D::Identity())) {
+        unionShape = new GeoShapeUnion(first, new GeoShapeShift(second,hepXf));
+    } else {
+        unionShape = new GeoShapeUnion(first, second);
+    }
+ 
+    return  const_cast<GeoShape*>(cacheShape(unionShape).get());
 }
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/src/MulticopyProcessor.cxx b/GeoModelTools/GeoModelXML/GeoModelXml/src/MulticopyProcessor.cxx
index ff5ad62d3c6a05e99ad55638552fb033d5fce5f2..c332ead2836195001543544202f5bbc5197069fa 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/src/MulticopyProcessor.cxx
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/src/MulticopyProcessor.cxx
@@ -151,7 +151,7 @@ DOMDocument *doc = element->getOwnerDocument();
             }
             GeoTrf::Transform3D hepXf=GeoTrf::Transform3D::Identity(); // Identity initially
             for (int i = 0; i < nCopies; ++i) {
-                xfList->push_back((GeoGraphNode *) new GeoTransform(hepXf));
+                xfList->push_back(makeTransform(hepXf));
                 hepXf = hepXf0 * hepXf;
             }
         }
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/src/ReplicaRPhiProcessor.cxx b/GeoModelTools/GeoModelXML/GeoModelXml/src/ReplicaRPhiProcessor.cxx
index 8a9596ba84ca48d56761681d9cab1eac336b3839..dc7601ee3419a5591a032b95e3b4be6c6b1e1609 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/src/ReplicaRPhiProcessor.cxx
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/src/ReplicaRPhiProcessor.cxx
@@ -162,14 +162,14 @@ DOMDocument *doc = element->getOwnerDocument();
                 hepXf0 = geoAXf->getTransform();
             }
             else {
-                geoXf = new GeoTransform (hepXf0);
+                geoXf = makeTransform (hepXf0);
                 hepXf0 = geoXf->getTransform();
             }
 	    double angle=offsetPhi;
             GeoTrf::Transform3D hepXf=hepXf0; 
             for (int i = 0; i < nCopies; ++i) {
 	    	hepXf=hepXf0*GeoTrf::TranslateX3D(radius*cos(angle))*GeoTrf::TranslateY3D(radius*sin(angle))*GeoTrf::RotateZ3D(angle);
-                xfList->push_back((GeoGraphNode *) new GeoTransform(hepXf));
+                xfList->push_back(makeTransform(hepXf));
                 // hepXf = hepXf * GeoTrf::RotateZ3D(stepPhi) ;
 		angle+=stepPhi;
             }
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/src/ReplicaXProcessor.cxx b/GeoModelTools/GeoModelXML/GeoModelXml/src/ReplicaXProcessor.cxx
index 56c8d3fdf71ddae96d2b490c0f85dd46079da258..1f2845d9f02fe58fc36121472ad8e87a2b9f03ac 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/src/ReplicaXProcessor.cxx
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/src/ReplicaXProcessor.cxx
@@ -141,12 +141,12 @@ DOMDocument *doc = element->getOwnerDocument();
                 hepXf0 = geoAXf->getTransform();
             }
             else {
-                geoXf = new GeoTransform (hepXf0);
+                geoXf = makeTransform (hepXf0);
                 hepXf0 = geoXf->getTransform();
             }
             GeoTrf::Transform3D hepXf=hepXf0; 
             for (int i = 0; i < nCopies; ++i) {
-                xfList->push_back((GeoGraphNode *) new GeoTransform(hepXf));
+                xfList->push_back(makeTransform(hepXf));
                 hepXf = hepXf * GeoTrf::TranslateX3D(stepX) ;
             }
     }
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/src/ReplicaXYArrays.cxx b/GeoModelTools/GeoModelXML/GeoModelXml/src/ReplicaXYArrays.cxx
index 3b52151ceed38590f9a4c877dccb8d2e731589f1..7918f27fe7c5ea6d6c986e2800d4a42aa507406c 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/src/ReplicaXYArrays.cxx
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/src/ReplicaXYArrays.cxx
@@ -160,13 +160,13 @@ DOMDocument *doc = element->getOwnerDocument();
                 hepXf0 = geoAXf->getTransform();
             }
             else {
-                geoXf = new GeoTransform (hepXf0);
+                geoXf = makeTransform (hepXf0);
                 hepXf0 = geoXf->getTransform();
             }
             GeoTrf::Transform3D hepXf=hepXf0; 
             for (int i = 0; i < nCopies; ++i) {
 	     	hepXf = GeoTrf::TranslateZ3D(zVal)*GeoTrf::TranslateX3D(xPos[i]) * GeoTrf::TranslateY3D(yPos[i]);
-                xfList->push_back((GeoGraphNode *) new GeoTransform(hepXf));
+                xfList->push_back(makeTransform(hepXf));
             }
     }
     else {
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/src/ReplicaYProcessor.cxx b/GeoModelTools/GeoModelXML/GeoModelXml/src/ReplicaYProcessor.cxx
index 438c0e620bfc54b6dea9042ac83d79c26603ce89..1004de4b2bc2ace794862d43b30afbb27213de71 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/src/ReplicaYProcessor.cxx
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/src/ReplicaYProcessor.cxx
@@ -142,12 +142,12 @@ DOMDocument *doc = element->getOwnerDocument();
                 hepXf0 = geoAXf->getTransform();
             }
             else {
-                geoXf = new GeoTransform (hepXf0);
+                geoXf = makeTransform (hepXf0);
                 hepXf0 = geoXf->getTransform();
             }
             GeoTrf::Transform3D hepXf=hepXf0; 
             for (int i = 0; i < nCopies; ++i) {
-                xfList->push_back((GeoGraphNode *) new GeoTransform(hepXf));
+                xfList->push_back(makeTransform(hepXf));
                 hepXf = hepXf * GeoTrf::TranslateY3D(stepY) ;
             }
     }
diff --git a/GeoModelTools/GeoModelXML/GeoModelXml/src/ReplicaZProcessor.cxx b/GeoModelTools/GeoModelXML/GeoModelXml/src/ReplicaZProcessor.cxx
index b9c247c21636bf2445c56a3e4e4c8a163383c428..b547f6e2ed006a89b7f840b5c97cbafbb64f23b2 100644
--- a/GeoModelTools/GeoModelXML/GeoModelXml/src/ReplicaZProcessor.cxx
+++ b/GeoModelTools/GeoModelXML/GeoModelXml/src/ReplicaZProcessor.cxx
@@ -143,12 +143,12 @@ DOMDocument *doc = element->getOwnerDocument();
                 hepXf0 = geoAXf->getTransform();
             }
             else {
-                geoXf = new GeoTransform (hepXf0);
+                geoXf = makeTransform (hepXf0);
                 hepXf0 = geoXf->getTransform();
             }
             GeoTrf::Transform3D hepXf=hepXf0; 
             for (int i = 0; i < nCopies; ++i) {
-                xfList->push_back((GeoGraphNode *) new GeoTransform(hepXf));
+                xfList->push_back(makeTransform(hepXf));
                 hepXf = hepXf * GeoTrf::TranslateZ3D(stepZ) ;
             }
     }