diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoAnnulusSurface.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoAnnulusSurface.h
index 2f85c16c96b264877eb20bb1bde33f2bdb2514e9..7dc31d3dc5e8e89b4be556755f90111b3d2f2a67 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoAnnulusSurface.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoAnnulusSurface.h
@@ -76,6 +76,7 @@ class GeoAnnulusSurface : public GeoVSurfaceShape
 
       void exec (GeoShapeAction *action) const override final;
 
+      virtual bool isOnSurface (const double Px, const double Py, const double Pz, const GeoTrf::Transform3D & trans) const override final;
     
     protected:
 
@@ -90,4 +91,4 @@ class GeoAnnulusSurface : public GeoVSurfaceShape
       double m_phi{2*M_PI};
 };
 
-#endif
\ No newline at end of file
+#endif
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoDiamondSurface.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoDiamondSurface.h
index a44801b50bed62ec1a8c5c2b941cef547d0904bb..bcf1b62cd6a63083768e28610e05b4aa76b9d9d3 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoDiamondSurface.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoDiamondSurface.h
@@ -65,6 +65,10 @@ class GeoDiamondSurface : public GeoVSurfaceShape
   //     Executes a GeoShapeAction
   void exec (GeoShapeAction *action) const override final;
 
+  //     Is the point (x,y,z) inside the shape?
+  virtual bool isOnSurface (const double Px, const double Py, const double Pz, const GeoTrf::Transform3D & trans) const override final;
+
+
  protected:
   virtual ~GeoDiamondSurface() = default;
   
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoRectSurface.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoRectSurface.h
index efc8e655b76836edbbdf029f94265981d5a154bc..c5a9935c961c2f8d412c0e43279f18995ba9d1ca 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoRectSurface.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoRectSurface.h
@@ -47,6 +47,9 @@ class GeoRectSurface : public GeoVSurfaceShape
   //     Executes a GeoShapeAction
   void exec (GeoShapeAction *action) const override final;
 
+  //     Is the point (x,y,z) inside the shape?
+  virtual bool isOnSurface (const double Px, const double Py, const double Pz, const GeoTrf::Transform3D & trans) const override final;
+
  protected:
   virtual ~GeoRectSurface() = default;
   
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTrapezoidSurface.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTrapezoidSurface.h
index 420966f1b6d4cd2efbff3ed141c17938321bf2b8..91a16c8dc2804c2fb2d84d0e07368693872999bc 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTrapezoidSurface.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoTrapezoidSurface.h
@@ -53,6 +53,9 @@ class GeoTrapezoidSurface : public GeoVSurfaceShape
   //     Executes a GeoShapeAction
   void exec (GeoShapeAction *action) const override final;
 
+  //     Is the point (x,y,z) inside the shape?
+  virtual bool isOnSurface (const double Px, const double Py, const double Pz, const GeoTrf::Transform3D & trans) const override final;
+
  protected:
   virtual ~GeoTrapezoidSurface() = default;
   
diff --git a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoVSurfaceShape.h b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoVSurfaceShape.h
index a378f9e99c537328b4ae8141e164bb430cdd4d5f..720047192946ca3374854da436f0fe5d5f36f67b 100644
--- a/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoVSurfaceShape.h
+++ b/GeoModelCore/GeoModelKernel/GeoModelKernel/GeoVSurfaceShape.h
@@ -43,6 +43,9 @@ class GeoVSurfaceShape : public RCBase
   //    TODO: Introduce GeoVSurfaceShapeAction in future
   virtual void exec (GeoShapeAction *action) const = 0;
 
+  //    Is the point (x,y,z) inside the shape?
+  virtual bool isOnSurface (const double Px, const double Py, const double Pz, const GeoTrf::Transform3D & trans) const = 0;
+
  protected:
   virtual ~GeoVSurfaceShape() = default;
 };
diff --git a/GeoModelCore/GeoModelKernel/src/GeoAnnulusSurface.cxx b/GeoModelCore/GeoModelKernel/src/GeoAnnulusSurface.cxx
index d230629d9feeed59b72aec393924445c888c107b..802a3cb071fc8b19fd9bcf6fff61144357cc8200 100644
--- a/GeoModelCore/GeoModelKernel/src/GeoAnnulusSurface.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoAnnulusSurface.cxx
@@ -36,6 +36,62 @@ void GeoAnnulusSurface::exec(GeoShapeAction *action) const{
     action->handleAnnulusSurface(this);
 }
 
