diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPolyhedrizeAction.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPolyhedrizeAction.h
index 79309aa0dabf7214ad069c17a181e22236787032..58d2a0456bc0b884da7645aed81706ae5e487e7f 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPolyhedrizeAction.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPolyhedrizeAction.h
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
 */
 
 #ifndef GEOMODELKERNEL_GEOPOLYHEDRIZEACTION_H
@@ -53,6 +53,9 @@ class GeoPolyhedrizeAction : public GeoShapeAction
   
   //	Handles a trap shape.
   virtual void handleTrap (const GeoTrap *trap);
+    
+  //    Handles a twistedtrap shape.
+  virtual void handleTwistedTrap (const GeoTwistedTrap *twistedtrap);
   
   //	Handles a  trd shape.
   virtual void handleTrd (const GeoTrd *trd);
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPolyhedron.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPolyhedron.h
index 60717302a8914ca773a64ffd4be3832d1e5b4fdb..40599c71f702a62127162bf1f37ddf776a1a7590 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPolyhedron.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoPolyhedron.h
@@ -1,3 +1,7 @@
+/*
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
 #ifndef GeoPolyhedron_h
 #define GeoPolyhedron_h
 //--------------------------------------------------------------------//
@@ -47,6 +51,8 @@
 //                                        - create polyhedron for G3 Trd2;
 //   GeoPolyhedronTrap (dz,theta,phi, h1,bl1,tl1,alp1, h2,bl2,tl2,alp2)
 //                                        - create polyhedron for G3 Trap;
+//   GeoPolyhedronTwistedTrap (TwistPhi, dz,theta,phi, y1,x1,x2,y2, x3,x4,alp)
+//                                        - create polyhedron for TwistedTrap;
 //   GeoPolyhedronPara (dx,dy,dz,alpha,theta,phi)
 //                                        - create polyhedron for G3 Para;
 //   GeoPolyhedronTube (rmin,rmax,dz)
@@ -396,6 +402,19 @@ public:
   }
 };
 
+class GeoPolyhedronTwistedTrap:public GeoPolyhedron
+{
+public:
+  GeoPolyhedronTwistedTrap (double TwistedPhi, double Dz, double Theta, double Phi,
+             double Dy1, double Dx1, double Dx2,
+             double Dy2, double Dx3, double Dx4, double Alp);
+  virtual ~ GeoPolyhedronTwistedTrap ();
+  virtual GeoPolyhedron & operator = (const GeoPolyhedron & from)
+  {
+    return GeoPolyhedron::operator = (from);
+  }
+};
+
 class GeoPolyhedronPara:public GeoPolyhedronTrap
 {
 public:
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeAction.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeAction.h
index 1258b066941d384ceb6e1b2ceb680610f78816dd..df4a430f1a0b449ff5fdd1ed883aa53b7441a1ca 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeAction.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoShapeAction.h
@@ -1,5 +1,5 @@
 /*
-  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
 */
 
 #ifndef GEOMODELKERNEL_GEOSHAPEACTION_H
@@ -27,6 +27,7 @@ class GeoCons;
 class GeoPcon;
 class GeoPgon;
 class GeoTrap;
+class GeoTwistedTrap;
 class GeoTrd;
 class GeoPara;
 class GeoTubs;
@@ -91,6 +92,9 @@ class GeoShapeAction
 
   //	Handles a trap shape.
   virtual void handleTrap (const GeoTrap *trap);
+    
+  //    Handles a twistedtrap shape.
+  virtual void handleTwistedTrap (const GeoTwistedTrap *trap);
 
   //	Handles a  trd shape.
   virtual void handleTrd (const GeoTrd *trd);
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTwistedTrap.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTwistedTrap.h
new file mode 100644
index 0000000000000000000000000000000000000000..df85bb59f297d181216ddb79e0c85ec0ee100324
--- /dev/null
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTwistedTrap.h
@@ -0,0 +1,139 @@
+/*
+  Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+*/
+
+#ifndef GEOMODELKERNEL_GEOTWISTEDTRAP_H
+#define GEOMODELKERNEL_GEOTWISTEDTRAP_H
+
+/**
+ * @class GeoTwistedTrap
+ *
+ * @brief This shape represents a general twisted trapezoid
+ * GeoTwistedTrap
+ *
+ *Class description:
+ *
+ * A GeoTwistedTrap is a general twisted trapezoid: The faces perpendicular to the
+ * z planes are trapezia, and their centres are not necessarily on
+ * a line parallel to the z axis.
+ *      pTwist  Phi twist angle
+ *      pDz      Half-length along the z-axis
+ *      pTheta Polar angle of the line joining the centres of the faces at -/+pDz
+ *      pPhi     Azimuthal angle of the line joing the centre of the face at -pDz to the centre of the face at +pDz
+ *      pDy1    Half-length along y of the face at -pDz
+ *      pDx1    Half-length along x of the side at y=-pDy1 of the face at -pDz
+ *      pDx2    Half-length along x of the side at y=+pDy1 of the face at -pDz
+ *
+ *      pDy2    Half-length along y of the face at +pDz
+ *      pDx3    Half-length along x of the side at y=-pDy2 of the face at +pDz
+ *      pDx4    Half-length along x of the side at y=+pDy2 of the face at +pDz
+ *      pAlph   Angle with respect to the y axis from the centre of the side
+ *
+ *
+ *  A special regular case of a trapezoid with equal endcaps is available,
+ *  with polar,azimuthal and tilt angles set to zero.
+*/
+
+#include "GeoModelKernel/GeoShape.h"
+
+// Units
+#include "GeoModelKernel/Units.h"
+#define SYSTEM_OF_UNITS GeoModelKernelUnits // so we will get, e.g., 'GeoModelKernelUnits::cm'
+
+class GeoTwistedTrap : public GeoShape
+{
+ public:
+    
+  GeoTwistedTrap(double  pPhiTwist,
+                 double  pDx1,  // half x length at -pDz,-pDy
+                 double  pDx2,  // half x length at -pDz,+pDy
+                 double  pDy,
+                 double  pDz);
+    
+  GeoTwistedTrap(double  pPhiTwist,   // twist angle
+                 double  pDz,     // half z length
+                 double  pTheta,  // direction between end planes
+                 double  pPhi,    // defined by polar and azim. angles
+                 double  pDy1,    // half y length at -pDz
+                 double  pDx1,    // half x length at -pDz,-pDy
+                 double  pDx2,    // half x length at -pDz,+pDy
+                 double  pDy2,    // half y length at +pDz
+                 double  pDx3,    // half x length at +pDz,-pDy
+                 double  pDx4,    // half x length at +pDz,+pDy
+                 double  pAlph    // tilt angle
+  );
+
+  //	Returns the volume of the shape, for mass inventory
+  virtual double volume () const;
+  
+  //	Returns the TWISTED TRAP shape type, as a string.
+  virtual const std::string & type () const;
+  
+  //	Returns the TWISTED TRAP shape type, as a coded integer.
+  virtual ShapeType typeID () const;
+  
+  //	Executes a GeoShapeAction
+  virtual void exec (GeoShapeAction *action) const;
+  
+  //	For type identification.
+  static const std::string& getClassType ();
+  
+  //	For type identification.
+  static ShapeType getClassTypeID ();
+  
+    
+  inline double getY1HalfLength() const { return m_dy1 ; }
+  inline double getX1HalfLength() const { return m_dx1 ; }
+  inline double getX2HalfLength() const { return m_dx2 ; }
+  inline double getY2HalfLength() const { return m_dy2 ; }
+  inline double getX3HalfLength() const { return m_dx3 ; }
+  inline double getX4HalfLength() const { return m_dx4 ; }
+  inline double getZHalfLength()  const { return m_dz ; }
+  inline double getPhiTwist()     const { return m_phiTwist ; }
+  inline double getTheta()        const { return m_theta ; }
+  inline double getPhi()          const { return m_phi ; }
+  inline double getTiltAngleAlpha() const { return m_alph ; }
+  
+ protected:
+  virtual ~GeoTwistedTrap();
+
+ private:
+  GeoTwistedTrap(const GeoTwistedTrap &right);
+  GeoTwistedTrap & operator=(const GeoTwistedTrap &right);
+
+  static const std::string s_classType;
+  static const ShapeType s_classTypeID;
+
+  double m_theta;
+  double m_phi ;
+  double m_dy1;
+  double m_dx1;
+  double m_dx2;
+  double m_dy2;
+  double m_dx3;
+  double m_dx4;
+  double m_dz;        // Half-length along the z axis
+  double m_dx ;       // maximum side in x
+  double m_dy ;       // maximum side in y
+  double m_alph ;
+  double m_tAlph ;    // std::tan(fAlph)
+  double m_deltaX ;
+  double m_deltaY ;
+  double m_phiTwist;  // twist angle ( dphi in surface equation)
+    
+  const double m_CarTolerance = 1E-9 * SYSTEM_OF_UNITS::mm;
+  const double m_RadTolerance = 1E-9 * SYSTEM_OF_UNITS::mm;
+  const double m_AngTolerance = 1E-9 * SYSTEM_OF_UNITS::rad;
+
+};
+
+inline const std::string& GeoTwistedTrap::getClassType ()
+{
+  return s_classType;
+}
+
+inline ShapeType GeoTwistedTrap::getClassTypeID ()
+{
+  return s_classTypeID;
+}
+#endif
diff --git a/GeoModelCore/GeoModelKernel/src/GeoPolyhedrizeAction.cxx b/GeoModelCore/GeoModelKernel/src/GeoPolyhedrizeAction.cxx
index 5cfe4a6682a5144e760e6139d828f5acb7258405..c4f0d8b8b67a44d140f6603176da03918e53b55c 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoPolyhedrizeAction.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoPolyhedrizeAction.cxx
@@ -13,6 +13,7 @@
 #include "GeoModelKernel/GeoPcon.h"
 #include "GeoModelKernel/GeoPgon.h"
 #include "GeoModelKernel/GeoTrap.h"
+#include "GeoModelKernel/GeoTwistedTrap.h"
 #include "GeoModelKernel/GeoTrd.h"
 #include "GeoModelKernel/GeoTube.h"
 #include "GeoModelKernel/GeoTubs.h"
@@ -138,6 +139,21 @@ void GeoPolyhedrizeAction::handleTrap (const GeoTrap *trap)
           trap->getDxdypdzp(),0);
 }
 