-// bool GeoAnnulusSurface::isInside(double x, double y){}
+bool GeoAnnulusSurface::isOnSurface(const double Px, const double Py, const double Pz, const GeoTrf::Transform3D & trans) const{
+    GeoTrf::Vector3D P_prime(Px, Py, Pz);    
+    Eigen::Vector4d P_prime_4d(Px, Py, Pz, 1.0);
+    P_prime_4d = trans.inverse() * P_prime_4d;
+    double Pp_x = P_prime_4d[0];
+    double Pp_y = P_prime_4d[1];
+    double Pp_z = P_prime_4d[2];
+
+    if(Pp_z != 0 && Pp_z > 1e-5 && Pp_z < -1e-5){
+        // now I take tolerance as 1e-5
+        return false;
+    }
+    
+    // The Annulus Shape starts from theta = 0
+    double real_theta = this -> getPhi();
+    int quotient = floor(real_theta/(2.0*M_PI));
+    double theta = real_theta - quotient*2.0*M_PI;
+
+    if (theta == 0.0 && real_theta > 0.0){
+        theta = 2.0*M_PI;
+    }
+    // std::cout << "theta: " << theta << std::endl;
+
+    double radius_in = this -> getRadiusIn();
+    double radius_out = this -> getRadiusOut();
+    double Ox = this -> getOx();
+    double Oy = this -> getOy();
+
+    // Make sure the Point is located in the general Annulus Ring
+    if (Pp_x * Pp_x + Pp_y * Pp_y - radius_in * radius_in < -1e-5 || Pp_x * Pp_x + Pp_y * Pp_y - radius_out * radius_out > 1e-5){
+        // std::cout << " RADIUS " << Pp_x * Pp_x + Pp_y * Pp_y <<  std::endl;
+        // std::cout << "FLAG 2" << std::endl;
+        return false;
+    }
+    // else if (Pp_x * Pp_x + Pp_y * Pp_y < radius_in * radius_in || Pp_x * Pp_x + Pp_y * Pp_y > radius_out * radius_out){
+    //   std::cout << " RADIUS QUESTIONED" << Pp_x * Pp_x + Pp_y * Pp_y <<  std::endl;
+    // }
+
+    // Make sure the Point is located in the Deviated Circular Sector
+    double Dx = Pp_x - Ox; // X direction distance from the deviation center to the point
+    double Dy = Pp_y - Oy; // Y direction distance from the deviation center to the point
+    // std::cout << "Dx: " << Dx << " Dy: " << Dy << std::endl;
+
+    double angle = atan2(Dy, Dx); // the range of atan2 is (-pi, pi]
+    if (angle < 0.0){
+        angle += 2*M_PI;
+    }
+    // std::cout << "angle: " << angle << " theta: " << theta << std::endl;
+    if (angle <= theta){
+        return true;
+    }
+    else{
+        // std::cout << "FLAG 3" << std::endl;
+        return false;
+    }
+
+}
 
 // double GeoAnnulusSurface::area() const{}
diff --git a/GeoModelCore/GeoModelKernel/src/GeoDiamondSurface.cxx b/GeoModelCore/GeoModelKernel/src/GeoDiamondSurface.cxx
index 4e8e221eb24dfa50015d2e6a1c9784dbc3d5031c..e4d04b77c6930e3c0b191ecf3be9d3092a407c1a 100644
--- a/GeoModelCore/GeoModelKernel/src/GeoDiamondSurface.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoDiamondSurface.cxx
@@ -19,3 +19,33 @@ void GeoDiamondSurface::exec (GeoShapeAction *action) const {
   action->handleDiamondSurface(this);
 }
 