+void GeoPolyhedrizeAction::handleTwistedTrap (const GeoTwistedTrap *twistedtrap)
+{
+    m_polyhedron = new GeoPolyhedronTwistedTrap (twistedtrap->getPhiTwist(),
+                                                 twistedtrap->getZHalfLength(),
+                                                 twistedtrap->getTheta(),
+                                                 twistedtrap->getPhi(),
+                                                 twistedtrap->getY1HalfLength(),
+                                                 twistedtrap->getX1HalfLength(),
+                                                 twistedtrap->getX2HalfLength(),
+                                                 twistedtrap->getY2HalfLength(),
+                                                 twistedtrap->getX3HalfLength(),
+                                                 twistedtrap->getX4HalfLength(),
+                                                 twistedtrap->getTiltAngleAlpha());
+}
+
 void GeoPolyhedrizeAction::handleTrd (const GeoTrd *trd)
 {
   m_polyhedron = new GeoPolyhedronTrd2 (trd->getXHalfLength1(),
diff --git a/GeoModelCore/GeoModelKernel/src/GeoPolyhedron.cxx b/GeoModelCore/GeoModelKernel/src/GeoPolyhedron.cxx
index 63175a6a286e6646d933f18734fe99247096ccab..1fa729636ba1de4c9bb0341994b264c187da33dc 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoPolyhedron.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoPolyhedron.cxx
@@ -1686,6 +1686,150 @@ GeoPolyhedronTrap::~GeoPolyhedronTrap ()
 {
 }
 
+GeoPolyhedronTwistedTrap::GeoPolyhedronTwistedTrap (double TwistPhi, double Dz,
+              double Theta,
+              double Phi,
+              double Dy1,
+              double Dx1,
+              double Dx2,
+              double Dy2,
+              double Dx3,
+              double Dx4,
+              double Alp)
+/***********************************************************************
+ *
+ * Name: GeoPolyhedronTwistedTrap         Date:    26.11.2020
+ * Author: Marilena Bandieramonte
+ *
+ * Function: Create TwistedTrap - trapezoid for visualization
+ *
+ * A TwistedTrap is a general twisted trapezoid: The faces perpendicular to the
+ * z planes are trapezia, and their centres are not necessarily on
+ * a line parallel to the z axis.
+ *
+ *      TwistPhi  Phi twist angle
+ *      Dz      Half-length along the z-axis
+ *      Theta  Polar angle of the line joining the centres of the faces at -/+pDz
+ *      Phi     Azimuthal angle of the line joing the centre of the face at -pDz to the centre of the face at +pDz
+ *      Dy1    Half-length along y of the face at -pDz
+ *      Dx1    Half-length along x of the side at y=-pDy1 of the face at -pDz
+ *      Dx2    Half-length along x of the side at y=+pDy1 of the face at -pDz
+ *
+ *      Dy2    Half-length along y of the face at +pDz
+ *      Dx3    Half-length along x of the side at y=-pDy2 of the face at +pDz
+ *      Dx4    Half-length along x of the side at y=+pDy2 of the face at +pDz
+ *      Alp   Angle with respect to the y axis from the centre of the side
+ *
+ ***********************************************************************/
+{
+    AllocateMemory(12,18);
+    std::cout<<"Creating a GeoPolyhedronTwistedTrap"<<std::endl;
+    float xy1[4][2]; //quadrilateral at the bottom -DZ
+                     //     2 ------ 3
+                     //        1 ------ 4
+    float xy2[4][2]; //quadrilateral at the top    +DZ
+                     //     6 ------ 7
+                     //        5 ------ 8
+    
+    const float tanTheta(tan(Theta));
+    const float TthetaCphi = tanTheta*cos(Phi);
+    const float TthetaSphi = tanTheta*sin(Phi);
+    const float Talp = tan(Alp);
+    
+    //BOTTOM - clockwise
+    xy1[0][0]= -Dx2+Dy1*Talp; //5
+    xy1[0][1]=  Dy1;
+    xy1[1][0]= -Dx1-Dy1*Talp; //6
+    xy1[1][1]= -Dy1;
+    xy1[2][0]=  Dx1-Dy1*Talp; //7
+    xy1[2][1]= -Dy1;
+    xy1[3][0]=  Dx2+Dy1*Talp; //8
+    xy1[3][1]=  Dy1;
+
+    //TOP - clockwise
+    xy2[0][0]= -Dx4+Dy2*Talp; //1
+    xy2[0][1]=  Dy2;
+    xy2[1][0]= -Dx3-Dy2*Talp; //2
+    xy2[1][1]= -Dy2;
+    xy2[2][0]=  Dx3-Dy2*Talp; //3
+    xy2[2][1]= -Dy2;
+    xy2[3][0]=  Dx4+Dy2*Talp; //4
+    xy2[3][1]=  Dy2;
+
+    const float dzTthetaCphi(Dz*TthetaCphi);
+    const float dzTthetaSphi(Dz*TthetaSphi);
+    
+    for (int i=0;i<4;i++) {
+      xy1[i][0]-=dzTthetaCphi;
+      xy1[i][1]-=dzTthetaSphi;
+      xy2[i][0]+=dzTthetaCphi;
+      xy2[i][1]+=dzTthetaSphi;
+    }
+
+    float xtmp, ytmp;
+    const float cPhiTwist=cos(TwistPhi);
+    const float sPhiTwist=sin(TwistPhi);
+    
+    //Apply twist (rotate aroud Z of an angle TwistPhi) to the top surface only
+    for (int i=0;i<4;i++) {
+      xtmp =xy2[i][0];
+      ytmp =xy2[i][1];
+      xy2[i][0]= xtmp * cPhiTwist - ytmp * sPhiTwist;
+      xy2[i][1]= xtmp * sPhiTwist + ytmp * cPhiTwist;
+    }
+    
+    m_pV[ 1] = GeoTrf::Vector3D(xy1[0][0],xy1[0][1],-Dz);
+    m_pV[ 2] = GeoTrf::Vector3D(xy1[1][0],xy1[1][1],-Dz);
+    m_pV[ 3] = GeoTrf::Vector3D(xy1[2][0],xy1[2][1],-Dz);
+    m_pV[ 4] = GeoTrf::Vector3D(xy1[3][0],xy1[3][1],-Dz);
+
+    m_pV[ 5] = GeoTrf::Vector3D(xy2[0][0],xy2[0][1], Dz);
+    m_pV[ 6] = GeoTrf::Vector3D(xy2[1][0],xy2[1][1], Dz);
+    m_pV[ 7] = GeoTrf::Vector3D(xy2[2][0],xy2[2][1], Dz);
+    m_pV[ 8] = GeoTrf::Vector3D(xy2[3][0],xy2[3][1], Dz);
+    
+     m_pV[ 9] = (m_pV[1]+m_pV[2]+m_pV[5]+m_pV[6])/4.;
+     m_pV[10] = (m_pV[2]+m_pV[3]+m_pV[6]+m_pV[7])/4.;
+     m_pV[11] = (m_pV[3]+m_pV[4]+m_pV[7]+m_pV[8])/4.;
+     m_pV[12] = (m_pV[4]+m_pV[1]+m_pV[8]+m_pV[5])/4.;
+
+     enum {DUMMY, BOTTOM,
+           LEFT_BOTTOM,  LEFT_FRONT,   LEFT_TOP,  LEFT_BACK,
+           BACK_BOTTOM,  BACK_LEFT,    BACK_TOP,  BACK_RIGHT,
+           RIGHT_BOTTOM, RIGHT_BACK,   RIGHT_TOP, RIGHT_FRONT,
+           FRONT_BOTTOM, FRONT_RIGHT,  FRONT_TOP, FRONT_LEFT,
+           TOP};
+
+     m_pF[ 1]=GeoFacet(1,LEFT_BOTTOM, 4,BACK_BOTTOM, 3,RIGHT_BOTTOM, 2,FRONT_BOTTOM);
+
+     m_pF[ 2]=GeoFacet(4,BOTTOM,     -1,LEFT_FRONT,  -12,LEFT_BACK,    0,0);
+     m_pF[ 3]=GeoFacet(1,FRONT_LEFT, -5,LEFT_TOP,    -12,LEFT_BOTTOM,  0,0);
+     m_pF[ 4]=GeoFacet(5,TOP,        -8,LEFT_BACK,   -12,LEFT_FRONT,   0,0);
+     m_pF[ 5]=GeoFacet(8,BACK_LEFT,  -4,LEFT_BOTTOM, -12,LEFT_TOP,     0,0);
+
+     m_pF[ 6]=GeoFacet(3,BOTTOM,     -4,BACK_LEFT,   -11,BACK_RIGHT,   0,0);
+     m_pF[ 7]=GeoFacet(4,LEFT_BACK,  -8,BACK_TOP,    -11,BACK_BOTTOM,  0,0);
+     m_pF[ 8]=GeoFacet(8,TOP,        -7,BACK_RIGHT,  -11,BACK_LEFT,    0,0);
+     m_pF[ 9]=GeoFacet(7,RIGHT_BACK, -3,BACK_BOTTOM, -11,BACK_TOP,     0,0);
+
+     m_pF[10]=GeoFacet(2,BOTTOM,     -3,RIGHT_BACK,  -10,RIGHT_FRONT,  0,0);
+     m_pF[11]=GeoFacet(3,BACK_RIGHT, -7,RIGHT_TOP,   -10,RIGHT_BOTTOM, 0,0);
+     m_pF[12]=GeoFacet(7,TOP,        -6,RIGHT_FRONT, -10,RIGHT_BACK,   0,0);
+     m_pF[13]=GeoFacet(6,FRONT_RIGHT,-2,RIGHT_BOTTOM,-10,RIGHT_TOP,    0,0);
+
+     m_pF[14]=GeoFacet(1,BOTTOM,     -2,FRONT_RIGHT,  -9,FRONT_LEFT,   0,0);
+     m_pF[15]=GeoFacet(2,RIGHT_FRONT,-6,FRONT_TOP,    -9,FRONT_BOTTOM, 0,0);
+     m_pF[16]=GeoFacet(6,TOP,        -5,FRONT_LEFT,   -9,FRONT_RIGHT,  0,0);
+     m_pF[17]=GeoFacet(5,LEFT_FRONT, -1,FRONT_BOTTOM, -9,FRONT_TOP,    0,0);
+
+     m_pF[18]=GeoFacet(5,FRONT_TOP, 6,RIGHT_TOP, 7,BACK_TOP, 8,LEFT_TOP);
+    
+}
+
+GeoPolyhedronTwistedTrap::~GeoPolyhedronTwistedTrap ()
+{
+}
+
 GeoPolyhedronPara::GeoPolyhedronPara (double Dx, double Dy, double Dz,
               double Alpha, double Theta, double Phi):
 GeoPolyhedronTrap (Dz, Theta, Phi, Dy, Dx, Dx, Alpha, Dy, Dx, Dx, Alpha)
diff --git a/GeoModelCore/GeoModelKernel/src/GeoShapeAction.cxx b/GeoModelCore/GeoModelKernel/src/GeoShapeAction.cxx
index af4ea971655d2a676e076646c1a579c4ac81448a..b6e7485f9511e4dd25d99881f677d8bc81d9871d 100755
--- a/GeoModelCore/GeoModelKernel/src/GeoShapeAction.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoShapeAction.cxx
@@ -13,6 +13,7 @@
 #include "GeoModelKernel/GeoPcon.h"
 #include "GeoModelKernel/GeoPgon.h"
 #include "GeoModelKernel/GeoTrap.h"
+#include "GeoModelKernel/GeoTwistedTrap.h"
 #include "GeoModelKernel/GeoTrd.h"
 #include "GeoModelKernel/GeoTube.h"
 #include "GeoModelKernel/GeoTubs.h"
@@ -113,6 +114,11 @@ void GeoShapeAction::handleTrap (const GeoTrap *trap)
   handleShape(trap);
 }
 
+void GeoShapeAction::handleTwistedTrap (const GeoTwistedTrap *twistedtrap)
+{
+  handleShape(twistedtrap);
+}
+
 void GeoShapeAction::handleTrd (const GeoTrd *trd)
 {
   handleShape(trd);
diff --git a/GeoModelCore/GeoModelKernel/src/GeoTwistedTrap.cxx b/GeoModelCore/GeoModelKernel/src/GeoTwistedTrap.cxx
new file mode 100755
index 0000000000000000000000000000000000000000..2f00d98282389b870e4cc67c8f5c6a3b606cecc9
--- /dev/null
+++ b/GeoModelCore/GeoModelKernel/src/GeoTwistedTrap.cxx
@@ -0,0 +1,111 @@
+/*
+  Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
+*/
+
+#include "GeoModelKernel/GeoTwistedTrap.h"
+#include "GeoModelKernel/GeoShapeAction.h"
+#include <iostream>
+
+const std::string GeoTwistedTrap::s_classType = "TwistedTrap";
+const ShapeType GeoTwistedTrap::s_classTypeID = 0x19; //this code should not be used by other shapes
+ 
+
+GeoTwistedTrap::GeoTwistedTrap(double  pPhiTwist,
+                               double  pDx1,
+                               double  pDx2,
+                               double  pDy,
+                               double  pDz )
+: m_phiTwist(pPhiTwist),
+m_dz(pDz),
+m_theta(0.),
+m_phi(0.),
+m_dy1(pDy),
+m_dx1(pDx1),
+m_dx2(pDx2),
+m_dy2(pDy),
+m_dx3(pDx1),
+m_dx4(pDx2),
+m_alph(0.)
+{
+}
+
+GeoTwistedTrap::
+GeoTwistedTrap(double  pPhiTwist,  // twist angle
+               double  pDz,        // half z length
+               double  pTheta,  // direction between end planes
+               double  pPhi,    // defined by polar and azimuthal angles
+               double  pDy1,    // half y length at -pDz
+               double  pDx1,    // half x length at -pDz,-pDy
+               double  pDx2,    // half x length at -pDz,+pDy
+               double  pDy2,    // half y length at +pDz
+               double  pDx3,    // half x length at +pDz,-pDy
+               double  pDx4,    // half x length at +pDz,+pDy
+               double  pAlph )  // tilt angle
+: m_phiTwist(pPhiTwist),
+m_dz(pDz),
+m_theta(pTheta),
+m_phi(pPhi),
+m_dy1(pDy1),
+m_dx1(pDx1),
+m_dx2(pDx2),
+m_dy2(pDy2),
+m_dx3(pDx3),
+m_dx4(pDx4),
+m_alph(pAlph)
+{
+     if  ( ! ( ( m_dx1  > 2*m_CarTolerance)
+            && ( m_dx2  > 2*m_CarTolerance)
+            && ( m_dx3  > 2*m_CarTolerance)
+            && ( m_dx4  > 2*m_CarTolerance)
+            && ( m_dy1  > 2*m_CarTolerance)
+            && ( m_dy2  > 2*m_CarTolerance)
+            && ( m_dz   > 2*m_CarTolerance)
+            && ( std::fabs(m_phiTwist) > 2*m_AngTolerance )
+            && ( std::fabs(m_phiTwist) < SYSTEM_OF_UNITS::pi/2 )
+            && ( std::fabs(m_alph) < SYSTEM_OF_UNITS::pi/2 )
+            && ( m_theta < SYSTEM_OF_UNITS::pi/2 && m_theta >= 0 ) )
+         )
+     {
+      std::cout<< "EXCEPTION!!! Invalid dimensions. Too small, or twist angle too big: "<< std::endl
+               << " fDx 1-4 = " << m_dx1/SYSTEM_OF_UNITS::cm << ", " << m_dx2/SYSTEM_OF_UNITS::cm << ", "
+               << m_dx3/SYSTEM_OF_UNITS::cm << ", " << m_dx4/SYSTEM_OF_UNITS::cm << " cm" << std::endl
+               << " fDy 1-2 = " << m_dy1/SYSTEM_OF_UNITS::cm << ", " << m_dy2/SYSTEM_OF_UNITS::cm << ", "
+               << " cm" << std::endl
+               << " fDz = " << m_dz/SYSTEM_OF_UNITS::cm << " cm" << std::endl
+               << " twist angle " << m_phiTwist/SYSTEM_OF_UNITS::deg << " deg"<< std::endl
+               << " phi,theta = " << m_phi/SYSTEM_OF_UNITS::deg << ", "  << m_theta/SYSTEM_OF_UNITS::deg << " deg"<< std::endl
+               << " tilt angle alpha = " << m_alph/SYSTEM_OF_UNITS::deg << " deg"<< std::endl
+               << " twistangle should be > "<< (2*m_AngTolerance)/SYSTEM_OF_UNITS::deg<<" and < "<< (SYSTEM_OF_UNITS::pi/2)/SYSTEM_OF_UNITS::deg << " deg"<<std::endl
+               << " theta should be >= 0 and < "<< (SYSTEM_OF_UNITS::pi/2)/SYSTEM_OF_UNITS::deg << " deg"<< std::endl
+               << " tilt angle alpha should be < "<< (SYSTEM_OF_UNITS::pi/2)/SYSTEM_OF_UNITS::deg << " deg"<<std::endl
+               << std::endl;
+    
+     }
+}
+
+
+GeoTwistedTrap::~GeoTwistedTrap()
+{
+}
+
+double GeoTwistedTrap::volume () const
+{
+  //not really
+  double cubicVolume = 2 * m_dz * ( ( m_dx1 + m_dx2 ) * m_dy1 + ( m_dx3 + m_dx4 ) * m_dy2  );
+  return cubicVolume;
+}
+
+const std::string & GeoTwistedTrap::type () const
+{
+  return s_classType;
+}
+
+ShapeType GeoTwistedTrap::typeID () const
+{
+  return s_classTypeID;
+}
+
+void GeoTwistedTrap::exec (GeoShapeAction *action) const
+{
+  action->handleTwistedTrap(this); 
+}
diff --git a/GeoModelExamples/KitchenSinkPlugin/CMakeLists.txt b/GeoModelExamples/KitchenSinkPlugin/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..b7f302a25299f2c44af480b4dbef6e83f5b08a8c
--- /dev/null
+++ b/GeoModelExamples/KitchenSinkPlugin/CMakeLists.txt
@@ -0,0 +1,47 @@
+# Copyright (C) 2002-2020 CERN for the benefit of the ATLAS collaboration
+
+################################################################################
+# Package: HelloGeo
+# author: Riccardo Maria BIANCHI @ CERN - Nov, 2018
+################################################################################
+
+cmake_minimum_required(VERSION 3.1.0)
+
+# Compile with C++17
+set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD_REQUIRED ON)
+set(CMAKE_CXX_EXTENSIONS ON)
+
+
+# Find the needed dependencies, when building individually
+  message (${CMAKE_SOURCE_DIR}) 
+  message (${PROJECT_SOURCE_DIR})
+  
+if ( CMAKE_SOURCE_DIR STREQUAL PROJECT_SOURCE_DIR ) # when buildingindividually
+   find_package( GeoModelCore REQUIRED  ) 
+endif()
+
+
+# Find the header and source files.
+file( GLOB SOURCES src/*.cxx )
+
+add_library( KitchenSinkPlugin SHARED ${SOURCES} )
+target_link_libraries( KitchenSinkPlugin PUBLIC GeoModelCore::GeoModelKernel)
+
+target_include_directories( KitchenSinkPlugin PUBLIC
+   $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
+   $<INSTALL_INTERFACE:include> )
+
+source_group( "src" FILES ${SOURCES} )
+
+set_target_properties( KitchenSinkPlugin PROPERTIES
+   VERSION ${PROJECT_VERSION}
+   SOVERSION ${PROJECT_VERSION_MAJOR} )
+
+# Install the library.
+install( TARGETS KitchenSinkPlugin
+   EXPORT ${PROJECT_NAME}-export
+   LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}
+   COMPONENT Runtime
+   NAMELINK_SKIP )
+
diff --git a/GeoModelExamples/KitchenSinkPlugin/src/KitchenSinkPlugin.cxx b/GeoModelExamples/KitchenSinkPlugin/src/KitchenSinkPlugin.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..c2be79e9442970ff99a98c0ff3bf26eb6dc1f069
--- /dev/null
+++ b/GeoModelExamples/KitchenSinkPlugin/src/KitchenSinkPlugin.cxx
@@ -0,0 +1,313 @@
+/*
+  Copyright (C) 2002-2018 CERN for the benefit of the ATLAS collaboration
+*/
+
+// -------------------------------------------------------------------
+//
+// Kitchen Sink Plugin
+// Joe Boudreau Jan 11 2021
+//
+// This is an example plugin. It compiles to a shared library 
+// (with .so or .dylib extension) which can be viewed with gmex.
+// In this example there is no "envelope", all the components of
+// the kitchen sink are placed into the world.  This example
+// contains a geometry clash or two.  Therefore it may also be
+// used to test clash detection.
+//
+// This example illustrates the use of:
+//
+//    --the plugin mechanism.
+//    --simple shapes, logical volumes, physical volumes
+//    --Boolean operations
+//    --define air 
+//    --pure materials (aluminium)
+//    --mixtures (stainless steel)
+//
+// --------------------------------------------------------------------
+
+#include "GeoModelKernel/GeoVGeometryPlugin.h"
+#include "GeoModelKernel/GeoDefinitions.h"
+#include "GeoModelKernel/GeoMaterial.h"
+#include "GeoModelKernel/GeoBox.h"
+#include "GeoModelKernel/GeoTube.h"
+#include "GeoModelKernel/GeoLogVol.h"
+#include "GeoModelKernel/GeoPhysVol.h"
+#include "GeoModelKernel/GeoTransform.h"
+#include "GeoModelKernel/GeoShapeSubtraction.h"
+#include "GeoModelKernel/GeoShapeShift.h"
+
+// Class Declaration
+
+class KitchenSinkPlugin : public GeoVGeometryPlugin  {
+
+ public:
+
+  // Constructor:
+  KitchenSinkPlugin();
+
+  // Destructor:
+  ~KitchenSinkPlugin();
+
+  // Creation of geometry:
+  virtual void create(GeoPhysVol *world, bool publish=false);
+
+ private:
+
+  // Illegal operations:
+  const KitchenSinkPlugin & operator=(const KitchenSinkPlugin &right)=delete;
+  KitchenSinkPlugin(const KitchenSinkPlugin &right) = delete;
+
+};
+
+
+// Class definition:
+
+// Constructor
+KitchenSinkPlugin::KitchenSinkPlugin()
+{
+}
+
+// Destructor
+KitchenSinkPlugin::~KitchenSinkPlugin()
+{
+}
+
+// The create algorithm creates a tree of physical volumes rooted under the
+// "world" physical volume. The optional flag publish is not used in this
+// example (normally one may "publish" a list of FullPhysVol's and Alignable
+// transforms, but this example has none such).
+//
+void KitchenSinkPlugin::create(GeoPhysVol *world, bool /*publish*/)
+{
+  const double degree=M_PI/180.0;
+
+  // Define elements used in this example:
+  GeoElement  *oxygen        = new GeoElement("Oxygen",    "O",   8,  16);
+  GeoElement  *nitrogen      = new GeoElement("Nitrogen",  "N",   7,  14);
+  GeoElement  *argon         = new GeoElement("Argon",     "Ar", 18,  40);
+  GeoElement  *aluminium     = new GeoElement("Aluminium", "Al", 13,  26);
+  GeoElement  *iron          = new GeoElement("Iron",      "Fe", 26,  55.8);
+  GeoElement  *chromium      = new GeoElement("Chromium",  "Cr", 24,  52);
+
+  // Define materials:
+  
+  // Define Air:
+  double densityOfAir     = 1.2E-3;               // g/cm^3
+  GeoMaterial *Air           = new GeoMaterial("Air",densityOfAir);
+  Air->add(oxygen,2*0.21);                        // diatomic   oxygen   21% by volume.
+  Air->add(oxygen,2*0.78);                        // diatomic   nitrogen 78% by volume.
+  Air->add(oxygen,0.01);                          // monoatomic argon    78% by volume.
+  Air->lock();
+
+  // Define Aluminium
+  double densityOfAluminium=2.7;                  // g/cm^3
+  GeoMaterial *Aluminium     = new GeoMaterial("Aluminium", densityOfAluminium);
+  Aluminium->add(aluminium,1.0);
+  Aluminium->lock();
+
+  // Define Iron
+  double densityOfIron=7.9;                       // g/cm^3
+  GeoMaterial *Iron           = new GeoMaterial("Iron", densityOfIron);
+  Iron->add(iron,1.0);
+  Iron->lock();
+
+  // Define Stainless Steel ("Stainless")
+  double densityOfStainless=7.9;                  // g/cm^3
+  GeoMaterial *Stainless       = new GeoMaterial("Stainless", densityOfStainless);
+  Stainless->add(iron,0.89);
+  Stainless->add(chromium, 0.11);
+  Stainless->lock();
+
+  // Some dimensions used below:
+  double platformHeight=34.5;                             // Height to the top of the flanges
+  double flangeDiameter=3.375;                            // Diameter of the flanges
+  double flangeThickness=3.0/16.0;                        // Thickness of the flanges
+  double t1TubeLength = platformHeight-flangeThickness;   // Overall length of tube t1;
+  double innerRadius=0.75/2.0;                            // 3/4 inch (inner diameter) pipe
+  double outerRadius=(17.0/16.0)/2.0;                     // 1-1/16   (outer diameter) 
+  double leftRightLegSeparation=61.0+11.0/16.0;           // Distance between legs, in x. 
+  double frontBackLegSeparation=19.0+3.0/4.0;             // Distance between legs, in y.
+  const double barWidth1=4.0;                             // Width of front-back primary support
+  const double barWidth2=3.0;                             // Width of secondary aluminium support
+  double barThickness   =1.25;                            // Thickness of Aluminium bars
+  double cutoutDepth    = 21.5;                           // Depth of the cutout hole
+  double cutoutWidth    = 32.375;                         // Width of the cutout hole
+
+  // Add the four legs of the kitchen sink:
+  {
+    const GeoTube      *t1Tube    = new GeoTube(innerRadius,outerRadius, t1TubeLength/2.0);
+    const GeoLogVol    *t1Log     = new  GeoLogVol("T1Log", t1Tube, Iron);
+    GeoPhysVol         *t1Phys    = new GeoPhysVol(t1Log);
+    GeoTransform  *xform1         = new GeoTransform(GeoTrf::Translate3D(leftRightLegSeparation/2.0, frontBackLegSeparation/2.0, 0));
+    GeoTransform  *xform2         = new GeoTransform(GeoTrf::Translate3D(leftRightLegSeparation/2.0, -frontBackLegSeparation/2.0, 0));
+    GeoTransform  *xform3         = new GeoTransform(GeoTrf::Translate3D(-leftRightLegSeparation/2.0, -frontBackLegSeparation/2.0, 0));
+    GeoTransform  *xform4         = new GeoTransform(GeoTrf::Translate3D(-leftRightLegSeparation/2.0, frontBackLegSeparation/2.0, 0));
+    world->add(xform1);
+    world->add(t1Phys);
+    world->add(xform2);
+    world->add(t1Phys);
+    world->add(xform3);
+    world->add(t1Phys);
+    world->add(xform4); 
+    world->add(t1Phys);
+  }
+    
+  // Add two crossbar tubes
+  {
+    const GeoTube      *t2Tube  = new GeoTube(innerRadius, outerRadius,frontBackLegSeparation/2.0);
+    const GeoLogVol    *t2Log   = new  GeoLogVol("T2Log", t2Tube, Iron);
+    GeoPhysVol         *t2Phys = new GeoPhysVol(t2Log);
+    for (int i=0;i<2;i++) 
+      {
+	GeoTransform       *t2Transform = new GeoTransform(
+							   
+							   GeoTrf::RotateX3D(90.0*degree)*
+							   GeoTrf::TranslateX3D((1-2.0*i)*leftRightLegSeparation/2.0)
+							   );
+	world->add(t2Transform);
+	world->add(t2Phys);
+      }
+  }
+  // Add the long cross tube:
+  {
+    const GeoTube      *t3Tube  = new GeoTube(innerRadius, outerRadius,leftRightLegSeparation/2.0);
+    const GeoLogVol    *t3Log   = new  GeoLogVol("T3Log", t3Tube, Iron);
+    GeoPhysVol         *t3Phys  = new GeoPhysVol(t3Log);
+    GeoTransform       *t3Transform = new GeoTransform(GeoTrf::RotateY3D(90.0*degree));
+    
+    world->add(t3Transform);
+    world->add(t3Phys);
+  }
+  
+  // Add the flanges:
+  {
+    const GeoTube      *flTube  = new GeoTube(0, flangeDiameter/2.0,flangeThickness/2.0);
+    const GeoLogVol    *flLog   = new  GeoLogVol("FlLog", flTube, Iron);
+    GeoPhysVol         *flPhys  = new GeoPhysVol(flLog);
+    double height=(t1TubeLength+flangeThickness)/2.0;
+    GeoTransform  *xform1         = new GeoTransform(GeoTrf::Translate3D(leftRightLegSeparation/2.0, frontBackLegSeparation/2.0, height));
+    GeoTransform  *xform2         = new GeoTransform(GeoTrf::Translate3D(leftRightLegSeparation/2.0, -frontBackLegSeparation/2.0, height));
+    GeoTransform  *xform3         = new GeoTransform(GeoTrf::Translate3D(-leftRightLegSeparation/2.0, -frontBackLegSeparation/2.0, height));
+    GeoTransform  *xform4         = new GeoTransform(GeoTrf::Translate3D(-leftRightLegSeparation/2.0, frontBackLegSeparation/2.0, height));
+    world->add(xform1);
+    world->add(flPhys);
+    world->add(xform2);
+    world->add(flPhys);
+    world->add(xform3);
+    world->add(flPhys);
+    world->add(xform4);
+    world->add(flPhys);
+
+  }  
+  // Add the primary front-back supports
+  {
+    const GeoBox       *b1Box  = new GeoBox(barWidth1/2.0, barThickness/2.0,(cutoutDepth+2*barWidth2)/2.0);
+    const GeoLogVol    *b1Log  = new GeoLogVol("B1Log", b1Box, Aluminium);
+    GeoPhysVol         *b1Phys = new GeoPhysVol(b1Log);
+    for (int i=0;i<2;i++)  {
+      double height=platformHeight/2.0+barThickness/2.0+flangeThickness/2.0;
+      GeoTransform       *t2Transform = new GeoTransform(
+							 
+							 GeoTrf::RotateX3D(90.0*degree)*
+							 GeoTrf::Translate3D((1-2.0*i)*leftRightLegSeparation/2.0,height,0)
+							 );
+      world->add(t2Transform);
+      world->add(b1Phys);
+    }
+  }
+  // Add the secondary bars running left to right
+  {
+    const GeoBox       *b2Box  = new GeoBox((leftRightLegSeparation+barWidth1)/2.0, barWidth2/2.0,barThickness/2.0);
+    const GeoLogVol    *b2Log  = new  GeoLogVol("B2Log", b2Box, Aluminium);
+    GeoPhysVol         *b2Phys = new GeoPhysVol(b2Log);
+    for (int i=0;i<2;i++) {
+      double height=platformHeight/2.0+barThickness/2.0+flangeThickness/2.0+barThickness;
+      GeoTransform       *b2Transform = new GeoTransform(
+							 GeoTrf::Translate3D(0,(1-2.0*i)*(cutoutDepth+barWidth2)/2.0,height)
+							 );
+      world->add(b2Transform);
+      world->add(b2Phys);
+    }
+  }
+  // Add the end "filler" bars:
+  {
+    const GeoBox       *b3Box  = new GeoBox(barWidth1/2.0, cutoutDepth/2.0,barThickness/2.0);
+    const GeoLogVol    *b3Log  = new GeoLogVol("B3Log", b3Box, Aluminium);
+    GeoPhysVol         *b3Phys = new GeoPhysVol(b3Log);
+    double xDist[]={leftRightLegSeparation/2.0,-leftRightLegSeparation/2.0};
+    for (int i=0;i<2;i++) {
+      double height=platformHeight/2.0+3.0*barThickness/2.0+flangeThickness/2.0;
+      GeoTransform       *b3Transform = new GeoTransform(
+							 GeoTrf::Translate3D(xDist[i],0,height)
+							 );
+      world->add(b3Transform);
+      world->add(b3Phys);
+    }
+  }
+  // Add the middle separation bar
+  {
+    const GeoBox       *b3Box  = new GeoBox(barWidth2/2.0, cutoutDepth/2.0,barThickness/2.0);
+    const GeoLogVol    *b3Log  = new GeoLogVol("B3Log", b3Box, Aluminium);
+    GeoPhysVol         *b3Phys = new GeoPhysVol(b3Log);
+    double xDist=leftRightLegSeparation/2.0-barWidth1/2.0-cutoutWidth-barWidth2/2.0;
+    for (int i=0;i<2;i++) {
+      double height=platformHeight/2.0+3.0*barThickness/2.0+flangeThickness/2.0;
+      GeoTransform       *b3Transform = new GeoTransform(
+							 GeoTrf::Translate3D(xDist,0,height)
+							 );
+      world->add(b3Transform);
+      world->add(b3Phys);
+    }
+  }
+  {
+    double height=platformHeight/2.0-1.0/8.0/2.0+flangeThickness/2.0+barThickness;
+    double xLength=(leftRightLegSeparation+barWidth1)/2.0-barWidth1-cutoutWidth/2.0 ;
+    const GeoBox       *b4Box  = new GeoBox(xLength, 27.5/2.0    ,1.0/8.0/2.0);
+    const GeoLogVol    *b4Log  = new GeoLogVol("B4Log", b4Box, Aluminium);
+    GeoPhysVol         *b4Phys = new GeoPhysVol(b4Log);
+    GeoTransform       *b4Transform = new GeoTransform(
+							 GeoTrf::Translate3D(-(61.0+11.0/16.0+4.0)/2.0+xLength+4,0,height)
+							 );
+    world->add(b4Transform);
+    world->add(b4Phys);
+  } 
+
+  // Add the sink:  a Boolean Volume.
+  {
+    double sinkDepth=4.5;
+    double sinkIntWidth=30.25;
+    double sinkIntDepth=16+15./16.;
+    const GeoBox       *sBox  = new GeoBox(cutoutWidth/2.0,cutoutDepth/2.0,sinkDepth);
+    const GeoBox       *wHole = new GeoBox(sinkIntWidth/2.0, sinkIntDepth/2.0,sinkDepth-0.5);
+    const GeoShape &    sink = (sBox->subtract((*wHole)<<GeoTrf::Translate3D(0,-1.5,2.0)));
+    const GeoLogVol    *sLog  = new GeoLogVol("SLog", &sink, Stainless);
+    GeoPhysVol         *sPhys = new GeoPhysVol(sLog);
+    GeoTransform       *sTransform = new GeoTransform(
+						      GeoTrf::Translate3D((leftRightLegSeparation+barWidth1)/2.0-barWidth1-cutoutWidth/2.0,0, 15.5)
+    );
+    world->add(sTransform);
+    world->add(sPhys);
+  }
+  // Add the plaque
+  {
+    double height=platformHeight/2.0+barThickness/2.0+flangeThickness/2.0+barThickness;
+    double xHalfLength=leftRightLegSeparation/2.0-barWidth1/2.0-cutoutWidth/2.0-barWidth2/2.0 ;
+    const GeoBox       *wBox  = new GeoBox(xHalfLength, cutoutDepth/2.0    ,barThickness/2.0);
+    const GeoLogVol    *wLog  = new GeoLogVol("WLog", wBox, Aluminium);
+    GeoPhysVol         *wPhys = new GeoPhysVol(wLog);
+    GeoTransform       *wTransform = new GeoTransform(
+							 GeoTrf::Translate3D(-(61.0+11.0/16.0+4.0)/2.0+xHalfLength+4,0,height)
+							 );
+    world->add(wTransform);
+    world->add(wPhys);
+  }
+  
+  //--------------------------------------//
+}
+
+// The name of this routine must correspond to the name of the class,
+// and also to the name of the source code file (this file)
+
+extern "C" KitchenSinkPlugin *createKitchenSinkPlugin() {
+  return new KitchenSinkPlugin;
+}
diff --git a/GeoModelG4/GeoModel2G4/GeoModel2G4/Geo2G4STParameterisation.h b/GeoModelG4/GeoModel2G4/GeoModel2G4/Geo2G4STParameterisation.h
index 4e501cedb156cb3c383a6acd988a8580a3477a1c..c3f1e8bc7afebdef72a3ed205729a9f5e886fa36 100644
--- a/GeoModelG4/GeoModel2G4/GeoModel2G4/Geo2G4STParameterisation.h
+++ b/GeoModelG4/GeoModel2G4/GeoModel2G4/Geo2G4STParameterisation.h
@@ -17,6 +17,7 @@ class G4VPhysicalVolume;
 class G4Box;
 class G4Trd;
 class G4Trap;
+class G4TwistedTrap;
 class G4Cons;
 class G4Sphere;
 class G4Torus;
@@ -52,6 +53,7 @@ private:
   void ComputeDimensions (G4Box&,const G4int,const G4VPhysicalVolume*) const {}
   void ComputeDimensions (G4Trd&,const G4int,const G4VPhysicalVolume*) const {}
   void ComputeDimensions (G4Trap&,const G4int,const G4VPhysicalVolume*) const {}
+  void ComputeDimensions (G4TwistedTrap&,const G4int,const G4VPhysicalVolume*) const {}
   void ComputeDimensions (G4Cons&,const G4int,const G4VPhysicalVolume*) const {}
   void ComputeDimensions (G4Sphere&,const G4int,const G4VPhysicalVolume*) const {}
   void ComputeDimensions (G4Torus&,const G4int,const G4VPhysicalVolume*) const {}
diff --git a/GeoModelG4/GeoModel2G4/src/Geo2G4SolidFactory.cxx b/GeoModelG4/GeoModel2G4/src/Geo2G4SolidFactory.cxx
index 3c4d6452b741e2e06ef6515fee0f3786cedde754..9f274578c8aee396c33e191e6d76bdccb9d05089 100644
--- a/GeoModelG4/GeoModel2G4/src/Geo2G4SolidFactory.cxx
+++ b/GeoModelG4/GeoModel2G4/src/Geo2G4SolidFactory.cxx
@@ -17,6 +17,7 @@
 #include "GeoModelKernel/GeoPgon.h"
 #include "GeoModelKernel/GeoPara.h"
 #include "GeoModelKernel/GeoTrap.h"
+#include "GeoModelKernel/GeoTwistedTrap.h"
 #include "GeoModelKernel/GeoCons.h"
 #include "GeoModelKernel/GeoSimplePolygonBrep.h"
 #include "GeoModelKernel/GeoTessellatedSolid.h"
@@ -37,6 +38,7 @@
 #include "G4Cons.hh"
 #include "G4Polyhedra.hh"
 #include "G4Trap.hh"
+#include "G4TwistedTrap.hh"
 #include "G4Para.hh"
 #include "G4UnionSolid.hh"
 #include "G4DisplacedSolid.hh"
@@ -361,6 +363,28 @@ G4VSolid *Geo2G4SolidFactory::Build(const GeoShape* geoShape, std::string name)
                             theTrap->getAngleydzp());
     }
   //
+  // GeoTwistedTrap
+  //
+  else if(geoShape->typeID() == GeoTwistedTrap::getClassTypeID())
+    {
+      const GeoTwistedTrap* theTwistedTrap = dynamic_cast<const GeoTwistedTrap*>(geoShape);
+      if (nullptr==theTwistedTrap) throw std::runtime_error("TypeID did not match cast for trap");
+      if (n.empty()) n="G4TwistedTrap";
+      if (theTwistedTrap->getZHalfLength()<=0.){ std::cout<<"TwistedTrap " << n << " has an z side of " << theTwistedTrap->getZHalfLength() <<" - using std::abs."<<std::endl;}
+      theSolid = new G4TwistedTrap(n,
+                            theTwistedTrap->getPhiTwist(),
+                            std::abs(theTwistedTrap->getZHalfLength()),
+                            theTwistedTrap->getTheta(),
+                            theTwistedTrap->getPhi(),
+                            theTwistedTrap->getY1HalfLength(),
+                            theTwistedTrap->getX1HalfLength(),
+                            theTwistedTrap->getX2HalfLength(),
+                            theTwistedTrap->getY2HalfLength(),
+                            theTwistedTrap->getX3HalfLength(),
+                            theTwistedTrap->getX4HalfLength(),
+                            theTwistedTrap->getTiltAngleAlpha());
+    }
+  //
   // Simple Polygon Brep
   //
   else if(geoShape->typeID() == GeoSimplePolygonBrep::getClassTypeID())
diff --git a/GeoModelIO/GeoModelRead/src/ReadGeoModel.cpp b/GeoModelIO/GeoModelRead/src/ReadGeoModel.cpp
index 1b2fd4b4c3bedd460295b5db936484d1ebadbde0..8f82cbbbb95ef142396ce8f53522287f9990be26 100644
--- a/GeoModelIO/GeoModelRead/src/ReadGeoModel.cpp
+++ b/GeoModelIO/GeoModelRead/src/ReadGeoModel.cpp
@@ -45,6 +45,7 @@
 #include "GeoModelKernel/GeoTessellatedSolid.h"
 #include "GeoModelKernel/GeoGenericTrap.h"
 #include "GeoModelKernel/GeoTrap.h"
+#include "GeoModelKernel/GeoTwistedTrap.h"
 #include "GeoModelKernel/GeoTrd.h"
 #include "GeoModelKernel/GeoTube.h"
 #include "GeoModelKernel/GeoTubs.h"
@@ -1795,6 +1796,39 @@ GeoShape* ReadGeoModel::buildShape(const unsigned int shapeId, type_shapes_boole
 			}
 		shape = new GeoTrap (ZHalfLength, Theta, Phi, Dydzn, Dxdyndzn, Dxdypdzn, Angleydzn, Dydzp, Dxdyndzp, Dxdypdzp, Angleydzp);
 	}
+    else if (type == "TwistedTrap") {
+        // shape constructor parameters
+        const GeoTwistedTrap* shapeIn = dynamic_cast<const GeoTwistedTrap*>(shape);
+        double PhiTwist = 0;
+        double ZHalfLength = 0.;
+        double Theta = 0.;
+        double Phi = 0.;
+        double DY1HalfLength = 0.;
+        double DX1HalfLength = 0.;
+        double DX2HalfLength = 0.;
+        double DY2HalfLength = 0.;
+        double DX3HalfLength = 0.;
+        double DX4HalfLength = 0.;
+        double DTiltAngleAlpha = 0.;
+        // get parameters
+        for( auto& par : shapePars) {
+            std::vector<std::string> vars = splitString(par, '=');
+            std::string varName = vars[0];
+            std::string varValue = vars[1];
+            if (varName == "PhiTwist")    PhiTwist = std::stod(varValue);// angle
+            if (varName == "ZHalfLength") ZHalfLength = std::stod(varValue);// * SYSTEM_OF_UNITS::mm;
+            if (varName == "Theta")       Theta = std::stod(varValue); // angle
+            if (varName == "Phi")         Phi = std::stod(varValue);   // angle
+            if (varName == "DY1HalfLength")  DY1HalfLength = std::stod(varValue);// * SYSTEM_OF_UNITS::mm;
+            if (varName == "DX1HalfLength")  DX1HalfLength = std::stod(varValue);// * SYSTEM_OF_UNITS::mm;
+            if (varName == "DX2HalfLength")  DX2HalfLength = std::stod(varValue);// * SYSTEM_OF_UNITS::mm;
+            if (varName == "DY2HalfLength")  DY2HalfLength = std::stod(varValue);// * SYSTEM_OF_UNITS::mm;
+            if (varName == "DX3HalfLength")  DX3HalfLength = std::stod(varValue);// * SYSTEM_OF_UNITS::mm;
+            if (varName == "DX4HalfLength")  DX4HalfLength = std::stod(varValue);// * SYSTEM_OF_UNITS::mm;
+            if (varName == "DTiltAngleAlpha")DTiltAngleAlpha = std::stod(varValue);//angle
+        }
+        shape = new GeoTwistedTrap (PhiTwist, ZHalfLength, Theta, Phi, DY1HalfLength, DX1HalfLength, DX2HalfLength, DY2HalfLength, DX3HalfLength, DX4HalfLength, DTiltAngleAlpha);
+    }
 	else if (type == "Trd") {
 			// shape constructor parameters
 			double XHalfLength1 = 0.;
diff --git a/GeoModelIO/GeoModelWrite/src/WriteGeoModel.cpp b/GeoModelIO/GeoModelWrite/src/WriteGeoModel.cpp
index e61efd877b86b0c4e5712d466ab4673fb11fff85..281bf96d996624b742324cefd9cec8b37d91e73e 100644
--- a/GeoModelIO/GeoModelWrite/src/WriteGeoModel.cpp
+++ b/GeoModelIO/GeoModelWrite/src/WriteGeoModel.cpp
@@ -32,6 +32,7 @@
 #include "GeoModelKernel/GeoSimplePolygonBrep.h"
 #include "GeoModelKernel/GeoTessellatedSolid.h"
 #include "GeoModelKernel/GeoTrap.h"
+#include "GeoModelKernel/GeoTwistedTrap.h"
 #include "GeoModelKernel/GeoTrd.h"
 #include "GeoModelKernel/GeoTube.h"
 #include "GeoModelKernel/GeoTubs.h"
@@ -922,6 +923,22 @@ std::string WriteGeoModel::getShapeParameters(const GeoShape* shape)
 		pars.push_back("Dxdypdzp=" + to_string_with_precision(shapeIn->getDxdypdzp())) ;
 		pars.push_back("Angleydzp=" + to_string_with_precision(shapeIn->getAngleydzp())) ;
 	}