+bool GeoDiamondSurface::isOnSurface (const double Px, const double Py, const double Pz, const GeoTrf::Transform3D & trans) const{
+    GeoTrf::Vector3D P_prime(Px, Py, Pz);    
+    Eigen::Vector4d P_prime_4d(Px, Py, Pz, 1.0);
+    P_prime_4d = trans.inverse() * P_prime_4d;
+    double Pp_x = P_prime_4d[0];
+    double Pp_y = P_prime_4d[1];
+    double Pp_z = P_prime_4d[2];
+
+    if(Pp_z != 0 && Pp_z > 1e-5 && Pp_z < -1e-5){
+        // now I take tolerance as 1e-5
+        return false;
+    }
+
+    double x_bot = this -> getXbottomHalf(); double y_bot = this -> getYbottomHalf();
+    double x_mid = this -> getXmidHalf();
+    double x_top = this -> getXtopHalf(); double y_top = this -> getYtopHalf();
+    double p1x = x_bot; double p1y = -y_bot;
+    double p2x = x_mid; double p2y = 0.0;
+    double p3x = x_top; double p3y = y_top;
+    double p4x = -x_top; double p4y = y_top;
+    double p5x = -x_mid; double p5y = 0.0;
+    double p6x = -x_bot; double p6y = -y_bot;
+    if( (p2x-p1x)*(Pp_y-p1y) - (p2y-p1y)*(Pp_x-p1x) < -1e-5 ) return false;
+    if( (p3x-p2x)*(Pp_y-p2y) - (p3y-p2y)*(Pp_x-p2x) < -1e-5 ) return false;
+    if( (p4x-p3x)*(Pp_y-p3y) - (p4y-p3y)*(Pp_x-p3x) < -1e-5 ) return false;
+    if( (p5x-p4x)*(Pp_y-p4y) - (p5y-p4y)*(Pp_x-p4x) < -1e-5 ) return false;
+    if( (p6x-p5x)*(Pp_y-p5y) - (p6y-p5y)*(Pp_x-p5x) < -1e-5 ) return false;
+    if( (p1x-p6x)*(Pp_y-p6y) - (p1y-p6y)*(Pp_x-p6x) < -1e-5 ) return false;
+    return true;
+}
diff --git a/GeoModelCore/GeoModelKernel/src/GeoRectSurface.cxx b/GeoModelCore/GeoModelKernel/src/GeoRectSurface.cxx
index af5b45cabeacf96268a13f50dedc396c67573dc9..c999cb2ef5c0e3e60e0acf8c075de8f99b3a97f2 100644
--- a/GeoModelCore/GeoModelKernel/src/GeoRectSurface.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoRectSurface.cxx
@@ -16,3 +16,28 @@ void GeoRectSurface::exec (GeoShapeAction *action) const {
   action->handleRectSurface(this);
 }
 
+bool GeoRectSurface::isOnSurface (const double Px, const double Py, const double Pz, const GeoTrf::Transform3D & trans) const{
+  GeoTrf::Vector3D P_prime(Px, Py, Pz);    
+  Eigen::Vector4d P_prime_4d(Px, Py, Pz, 1.0);
+  P_prime_4d = trans.inverse() * P_prime_4d;
+  double Pp_x = P_prime_4d[0];
+  double Pp_y = P_prime_4d[1];
+  double Pp_z = P_prime_4d[2];
+
+  if(Pp_z != 0 && Pp_z > 1e-5 && Pp_z < -1e-5){
+    // now I take tolerance as 1e-5
+    return false;
+  }
+
+  double half_x = this -> getXHalfLength();
+  double half_y = this -> getYHalfLength();
+  double p1x = half_x; double p1y = -half_y;
+  double p2x = half_x; double p2y = half_y;
+  double p3x = -half_x; double p3y = half_y;
+  double p4x = -half_x; double p4y = -half_y;
+  if( (p2x-p1x)*(Pp_y-p1y) - (p2y-p1y)*(Pp_x-p1x) < -1e-5 ) return false;
+  if( (p3x-p2x)*(Pp_y-p2y) - (p3y-p2y)*(Pp_x-p2x) < -1e-5 ) return false;
+  if( (p4x-p3x)*(Pp_y-p3y) - (p4y-p3y)*(Pp_x-p3x) < -1e-5 ) return false;
+  if( (p1x-p4x)*(Pp_y-p4y) - (p1y-p4y)*(Pp_x-p4x) < -1e-5 ) return false;
+  return true;  
+}
diff --git a/GeoModelCore/GeoModelKernel/src/GeoTrapezoidSurface.cxx b/GeoModelCore/GeoModelKernel/src/GeoTrapezoidSurface.cxx
index 765d70a89d2e7b9b38cde09ebcdfafa9128e1810..0ce90f5bc7a2f7d8a82256a15e8b9d9d81098398 100644
--- a/GeoModelCore/GeoModelKernel/src/GeoTrapezoidSurface.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoTrapezoidSurface.cxx
@@ -17,3 +17,29 @@ void GeoTrapezoidSurface::exec (GeoShapeAction *action) const {
   action->handleTrapezoidSurface(this);
 }
 