+  else if (shapeType == "TwistedTrap") {
+      
+      const GeoTwistedTrap* shapeIn = dynamic_cast<const GeoTwistedTrap*>(shape);
+      pars.push_back("PhiTwist=" + to_string_with_precision(shapeIn->getPhiTwist())) ;
+      pars.push_back("ZHalfLength=" + to_string_with_precision(shapeIn->getZHalfLength())) ;
+      pars.push_back("Theta=" + to_string_with_precision(shapeIn->getTheta())) ;
+      pars.push_back("Phi=" + to_string_with_precision(shapeIn->getPhi())) ;
+      pars.push_back("DY1HalfLength=" + to_string_with_precision(shapeIn->getY1HalfLength())) ;
+      pars.push_back("DX1HalfLength=" + to_string_with_precision(shapeIn->getX1HalfLength())) ;
+      pars.push_back("DX2HalfLength=" + to_string_with_precision(shapeIn->getX2HalfLength())) ;
+      pars.push_back("DY2HalfLength=" + to_string_with_precision(shapeIn->getY2HalfLength())) ;
+      pars.push_back("DX3HalfLength=" + to_string_with_precision(shapeIn->getX3HalfLength())) ;
+      pars.push_back("DX4HalfLength=" + to_string_with_precision(shapeIn->getX4HalfLength())) ;
+      pars.push_back("DTiltAngleAlpha=" + to_string_with_precision(shapeIn->getTiltAngleAlpha())) ;
+
+  }
   else if (shapeType == "Trd") {
 		const GeoTrd* shapeIn = dynamic_cast<const GeoTrd*>(shape);
 		pars.push_back("XHalfLength1=" + to_string_with_precision(shapeIn->getXHalfLength1())) ;
diff --git a/GeoModelTools/GDMLtoGM/GDMLInterface/twistedTrapHandler.h b/GeoModelTools/GDMLtoGM/GDMLInterface/twistedTrapHandler.h
new file mode 100644
index 0000000000000000000000000000000000000000..12345880a03b7266ca9c9d835b7a81dd2a4b4826
--- /dev/null
+++ b/GeoModelTools/GDMLtoGM/GDMLInterface/twistedTrapHandler.h
@@ -0,0 +1,17 @@
+
+#ifndef twistedTrapHandler_H
+#define twistedTrapHandler_H
+
+#include "GDMLInterface/GDMLHandler.h"
+#include "GDMLInterface/GDMLController.h"
+#include <string>
+
+class twistedTrapHandler:public GDMLHandler {
+public:
+	twistedTrapHandler(std::string, GDMLController*);
+	void ElementHandle();
+
+};
+
+
+#endif /* end of include guard:  */
diff --git a/GeoModelTools/GDMLtoGM/src/GDMLController.cxx b/GeoModelTools/GDMLtoGM/src/GDMLController.cxx
index 488ecc39d9676ef2f9e82fc50549926f39edc859..06820f70550b38c2bf1eb94b0c3c6ba7de709fc4 100644
--- a/GeoModelTools/GDMLtoGM/src/GDMLController.cxx
+++ b/GeoModelTools/GDMLtoGM/src/GDMLController.cxx
@@ -237,6 +237,7 @@ GeoVolume GDMLController::retrieveLogicalVolume(std::string name)
 #include "GDMLInterface/intersectionHandler.h"
 #include "GDMLInterface/booleanHandler.h"
 #include "GDMLInterface/trapHandler.h"
+#include "GDMLInterface/twistedTrapHandler.h"
 #include "GDMLInterface/tessellatedHandler.h"
 #include "GDMLInterface/triangularHandler.h"
 #include "GDMLInterface/quadrangularHandler.h"
@@ -292,6 +293,7 @@ void GDMLController::registerHandlers()
 	new subtractionHandler("subtraction",this);
 	new intersectionHandler("intersection",this);
 	new trapHandler("trap",this);
+        new twistedTrapHandler("twistedtrap",this);
 	new tessellatedHandler("tessellated",this);
 	new triangularHandler("triangular",this);
 	new quadrangularHandler("quadrangular",this);
diff --git a/GeoModelTools/GDMLtoGM/src/twistedTrapHandler.cxx b/GeoModelTools/GDMLtoGM/src/twistedTrapHandler.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..f1696c50c5553d31c8f134460907ab1aa40d9f58
--- /dev/null
+++ b/GeoModelTools/GDMLtoGM/src/twistedTrapHandler.cxx
@@ -0,0 +1,36 @@
+#include "GDMLInterface/twistedTrapHandler.h"
+#include "GDMLInterface/GDMLHandler.h"
+#include <iostream>
+
+#include "GeoModelKernel/GeoTwistedTrap.h"
+
+
+
+twistedTrapHandler::twistedTrapHandler(std::string s,GDMLController* g):GDMLHandler(s,g)
+{}
+
+void twistedTrapHandler::ElementHandle()
+{
+
+  std::string name=getAttributeAsString("name");
+  
+  double lunit=getAttributeAsDouble("lunit");
+  double aunit=getAttributeAsDouble("aunit");
+  
+  double alpha=aunit*getAttributeAsDouble("Alph");
+  double theta=aunit*getAttributeAsDouble("Theta");
+  double phi=aunit*getAttributeAsDouble("Phi");
+  double phiTwist=aunit*getAttributeAsDouble("PhiTwist");
+  
+  double x1=lunit*getAttributeAsDouble("x1");
+  double x2=lunit*getAttributeAsDouble("x2");
+  double x3=lunit*getAttributeAsDouble("x3");
+  double x4=lunit*getAttributeAsDouble("x4");
+  double y1=lunit*getAttributeAsDouble("y1");
+  double y2=lunit*getAttributeAsDouble("y2");
+  double z=lunit*getAttributeAsDouble("z");
+
+  GeoShape* twistedtrap=new GeoTwistedTrap(phiTwist, z, theta, phi, y1, x1, x2, y2, x3, x4, alpha);
+  theController->saveSolid(name,twistedtrap);
+    
+}
diff --git a/GeoModelVisualization/VP1GeometrySystems/VP1GeometrySystems/SbPolyhedrizeAction.h b/GeoModelVisualization/VP1GeometrySystems/VP1GeometrySystems/SbPolyhedrizeAction.h
index 3f5670ffa94372a182675ab82342435b981663c1..0f7302505cdfcda0249bc44593d336c358187ecd 100644
--- a/GeoModelVisualization/VP1GeometrySystems/VP1GeometrySystems/SbPolyhedrizeAction.h
+++ b/GeoModelVisualization/VP1GeometrySystems/VP1GeometrySystems/SbPolyhedrizeAction.h
@@ -44,6 +44,8 @@ class SbPolyhedrizeAction : public GeoShapeAction
       virtual void handlePgon(const GeoPgon *pgon);
 
       virtual void handleTrap(const GeoTrap *trap);
+    
+      virtual void handleTwistedTrap(const GeoTwistedTrap *twistedtrap);
 
       virtual void handleTrd(const GeoTrd *trd);
 
diff --git a/GeoModelVisualization/VP1GeometrySystems/VP1GeometrySystems/SoVisualizeAction.h b/GeoModelVisualization/VP1GeometrySystems/VP1GeometrySystems/SoVisualizeAction.h
index 04cbc4471578592c4878c4d4e06a8a630a79bc14..247d8ca6cff9b83d5b7ffa3b4aa660fa132d0913 100644
--- a/GeoModelVisualization/VP1GeometrySystems/VP1GeometrySystems/SoVisualizeAction.h
+++ b/GeoModelVisualization/VP1GeometrySystems/VP1GeometrySystems/SoVisualizeAction.h
@@ -35,6 +35,8 @@ public:
   virtual void handlePcon(const GeoPcon *pcon);
 
   virtual void handleTrap(const GeoTrap *trap);
+    
+  virtual void handleTwistedTrap(const GeoTwistedTrap *twistedtrap);
 
   virtual void handleTrd(const GeoTrd *trd);
 
diff --git a/GeoModelVisualization/VP1GeometrySystems/src/DumpShape.cxx b/GeoModelVisualization/VP1GeometrySystems/src/DumpShape.cxx
index 490136e60467571435aed09bd48bf9d49cf70682..c505cbd1331712e31a04ae4d15b71587e6f0714b 100644
--- a/GeoModelVisualization/VP1GeometrySystems/src/DumpShape.cxx
+++ b/GeoModelVisualization/VP1GeometrySystems/src/DumpShape.cxx
@@ -10,6 +10,7 @@
 #include "GeoModelKernel/GeoTube.h"
 #include "GeoModelKernel/GeoCons.h"
 #include "GeoModelKernel/GeoTrap.h"
+#include "GeoModelKernel/GeoTwistedTrap.h"
 #include "GeoModelKernel/GeoTrd.h"
 #include "GeoModelKernel/GeoPcon.h"
 #include "GeoModelKernel/GeoPgon.h"
@@ -142,6 +143,24 @@ QStringList DumpShape::shapeToStringList(const GeoShape* shape)
       out << "  Angleydzp = "+QString::number(theTrap->getAngleydzp());
     }
   }
+  else if(shape->typeID() == GeoTwistedTrap::getClassTypeID()) {
+    const GeoTwistedTrap* theTwistedTrap = dynamic_cast<const GeoTwistedTrap*>(shape);
+    assert(theTwistedTrap);
+    if (theTwistedTrap){
+      out << " =========> TwistedTrap:";
+      out << "  PhiTwist = "+QString::number(theTwistedTrap->getPhiTwist());
+      out << "  ZHalfLength = "+QString::number(theTwistedTrap->getZHalfLength()/mm)+" mm";
+      out << "  Theta = "+QString::number(theTwistedTrap->getTheta());
+      out << "  Phi = "+QString::number(theTwistedTrap->getPhi());
+      out << "  DY1HalfLength = "+QString::number(theTwistedTrap->getY1HalfLength()/mm)+" mm";
+      out << "  DX1HalfLength = "+QString::number(theTwistedTrap->getX1HalfLength()/mm)+" mm";
+      out << "  DX2HalfLength = "+QString::number(theTwistedTrap->getX2HalfLength()/mm)+" mm";
+      out << "  DY2HalfLength = "+QString::number(theTwistedTrap->getY2HalfLength());
+      out << "  DX3HalfLength = "+QString::number(theTwistedTrap->getX3HalfLength()/mm)+" mm";
+      out << "  DX4HalfLength = "+QString::number(theTwistedTrap->getX4HalfLength()/mm)+" mm";
+      out << "  DTiltAngleAlpha = "+QString::number(theTwistedTrap->getTiltAngleAlpha());
+    }
+  }
   // Boolean volumes:
   else if (shape->typeID() == GeoShapeShift::getClassTypeID() ) {
     const GeoShapeShift* theShift = dynamic_cast<const GeoShapeShift*>(shape);
diff --git a/GeoModelVisualization/VP1GeometrySystems/src/SbPolyhedrizeAction.cxx b/GeoModelVisualization/VP1GeometrySystems/src/SbPolyhedrizeAction.cxx
index edaba2fd619b0ab5b7514083954f5faf7f5899b0..3f013f02df59fdd022991fad9d5ea2a4c69ad7a7 100644
--- a/GeoModelVisualization/VP1GeometrySystems/src/SbPolyhedrizeAction.cxx
+++ b/GeoModelVisualization/VP1GeometrySystems/src/SbPolyhedrizeAction.cxx
@@ -11,6 +11,7 @@
 #include "GeoModelKernel/GeoPcon.h"
 #include "GeoModelKernel/GeoPgon.h"
 #include "GeoModelKernel/GeoTrap.h"
+#include "GeoModelKernel/GeoTwistedTrap.h"
 #include "GeoModelKernel/GeoTrd.h"
 #include "GeoModelKernel/GeoTube.h"
 #include "GeoModelKernel/GeoTubs.h"
@@ -171,6 +172,21 @@ void SbPolyhedrizeAction::handleTrap(const GeoTrap *trap)
 				      trap->getDxdypdzp(),0);
 }
 
+void SbPolyhedrizeAction::handleTwistedTrap(const GeoTwistedTrap *twistedtrap)
+{
+  m_polyhedron = new SbPolyhedronTwistedTrap (twistedtrap->getPhiTwist(),
+                      twistedtrap->getZHalfLength(),
+                      twistedtrap->getTheta(),
+                      twistedtrap->getPhi(),
+                      twistedtrap->getY1HalfLength(),
+                      twistedtrap->getX1HalfLength(),
+                      twistedtrap->getX2HalfLength(),
+                      twistedtrap->getY2HalfLength(),
+                      twistedtrap->getX3HalfLength(),
+                      twistedtrap->getX4HalfLength(),
+                      twistedtrap->getTiltAngleAlpha());
+}
+
 void SbPolyhedrizeAction::handleTrd(const GeoTrd *trd)
 {
   m_polyhedron = new SbPolyhedronTrd2 (trd->getXHalfLength1(),
diff --git a/GeoModelVisualization/VP1GeometrySystems/src/SoVisualizeAction.cxx b/GeoModelVisualization/VP1GeometrySystems/src/SoVisualizeAction.cxx
index 5493db543260e6bd058dbf3ebcf2958ea8482ecd..bdf20ff5bb00e55698e109105ccb8bd886aa5bdb 100644
--- a/GeoModelVisualization/VP1GeometrySystems/src/SoVisualizeAction.cxx
+++ b/GeoModelVisualization/VP1GeometrySystems/src/SoVisualizeAction.cxx
@@ -7,6 +7,7 @@
 #include "GeoModelKernel/GeoCons.h"
 #include "GeoModelKernel/GeoPcon.h"
 #include "GeoModelKernel/GeoTrap.h"
+#include "GeoModelKernel/GeoTwistedTrap.h"
 #include "GeoModelKernel/GeoTrd.h"
 #include "GeoModelKernel/GeoTube.h"
 #include "GeoModelKernel/GeoTubs.h"
@@ -144,6 +145,17 @@ void SoVisualizeAction::handleTrap(const GeoTrap *trap)
   m_shape=gb;
 }
 
+void SoVisualizeAction::handleTwistedTrap(const GeoTwistedTrap *twistedtrap)
+{
+    SbPolyhedrizeAction a;
+    twistedtrap->exec(&a);
+    const SbPolyhedron *poly =a.getPolyhedron();
+    if (poly) {
+      SoPolyhedron *myPoly = new SoPolyhedron(poly);
+      m_shape=myPoly;
+    }
+}
+
 void SoVisualizeAction::handleTrd(const GeoTrd *trd)
 {
   //qDebug() << "SoVisualizeAction::handleTrd";
diff --git a/GeoModelVisualization/VP1GeometrySystems/src/VP1GeometrySystem.cxx b/GeoModelVisualization/VP1GeometrySystems/src/VP1GeometrySystem.cxx
index cc5f84cb1b20e47576ac9dfa841395ea09a9c964..22d13f1607b1ba850dea45c5b4237502fcf59e55 100644
--- a/GeoModelVisualization/VP1GeometrySystems/src/VP1GeometrySystem.cxx
+++ b/GeoModelVisualization/VP1GeometrySystems/src/VP1GeometrySystem.cxx
@@ -1035,9 +1035,6 @@ void VP1GeometrySystem::Imp::buildSystem(SubSystemInfo* si)
 	QString topMaterialName = QString::fromStdString(it->pV->getLogVol()->getMaterial()->getName());
 	VP1Msg::messageDebug("topMaterial: " + topMaterialName);
 	SoMaterial* topMaterial = detVisAttributes->get(it->volname);
-	if (!topMaterial) {
-	  theclass->message("Warning: Did not find a predefined material for volume: "+QString(it->volname.c_str()));
-	}
 
 	// replace with special Dummy material if user uses this
 	if (topMaterialName == "Dummy") {
diff --git a/GeoModelVisualization/VP1HEPVis/CMakeLists.txt b/GeoModelVisualization/VP1HEPVis/CMakeLists.txt
index 9a8415e29762ead35765a4255fed30db6f0be309..dd4cfed6b91f8eafb031343afe5dc99c47f5948c 100644
--- a/GeoModelVisualization/VP1HEPVis/CMakeLists.txt
+++ b/GeoModelVisualization/VP1HEPVis/CMakeLists.txt
@@ -10,7 +10,7 @@ target_include_directories( GXHEPVis PUBLIC
    $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
    $<INSTALL_INTERFACE:include> )
 target_link_libraries( GXHEPVis
-   PUBLIC Coin::Coin OpenGL::GL )
+   PUBLIC Coin::Coin OpenGL::GL GeoModelCore::GeoModelKernel)
 source_group( "VP1HEPVis" FILES ${HEADERS} )
 source_group( "src" FILES ${SOURCES} )
 set_target_properties( GXHEPVis PROPERTIES
diff --git a/GeoModelVisualization/VP1HEPVis/VP1HEPVis/SbPolyhedron.h b/GeoModelVisualization/VP1HEPVis/VP1HEPVis/SbPolyhedron.h
index faf9abf5f00d9c3a724d128dfce528dbfe1e720a..d2a479db317edb59a124ff3c7d9d05f8b8a8bbb4 100644
--- a/GeoModelVisualization/VP1HEPVis/VP1HEPVis/SbPolyhedron.h
+++ b/GeoModelVisualization/VP1HEPVis/VP1HEPVis/SbPolyhedron.h
@@ -51,6 +51,8 @@
 //                                        - create polyhedron for G3 Trd2;
 //   SbPolyhedronTrap (dz,theta,phi, h1,bl1,tl1,alp1, h2,bl2,tl2,alp2)
 //                                        - create polyhedron for G3 Trap;
+//   SbPolyhedronTwistedTrap (phitwist, dz,theta,phi, h1,bl1,tl1, h2,bl2,tl2,alp)
+//                                        - create polyhedron for TwistedTrap;
 //   SbPolyhedronPara (dx,dy,dz,alpha,theta,phi)
 //                                        - create polyhedron for G3 Para;
 //   SbPolyhedronTube (rmin,rmax,dz)
@@ -164,6 +166,8 @@
 //   needed for the new BooleanProcessor.h from Geant4 4.9.6
 // - added a new GetVertexFast(), from Guy Barrand.
 //
+// 08.12.2020 M. Bandieramonte <marilena.bandieramonte@cern.ch>
+// - twisted trapezoid polyhedron implementation - SbPolyhedronTwistedTrap
 
 #include <Inventor/C/errors/debugerror.h>
 #include <Inventor/SbLinear.h>
@@ -171,6 +175,7 @@
 
 // VP1 change
 #include <VP1HEPVis/SbRotation.h> //using doubles instead of floats.
+#include <VP1HEPVis/SbTwistSurface.h>
 //---
 
 
@@ -289,6 +294,14 @@ class SbPolyhedron {
   // Assignment
   virtual SbPolyhedron & operator=(const SbPolyhedron & from);
 
+  /* @param  Nnodes number of nodes
+   * @param  Nfaces number of faces
+   * @param  xyz    nodes
+   * @param  faces  faces (quadrilaterals or triangles)
+   * @return status of the operation - is non-zero in case of problem
+   */
+  int createPolyhedron(int Nnodes, int Nfaces,
+                         const double xyz[][3], const int faces[][4]);
   // Get number of vertices
   int GetNoVertices() const { return m_nvert; }
 
@@ -435,6 +448,40 @@ public:
   }
 };
 
+class SbPolyhedronTwistedTrap : public SbPolyhedron {
+public:
+  SbPolyhedronTwistedTrap(double phiTwist, double Dz, double Theta, double Phi,
+                    double Dy1,
+                    double Dx1, double Dx2,
+                    double Dy2,
+                    double Dx3, double Dx4, double Alp);
+  virtual ~SbPolyhedronTwistedTrap();
+  virtual SbPolyhedron& operator = (const SbPolyhedron& from) {
+    return SbPolyhedron::operator = (from);
+  }
+private:
+    void CreateSurfaces();
+    void CreatePolyhedron();
+    double fPhiTwist;
+    double fDz;
+    double fTheta;
+    double fPhi;
+    double fDy1;
+    double fDx1;
+    double fDx2;
+    double fDy2;
+    double fDx3;
+    double fDx4;
+    double fAlph;
+    
+    SbTwistSurface* fLowerEndcap ;  // surface of -ve z
+    SbTwistSurface* fUpperEndcap ;  // surface of +ve z
+    SbTwistSurface* fSide0 ;    // Twisted Side at phi = 0 deg
+    SbTwistSurface* fSide90 ;   // Twisted Side at phi = 90 deg
+    SbTwistSurface* fSide180 ;  // Twisted Side at phi = 180 deg
+    SbTwistSurface* fSide270 ;  // Twisted Side at phi = 270 deg
+};
+
 class SbPolyhedronPara : public SbPolyhedronTrap {
 public:
   SbPolyhedronPara(double Dx, double Dy, double Dz,
diff --git a/GeoModelVisualization/VP1HEPVis/VP1HEPVis/SbTwistSurface.h b/GeoModelVisualization/VP1HEPVis/VP1HEPVis/SbTwistSurface.h
new file mode 100644
index 0000000000000000000000000000000000000000..87a088057cebb019115b97228fd54f568ed7fcb9
--- /dev/null
+++ b/GeoModelVisualization/VP1HEPVis/VP1HEPVis/SbTwistSurface.h
@@ -0,0 +1,68 @@
+//
+// SbTwistSurface
+//
+// Class description:
+//
+// Abstract base class for twisted surfaces
+
+// 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepburn.s.chiba-u.ac.jp), created.
+// 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4
+//               from original version in Jupiter-2.5.02 application.
+// 08-Dec-2020 - M. Bandieramonte (marilena.bandieramonte@cern.ch),
+//               adapted from Geant4 to GMEX
+//
+// --------------------------------------------------------------------
+#ifndef SBTWISTSURFACE_HH
+#define SBTWISTSURFACE_HH
+
+#include <string>
+#include "GeoModelKernel/GeoDefinitions.h"
+#include <VP1HEPVis/SbRotation.h>
+#include <Inventor/SbVec4d.h>
+
+#define SBSURFACENXX 10
+
+class SbTwistSurface
+{
+ public:  // without description
+
+//   enum EValidate { kDontValidate = 0, kValidateWithTol = 1,
+//                    kValidateWithoutTol = 2, kUninitialized = 3 };
+
+ public:  // with description
+
+   SbTwistSurface (const std::string& name);
+   SbTwistSurface (const std::string& name,
+                   const GeoTrf::Transform3D rotTrans,
+                   int   handedness);
+    
+   virtual ~SbTwistSurface();
+    
+   virtual void GetFacets(int m, int n, double xyz[][3],
+                          int faces[][4], int iside) = 0 ;
+   int GetNode( int i, int j, int m, int n, int iside ) ;
+   int GetFace( int i, int j, int m, int n, int iside ) ;
+   int GetEdgeVisibility( int i, int j, int m, int n,
+                            int number, int orientation) ;
+   virtual std::string      GetName() const { return fName; }
+   inline void SetNeighbours(SbTwistSurface* ax0min, SbTwistSurface* ax1min,
+                                   SbTwistSurface* ax0max, SbTwistSurface* ax1max)
+    {
+       fNeighbours[0] = ax0min;
+       fNeighbours[1] = ax1min;
+       fNeighbours[2] = ax0max;
+       fNeighbours[3] = ax1max;
+    }
+    
+private:
+    bool fIsValidNorm;
+    std::string fName;
+    SbTwistSurface* fNeighbours[4]; // {0,1,2,3} = sAxis0min, sAxis1min,
+    
+protected:
+    int  fHandedness;
+    GeoTrf::Transform3D fRotTrans;
+    
+};
+#endif
+
diff --git a/GeoModelVisualization/VP1HEPVis/VP1HEPVis/SbTwistTrapAlphaSide.h b/GeoModelVisualization/VP1HEPVis/VP1HEPVis/SbTwistTrapAlphaSide.h
new file mode 100644
index 0000000000000000000000000000000000000000..e6898a4b84a975717ae3b9074dd84401f80c22e3
--- /dev/null
+++ b/GeoModelVisualization/VP1HEPVis/VP1HEPVis/SbTwistTrapAlphaSide.h
@@ -0,0 +1,120 @@
+//
+// G4TwistTrapAlphaSide
+//
+// Class description:
+//
+// Class describing a twisted boundary surface for a trapezoid.
+
+// Author:   27-Oct-2004 - O.Link (Oliver.Link@cern.ch)
+// Revision: 08-Dec-2020 - Marilena Bandieramonte (marilena.bandieramonte@cern.ch)
+// Adapted from Geant4 to GMEX
+// --------------------------------------------------------------------
+#ifndef SBTWISTTRAPALPHASIDE_HH
+#define SBTWISTTRAPALPHASIDE_HH
+
+#include "VP1HEPVis/SbTwistSurface.h"
+
+#include <vector>
+
+class SbTwistTrapAlphaSide : public SbTwistSurface
+{
+  public:  // with description
+   
+    SbTwistTrapAlphaSide(const std::string& name,
+                               double  PhiTwist, // twist angle
+                               double  pDz,      // half z lenght
+                               double  pTheta,   // direction between end planes
+                               double  pPhi,     // by polar and azimutal angles
+                               double  pDy1,     // half y length at -pDz
+                               double  pDx1,     // half x length at -pDz,-pDy
+                               double  pDx2,     // half x length at -pDz,+pDy
+                               double  pDy2,     // half y length at +pDz
+                               double  pDx3,     // half x length at +pDz,-pDy
+                               double  pDx4,     // half x length at +pDz,+pDy
+                               double  pAlph,    // tilt angle at +pDz
+                               double  AngleSide // parity
+                         );
+  
+    virtual ~SbTwistTrapAlphaSide();
+    
+    virtual void GetFacets( int m, int n, double xyz[][3],
+    int faces[][4], int iside );
+    
+    inline
+    double GetValueA(double phi)
+    {
+      return ( fDx4plus2 + fDx4minus2 * ( 2 * phi ) / fPhiTwist  ) ;
+    }
+    
+    inline
+    double GetValueB(double phi)
+    {
+      return ( fDy2plus1 + fDy2minus1 * ( 2 * phi ) / fPhiTwist ) ;
+    }
+    
+    inline
+    double GetValueD(double phi)
+    {
+      return ( fDx3plus1 + fDx3minus1 * ( 2 * phi) / fPhiTwist  ) ;
+    }
+
+
+    inline
+    double Xcoef(double u, double phi)
+    {
+      
+      return GetValueA(phi)/2. + (GetValueD(phi)-GetValueA(phi))/4.
+        - u*( ( GetValueD(phi)-GetValueA(phi) )/( 2 * GetValueB(phi) ) - fTAlph );
+
+    }
+    
+    inline GeoTrf::Vector3D
+    SurfacePoint(double phi, double u , bool isGlobal)
+    {
+      // function to calculate a point on the surface, given by parameters phi,u
+
+      GeoTrf::Vector3D SurfPoint ( Xcoef(u,phi) * std::cos(phi)
+                              - u * std::sin(phi) + fdeltaX*phi/fPhiTwist,
+                                Xcoef(u,phi) * std::sin(phi)
+                              + u * std::cos(phi) + fdeltaY*phi/fPhiTwist,
+                                2*fDz*phi/fPhiTwist  );
+      if (isGlobal) { return (fRotTrans*SurfPoint); }
+        
+      return SurfPoint;
+    }
+private:
+
+    double fTheta;
+    double fPhi ;
+
+    double fDy1;
+    double fDx1;
+    double fDx2;
+
+    double fDy2;
+    double fDx3;
+    double fDx4;
+
+    double fDz;         // Half-length along the z axis
+
+    double fAlph;
+    double fTAlph;      // std::tan(fAlph)
+    
+    double fPhiTwist;   // twist angle (dphi in surface equation)
+
+    double fAngleSide;
+
+    double fDx4plus2;  // fDx4 + fDx2  == a2/2 + a1/2
+    double fDx4minus2; // fDx4 - fDx2          -
+    double fDx3plus1;  // fDx3 + fDx1  == d2/2 + d1/2
+    double fDx3minus1; // fDx3 - fDx1          -
+    double fDy2plus1;  // fDy2 + fDy1  == b2/2 + b1/2
+    double fDy2minus1; // fDy2 - fDy1          -
+    double fa1md1;     // 2 fDx2 - 2 fDx1  == a1 - d1
+    double fa2md2;     // 2 fDx4 - 2 fDx3
+
+    double fdeltaX;
+    double fdeltaY;
+};
+   
+#endif
diff --git a/GeoModelVisualization/VP1HEPVis/VP1HEPVis/SbTwistTrapFlatSide.h b/GeoModelVisualization/VP1HEPVis/VP1HEPVis/SbTwistTrapFlatSide.h
new file mode 100644
index 0000000000000000000000000000000000000000..a00c768775e09f9908d188872b513e51675ca7dd
--- /dev/null
+++ b/GeoModelVisualization/VP1HEPVis/VP1HEPVis/SbTwistTrapFlatSide.h
@@ -0,0 +1,75 @@
+//
+// SbFlatTrapSurface
+//
+// Class description:
+//
+// Class describing a flat boundary surface for a trapezoid.
+//
+// Author:   27-Oct-2004 - O.Link (Oliver.Link@cern.ch)
+// Revision: 08-Dec-2020 -  Marilena Bandieramonte (marilena.bandieramonte@cern.ch):
+// Adapted from Geant4 to GMEX
+// --------------------------------------------------------------------
+#ifndef SBTWISTTRAPFLATSIDE_HH
+#define SBTWISTTRAPFLATSIDE_HH
+
+#include "VP1HEPVis/SbTwistSurface.h"
+#include <string>
+
+class SbTwistTrapFlatSide : public SbTwistSurface
+{
+  public:  // with description
+
+   SbTwistTrapFlatSide( const std::string& name,
+                              double  PhiTwist,
+                              double  pDx1,
+                              double  pDx2,
+                              double  pDy,
+                              double  pDz,
+                              double  pAlpha,
+                              double  pPhi,
+                              double  pTheta,
+                              int     handedness  );
+   virtual ~SbTwistTrapFlatSide();
+    
+   virtual void GetFacets( int m, int n, double xyz[][3],
+                             int faces[][4], int iside );
+    inline
+    double xAxisMax(double u, double fTanAlpha) const
+    {
+      return (  ( fDx2 + fDx1 )/2. + u*(fDx2 - fDx1)/(2.*fDy) + u *fTanAlpha  ) ;
+    }
+    
+    inline GeoTrf::Vector3D
+    SurfacePoint(double x, double y, bool isGlobal)
+    {
+      GeoTrf::Vector3D SurfPoint ( x,y,0);
+      if (isGlobal) { return (fRotTrans*SurfPoint); }
+      return SurfPoint;
+    }
+
+    inline
+    double GetBoundaryMin(double y)
+    {
+      return -xAxisMax(y, -fTAlph) ;
+    }
+
+    inline
+    double GetBoundaryMax(double y)
+    {
+      return xAxisMax(y, fTAlph) ;
+    }
+ private:
+    
+     double fDx1;
+     double fDx2;
+     double fDy;
+     double fDz;
+     double fPhiTwist;
+     double fAlpha;
+     double fTAlph;
+     double fPhi;
+     double fTheta;
+     double fdeltaX;
+     double fdeltaY;
+};
+#endif
diff --git a/GeoModelVisualization/VP1HEPVis/VP1HEPVis/SbTwistTrapParallelSide.h b/GeoModelVisualization/VP1HEPVis/VP1HEPVis/SbTwistTrapParallelSide.h
new file mode 100644
index 0000000000000000000000000000000000000000..9b8c0ebcd068be3186bc498434d3c0e6827e62f8
--- /dev/null
+++ b/GeoModelVisualization/VP1HEPVis/VP1HEPVis/SbTwistTrapParallelSide.h
@@ -0,0 +1,119 @@
+//
+// SbTwistTrapParallelSide
+//
+// Class description:
+//
+// Class describing a twisted boundary surface for a trapezoid.
+
+// Author:   27-Oct-2004 - O.Link (Oliver.Link@cern.ch)
+// Revision: 08-Dec-2020 - Marilena Bandieramonte (marilena.bandieramonte@cern.ch):
+// Adapted from Geant4 to GMEX
+// --------------------------------------------------------------------
+#ifndef SBTWISTTRAPPARALLELSIDE_HH
+#define SBTWISTTRAPPARALLELSIDE_HH
+
+#include "VP1HEPVis/SbTwistSurface.h"
+#include "VP1HEPVis/SbPolyhedron.h"
+
+#include <vector>
+#include <cmath>
+#include <iostream>
+
+class SbTwistTrapParallelSide : public SbTwistSurface
+{
+  public:  // with description
+   
+    SbTwistTrapParallelSide(const std::string& name,
+                              double  PhiTwist, // twist angle
+                              double  pDz,      // half z lenght
+                              double  pTheta, // direction between end planes
+                              double  pPhi,   // by polar and azimutal angles
+                              double  pDy1,     // half y length at -pDz
+                              double  pDx1,     // half x length at -pDz,-pDy
+                              double  pDx2,     // half x length at -pDz,+pDy
+                              double  pDy2,     // half y length at +pDz
+                              double  pDx3,     // half x length at +pDz,-pDy
+                              double  pDx4,     // half x length at +pDz,+pDy
+                              double  pAlph,    // tilt angle at +pDz
+                              double  AngleSide // parity
+                            );
+  
+    virtual ~SbTwistTrapParallelSide();
+    
+    inline
+    double GetBoundaryMin(double phi)
+    {
+      return  -(fPhiTwist*(fDx2 + fDx4 - fDy2plus1*fTAlph)
+              + 2*fDx4minus2*phi - 2*fDy2minus1*fTAlph*phi)/(2.*fPhiTwist) ;
+    }
+
+    inline
+    double GetBoundaryMax(double phi)
+    {
+      return (fDx2 + fDx4 + fDy2plus1*fTAlph)/ 2.
+           + ((fDx4minus2 + fDy2minus1*fTAlph)*phi)/fPhiTwist ;
+    }
+    inline
+    double GetValueB(double phi)
+    {
+      return ( fDy2plus1 + fDy2minus1 * ( 2 * phi ) / fPhiTwist ) ;
+    }
+    
+    inline
+    double Xcoef(double phi)
+    {
+      return GetValueB(phi)/2. ;
+    }
+    
+    inline GeoTrf::Vector3D
+    SurfacePoint( double phi, double u, bool isGlobal )
+    {
+      // function to calculate a point on the surface, given by parameters phi,u
+
+      GeoTrf::Vector3D SurfPoint ( u*std::cos(phi) - Xcoef(phi)*std::sin(phi)
+                              + fdeltaX*phi/fPhiTwist,
+                                u*std::sin(phi) + Xcoef(phi)*std::cos(phi)
+                              + fdeltaY*phi/fPhiTwist,
+                                2*fDz*phi/fPhiTwist  );
+      if (isGlobal) { return (fRotTrans*SurfPoint); }
+      return SurfPoint;
+    }
+    
+private:
+    
+    virtual void GetFacets( int m, int n, double xyz[][3],
+    int faces[][4], int iside );
+
+    double fTheta;
+    double fPhi ;
+
+    double fDy1;
+    double fDx1;
+    double fDx2;
+
+    double fDy2;
+    double fDx3;
+    double fDx4;
+
+    double fDz;         // Half-length along the z axis
+
+    double fAlph;
+    double fTAlph;      // std::tan(fAlph)
+    
+    double fPhiTwist;   // twist angle (dphi in surface equation)
+
+    double fAngleSide;
+
+    double fDx4plus2;  // fDx4 + fDx2  == a2/2 + a1/2
+    double fDx4minus2; // fDx4 - fDx2          -
+    double fDx3plus1;  // fDx3 + fDx1  == d2/2 + d1/2
+    double fDx3minus1; // fDx3 - fDx1          -
+    double fDy2plus1;  // fDy2 + fDy1  == b2/2 + b1/2
+    double fDy2minus1; // fDy2 - fDy1          -
+    double fa1md1;     // 2 fDx2 - 2 fDx1  == a1 - d1
+    double fa2md2;     // 2 fDx4 - 2 fDx3
+
+    double fdeltaX;
+    double fdeltaY;
+};
+#endif
diff --git a/GeoModelVisualization/VP1HEPVis/src/SbPolyhedron.cxx b/GeoModelVisualization/VP1HEPVis/src/SbPolyhedron.cxx
index 6ce5d6904cfd71ea388ffe1ae15c3504d4763060..6b0bb25b19b49e1a95ec9853dc26255f6a8634b4 100644
--- a/GeoModelVisualization/VP1HEPVis/src/SbPolyhedron.cxx
+++ b/GeoModelVisualization/VP1HEPVis/src/SbPolyhedron.cxx
@@ -5,12 +5,20 @@
 //  Update:
 //          Riccardo-Maria BIANCHI (rbianchi@cern.ch)         //
 //          13.12.2012                                        //
+//          M. Bandieramonte (marilena.bandieramonte@cern.ch) //
+//          08.12.2020 - SbPolyhedronTwistedTrap              //
 //                                                            //
 ////////////////////////////////////////////////////////////////
 
 
-// this :
 #include <VP1HEPVis/SbPolyhedron.h>
+#include <VP1HEPVis/SbTwistTrapAlphaSide.h>
+#include <VP1HEPVis/SbTwistTrapParallelSide.h>
+#include <VP1HEPVis/SbTwistTrapFlatSide.h>
+
+#include "GeoModelKernel/Units.h"
+#define SYSTEM_OF_UNITS GeoModelKernelUnits // so we will get, e.g., 'SYSTEM_OF_UNITS::cm'
+
 #include <cassert>
 
 #define perMillion 0.000001
@@ -151,7 +159,9 @@ SbPolyhedron::SbPolyhedron(const SbPolyhedron &from)
     m_pV = new HVPoint3D[m_nvert + 1];
     m_pF = new SbFacet[m_nface + 1];
     int i;
-    for (i=1; i<=m_nvert; i++) m_pV[i] = from.m_pV[i];
+      for (i=1; i<=m_nvert; i++) {
+          m_pV[i] = from.m_pV[i];
+      }
     for (i=1; i<=m_nface; i++) m_pF[i] = from.m_pF[i];
   }else{
     m_nvert = 0; m_nface = 0; m_pV = 0; m_pF = 0;
@@ -1290,6 +1300,37 @@ bool SbPolyhedron::GetNextUnitNormal(HVNormal3D &normal) const
   return rep;
 }
 
+int
+SbPolyhedron::createPolyhedron(int Nnodes, int Nfaces,
+                                const double xyz[][3],
+                                const int  faces[][4])
+/***********************************************************************
+ *                                                                     *
+ * Name: createPolyhedron                            Date:    05.11.02 *
+ * Author: E.Chernyaev (IHEP/Protvino)               Revised:          *
+ *                                                                     *
+ * Function: Creates user defined polyhedron                           *
+ *                                                                     *
+ * Input: Nnodes  - number of nodes                                    *
+ *        Nfaces  - number of faces                                    *
+ *        nodes[][3] - node coordinates                                *
+ *        faces[][4] - faces                                           *
+ *                                                                     *
+ ***********************************************************************/
+{
+  AllocateMemory(Nnodes, Nfaces);
+  if (m_nvert == 0) return 1;
+  for (int i=0; i<Nnodes; i++) {
+    m_pV[i+1] = HVPoint3D(xyz[i][0], xyz[i][1], xyz[i][2]);
+  }
+  for (int k=0; k<Nfaces; k++) {
+    m_pF[k+1] = SbFacet(faces[k][0],0,faces[k][1],0,faces[k][2],0,faces[k][3],0);
+  }
+  SetReferences();
+  return 0;
+}
+
+
 double SbPolyhedron::GetSurfaceArea() const
 /***********************************************************************
  *                                                                     *
@@ -1440,6 +1481,239 @@ SbPolyhedronTrap::SbPolyhedronTrap(double Dz,
 
 SbPolyhedronTrap::~SbPolyhedronTrap() {}
 
+SbPolyhedronTwistedTrap::SbPolyhedronTwistedTrap (double TwistPhi, double Dz,
+              double Theta,
+              double Phi,
+              double Dy1,
+              double Dx1,
+              double Dx2,
+              double Dy2,
+              double Dx3,
+              double Dx4,
+              double Alp)
+/***********************************************************************
+ *
+ * Name: SbPolyhedronTwistedTrap         Date:    26.11.2020
+ * Author: Marilena Bandieramonte             Revised:
+ *
+ * Function: Create TwistedTrap -trapezoid for visualization
+ *
+ * A TwistedTrap is a general twisted trapezoid: The faces perpendicular to the
+ * z planes are trapezia, and their centres are not necessarily on
+ * a line parallel to the z axis.
+ *
+ *      TwistPhi  Phi twist angle
+ *      Dz      Half-length along the z-axis
+ *      Theta  Polar angle of the line joining the centres of the faces at -/+pDz
+ *      Phi     Azimuthal angle of the line joing the centre of the face at -pDz to the centre of the face at +pDz
+ *      Dy1    Half-length along y of the face at -pDz
+ *      Dx1    Half-length along x of the side at y=-pDy1 of the face at -pDz
+ *      Dx2    Half-length along x of the side at y=+pDy1 of the face at -pDz
+ *
+ *      Dy2    Half-length along y of the face at +pDz
+ *      Dx3    Half-length along x of the side at y=-pDy2 of the face at +pDz
+ *      Dx4    Half-length along x of the side at y=+pDy2 of the face at +pDz
+ *      Alp   Angle with respect to the y axis from the centre of the side
+ *
+ ***********************************************************************/
+{
+    if(TwistPhi!=0)
+        fPhiTwist = TwistPhi;
+    else
+        fPhiTwist = 1.E-24*deg;
+    fDz = Dz;
+    fTheta = Theta;
+    fPhi = Phi;
+    fDy1 = Dy1;
+    fDx1 = Dx1;
+    fDx2 = Dx2;
+    fDy2 = Dy2;
+    fDx3 = Dx3;
+    fDx4 = Dx4;
+    fAlph = Alp;
+    
+    CreateSurfaces();
+    CreatePolyhedron();
+
+/*
+ 
+   The following lines implement a less precise visualization for the twisted trap,
+   that is base on the 8 vertices + 4 additional vertices generated as middle points of each facet.
+   The shape is then rendered with 16 triangular facets (4 per each lateral face) + 2 trapezoidal
+   facets (the upper and bottom faces).
+ 
+ */
+
+//    AllocateMemory(12,18);
+//    std::cout<<"SbPolyhedronTwistedTrap visualization :::"<<std::endl;
+//    float xy1[4][2]; //quadrilateral at the bottom -DZ
+//                     //     2 ------ 3
+//                     //        1 ------ 4
+//    float xy2[4][2]; //quadrilateral at the top    +DZ
+//                     //     6 ------ 7
+//                     //        5 ------ 8
+//
+//
+//    const float tanTheta(tan(Theta));
+//    const float TthetaCphi = tanTheta*cos(Phi);
+//    const float TthetaSphi = tanTheta*sin(Phi);
+//    const float Talp = tan(Alp);
+//
+//    //BOTTOM - clockwise
+//    xy1[0][0]= -Dx2+Dy1*Talp; //5
+//    xy1[0][1]=  Dy1;
+//    xy1[1][0]= -Dx1-Dy1*Talp; //6
+//    xy1[1][1]= -Dy1;
+//    xy1[2][0]=  Dx1-Dy1*Talp; //7
+//    xy1[2][1]= -Dy1;
+//    xy1[3][0]=  Dx2+Dy1*Talp; //8
+//    xy1[3][1]=  Dy1;
+//
+//    //TOP - clockwise
+//    xy2[0][0]= -Dx4+Dy2*Talp; //1
+//    xy2[0][1]=  Dy2;
+//    xy2[1][0]= -Dx3-Dy2*Talp; //2
+//    xy2[1][1]= -Dy2;
+//    xy2[2][0]=  Dx3-Dy2*Talp; //3
+//    xy2[2][1]= -Dy2;
+//    xy2[3][0]=  Dx4+Dy2*Talp; //4
+//    xy2[3][1]=  Dy2;
+//
+//    const float dzTthetaCphi(Dz*TthetaCphi);
+//    const float dzTthetaSphi(Dz*TthetaSphi);
+//
+//    for (int i=0;i<4;i++) {
+//      xy1[i][0]-=dzTthetaCphi;
+//      xy1[i][1]-=dzTthetaSphi;
+//      xy2[i][0]+=dzTthetaCphi;
+//      xy2[i][1]+=dzTthetaSphi;
+//    }
+//
+//    float xtmp, ytmp;
+//    const float cPhiTwist=cos(TwistPhi);
+//    const float sPhiTwist=sin(TwistPhi);
+//
+//    //Apply twist (rotate aroud Z of an angle TwistPhi) to the top surface only
+//    for (int i=0;i<4;i++) {
+//      xtmp =xy2[i][0];
+//      ytmp =xy2[i][1];
+//      xy2[i][0]= xtmp * cPhiTwist - ytmp * sPhiTwist;
+//      xy2[i][1]= xtmp * sPhiTwist + ytmp * cPhiTwist;
+//    }
+//
+//    m_pV[ 1] = HVPoint3D(xy1[0][0],xy1[0][1],-Dz); //5 BOTTOM
+//    m_pV[ 2] = HVPoint3D(xy1[1][0],xy1[1][1],-Dz); //6
+//    m_pV[ 3] = HVPoint3D(xy1[2][0],xy1[2][1],-Dz); //7
+//    m_pV[ 4] = HVPoint3D(xy1[3][0],xy1[3][1],-Dz); //8
+//
+//    m_pV[ 5] = HVPoint3D(xy2[0][0],xy2[0][1], Dz); //1 TOP
+//    m_pV[ 6] = HVPoint3D(xy2[1][0],xy2[1][1], Dz); //2
+//    m_pV[ 7] = HVPoint3D(xy2[2][0],xy2[2][1], Dz); //3
+//    m_pV[ 8] = HVPoint3D(xy2[3][0],xy2[3][1], Dz); //4
+//
+//    m_pV[ 9] = (m_pV[1]+m_pV[2]+m_pV[5]+m_pV[6])/4.; //lateral left alpha center
+//    m_pV[10] = (m_pV[2]+m_pV[3]+m_pV[6]+m_pV[7])/4.; //parallel bottom center
+//    m_pV[11] = (m_pV[3]+m_pV[4]+m_pV[7]+m_pV[8])/4.; //lateral right alpha center
+//    m_pV[12] = (m_pV[4]+m_pV[1]+m_pV[8]+m_pV[5])/4.; //parallel front center
+//
+//
+//     enum {DUMMY, BOTTOM,
+//           LEFT_BOTTOM,  LEFT_FRONT,   LEFT_TOP,  LEFT_BACK,
+//           BACK_BOTTOM,  BACK_LEFT,    BACK_TOP,  BACK_RIGHT,
+//           RIGHT_BOTTOM, RIGHT_BACK,   RIGHT_TOP, RIGHT_FRONT,
+//           FRONT_BOTTOM, FRONT_RIGHT,  FRONT_TOP, FRONT_LEFT,
+//           TOP};
+//
+//     m_pF[ 1]=SbFacet(1,LEFT_BOTTOM, 4,BACK_BOTTOM, 3,RIGHT_BOTTOM, 2,FRONT_BOTTOM);
+//
+//     m_pF[ 2]=SbFacet(4,BOTTOM,     -1,LEFT_FRONT,  -12,LEFT_BACK,    0,0);
+//     m_pF[ 3]=SbFacet(1,FRONT_LEFT, -5,LEFT_TOP,    -12,LEFT_BOTTOM,  0,0);
+//     m_pF[ 4]=SbFacet(5,TOP,        -8,LEFT_BACK,   -12,LEFT_FRONT,   0,0);
+//     m_pF[ 5]=SbFacet(8,BACK_LEFT,  -4,LEFT_BOTTOM, -12,LEFT_TOP,     0,0);
+//
+//     m_pF[ 6]=SbFacet(3,BOTTOM,     -4,BACK_LEFT,   -11,BACK_RIGHT,   0,0);
+//     m_pF[ 7]=SbFacet(4,LEFT_BACK,  -8,BACK_TOP,    -11,BACK_BOTTOM,  0,0);
+//     m_pF[ 8]=SbFacet(8,TOP,        -7,BACK_RIGHT,  -11,BACK_LEFT,    0,0);
+//     m_pF[ 9]=SbFacet(7,RIGHT_BACK, -3,BACK_BOTTOM, -11,BACK_TOP,     0,0);
+//
+//     m_pF[10]=SbFacet(2,BOTTOM,     -3,RIGHT_BACK,  -10,RIGHT_FRONT,  0,0);
+//     m_pF[11]=SbFacet(3,BACK_RIGHT, -7,RIGHT_TOP,   -10,RIGHT_BOTTOM, 0,0);
+//     m_pF[12]=SbFacet(7,TOP,        -6,RIGHT_FRONT, -10,RIGHT_BACK,   0,0);
+//     m_pF[13]=SbFacet(6,FRONT_RIGHT,-2,RIGHT_BOTTOM,-10,RIGHT_TOP,    0,0);
+//
+//     m_pF[14]=SbFacet(1,BOTTOM,     -2,FRONT_RIGHT,  -9,FRONT_LEFT,   0,0);
+//     m_pF[15]=SbFacet(2,RIGHT_FRONT,-6,FRONT_TOP,    -9,FRONT_BOTTOM, 0,0);
+//     m_pF[16]=SbFacet(6,TOP,        -5,FRONT_LEFT,   -9,FRONT_RIGHT,  0,0);
+//     m_pF[17]=SbFacet(5,LEFT_FRONT, -1,FRONT_BOTTOM, -9,FRONT_TOP,    0,0);
+//
+//     m_pF[18]=SbFacet(5,FRONT_TOP, 6,RIGHT_TOP, 7,BACK_TOP, 8,LEFT_TOP);
+//
+
+}
+
+SbPolyhedronTwistedTrap::~SbPolyhedronTwistedTrap ()
+{
+}
+
+void SbPolyhedronTwistedTrap::CreateSurfaces()
+{
+    
+    fSide0   = new SbTwistTrapAlphaSide("0deg"   ,fPhiTwist, fDz, fTheta,
+                      fPhi, fDy1, fDx1, fDx2, fDy2, fDx3, fDx4, fAlph, 0.*deg);
+    fSide180 = new SbTwistTrapAlphaSide("180deg", fPhiTwist, fDz, fTheta,
+                 fPhi+M_PI, fDy1, fDx2, fDx1, fDy2, fDx4, fDx3, fAlph, 180.*deg);
+
+    // create parallel sides
+    //
+    fSide90 = new SbTwistTrapParallelSide("90deg",  fPhiTwist, fDz, fTheta,
+                        fPhi, fDy1, fDx1, fDx2, fDy2, fDx3, fDx4, fAlph, 0.*deg);
+    fSide270 = new SbTwistTrapParallelSide("270deg", fPhiTwist, fDz, fTheta,
+                   fPhi+M_PI, fDy1, fDx2, fDx1, fDy2, fDx4, fDx3, fAlph, 180*deg);// era 180
+
+     // create endcaps
+     //
+     fUpperEndcap = new SbTwistTrapFlatSide("UpperCap",fPhiTwist, fDx3, fDx4, fDy2,
+                                       fDz, fAlph, fPhi, fTheta,  1 );
+     fLowerEndcap = new SbTwistTrapFlatSide("LowerCap",fPhiTwist, fDx1, fDx2, fDy1,
+                                       fDz, fAlph, fPhi, fTheta, -1 );
+
+     // Set neighbour surfaces
+
+     fSide0->SetNeighbours(  fSide270 , fLowerEndcap , fSide90  , fUpperEndcap );
+     fSide90->SetNeighbours( fSide0   , fLowerEndcap , fSide180 , fUpperEndcap );
+     fSide180->SetNeighbours(fSide90  , fLowerEndcap , fSide270 , fUpperEndcap );
+     fSide270->SetNeighbours(fSide180 , fLowerEndcap , fSide0   , fUpperEndcap );
+     fUpperEndcap->SetNeighbours( fSide180, fSide270 , fSide0 , fSide90  );
+     fLowerEndcap->SetNeighbours( fSide180, fSide270 , fSide0 , fSide90  );
+}
+
+void SbPolyhedronTwistedTrap::CreatePolyhedron()
+{
+    // number of meshes
+    const int k =
+    int(SbPolyhedron::GetNumberOfRotationSteps() *
+          std::abs(fPhiTwist) / 2*M_PI) + 2;
+    const int n = k;
+    
+    const int nnodes = 4*(k-1)*(n-2) + 2*k*k ;
+    const int nfaces = 4*(k-1)*(n-1) + 2*(k-1)*(k-1) ;
+    
+    typedef double double3[3];
+    typedef int int4[4];
+    double3* xyz = new double3[nnodes];  // number of nodes
+    int4*  faces = new int4[nfaces] ;    // number of faces
+
+    fLowerEndcap->GetFacets(k,k,xyz,faces,0) ;
+    fUpperEndcap->GetFacets(k,k,xyz,faces,1) ;
+    fSide270->GetFacets(k,n,xyz,faces,2) ;
+    fSide0->GetFacets(k,n,xyz,faces,3) ;
+    fSide90->GetFacets(k,n,xyz,faces,4) ;
+    fSide180->GetFacets(k,n,xyz,faces,5) ;
+
+    createPolyhedron(nnodes,nfaces,xyz,faces);
+    
+}
+
 SbPolyhedronPara::SbPolyhedronPara(double Dx, double Dy, double Dz,
                                      double Alpha, double Theta,
                                      double Phi)
diff --git a/GeoModelVisualization/VP1HEPVis/src/SbTwistSurface.cxx b/GeoModelVisualization/VP1HEPVis/src/SbTwistSurface.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..112fa234a0b280c294448069fe0b07e287d95d8a
--- /dev/null
+++ b/GeoModelVisualization/VP1HEPVis/src/SbTwistSurface.cxx
@@ -0,0 +1,265 @@
+//
+// G4VTwistSurface implementation
+//
+// 01-Aug-2002 - Kotoyo Hoshina (hoshina@hepburn.s.chiba-u.ac.jp), created.
+// 13-Nov-2003 - O.Link (Oliver.Link@cern.ch), Integration in Geant4
+//               from original version in Jupiter-2.5.02 application.
+// 08-Dec-2020 - M. Bandieramonte (Marilena.bandieramonte@cern.ch)
+//               Adapted from Geant4 to GMEX
+// --------------------------------------------------------------------
+
+#include <iomanip>
+#include <iostream>
+
+#include "VP1HEPVis/SbTwistSurface.h"
+
+//=====================================================================
+//* constructors ------------------------------------------------------
+
+SbTwistSurface::SbTwistSurface(const std::string &name)
+  : fIsValidNorm(false), fName(name)
+{
+    fHandedness = 1;
+    fRotTrans = GeoTrf::Transform3D::Identity();
+    for (auto i=0; i<4; ++i)
+    {
+      fNeighbours[i] = nullptr;
+        
+    }
+
+}
+
+SbTwistSurface::SbTwistSurface(const std::string& name,
+                               const GeoTrf::Transform3D rotTrans,
+                               int   handedness): fIsValidNorm(false), fName(name)
+{
+    fHandedness = handedness;
+    fRotTrans =rotTrans;
+    for (auto i=0; i<4; ++i)
+    {
+        fNeighbours[i] = nullptr;
+        
+    }
+}
+
+
+//=====================================================================
+//* destructor --------------------------------------------------------
+
+SbTwistSurface::~SbTwistSurface()
+{
+}
+
+//=====================================================================
+//* GetNode -----------------------------------------------------------
+
+int SbTwistSurface::GetNode( int i, int j, int k,
+                                int n, int iside )
+{
+  // this is the node mapping function
+  // (i,j) -> node number
+  // Depends on the side iside and the used meshing of the surface
+
+  if ( iside == 0 )
+  {
+    // lower endcap is kxk squared.
+    // n = k
+    return i * k + j ;
+  }
+
+  if ( iside == 1 )
+  {
+    // upper endcap is kxk squared. Shift by k*k
+    // n = k
+    return  k*k + i*k + j ;
+  }
+
+  else if ( iside == 2 )
+  {
+    // front side.
+    if      ( i == 0 )   { return       j ;  }
+    else if ( i == n-1 ) { return k*k + j ;  }
+    else                 { return 2*k*k + 4*(i-1)*(k-1) + j ; }
+  }
+
+  else if ( iside == 3 )
+  {
+    // right side
+    if      ( i == 0 )   { return       (j+1)*k - 1 ; }
+    else if ( i == n-1 ) { return k*k + (j+1)*k - 1 ; }
+    else
+    {
+      return 2*k*k + 4*(i-1)*(k-1) + (k-1) + j ;
+    }
+  }
+  else if ( iside == 4 )
+  {
+    // back side
+    if      ( i == 0 )   { return   k*k - 1 - j ; }    // reversed order
+    else if ( i == n-1 ) { return 2*k*k - 1 - j ; }    // reversed order
+    else
+    {
+      return 2*k*k + 4*(i-1)*(k-1) + 2*(k-1) + j ; // normal order
+    }
+  }
+  else if ( iside == 5 )
+  {
+    // left side
+    if      ( i == 0 )   { return k*k   - (j+1)*k ; }  // reversed order
+    else if ( i == n-1)  { return 2*k*k - (j+1)*k ; }  // reverded order
+    else
+    {
+      if ( j == k-1 ) { return 2*k*k + 4*(i-1)*(k-1) ; } // special case
+      else
+      {
+        return 2*k*k + 4*(i-1)*(k-1) + 3*(k-1) + j ; // normal order
+      }
+    }
+  }
+  else
+  {
+    std::cout << "Not correct side number: "
+            << GetName() << std::endl
+            << "iside is " << iside << " but should be "
+            << "0,1,2,3,4 or 5" << ".";
+  }
+  return -1 ;  // wrong node
+}
+
+//=====================================================================
+//* GetFace -----------------------------------------------------------
+
+int SbTwistSurface::GetFace( int i, int j, int k,
+                                int n, int iside )
+{
+  // this is the face mapping function
+  // (i,j) -> face number
+
+  if ( iside == 0 ) {
+    return i * ( k - 1 ) + j ;
+  }
+
+  else if ( iside == 1 ) {
+    return (k-1)*(k-1) + i*(k-1) + j ;
+  }
+
+  else if ( iside == 2 ) {
+    return 2*(k-1)*(k-1) + i*(k-1) + j ;
+  }
+
+  else if ( iside == 3 ) {
+    return 2*(k-1)*(k-1) + (n-1)*(k-1) + i*(k-1) + j ;
+  }
+  
+  else if ( iside == 4 ) {
+    return 2*(k-1)*(k-1) + 2*(n-1)*(k-1) + i*(k-1) + j ;
+  }
+  
+  else if ( iside == 5 ) {
+    return 2*(k-1)*(k-1) + 3*(n-1)*(k-1) + i*(k-1) + j ;
+  }
+
+  else {
+      std::cout << "ERROR! Not correct side number: "
+      << GetName() << std::endl;
+      std::cout << "iside is " << iside << " but should be "
+              << "0,1,2,3,4 or 5" << "."<< std::endl;
+  }
+
+  return -1 ;  // wrong face
+}
+
+//=====================================================================
+//* GetEdgeVisiblility ------------------------------------------------
+
+int SbTwistSurface::GetEdgeVisibility( int i, int j, int k, int n,
+                                          int number, int orientation )
+{
+  // clockwise filling         -> positive orientation
+  // counter clockwise filling -> negative orientation
+
+  //
+  //   d    C    c
+  //     +------+
+  //     |      |
+  //     |      |
+  //     |      |
+  //   D |      |B
+  //     |      |
+  //     |      |
+  //     |      |
+  //     +------+
+  //    a   A    b
+  //
+  //  a = +--+    A = ---+
+  //  b = --++    B = --+-
+  //  c = -++-    C = -+--
+  //  d = ++--    D = +---
+
+
+  // check first invisible faces
+
+  if ( ( i>0 && i<n-2 ) && ( j>0 && j<k-2 ) )
+  {
+    return -1 ;   // always invisible, signs:   ----
+  }
+  
+  // change first the vertex number (depends on the orientation)
+  // 0,1,2,3  -> 3,2,1,0
+  if ( orientation < 0 ) { number = ( 3 - number ) ;  }
+  
+  // check true edges
+  if ( ( j>=1 && j<=k-3 ) )
+  {
+    if ( i == 0 )  {        // signs (A):  ---+
+      return ( number == 3 ) ? 1 : -1 ;
+    }
+    
+    else if ( i == n-2 ) {  // signs (C):  -+--
+      return  ( number == 1 ) ? 1 : -1 ;
+    }
+    
+    else
+    {
+      std::cout << "Not correct face number: " << GetName() << " !"<<std::endl;
+    
+    }
+  }
+  
+  if ( ( i>=1 && i<=n-3 ) )
+  {
+    if ( j == 0 )  {        // signs (D):  +---
+      return ( number == 0 ) ? 1 : -1 ;
+    }
+
+    else if ( j == k-2 ) {  // signs (B):  --+-
+      return  ( number == 2 ) ? 1 : -1 ;
+    }
+
+    else
+    {
+      std::cout << "Not correct face number: " << GetName() << " !"<<std::endl;
+    }
+  }
+  
+  // now the corners
+  if ( i == 0 && j == 0 ) {          // signs (a) : +--+
+    return ( number == 0 || number == 3 ) ? 1 : -1 ;
+  }
+  else if ( i == 0 && j == k-2 ) {   // signs (b) : --++
+    return ( number == 2 || number == 3 ) ? 1 : -1 ;
+  }
+  else if ( i == n-2 && j == k-2 ) { // signs (c) : -++-
+    return ( number == 1 || number == 2 ) ? 1 : -1 ;
+  }
+  else if ( i == n-2 && j == 0 ) {   // signs (d) : ++--
+    return ( number == 0 || number == 1 ) ? 1 : -1 ;
+  }
+  else
+  {
+      std::cout << "Not correct face number: " << GetName() << " !"<<std::endl;
+  }
+  std::cout << "Not correct face number: " << GetName() << " !"<<std::endl;
+
+  return 0 ;
+}
diff --git a/GeoModelVisualization/VP1HEPVis/src/SbTwistTrapAlphaSide.cxx b/GeoModelVisualization/VP1HEPVis/src/SbTwistTrapAlphaSide.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..a88642542d2ea1ee93f8b280955f78b4551629ac
--- /dev/null
+++ b/GeoModelVisualization/VP1HEPVis/src/SbTwistTrapAlphaSide.cxx
@@ -0,0 +1,134 @@
+//
+// G4TwistTrapAlphaSide implementation
+//
+// Author:   18/03/2005 - O.Link (Oliver.Link@cern.ch)
+// Revision: 08/Dec/2020 - Marilena Bandieramonte (marilena.bandieramonte@cern.ch):
+// Adapted from Geant4 to GMEX
+// --------------------------------------------------------------------
+
+#include <cmath>
+#include <iostream>
+
+#include "VP1HEPVis/SbTwistTrapAlphaSide.h"
+
+//=====================================================================
+//* constructors ------------------------------------------------------
+
+SbTwistTrapAlphaSide::
+SbTwistTrapAlphaSide(const std::string& name,
+                     double      PhiTwist,  // twist angle
+                     double      pDz,       // half z lenght
+                     double      pTheta,    // direction between end planes
+                     double      pPhi,      // by polar and azimutal angles
+                     double      pDy1,      // half y length at -pDz
+                     double      pDx1,      // half x length at -pDz,-pDy
+                     double      pDx2,      // half x length at -pDz,+pDy
+                     double      pDy2,      // half y length at +pDz
+                     double      pDx3,      // half x length at +pDz,-pDy
+                     double      pDx4,      // half x length at +pDz,+pDy
+                     double      pAlph,     // tilt angle at +pDz
+                     double      AngleSide  // parity
+                                             )
+  : SbTwistSurface(name)
+{
+  
+  fDx1  = pDx1 ;
+  fDx2  = pDx2 ;
+  fDx3  = pDx3 ;
+  fDx4  = pDx4 ;
+
+  fDy1   = pDy1 ;
+  fDy2   = pDy2 ;
+
+  fDz   = pDz ;
+
+  fAlph = pAlph  ;
+  fTAlph = std::tan(fAlph) ;
+
+  fTheta = pTheta ;
+  fPhi   = pPhi ;
+
+  // precalculate frequently used parameters
+  fDx4plus2  = fDx4 + fDx2 ;
+  fDx4minus2 = fDx4 - fDx2 ;
+  fDx3plus1  = fDx3 + fDx1 ;
+  fDx3minus1 = fDx3 - fDx1 ;
+  fDy2plus1  = fDy2 + fDy1 ;
+  fDy2minus1 = fDy2 - fDy1 ;
+
+  fa1md1 = 2*fDx2 - 2*fDx1  ;
+  fa2md2 = 2*fDx4 - 2*fDx3 ;
+
+  fPhiTwist = PhiTwist ;     // dphi
+  fAngleSide = AngleSide ;  // 0,90,180,270 deg
+
+  fdeltaX = 2 * fDz * std::tan(fTheta) * std::cos(fPhi);
+    // dx in surface equation
+  fdeltaY = 2 * fDz * std::tan(fTheta) * std::sin(fPhi);
+    // dy in surface equation
+    
+  GeoTrf::Transform3D rotation = GeoTrf::RotateZ3D(AngleSide);
+  GeoTrf::Vector3D translation=GeoTrf::Vector3D(0,0,0);
+  fRotTrans = GeoTrf::Translate3D(translation.x(),translation.y(),translation.z())*rotation;
+    
+}
+
+
+//=====================================================================
+//* destructor --------------------------------------------------------
+
+SbTwistTrapAlphaSide::~SbTwistTrapAlphaSide()
+{
+}
+
+//=====================================================================
+//* GetFacets ---------------------------------------------------------
+
+void
+SbTwistTrapAlphaSide::GetFacets( int k, int n, double xyz[][3],
+                                 int faces[][4], int iside )
+{
+
+  double phi ;
+  double b ;
+
+  double z, u ;     // the two parameters for the surface equation
+  GeoTrf::Vector3D  p ;  // a point on the surface, given by (z,u)
+
+  int nnode ;
+  int nface ;
+
+  // calculate the (n-1)*(k-1) vertices
+
+  for ( int i = 0 ; i<n ; ++i )
+  {
+    z = -fDz+i*(2.*fDz)/(n-1) ;
+    phi = z*fPhiTwist/(2*fDz) ;
+    b = GetValueB(phi) ;
+
+    for ( int j = 0 ; j<k ; ++j )
+    {
+      nnode = GetNode(i,j,k,n,iside) ;
+      u = -b/2 +j*b/(k-1) ;
+      p = SurfacePoint(phi,u,true) ;  // surface point in global coordinates
+
+      xyz[nnode][0] = p.x() ;
+      xyz[nnode][1] = p.y() ;
+      xyz[nnode][2] = p.z() ;
+
+      if ( i<n-1 && j<k-1 )  // conterclock wise filling
+      {
+        nface = GetFace(i,j,k,n,iside) ;
+        faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,-1)
+                        * (GetNode(i  ,j  ,k,n,iside)+1) ;  // f77 numbering
+        faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,-1)
+                        * (GetNode(i  ,j+1,k,n,iside)+1) ;
+        faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,-1)
+                        * (GetNode(i+1,j+1,k,n,iside)+1) ;
+        faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,-1)
+                        * (GetNode(i+1,j  ,k,n,iside)+1) ;
+      }
+    }
+  }
+}
+
diff --git a/GeoModelVisualization/VP1HEPVis/src/SbTwistTrapFlatSide.cxx b/GeoModelVisualization/VP1HEPVis/src/SbTwistTrapFlatSide.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..7f63c25c1fefcf71e5b41575f1b1d01b60354b99
--- /dev/null
+++ b/GeoModelVisualization/VP1HEPVis/src/SbTwistTrapFlatSide.cxx
@@ -0,0 +1,123 @@
+//
+// SbTwistTrapFlatSide implementation
+//
+// Author:   30-Aug-2002 - O.Link (Oliver.Link@cern.ch)
+// Revision: 08-Dec-2020 - Marilena Bandieramonte (marilena.bandieramonte@cern.ch):
+// Adapted from Geant4 to GMEX
+// --------------------------------------------------------------------
+
+#include "VP1HEPVis/SbTwistTrapFlatSide.h"
+#include <cmath>
+#include <iostream> 
+
+//=====================================================================
+//* constructors ------------------------------------------------------
+
+SbTwistTrapFlatSide::SbTwistTrapFlatSide( const std::string& name,
+                                                double  PhiTwist,
+                                                double  pDx1,
+                                                double  pDx2,
+                                                double  pDy,
+                                                double  pDz,
+                                                double  pAlpha,
+                                                double  pPhi,
+                                                double  pTheta,
+                                                int     handedness)
+  : SbTwistSurface(name)
+{
+   fHandedness = handedness;   // +z = +ve, -z = -ve
+
+   fDx1 = pDx1 ;
+   fDx2 = pDx2 ;
+   fDy = pDy ;
+   fDz = pDz ;
+   fAlpha = pAlpha ;
+   fTAlph = std::tan(fAlpha) ;
+   fPhi  = pPhi ;
+   fTheta = pTheta ;
+
+   fdeltaX = 2 * fDz * std::tan(fTheta) * std::cos(fPhi)  ;
+     // dx in surface equation
+   fdeltaY = 2 * fDz * std::tan(fTheta) * std::sin(fPhi)  ;
+     // dy in surface equation
+
+   fPhiTwist = PhiTwist ;
+    
+   GeoTrf::Transform3D rotation = GeoTrf::RotateZ3D(fHandedness > 0 ? 0.5 * fPhiTwist : -0.5 * fPhiTwist );
+   GeoTrf::Vector3D translation=GeoTrf::Vector3D(fHandedness > 0 ? 0.5*fdeltaX : -0.5*fdeltaX,fHandedness > 0 ? 0.5*fdeltaY : -0.5*fdeltaY,fHandedness > 0 ? fDz : -fDz);
+   fRotTrans = GeoTrf::Translate3D(translation.x(),translation.y(),translation.z())*rotation;
+    
+}
+
+//=====================================================================
+//* destructor --------------------------------------------------------
+
+SbTwistTrapFlatSide::~SbTwistTrapFlatSide()
+{
+}
+
+
+//=====================================================================
+//* GetFacets() -------------------------------------------------------
+
+void SbTwistTrapFlatSide::GetFacets( int k, int n, double xyz[][3],
+                                     int faces[][4], int iside )
+{
+  double x,y    ;     // the two parameters for the surface equation
+  GeoTrf::Vector3D  p ;  // a point on the surface, given by (z,u)
+
+  int nnode ;
+  int nface ;
+
+  double xmin,xmax ;
+
+  // calculate the (n-1)*(k-1) vertices
+
+  for ( int i = 0 ; i<n ; ++i )
+  {
+    y = -fDy + i*(2*fDy)/(n-1) ;
+
+    for ( int j = 0 ; j<k ; ++j )
+    {
+      xmin = GetBoundaryMin(y) ;
+      xmax = GetBoundaryMax(y) ;
+      x = xmin + j*(xmax-xmin)/(k-1) ;
+
+      nnode = GetNode(i,j,k,n,iside) ;
+      p = SurfacePoint(x,y,true) ;  // surface point in global coordinate system
+      
+      xyz[nnode][0] = p.x() ;
+      xyz[nnode][1] = p.y() ;
+      xyz[nnode][2] = p.z() ;
+
+      if ( i<n-1 && j<k-1 )
+      {
+        nface = GetFace(i,j,k,n,iside) ;
+
+        if (fHandedness < 0)  // lower side
+        {
+          faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,1)
+                          * ( GetNode(i  ,j  ,k,n,iside)+1) ;
+          faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,1)
+                          * ( GetNode(i+1,j  ,k,n,iside)+1) ;
+          faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,1)
+                          * ( GetNode(i+1,j+1,k,n,iside)+1) ;
+          faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,1)
+                          * ( GetNode(i  ,j+1,k,n,iside)+1) ;
+        }
+        else                  // upper side
+        {
+          faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,-1)
+                          * ( GetNode(i  ,j  ,k,n,iside)+1) ;
+          faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,-1)
+                          * ( GetNode(i  ,j+1,k,n,iside)+1) ;
+          faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,-1)
+                          * ( GetNode(i+1,j+1,k,n,iside)+1) ;
+          faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,-1)
+                          * ( GetNode(i+1,j  ,k,n,iside)+1) ;
+        }
+      }
+    }
+  }
+}
+
diff --git a/GeoModelVisualization/VP1HEPVis/src/SbTwistTrapParallelSide.cxx b/GeoModelVisualization/VP1HEPVis/src/SbTwistTrapParallelSide.cxx
new file mode 100644
index 0000000000000000000000000000000000000000..bc20ad016e99a9ebc8c4e462af6ea5f2334327e0
--- /dev/null
+++ b/GeoModelVisualization/VP1HEPVis/src/SbTwistTrapParallelSide.cxx
@@ -0,0 +1,126 @@
+//
+// SbTwistTrapParallelSide implementation
+//
+// Author: Oliver Link (Oliver.Link@cern.ch)
+// Revision: Marilena Bandieramonte (marilena.bandieramonte@cern.ch):
+// Adapted from Geant4 to GMEX
+// --------------------------------------------------------------------
+
+#include <cmath>
+#include <iostream>
+
+#include "VP1HEPVis/SbTwistTrapParallelSide.h"
+
+//=====================================================================
+//* constructors ------------------------------------------------------
+
+SbTwistTrapParallelSide::SbTwistTrapParallelSide(const std::string& name,
+                           double PhiTwist,  // twist angle
+                           double pDz,       // half z lenght
+                           double pTheta,    // direction between end planes
+                           double pPhi,      // by polar and azimutal angles
+                           double pDy1,      // half y length at -pDz
+                           double pDx1,      // half x length at -pDz,-pDy
+                           double pDx2,      // half x length at -pDz,+pDy
+                           double pDy2,      // half y length at +pDz
+                           double pDx3,      // half x length at +pDz,-pDy
+                           double pDx4,      // half x length at +pDz,+pDy
+                           double pAlph,     // tilt angle at +pDz
+                           double AngleSide  // parity
+                                               )
+  : SbTwistSurface(name)
+{
+  
+  fDx1  = pDx1 ;
+  fDx2  = pDx2 ;
+  fDx3  = pDx3 ;
+  fDx4  = pDx4 ;
+
+  fDy1   = pDy1 ;
+  fDy2   = pDy2 ;
+
+  fDz   = pDz ;
+
+  fAlph = pAlph  ;
+  fTAlph = std::tan(fAlph) ;
+
+  fTheta = pTheta ;
+  fPhi   = pPhi ;
+
+  // precalculate frequently used parameters
+  //
+  fDx4plus2  = fDx4 + fDx2 ;
+  fDx4minus2 = fDx4 - fDx2 ;
+  fDx3plus1  = fDx3 + fDx1 ;
+  fDx3minus1 = fDx3 - fDx1 ;
+  fDy2plus1  = fDy2 + fDy1 ;
+  fDy2minus1 = fDy2 - fDy1 ;
+
+  fa1md1 = 2*fDx2 - 2*fDx1  ;
+  fa2md2 = 2*fDx4 - 2*fDx3 ;
+
+  fPhiTwist = PhiTwist ;    // dphi
+  fAngleSide = AngleSide ;  // 0,90,180,270 deg
+
+  fdeltaX = 2*fDz*std::tan(fTheta)*std::cos(fPhi); // dx in surface equation
+  fdeltaY = 2*fDz*std::tan(fTheta)*std::sin(fPhi); // dy in surface equation
+
+  GeoTrf::Transform3D rotation = GeoTrf::RotateZ3D(AngleSide);
+  GeoTrf::Vector3D translation=GeoTrf::Vector3D(0,0,0);
+  fRotTrans = GeoTrf::Translate3D(translation.x(),translation.y(),translation.z())*rotation;
+}
+
+
+SbTwistTrapParallelSide::~SbTwistTrapParallelSide()
+{
+}
+
+//=====================================================================
+//* GetFacets() -------------------------------------------------------
+
+void SbTwistTrapParallelSide::GetFacets( int k, int n, double xyz[][3],
+                                         int faces[][4], int iside )
+{
+  double phi ;
+  double z, u ;     // the two parameters for the surface equation
+  GeoTrf::Vector3D  p ;  //a point on the surface, given by (z,u)
+  
+  int nnode ;
+  int nface ;
+
+  double umin, umax ;
+
+  // calculate the (n-1)*(k-1) vertices
+
+  for ( int i = 0 ; i<n ; ++i )
+  {
+    z = -fDz+i*(2.*fDz)/(n-1) ;
+    phi = z*fPhiTwist/(2*fDz) ;
+    umin = GetBoundaryMin(phi) ;
+    umax = GetBoundaryMax(phi) ;
+
+    for ( int j = 0 ; j<k ; ++j )
+    {
+      nnode = GetNode(i,j,k,n,iside) ;
+      u = umax - j*(umax-umin)/(k-1) ;
+      p = SurfacePoint(phi,u,true) ;  // surface point in global coords
+
+      xyz[nnode][0] = p.x() ;
+      xyz[nnode][1] = p.y() ;
+      xyz[nnode][2] = p.z() ;
+
+      if ( i<n-1 && j<k-1 )    // conterclock wise filling
+      {
+        nface = GetFace(i,j,k,n,iside) ;
+        faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,-1)
+                        * (GetNode(i  ,j  ,k,n,iside)+1) ; // fortran numbering
+        faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,-1)
+                        * (GetNode(i  ,j+1,k,n,iside)+1) ;
+        faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,-1)
+                        * (GetNode(i+1,j+1,k,n,iside)+1) ;
+        faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,-1)
+                        * (GetNode(i+1,j  ,k,n,iside)+1) ;
+      }
+    }
+  }
+}
diff --git a/documentation/.gitignore b/documentation/.gitignore
index 5bc1a7f959be987ea485686199204789bb4e8c79..0b7f4b423869a7555f41bcb79571c98fec6e222b 100755
--- a/documentation/.gitignore
+++ b/documentation/.gitignore
@@ -1,5 +1,8 @@
 # MkDocs output
 _build
 