+bool GeoTrapezoidSurface::isOnSurface (const double Px, const double Py, const double Pz, const GeoTrf::Transform3D & trans) const{
+    GeoTrf::Vector3D P_prime(Px, Py, Pz);    
+    Eigen::Vector4d P_prime_4d(Px, Py, Pz, 1.0);
+    P_prime_4d = trans.inverse() * P_prime_4d;
+    double Pp_x = P_prime_4d[0];
+    double Pp_y = P_prime_4d[1];
+    double Pp_z = P_prime_4d[2];
+
+    if(Pp_z != 0 && Pp_z > 1e-5 && Pp_z < -1e-5){
+        // now I take tolerance as 1e-5
+        return false;
+    }
+
+    double half_x_max = this -> getXHalfLengthMax();
+    double half_x_min = this -> getXHalfLengthMin();
+    double half_y = this -> getYHalfLength();
+    double p1x = half_x_min; double p1y = -half_y;
+    double p2x = half_x_max; double p2y = half_y;
+    double p3x = -half_x_max; double p3y = half_y;
+    double p4x = -half_x_min; double p4y = -half_y;
+    if( (p2x-p1x)*(Pp_y-p1y) - (p2y-p1y)*(Pp_x-p1x) < -1e-5 ) return false;
+    if( (p3x-p2x)*(Pp_y-p2y) - (p3y-p2y)*(Pp_x-p2x) < -1e-5 ) return false;
+    if( (p4x-p3x)*(Pp_y-p3y) - (p4y-p3y)*(Pp_x-p3x) < -1e-5 ) return false;
+    if( (p1x-p4x)*(Pp_y-p4y) - (p1y-p4y)*(Pp_x-p4x) < -1e-5 ) return false;
+    return true;
+}
diff --git a/GeoModelCore/GeoModelKernel/src/GeoVSurface.cxx b/GeoModelCore/GeoModelKernel/src/GeoVSurface.cxx
index 332646b31bb2eccc9b005e58a3f529b5cf4259b8..29a918b87fc746cb70b86fea6817e7d139a7c752 100644
--- a/GeoModelCore/GeoModelKernel/src/GeoVSurface.cxx
+++ b/GeoModelCore/GeoModelKernel/src/GeoVSurface.cxx
@@ -6,120 +6,7 @@ void GeoVSurface::exec(GeoNodeAction *action) const{
   action->handleVSurface(this);
 }
 