+# Python temp files
+__pycache__
+
 # Python3 virtual env for MkDocs packages
 venvp3
diff --git a/documentation/create_activate_virtualenv_p3.sh b/documentation/create_activate_virtualenv_p3.sh
index 98bdc8c521544efdf97bc34bc48913c84cccf098..2b30df617fd4ad9a3314cd3b72399410a3cdbee8 100644
--- a/documentation/create_activate_virtualenv_p3.sh
+++ b/documentation/create_activate_virtualenv_p3.sh
@@ -22,7 +22,7 @@ pip install mkdocs-mermaid2-plugin # adds support for mermaid diagrams / flowcha
 
 
 
-echo "\n\n\nSuccess!!\nYou can now build and serve the web iste locally by running: 'brew serve'"
+echo "\n\n\nSuccess!!\nYou can now build and serve the web iste locally by running: 'mkdocs serve'"
 echo "\n[Note: if you have mkdocs installed on your system folders already, to avoid errors about missing packages, you might have to run the one in the virtualenv explicitely, by running: './venvp3/bin/mkdocs serve']"
 
 echo "\nThen, when you finished working with mkdocs, you can close the virtual environment by running: 'deactivate'\n\n"
diff --git a/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoBox.md b/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoBox.md
index 51c6a891395c1ee4489159d7504dd0339146dadb..cf8fea79a2188c02dd94f61e99af2c0fb5f121bc 100644
--- a/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoBox.md
+++ b/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoBox.md
@@ -16,10 +16,12 @@
 The constructor fills the box with the $x$, $y$, and $z$ *half-lengths*, and the accessors return those quantities.
 
 
-<figure>
-  <img src="/components/kernel/reference/RCBase/GeoShape/GeoBox.png" width="400" />
-  <figcaption>Figure 2: A GeoBox object, representing a rectangular prism or “box”.</figcaption>
-</figure>
+{{ imgutils_image_caption('RCBase/GeoShape/GeoBox.png', 
+   alt='The GeoBox shape', 
+   cap='Figure 2: A GeoBox object, representing a rectangular prism or “box”.',
+   urlFix=False) 
+}}
+
 
 
 
diff --git a/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoCons.md b/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoCons.md
index c8834291b72fc639b8c8b0f3cb561e1322b3432e..0d59b1ab6d6743cac70852577793ab37b26b439b 100644
--- a/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoCons.md
+++ b/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoCons.md
@@ -1,5 +1,6 @@
 ### GeoCons
 
+
 ```cpp
 // === GeoCons ===
 
@@ -20,9 +21,12 @@
 
 A `GeoCons` represents a cone section positioned with its axis along the $z$ direction, and is specified by a starting value of $\phi$, a total subtended angle in $\phi$, a *half-width* in $z$, and minimum and maximum values of radius at both extremities in $z$.  The constructor specifies the values of these constants, and the accessors return them.
 
-<figure>
-  <img src="/components/kernel/reference/RCBase/GeoShape/GeoCons.png" width="400" />
-  <figcaption>Figure 3:  A GeoCons Object, representing a cone section.</figcaption>
-</figure>
+
+{{ imgutils_image_caption('RCBase/GeoShape/GeoCons.png', 
+   alt='The GeoCons shape', 
+   cap='Figure: A GeoCons Object, representing a cone section.',
+   urlFix=False) 
+}}
+
 
 
diff --git a/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoPara.md b/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoPara.md
index 8406fbc7a83fcea71bdc33a7780e8b0fb663c819..8ad30b46cc75226906ed45751d41ac261c372025 100644
--- a/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoPara.md
+++ b/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoPara.md
@@ -19,11 +19,18 @@
 
 The `GeoPara` class represents a parallelepiped. Faces at $\pm z$ are parallel to the $x-y$ plane. One edge of each of these faces is parallel to the $x$-axis, while the other edge makes an angle $\alpha$ with respect to the $y$-axis.  The remaining edge of the parallelepiped is oriented along a vector whose polar angle is $\theta$ and whose azimuthal angle is $\phi$.  *Half-lengths* in $x$, $y$, and $z$ describe the projections of the sides of the parallelepiped project onto the coordinate axes.  The constructor fills these data, while the accessors return them.
 
-*Note:  Visualization of GeoPara is on the to-do list.*
 
-<figure>
-  <img src="https://dummyimage.com/600x400/eee/aaa" width="400" />
-  <figcaption>Figure 4: GeoPara object, representing a parallelepiped.</figcaption>
-</figure>
+
+!!! warning
+
+    Visualization of GeoPgon is on the to-do list.
+
+
+
+{{ imgutils_image_caption('https://dummyimage.com/600x400/eee/aaa', 
+   alt='The GeoPara shape', 
+   cap='Figure: GeoPara object, representing a parallelepiped.',
+   urlFix=False) 
+}}
 
 
diff --git a/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoPcon.md b/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoPcon.md
index 31264dda200ac6f3aa193286d70e639500f98b78..7d2810c1a0a2d4c64006bc28289d990ca3561cbf 100644
--- a/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoPcon.md
+++ b/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoPcon.md
@@ -22,14 +22,17 @@
 
 When the polycone is constructed, only $\phi_0$ and $\phi$ are given; then, the information at each $z$ location is given, one plane at a time, by using the `addPlane()` method.  At least two planes must be given, otherwise the polycone is not valid and methods such as `volume()` will fail and throw an exception.  The `isValid()` method can be used to determine whether the polycone has at least two planes.
 
+
 !!! note
 
     A polycone (`GeoPcon`) with exactly two planes is equivalent geometrically to a cone section (`GeoCons`).
 
 
-<figure>
-  <img src="/components/kernel/reference/RCBase/GeoShape/GeoPcon.png" width="400" />
-  <figcaption>Figure 2: A GeoPcon object, representing a "polycone", which is a union of cone sections.</figcaption>
-</figure>
+
+{{ imgutils_image_caption('RCBase/GeoShape/GeoPcon.png', 
+   alt='The GeoPcon shape', 
+   cap='Figure: A GeoPcon object, representing a "polycone", which is a union of cone sections.',
+   urlFix=False) 
+}}
 
 
diff --git a/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoPgon.md b/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoPgon.md
index b6bcc92c9044c252ba3c53e073f8393e82257de7..2813c4b782a4f478eec8f849c22a4dda14eceff2 100644
--- a/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoPgon.md
+++ b/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoPgon.md
@@ -23,12 +23,19 @@
 
 The constructor takes $\phi$ , $\phi_0$, and the number of sides in the cross-section as arguments; then, the information at each $z$ location is given, one plane at a time, by using the `addPlane()` method.  At least two planes must be given, otherwise the polygon is not valid and methods such as `volume()` will fail and throw an exception.  The `isValid()` method can be used to determine whether the polygon has at least two planes.
 
-*Note:  Visualization of GeoPgon is on the to-do list.*
 
-<figure>
-  <img src="https://dummyimage.com/600x400/eee/aaa" width="400" />
-  <figcaption>Figure 4: GeoPgon object, representing a polycone with a polygonal cross section.</figcaption>
-</figure>
+
+!!! warning
+
+    Visualization of GeoPgon is on the to-do list.
+
+
+{{ imgutils_image_caption('https://dummyimage.com/600x400/eee/aaa', 
+   alt='The GeoPgon shape', 
+   width='400', 
+   cap='Figure: GeoPgon object, representing a polycone with a polygonal cross section.',
+   urlFix=False) 
+}}
 
 
 
diff --git a/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoTrap.md b/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoTrap.md
index 791e5fac021e6b1ff9cf1d20fde223c345a15dad..ffb02570506840f55f50b02b0a0012c44e3b2181 100644
--- a/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoTrap.md
+++ b/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoTrap.md
@@ -30,7 +30,13 @@ The two edges not parallel to the $x$-axis make an angle of $\alpha_n$ (`Angleyd
 
 The constructor fills the `GeoTrap` with these values and the accessors return them.
 