-bool GeoVSurface::isOnSurface(const double Px, const double Py, const double Pz) const{
-  GeoTrf::Vector3D P_prime(Px, Py, Pz);    
-  Eigen::Vector4d P_prime_4d(Px, Py, Pz, 1.0);
-  P_prime_4d = this->getX().inverse() * P_prime_4d;
-  double Pp_x = P_prime_4d[0];
-  double Pp_y = P_prime_4d[1];
-  double Pp_z = P_prime_4d[2];
-  // std::cout << " " << std::endl;
-  // std::cout << " AFTER INVERSE TRANSFORMATION " << std::endl;
-  // std::cout << "Pp_x: " << Pp_x << " Pp_y: " << Pp_y << " Pp_z: " << Pp_z << std::endl;
-// The Point is not even on the same plane as the Annulus Surface
-  if(Pp_z != 0 && Pp_z > 1e-5 && Pp_z < -1e-5){
-    // what is the tolerance for the point to be on the surface?
-    // now I take it as 1e-5
-    // std::cout << "FLAG 1" << std::endl;
-    return false;
-  }
-  if (dynamic_cast<const GeoAnnulusSurface*>(m_surfaceshape.get())){
-    const GeoAnnulusSurface *annulus_shape = dynamic_cast<const GeoAnnulusSurface*>(m_surfaceshape.get());
-
-    // The Annulus Shape starts from theta = 0
-    double real_theta = annulus_shape -> getPhi();
-    int quotient = floor(real_theta/(2.0*M_PI));
-    double theta = real_theta - quotient*2.0*M_PI;
-    
-    if (theta == 0.0 && real_theta > 0.0){
-      theta = 2.0*M_PI;
-    }
-    // std::cout << "theta: " << theta << std::endl;
-
-    double radius_in = annulus_shape -> getRadiusIn();
-    double radius_out = annulus_shape -> getRadiusOut();
-    double Ox = annulus_shape -> getOx();
-    double Oy = annulus_shape -> getOy();
-    
-  // Make sure the Point is located in the general Annulus Ring
-    if (Pp_x * Pp_x + Pp_y * Pp_y - radius_in * radius_in < -1e-5 || Pp_x * Pp_x + Pp_y * Pp_y - radius_out * radius_out > 1e-5){
-      // std::cout << " RADIUS " << Pp_x * Pp_x + Pp_y * Pp_y <<  std::endl;
-      // std::cout << "FLAG 2" << std::endl;
-      return false;
-    }
-    // else if (Pp_x * Pp_x + Pp_y * Pp_y < radius_in * radius_in || Pp_x * Pp_x + Pp_y * Pp_y > radius_out * radius_out){
-    //   std::cout << " RADIUS QUESTIONED" << Pp_x * Pp_x + Pp_y * Pp_y <<  std::endl;
-    // }
-
-  // Make sure the Point is located in the Deviated Circular Sector
-    double Dx = Pp_x - Ox; // X direction distance from the deviation center to the point
-    double Dy = Pp_y - Oy; // Y direction distance from the deviation center to the point
-    // std::cout << "Dx: " << Dx << " Dy: " << Dy << std::endl;
- 
-    double angle = atan2(Dy, Dx); // the range of atan2 is (-pi, pi]
-    if (angle < 0.0){
-      angle += 2*M_PI;
-    }
-    // std::cout << "angle: " << angle << " theta: " << theta << std::endl;
-    if (angle <= theta){
-      return true;
-    }
-    else{
-      // std::cout << "FLAG 3" << std::endl;
-      return false;
-    }
-  }
-  else if (dynamic_cast<const GeoRectSurface*>(m_surfaceshape.get())){
-    const GeoRectSurface* rect_shape = dynamic_cast<const GeoRectSurface*>(m_surfaceshape.get());
-    double half_x = rect_shape -> getXHalfLength();
-    double half_y = rect_shape -> getYHalfLength();
-    double p1x = half_x; double p1y = -half_y;
-    double p2x = half_x; double p2y = half_y;
-    double p3x = -half_x; double p3y = half_y;
-    double p4x = -half_x; double p4y = -half_y;
-    if( (p2x-p1x)*(Pp_y-p1y) - (p2y-p1y)*(Pp_x-p1x) < -1e-5 ) return false;
-    if( (p3x-p2x)*(Pp_y-p2y) - (p3y-p2y)*(Pp_x-p2x) < -1e-5 ) return false;
-    if( (p4x-p3x)*(Pp_y-p3y) - (p4y-p3y)*(Pp_x-p3x) < -1e-5 ) return false;
-    if( (p1x-p4x)*(Pp_y-p4y) - (p1y-p4y)*(Pp_x-p4x) < -1e-5 ) return false;
-    return true;
-  }
-  else if (dynamic_cast<const GeoTrapezoidSurface*>(m_surfaceshape.get())){
-    const GeoTrapezoidSurface* trapezoid_shape = dynamic_cast<const GeoTrapezoidSurface*>(m_surfaceshape.get());
-    double half_x_max = trapezoid_shape -> getXHalfLengthMax();
-    double half_x_min = trapezoid_shape -> getXHalfLengthMin();
-    double half_y = trapezoid_shape -> getYHalfLength();
-    double p1x = half_x_min; double p1y = -half_y;
-    double p2x = half_x_max; double p2y = half_y;
-    double p3x = -half_x_max; double p3y = half_y;
-    double p4x = -half_x_min; double p4y = -half_y;
-    if( (p2x-p1x)*(Pp_y-p1y) - (p2y-p1y)*(Pp_x-p1x) < -1e-5 ) return false;
-    if( (p3x-p2x)*(Pp_y-p2y) - (p3y-p2y)*(Pp_x-p2x) < -1e-5 ) return false;
-    if( (p4x-p3x)*(Pp_y-p3y) - (p4y-p3y)*(Pp_x-p3x) < -1e-5 ) return false;
-    if( (p1x-p4x)*(Pp_y-p4y) - (p1y-p4y)*(Pp_x-p4x) < -1e-5 ) return false;
-    return true;
-  }
-  else if (dynamic_cast<const GeoDiamondSurface*>(m_surfaceshape.get())){
-    const GeoDiamondSurface* diamond_shape = dynamic_cast<const GeoDiamondSurface*>(m_surfaceshape.get());
-    double x_bot = diamond_shape -> getXbottomHalf(); double y_bot = diamond_shape -> getYbottomHalf();
-    double x_mid = diamond_shape -> getXmidHalf();
-    double x_top = diamond_shape -> getXtopHalf(); double y_top = diamond_shape -> getYtopHalf();
-    double p1x = x_bot; double p1y = -y_bot;
-    double p2x = x_mid; double p2y = 0.0;
-    double p3x = x_top; double p3y = y_top;
-    double p4x = -x_top; double p4y = y_top;
-    double p5x = -x_mid; double p5y = 0.0;
-    double p6x = -x_bot; double p6y = -y_bot;
-    if( (p2x-p1x)*(Pp_y-p1y) - (p2y-p1y)*(Pp_x-p1x) < -1e-5 ) return false;
-    if( (p3x-p2x)*(Pp_y-p2y) - (p3y-p2y)*(Pp_x-p2x) < -1e-5 ) return false;
-    if( (p4x-p3x)*(Pp_y-p3y) - (p4y-p3y)*(Pp_x-p3x) < -1e-5 ) return false;
-    if( (p5x-p4x)*(Pp_y-p4y) - (p5y-p4y)*(Pp_x-p4x) < -1e-5 ) return false;
-    if( (p6x-p5x)*(Pp_y-p5y) - (p6y-p5y)*(Pp_x-p5x) < -1e-5 ) return false;
-    if( (p1x-p6x)*(Pp_y-p6y) - (p1y-p6y)*(Pp_x-p6x) < -1e-5 ) return false;
-    return true;
-  }
-  else{
-    THROW_EXCEPTION("The Surface Shape is Unknown!");
-  }
-  // std::cout << "FLAG 10" << std::endl;
-  return false;
-}
+bool GeoVSurface::isOnSurface (const double Px, const double Py, const double Pz) const{
+  GeoTrf::Transform3D trans = this->getX();
+  return m_surfaceshape.get()->isOnSurface(Px, Py, Pz, trans);
+}
\ No newline at end of file
diff --git a/GeoModelExamples/SurfAnnulusDemo/src/SurfAnnulusDemo.cxx b/GeoModelExamples/SurfAnnulusDemo/src/SurfAnnulusDemo.cxx
index 1a21ac27f564d45c9c37f8bc755ef0b22c8d6ae8..1cb6cd48c6dad55bc34504f2a07cebc6ce8f99ee 100644
--- a/GeoModelExamples/SurfAnnulusDemo/src/SurfAnnulusDemo.cxx
+++ b/GeoModelExamples/SurfAnnulusDemo/src/SurfAnnulusDemo.cxx
@@ -143,14 +143,11 @@ void SurfAnnulusDemo::create(GeoVPhysVol *world, bool /*publish*/) {
     GeoLogVol  *boxLog = new GeoLogVol("BoxLog",boxShape,Air);
     GeoFullPhysVol *boxPhys3 = new GeoFullPhysVol(boxLog);
 
-    GeoAnnulusSurface* annulus = new GeoAnnulusSurface(-4.0, -5.0, 7.0, 10.0, 2*M_PI); // Ox, Oy, Rin, Rout, Phi
+    GeoAnnulusSurface* annulus = new GeoAnnulusSurface(-4.0, -5.0, 7.0, 10.0, 0.3*M_PI); // Ox, Oy, Rin, Rout, Phi
     GeoVSurface* surf3 = new GeoVSurface(annulus);
     GeoAlignableTransform *move_together3=new GeoAlignableTransform(GeoTrf::RotateY3D(0.4*M_PI)*GeoTrf::TranslateX3D(20));
     //GeoTransform * point_trans = new GeoTransform(GeoTrf::RotateY3D(0.4*M_PI)*GeoTrf::TranslateX3D(20));
 
-
-    boxPhys->add(move_together3);
-    boxPhys->add(boxPhys3);
     boxPhys->add(move_together3);   
     boxPhys->add(surf3);
     // boxPhys->add(transform_point); // move the point (initially at (0,0,0)) to the starting position (5.66, -5.66, 0)
diff --git a/cmake/GeoModel-version.cmake b/cmake/GeoModel-version.cmake
index e45a1f3b5f79205afc9978fa6de61a00ba626cf7..80eb7d554e0435b1d28188d4f45cec1c832493ce 100644
--- a/cmake/GeoModel-version.cmake
+++ b/cmake/GeoModel-version.cmake
@@ -1,5 +1,5 @@
 # Set up the version of GeoModel as a cache variable, so that other
 # sub-projects could use this value.
-set( GeoModel_VERSION "7.0.0" CACHE STRING
+set( GeoModel_VERSION "6.6.0" CACHE STRING
     "Version of the GeoModel project" )