-<figure>
-  <img src="/components/kernel/reference/RCBase/GeoShape/GeoTrap.png" width="400" />
-  <figcaption>Figure 2: A GeoTrap object, representing a very general kind of trapezoid.</figcaption>
-</figure>
+
+{{ imgutils_image_caption('RCBase/GeoShape/GeoTrap.png', 
+   alt='The GeoTrap shape', 
+   cap='Figure: A GeoTrap object, representing a very general kind of trapezoid.',
+   urlFix=False) 
+}}
+
+
+
+
diff --git a/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoTrd.md b/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoTrd.md
index 9756a9fc3d41f7889afb8f60afde7592649e16e9..ce4caab3e7021e25346a461ba1c75edd5a310b92 100644
--- a/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoTrd.md
+++ b/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoTrd.md
@@ -22,9 +22,11 @@ A `GeoTrd` is a simple trapezoid.  Two faces at $\pm \Delta z$ are parallel to e
 The constructor fills the object with these values and the accessors return them.
 
 
-<figure>
-  <img src="/components/kernel/reference/RCBase/GeoShape/GeoTrd.png" width="400" />
-  <figcaption>Figure 8: A GeoTrd object, representing a simple trapezoid.</figcaption>
-</figure>
+
+{{ imgutils_image_caption('RCBase/GeoShape/GeoTrd.png', 
+   alt='The GeoTrd shape', 
+   cap='Figure: A GeoTrd object, representing a simple trapezoid.',
+   urlFix=False) 
+}}
 
 
diff --git a/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoTube.md b/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoTube.md
index aac89285a54f3725ac0676e0dc64efc809bb2ecf..a3968604ba4bbaa1efc46ce618be7c6ff595c7c4 100644
--- a/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoTube.md
+++ b/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoTube.md
@@ -18,10 +18,13 @@ The `GeoTube` class represents a tube, specified by an inner radius (`RMin`), an
 The constructor fills these quantities and the accessors return them.
 
 
-<figure>
-  <img src="/components/kernel/reference/RCBase/GeoShape/GeoTube.png" width="400" />
-  <figcaption>Figure 2: A GeoTube object, representing a tube or "pipe".</figcaption>
-</figure>
+
+{{ imgutils_image_caption('RCBase/GeoShape/GeoTube.png', 
+   alt='The GeoTube shape', 
+   cap='Figure N: A GeoTube object, representing a tube or "pipe".',
+   urlFix=False) 
+}}
+
 
 
 
diff --git a/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoTubs.md b/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoTubs.md
index 7a00f53b8a299abb89668232aff1b7ec8577560c..a71eb394dd59b8347bae22d60cc0ad983f9254a2 100644
--- a/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoTubs.md
+++ b/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoTubs.md
@@ -20,8 +20,13 @@ A `GeoTubs` is a tube section; a tube that subtends some plane angle (less than
 
 Member functions provide access to these quantities.
 
-<figure>
-  <img src="/components/kernel/reference/RCBase/GeoShape/GeoTubs.png" width="400" />
-  <figcaption>Figure 7: A GeoTubs object, representing a tube section.</figcaption>
-</figure>
+
+{{ imgutils_image_caption('RCBase/GeoShape/GeoTubs.png', 
+   alt='The GeoTubs shape', 
+   cap='Figure: A GeoTubs object, representing a tube section.',
+   urlFix=False) 
+}}
+
+
+
 
diff --git a/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoTwistedTrap.md b/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoTwistedTrap.md
new file mode 100644
index 0000000000000000000000000000000000000000..bcba627d7812a083341c698be596a2b4e2f7114c
--- /dev/null
+++ b/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoTwistedTrap.md
@@ -0,0 +1,44 @@
+
+### GeoTwistedTrap
+
+```cpp
+//=== GeoTwistedTrap ===
+
+  // Constructor:
+  GeoTwistedTrap(double  PhiTwist,   // twist angle
+                 double  Dz,     // half z length
+                 double  Theta,  // direction between end planes
+                 double  Phi,    // defined by polar and azim. angles
+                 double  Dy1,    // half y length at -pDz
+                 double  Dx1,    // half x length at -pDz,-pDy
+                 double  Dx2,    // half x length at -pDz,+pDy
+                 double  Dy2,    // half y length at +pDz
+                 double  Dx3,    // half x length at +pDz,-pDy
+                 double  Dx4,    // half x length at +pDz,+pDy
+                 double  Alph    // tilt angle
+                 )
+
+  // Public Methods:
+ double getY1HalfLength() const 
+ double getX1HalfLength() const 
+ double getX2HalfLength() const 
+ double getY2HalfLength() const 
+ double getX3HalfLength() const 
+ double getX4HalfLength() const 
+ double getZHalfLength()  const 
+ double getPhiTwist()     const 
+ double getTheta()        const 
+ double getPhi()          const 
+ double getTiltAngleAlpha() const 
+```
+
+The `GeoTwistedTrap` class represents a twisted trapezoid.  Two faces at $\pm \Delta z$ (`Dz`)  are parallel to each other and to the $x-y$ plane.  The centers of the faces are offset by a vector whose polar and azimuthal angles respectively are $\theta$ and $\phi$.  At $-\Delta z$, two edges parallel to the $x$-axis are offset by $\pm \Delta y_1$ (`Dy1`) from the face’s center, and these two faces have *half-lengths* of $\Delta x_{1}$ (`Dx1`) and $\Delta x_{2}$ (`Dx2`). The face at $+\Delta z$ are similar:  two edges parallel to the $x$-axis are offset by $\pm \Delta y_2$ (`Dy2`) from the face’s center, and these two faces have *half-lengths* of $\Delta x_{3}$ (`Dx3`) and $\Delta x_{4}$ (`Dx4`).
+
+The face at $+ \Delta z$ is twisted with respect to the one at $- \Delta z$  of an angle defined by $\phi Twist$. The tilt angle $\alpha$ defines the angle with respect to the y axis from the centre of the top and bottom trapezoids. It transforms the two trapezoids from isosceles to scalene. 
+
+The constructor fills the `GeoTwistedTrap` with these values and the accessors return them.
+
+<figure>
+  <img src="/reference/RCBase/GeoShape/GeoTwistedTrap.png" width="400" />
+  <figcaption>Figure 2: A GeoTwistedTrap object, representing a twisted trapezoid.</figcaption>
+</figure>
diff --git a/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoTwistedTrap.png b/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoTwistedTrap.png
new file mode 100644
index 0000000000000000000000000000000000000000..c99b7a5022e19342d57914d32eabc943364237e4
Binary files /dev/null and b/documentation/docs/components/kernel/reference/RCBase/GeoShape/GeoTwistedTrap.png differ
diff --git a/documentation/docs/components/kernel/reference/Shapes.md b/documentation/docs/components/kernel/reference/Shapes.md
index a42d3b0d07cb75171706d005b17f4e91eb245e90..e3c1ae32f1921dd76ad468e0746179192258d501 100644
--- a/documentation/docs/components/kernel/reference/Shapes.md
+++ b/documentation/docs/components/kernel/reference/Shapes.md
@@ -2,9 +2,10 @@
 
 ##	Shapes
 
+
 ### Introduction
 
-The shape classes in the geometry kernel are data structures designed to describe several geometrical primitives. Table 1 describes the different shapes presently provided within the geometry kernel.  This set is extensible; one only needs to derive a new shape from the base class and insure that it fits the pattern described below. Shapes are reference- counted objects, as described in [1.3.2.6][#].
+The shape classes in the geometry kernel are data structures designed to describe several geometrical primitives. Table 1 describes the different shapes presently provided within the geometry kernel.  This set is extensible; one only needs to derive a new shape from the base class and insure that it fits the pattern described below. Shapes are reference-counted objects, as described in [Reference Counting](#reference-counting).
 
 
 | Class   | Shape |
@@ -86,23 +87,23 @@ The classes `GeoShapeShift`, `GeoShapeUnion`, `GeoShapeSubtraction`, and `GeoSha
 
 We now present the interfaces to specific shapes.  In general these shapes are by default constructed as symmetrically around the origin as possible.  
 
-{!components/kernel/reference/RCBase/GeoShape/GeoBox.md!}
+{% include 'components/kernel/reference/RCBase/GeoShape/GeoBox.md' %}
 
-{!components/kernel/reference/RCBase/GeoShape/GeoCons.md!}
+{% include 'components/kernel/reference/RCBase/GeoShape/GeoCons.md' %}
 
-{!components/kernel/reference/RCBase/GeoShape/GeoPara.md!}
+{% include 'components/kernel/reference/RCBase/GeoShape/GeoPara.md' %}
 
-{!components/kernel/reference/RCBase/GeoShape/GeoPcon.md!}
+{% include 'components/kernel/reference/RCBase/GeoShape/GeoPcon.md' %}
 
-{!components/kernel/reference/RCBase/GeoShape/GeoPgon.md!}
+{% include 'components/kernel/reference/RCBase/GeoShape/GeoPgon.md' %}
 
-{!components/kernel/reference/RCBase/GeoShape/GeoTrap.md!}
+{% include 'components/kernel/reference/RCBase/GeoShape/GeoTrap.md' %}
 
-{!components/kernel/reference/RCBase/GeoShape/GeoTube.md!}
+{% include 'components/kernel/reference/RCBase/GeoShape/GeoTube.md' %}
 
-{!components/kernel/reference/RCBase/GeoShape/GeoTubs.md!}
+{% include 'components/kernel/reference/RCBase/GeoShape/GeoTubs.md' %}
 
-{!components/kernel/reference/RCBase/GeoShape/GeoTrd.md!}
+{% include 'components/kernel/reference/RCBase/GeoShape/GeoTrd.md' %}
 
 
 
diff --git a/documentation/docs/components/kernel/reference/index.md b/documentation/docs/components/kernel/reference/index.md
index 65add170724dfa8d5f37d0041f52577cdd82a6e7..8b8638701fadfba01caf388de80080ba2bdfa6aa 100644
--- a/documentation/docs/components/kernel/reference/index.md
+++ b/documentation/docs/components/kernel/reference/index.md
@@ -1,24 +1,25 @@
 # GeoModel Kernel Class Reference
 
+
 This section describes in more detail the classes in the geometry kernel. In most cases we provide the class interface. In cases where part of the interface is used only internally by the geometry kernel itself and not by other users. In such cases we present a simplified picture of the interfaces.
 
 Detailed descriptions of the geometry kernel classes follow.
 
 
 
-{!components/kernel/reference/ReferenceCounting.md!}
+{% include 'components/kernel/reference/ReferenceCounting.md' %}
 
-{!components/kernel/reference/Materials_Elements.md!}
+{% include 'components/kernel/reference/Materials_Elements.md' %}
 
-{!components/kernel/reference/Shapes.md!}
+{% include 'components/kernel/reference/Shapes.md' %} 
 
-{!components/kernel/reference/LogicalVolumes.md!}
+{% include 'components/kernel/reference/LogicalVolumes.md' %}
 
-{!components/kernel/reference/PhysicalVolumes.md!}
+{% include 'components/kernel/reference/PhysicalVolumes.md' %}
 
-{!components/kernel/reference/Transformations.md!}
+{% include 'components/kernel/reference/Transformations.md' %}
 
-{!components/kernel/reference/NameTags.md!}
+{% include 'components/kernel/reference/NameTags.md' %}
 
 
 
diff --git a/documentation/docs/index.md b/documentation/docs/index.md
index d53f9f6673569f281b1a009b25a004369533434c..167cff67b423618f63b42bbc1223dc0015949fe8 100755
--- a/documentation/docs/index.md
+++ b/documentation/docs/index.md
@@ -15,4 +15,12 @@ The external dependencies are minimal:
 
 ----
 
-*This Documentation is Work In Progress*
+<i>Last update:</i> 
+{% if git.status %}
+  <i> {{ git.date.strftime("%b %d, %Y -- %H:%M") or now() }} </i>
+{% else %}
+  <i> {{ now() }} </i> 
+{% endif %}
+
+
+
diff --git a/documentation/imgutils.py b/documentation/imgutils.py
new file mode 100644
index 0000000000000000000000000000000000000000..e53637daa46f87863eb3c56fc86a37a91d607470
--- /dev/null
+++ b/documentation/imgutils.py
@@ -0,0 +1,56 @@
+"""
+
+'imgutils'
+
+Utility macros to handle images, captions, and related things.
+
+Author: Riccardo Maria BIANCHI <riccardo.maria.bianchi@cern.ch>
+
+12 Jan 2021
+"""
+
+import math
+
+def define_env(env):
+    """
+    Utility macros to handle images, captions, and related things.
+    """
+
+    # import the predefined macro
+    fix_url = env.variables.fix_url # make relative urls point to root
+    
+    # activate trace
+    chatter = env.start_chatting("imgutils")
+
+
+    @env.macro
+    def imgutils_image(url:str, alt:str=''):
+        "Add an image within <figure> tags. Optionally, you can set an alternate-text" 
+        url2:str = url
+        if urlFix : 
+            url2 = fix_url(url)
+            chatter("Fixing image's URL: %s --> %s" % (url, url2))
+        return '<img src="%s", alt="%s">' % (url2, alt)
+
+    @env.macro
+    def imgutils_image_caption(url:str, alt:str='', width:int=400, cap:str='', urlFix:bool=True):
+        "Add an image within <figure> tags. Optionally, you can set a caption, an alternate-text, and a width value." 
+        url2:str = url
+        if urlFix : 
+            url2:str = fix_url(url)
+            chatter("Fixing image's URL: %s --> %s" % (url, url2))
+        return '<figure> <img src="%s", alt="%s", width="%s"> <figcaption>%s<figcaption> </figure>' % ( url2, alt, width, cap)
+
+
+    @env.macro
+    def imgutils_doc_env():
+        """
+        Document the environment, for development/debugging. 
+        (Hint: Use it in your .md file within a code block, i.e. as: 
+           ```
+           {{ doc_env() | pprint }}
+           ```
+        """
+        return {name:getattr(env, name) for name in dir(env) if not name.startswith('_')}
+
+
diff --git a/documentation/mkdocs.yml b/documentation/mkdocs.yml
index 3b15eaadbc6c760ea2e45cadcaf46b4aa31021e3..2d3ca16c7a6cbd695a794ddcd8250d37d939cc6b 100755
--- a/documentation/mkdocs.yml
+++ b/documentation/mkdocs.yml
@@ -4,6 +4,7 @@ site_author: ATLAS Collaboration & Contributors
 site_url: https://cern.ch/geomodel/home
 repo_name: GeoModel
 repo_url: https://gitlab.cern.ch/GeoModelDev/GeoModel
+site_author: geomodel-developers@cern.ch
 
 theme:
     name: material
@@ -29,6 +30,7 @@ nav:
         - Kernel:
             - 'Overview': 'components/kernel/overview/index.md'
             - 'Class Reference': 'components/kernel/reference/index.md'
+              #- 'Box': 'components/kernel/reference/RCBase/GeoShape/GeoBox.md'
     - For Developers:
         - 'Build the Libraries and Tools': 'dev/index.md'
         - 'Install single packages': 'dev/install-single-packages.md'
@@ -60,9 +62,11 @@ markdown_extensions:
 
 plugins:
   - search
-  - macros
+  - macros:
+        module_name: imgutils
+        verbose: true # enables debug through the 'chatter' method
   - mermaid2 # for diagrams
-  #- git-revision-date #TODO: teh mkdocs docker image apparently lacks the git exec, so this cannot be used. Update the image!
+  #- git-revision-date #TODO: the mkdocs docker image apparently lacks the git exec, so this cannot be used. Update the image!
 
 extra_javascript:
   - javascripts/config.js